In [1]:

    

# Definimos el cuerpo donde vamos a trabajar
m = 4;rango = 3;N = 2^m
K_.<a> = GF(2)
F.<a> = GF(2^m) #creado el campo donde "vivirá" el código


    

In [2]:

    

# Creamos el anillo de polinomios
PR = PolynomialRing(F,'X')
X = PR.gen()
g = X^3+X+1                     # Polinomio de Goppa
L = [a^i for i in range(N)]     # Soporte del codigo


    

In [3]:

    

# -------------------
# Funcion auxiliar que permite descomponer un polinomio en irreducibles
# -------------------
def descomponer_polinomio(p):
    # El siguiente metodo permite descomponer un polinomio p en factores irreducibles p(z) = p0 (z) + z p1 (z)
    # Entrada: Polinomio p
    Phi1 = p.parent()
    p0 = Phi1([sqrt(c) for c in p.list()[0::2]])
    p1 = Phi1([sqrt(c) for c in p.list()[1::2]])
    return (p0,p1)

# -------------------
# Algoritmo de euclides extendido: Obtener MCD y los s,t que lo generan.
# -------------------
def algoritmo_euclides_extendido(self, other):
    delta = self.degree() #grado de polinomio 1
    if other.is_zero(): # si el polinomio introducido es
        ring = self.parent() #comprobamos el cuerpo en el que trabajamos
        return self, R.one(), R.zero() #mcd = mismo polinomio y devuelve un uno (s) y un cero (t) en el cuerpo que trabajamos.

    # mcd (a,b) = as+bt

    ring = self.parent() #comprobamos el cuerpo en el que trabajamos
    a = self # guardamos una copia del primer polinomio 1 (self)
    b = other # guardamos una copia del segundo polinomio (other)

    s = ring.one() # guardamos en s el uno del anillo
    t = ring.zero() # guardamos en t el cero del anillo

    resto0 = a
    resto1 = b

    while true:
        cociente,resto_auxiliar = resto0.quo_rem(resto1) # La funcion quo_rem de Sage devuelve el cociente y el resto. Que guardamos en Q y ring.
        resto0 = resto1
        resto1 = resto_auxiliar

        s = t
        t = s - t*cociente

        if resto1.degree() <= floor((delta-1)/2) and resto0.degree() <= floor((delta)/2):
             break

    V = (resto0-a*s)//b
    coeficiente_lider = resto0.leading_coefficient() # guardamos el coeficiente lider del resto 0

    return resto0/coeficiente_lider, s/coeficiente_lider, V/coeficiente_lider

# -------------------
# Funcion que calcula la inversa de un polinomio utilizando el algoritmo de euclides de Sage
# -------------------
def inversa_g(p,g):
    (d,u,v) = xgcd(p,g)
    return u.mod(g)

# -------------------
# Funcion de decodificacion de Patterson
# -------------------
def decodePatterson(y):
    # Calculamos primero el vector alpha con los elementos primitivos.
    alpha = vector(H*y)

    # Consideramos nuestras matrices T,Y,Z definidas así como nuestro polinomio irreducible g

    polinomioS = PR(0) # Inicializa como el polinomio 0 del anillo
    for i in range(len(alpha)):
        polinomioS = polinomioS + alpha[i]*(X^(len(alpha)-i-1)) # Lo vamos rellenando con los alpha

    vector_g = descomponer_polinomio(g) # Guardamos en vector_g el par de polinomios irreducibles
    w = ((vector_g[0])*inversa_g(vector_g[1],g)).mod(g)
    vector_t = descomponer_polinomio(inversa_g(polinomioS,g) + X)

    R = (vector_t[0]+(w)*(vector_t[1])).mod(g)

    (a11,b11,c11) = algoritmo_euclides_extendido(g,R)

    # Definimos el polinomio sigma
    sigma = a11^2+X*(c11^2)

    # Vamos comprobando uno a uno los coeficientes de sigma
    # para asi determinar el conjunto de posiciones de error E - {i tal que e_i es distinto de 0}
    for i in range(N):
        if (sigma(a^i)==0): # an error occured
            print ("Error encontrado en la posición: " + str(i))
            y[i] = y[i] + 1
    return y


    

