TP1 : codes de Reed-Solomon

Ce TP a pour but d'implanter effectivement un code de Reed-Solomon et son algorithme de décodage.

Exercice 1 : Fonctions de base

Question : Implanter une fonction random_message(F, k) qui retourne un vecteur aléatoire tiré uniformément dans l'esspace vectoriel sur F de dimentsion k.

In [286]:
def random_message(F, k):
    return vector(F, [F.random_element() for i in range(k)])

Question : Implanter une fonction random_error(F, n, t) qui retourne un vecteur d'erreur de poids t et de longueur n sur le corps fini F. (on pourra supposer que le nombre d'erreurs est assez faible comparé à la longueur, typiquement $t \le n$)

In [287]:
def random_error(F, n, t):
    e = [F.zero() for j in range(n)]
    i = 0
    while i < t:
        j = randint(0, n-1)
        if e[j] == 0:
            a = F.random_element()
            while a == 0:
                a = F.random_element()
            e[j] = a
        i += 1
    return vector(e)

Question : En déduire une fonction corrupt(c, t) qui ajoute une erreur de poids t à un mot c.

In [288]:
def corrupt(c, t):
    F = c[0].parent()
    return c + random_error(F, len(c), t)

Question : Implanter une fonction parity_check_matrix(G) qui retourne une matrice de parité du code engendré par la matrice génératrice G.

In [289]:
def parity_check_matrix(G):
    return G.right_kernel_matrix()
In [290]:
def TEST_PARITY_CHECK_MATRIX():
    r, k, n = 0, 6, 10
    G = random_matrix(GF(7), k, n)
    print(G * parity_check_matrix(G).transpose() == 0)

Question : Implanter une fonction systematic_form(G) qui retourne :

  • une matrice génératrice $G'$ sous forme systématique du code engendré par la matrice G
  • les indices des colonnes de $G'$ sur lequelles $G'$ est égale à l'identité
In [291]:
def systematic_form(G):
    return G.rref(), G.pivots()
In [292]:
def TEST_SYSTEMATIC_FORM():
    r, k, n = 0, 6, 10
    while r != k:
        G = random_matrix(GF(7), k, n)
        r = G.rank()
    Gp, I = systematic_form(G)
    print(G.row_space() == Gp.row_space() and Gp[:,I] == identity_matrix(GF(7), k))

Question : Implanter une fonction uncode(G, c) qui prend en entrée une matrice génératrice G d'un code $\mathcal{C}$ et un mot c de $\mathcal{C}$, et qui retourne le message associé au mot c.

In [293]:
def uncode(G, c):
    return c/G
In [294]:
def TEST_UNCODE():
    r, k, n = 0, 6, 10
    while r != k:
        G = random_matrix(GF(7), k, n)
        r = G.rank()
    m = random_message(GF(7), k)
    c = m * G
    print(uncode(G, c) == m)
In [295]:
TEST_PARITY_CHECK_MATRIX()
TEST_SYSTEMATIC_FORM()
TEST_UNCODE()

True
True
True

Exercice 2 : Codes de Reed--Solomon

Question : Implanter une fonction generator_matrix_reed_solomon(x, k) qui prend en entrée un vecteur x de points d'évaluation et un entier k, et qui retourne une matrice génératrice du code de Reed--Solomon $\mathrm{RS}_k(\mathbf{x})$.

In [296]:
def generator_matrix_reed_solomon(x, k):
    return matrix([[a**j for a in x] for j in range(k)])

Question : Implanter une fonction random_evaluation_points(F, n) qui prend en entrée un corps fini F (de cardinal $q$) et un entier n ($\le q$), et qui retourne un vecteur aléatoire de longueur n formé d'éléments de F deux à deux distincts.

In [297]:
def random_evaluation_points(F, n):
    L = F.list()
    q = len(L)
    x = []
    for j in range(n):
        i = randint(0, q-1)
        a = L.pop(i)
        x.append(a)
        q -= 1
    return vector(x)

Question : Implanter une fonction get_evaluation_vector(P, x) qui prend en entrée un polynôme P et des points d'évaluation xet qui retourne le vecteur d'évaluation de P sur x.

In [298]:
def get_evaluation_vector(P, x):
    return vector([P(a) for a in x])

Exercice 3 : Décodage de Gao

Question : Implanter une fonction interpolate(x, y) qui prend en entrée un vecteur de points d'évaluation x et un vecteur y de même longueur, et qui retourne le polynôme $P$ de plus petit degré tel que $P(x_i) = y_i$ pour tout $i$.

In [299]:
def interpolate(x, y):
    n = len(x)
    R = PolynomialRing(x[0].parent(), "X")
    X = R.gen()
    res = R.zero()
    for i in range(n):
        P = R.one()
        for j in range(n):
            if i != j:
                P *= (X-x[j])/(x[i]-x[j])
        res += (y[i] * P)    
    return res 
    
    # on peut aussi utiliser la fonction sage :
    # return R.lagrange_polynomial([[x[i], y[i]] for i in range(n)])

Question : Implanter une fonction partial_xgcd(A, B, s) qui prend en entrée deux polynômes A et B, qui effectue l'algorithme d'Euclide étendu jusque ce que le degré du reste soit $<s$, puis qui retourne les coefficients de Bezout et le reste associé.

In [300]:
def partial_xgcd(A, B, s):
    R0, R1, U0, U1, V0, V1 = A, B, 1, 0, 0, 1
    while R1.degree() >= s:
        Q = R0.quo_rem(R1)[0]
        R0, R1 = R1, R0 - Q*R1
        U0, U1 = U1, U0 - Q*U1
        V0, V1 = V1, V0 - Q*V1
    return R1, U1, V1

Question : Implanter l'algorithme de Gao tel que vu en TD. La fonction à implanter sera notée gao_decoding(y, x, k), et retournera un message sous la forme d'un vecteur $\mathbf{m} = (m_0, \dots, m_{k-1})$ tel que y est le mot le plus proche de get_evaluation_vector(P, x) où $P(X) = \sum_{i=0}^{k-1} m_i X^i$.

In [301]:
def gao_decoding(y, x, k):
    n = len(y)
    F = y[0].parent()
    PolyRing = PolynomialRing(F, "X")
    X = PolyRing.gen()
    A = PolyRing.one()
    for i in range(n):
        A *= (X - x[i])
    Y = interpolate(x, y)
    s = (k+n)//2
    R, U, V = partial_xgcd(A, Y, s)
    C = R//V
    return vector(C.coefficients(sparse=False))

Question : Tester votre algorithme avec la fonction suivante :

In [302]:
def TEST_GAO():
    q, n, k = 31, 26, 12
    t = (n-k)//2
    F = GF(q)
    x = random_evaluation_points(F, n)
    G = generator_matrix_reed_solomon(x, k)
    m = random_message(F, k)
    c = m * G
    y = corrupt(c, t)
    m_dec = gao_decoding(y, x, k)
    print(m == m_dec)

TEST_GAO()

True