In [4]:

    

T = matrix(F,rango,rango)
for i in range(rango):
    count = rango - i
    for j in range(rango):
        if i > j:
            T[i,j]=g.list()[count]
            count = count + 1
        if i < j:
            T[i,j] = 0
        if i == j:
            T[i,j] = 1


    

In [5]:

    

print ("Matriz T: ")
show(T)


    

    




Matriz T: 






    
Out[5]:




 $$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr}
1 & 0 & 0 \\
0 & 1 & 0 \\
1 & 0 & 1
\end{array}\right)$$ 

In [6]:

    

V = matrix([[L[j]^i for j in range(N)] for i in range(rango)])
D = diagonal_matrix([1/g(L[i]) for i in range(N)])
H = T*V*D
H_Goppa_K = matrix(K_, m*H.nrows(),H.ncols())
for i in range(H.nrows()):
    for j in range(H.ncols()):
        be = bin(eval(H[i,j]._int_repr()))[2:];
        be = '0'*(m-len(be))+be; be = list(be);
        H_Goppa_K[m*i:m*(i+1),j]=vector(map(int,be));
show(H_Goppa_K)


    

    
Out[6]:




 $$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrrrrrrrrrr}
0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 1 & 0 & 1 & 1 & 1 & 0 & 1 & 0 & 1 & 1 & 1 & 0 & 1 & 1 & 0 \\
0 & 0 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 0 \\
1 & 1 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 1 & 1 & 0 & 1 \\
0 & 1 & 1 & 1 & 1 & 0 & 1 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 1 & 1 & 0 & 1 & 0 \\
1 & 0 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 1 & 1 \\
0 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 0 & 1 & 0 & 1 & 1 & 1 & 1 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 \\
0 & 1 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 1 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & 1 & 0 & 0
\end{array}\right)$$ 

In [7]:

    

Krnl = H_Goppa_K.right_kernel();
G = Krnl.basis_matrix();
S = random_matrix(GF(2), N-m*rango)

while (S.determinant()==0):
    S = random_matrix(GF(2), N-m*rango)

rng = range(N)

P = matrix(GF(2),N);

for i in range(N):
    p = floor(len(rng)*random());
    P[i,rng[p]]=1; rng=rng[:p]+rng[p+1:];

G_prima = S*G*P
show(G_prima)


    

    
Out[7]:




 $$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrrrrrrrrrr}
1 & 1 & 0 & 1 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 1 & 1 & 0 & 1 & 0 \\
0 & 0 & 1 & 0 & 1 & 1 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 1 & 0 & 1 \\
0 & 0 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 1 & 1 \\
1 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 1 & 0
\end{array}\right)$$ 

In [8]:

    

u = vector(K_,[randint(0,1) for _ in range(G_prima.nrows())])
c = u*G_prima
print 'Vector u' ; show(u)
print 'Vector c' ; show(c)
e = vector(K_,N)
e[8] = 1
e[9] = 1
print 'Vector de errores e' ; show(e)
y = c + e
print "Vector codificado y" ; show(y)


    

    




Vector u






    
Out[8]:




 $$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0,\,1,\,0,\,1\right)$$ 






    




Vector c






    
Out[8]:




 $$\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1\right)$$ 






    




Vector de errores e






    
Out[8]:




 $$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,1,\,1,\,0,\,0,\,0,\,0,\,0,\,0\right)$$ 






    




Vector codificado y






    
Out[8]:




 $$\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,0,\,0,\,1,\,1,\,1,\,1,\,1,\,1\right)$$ 

In [9]:

    

yP = y*(P.inverse())
yd = decodePatterson(yP)
corregido = (G.transpose()\yd)*S.inverse()
show(corregido)


    

    




Error encontrado en la posición: 1
Error encontrado en la posición: 14






    
Out[9]:




 $$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0,\,1,\,0,\,1\right)$$