El algoritmo de Berlekamp-Massey

Los síndromes de un código BCH dan un sistema lineal para los coeficientes del polinomio localizador. El método de Peterson lo resuelve con determinantes, y su costo crece como $t^3$. El algoritmo iterativo llega al mismo $\sigma (X)$ en del orden de $t^2$ operaciones y sin un solo determinante, porque mira el problema de otra manera: busca el registro de desplazamiento más corto que genera los síndromes. E. R. Berlekamp lo presentó en 1965 para los BCH binarios, y J. L. Massey lo leyó como síntesis de registros en 1969.

01 · El localizador es un registro

Sea un código BCH binario de sentido estricto que corrige $t$ errores, con síndromes $S_1,\dots ,S_{2t}$. Si hubo $\nu \le t$ errores con números localizadores $X_1,\dots ,X_{\nu }$, el localizador es

$$\sigma (X)=\prod _{l=1}^{\nu }(1+X_{l}X)=1+{\sigma }_{1}X+\cdots +{\sigma }_{\nu }{X}^{\nu }$$

y las identidades de Newton, para $i>\nu$, se escriben como una recurrencia:

$$S_i={\sigma }_1{S}_{i-1}+{\sigma }_2{S}_{i-2}+\cdots +{\sigma }_{\nu }{S}_{i-\nu }$$

Cada síndrome es una combinación de los $\nu$ anteriores, con coeficientes fijos. Es lo que hace un LFSR de $\nu$ etapas con polinomio de conexión $\sigma (X)$: cargado con $S_1,\dots ,S_{\nu }$, saca el resto de la sucesión. La única diferencia con los de la Unidad 5 es que las etapas guardan elementos de $GF(2^m)$, no bits.

Encontrar $\sigma (X)$ pasa a ser un problema de síntesis: el registro más corto que genera $S_1,\dots ,S_{2t}$. Más corto, porque el patrón más probable es el de menos errores, y el grado de $\sigma (X)$ es el número de errores. El algoritmo lo arma de a una identidad por vez: arranca con un registro que genera el primer síndrome y, en cada paso, lo corrige lo mínimo indispensable para que genere también el siguiente.

02 · La discrepancia y la elección de ρ

El algoritmo produce una sucesión de polinomios ${\sigma }^{(\mu )}(X)$, cada uno el de menor grado que satisface las primeras $\mu$ identidades. Su grado es $l_{\mu }$. Para pasar al siguiente se mide cuánto falla ${\sigma }^{(\mu )}$ en la identidad número $\mu +1$. Esa medida es la discrepancia:

$$d_{\mu }=S_{\mu +1}+{\sigma }_1^{(\mu )}{S}_{\mu }+\cdots +{\sigma }_{l_{\mu }}^{(\mu )}{S}_{\mu +1-l_{\mu }}$$

Es lo que el registro predice para $S_{\mu +1}$, sumado a lo que vale. Si $d_{\mu }=0$, el registro acierta y pasa igual: ${\sigma }^{(\mu +1)}={\sigma }^{(\mu )}$, $l_{\mu +1}=l_{\mu }$. Si no, se busca un paso anterior $\rho$ con $d_{\rho }\ne 0$ y $\rho -l_{\rho }$ lo más grande posible, y se corrige con él:

$${\sigma }^{(\mu +1)}(X)={\sigma }^{(\mu )}(X)+d_{\mu }{d}_{\rho }^{-1}{X}^{\mu -\rho }{\sigma }^{(\rho )}(X)$$

y el grado pasa a ser

$$l_{\mu +1}=\max (l_{\mu },l_{\rho }+\mu -\rho )$$

Por qué el término agregado arregla la identidad $\mu +1$ sin romper las anteriores:

Paso 1: correr es atrasar la identidad. Multiplicar ${\sigma }^{(\rho )}$ por $X^{\mu -\rho }$ corre sus coeficientes $\mu -\rho$ lugares. Lo que el polinomio corrido aporta a la identidad $j$ es lo que ${\sigma }^{(\rho )}$ aportaba a la identidad $j-(\mu -\rho )$.

Paso 2: en las identidades viejas aporta cero. Para $j\le \mu$, la identidad atrasada tiene número $\le \rho$, y ${\sigma }^{(\rho )}$ las cumplía todas. Lo que ${\sigma }^{(\mu )}$ ya acertaba, lo sigue acertando.

Paso 3: en la nueva aporta exactamente la discrepancia. Para $j=\mu +1$, la identidad atrasada es la $\rho +1$, donde ${\sigma }^{(\rho )}$ fallaba por $d_{\rho }$. Escalado por $d_{\mu }{d}_{\rho }^{-1}$, aporta $d_{\mu }$, que se cancela con la discrepancia de ${\sigma }^{(\mu )}$. En característica dos, sumar es restar.

Cualquier $\rho$ con $d_{\rho }\ne 0$ arregla la identidad. La elección decide el grado: $l_{\rho }+\mu -\rho$ es más chico cuanto más grande es $\rho -l_{\rho }$. Massey probó que con esa elección el registro es, en cada paso, el más corto posible: su largo cumple $l_{\mu +1}=\max (l_{\mu },\mu +1-l_{\mu })$, y no hay registro más corto que genere esos $\mu +1$ términos. Al cabo de $2t$ pasos, si hubo a lo sumo $t$ errores, ${\sigma }^{(2t)}(X)$ es el localizador.

El cálculo se ordena en una tabla con una fila por paso: $\mu$, ${\sigma }^{(\mu )}$, $d_{\mu }$, $l_{\mu }$ y $\mu -l_{\mu }$, la columna que sirve para elegir $\rho$. Las dos primeras filas no son iteraciones, son el arranque: $\mu =-1$ con ${\sigma }^{(-1)}=1$ y $d_{-1}=1$, que existe para que el primer paso tenga un $\rho$ del cual tirar; y $\mu =0$ con ${\sigma }^{(0)}=1$ y $d_0=S_1$.

03 · En binario, la mitad de los pasos no hace nada

En un código binario los síndromes pares son cuadrados: $S_{2i}=S_i^2$. Con eso se prueba que, si ${\sigma }^{(\mu )}$ cumple las identidades hasta un orden $\mu$ par, cumple también la siguiente: la discrepancia de todo paso de índice impar es cero, y ese paso copia la fila. Se ve en la forma general aplicada al (15,5) con errores en las posiciones 3, 5 y 12:

La columna $\rho$ dice con qué fila se armó cada polinomio. Las filas 1, 3 y 5 tienen discrepancia nula, y no es suerte de este patrón: con la forma general corrida sobre los 575 patrones de uno, dos y tres errores del (15,5), las discrepancias de índice impar dieron cero en todos, y el $\sigma (X)$ final fue el mismo que el de la versión binaria.

La versión binaria recorre sólo los pasos pares. Se los reindexa con $\mu =-\frac{1}{2},0,1,\dots ,t$, y el algoritmo termina en $t$ iteraciones en lugar de $2t$. Cambian dos cosas. La discrepancia se calcula contra los síndromes impares:

$$d_{\mu }=S_{2\mu +1}+{\sigma }_1^{(\mu )}{S}_{2\mu }+\cdots +{\sigma }_{l_{\mu }}^{(\mu )}{S}_{2\mu +1-l_{\mu }}$$

Y la corrección se corre el doble, con $\rho$ elegido por el mayor $2\rho -l_{\rho }$:

$${\sigma }^{(\mu +1)}(X)={\sigma }^{(\mu )}(X)+d_{\mu }{d}_{\rho }^{-1}{X}^{2(\mu -\rho )}{\sigma }^{(\rho )}(X)$$

El arranque es $\mu =-\frac{1}{2}$, con $\sigma =1$, $d=1$, $l=0$, y $\mu =0$, con $\sigma =1$ y $d_0=S_1$. El medio del primero es para que el corrimiento $2(\mu -\rho )$ desde él dé un número impar: es la fila $-1$ de la forma general, contada en pasos dobles. La simplificación vale sólo para códigos binarios. Para uno con símbolos en $GF(2^m)$, como un Reed-Solomon, $S_{2i}=S_i^2$ deja de valer y hay que dar los $2t$ pasos.

04 · La Tabla 6, fila por fila

El (15,5) corrige $t=3$ errores sobre $GF(2^4)$, con $1+X+X^4$. Se transmite la palabra nula y llega $r(X)=X^3+X^5+X^{12}$. Los seis síndromes son

$$S_1=1,S_2=1,S_3={\alpha }^{10},S_4=1,S_5={\alpha }^{10},S_6={\alpha }^5$$

con los pares como cuadrados de los impares: $S_6={({\alpha }^{10})}^2={\alpha }^{20}={\alpha }^5$. El algoritmo binario da tres iteraciones:

Fila 0. $d_0=S_1=1$. El único paso previo con discrepancia no nula es $\rho =-\frac{1}{2}$, y el corrimiento es $2(0+\frac{1}{2})=1$: ${\sigma }^{(1)}=1+X$, de grado 1.

Fila 1. $d_1=S_3+{\sigma }_1{S}_2={\alpha }^{10}+1={\alpha }^5$. Los candidatos son $\rho =-\frac{1}{2}$, con $2\rho -l_{\rho }=-1$, y $\rho =0$, con $0$. Gana el segundo, corrimiento 2: ${\sigma }^{(2)}=1+X+{\alpha }^5{X}^2$.

Fila 2. $d_2=S_5+{\sigma }_1{S}_4+{\sigma }_2{S}_3={\alpha }^{10}+1+{\alpha }^5{\alpha }^{10}={\alpha }^{10}$, porque ${\alpha }^{15}=1$ cancela el 1. Ahora gana $\rho =1$, con $2\rho -l_{\rho }=1$; el factor es ${\alpha }^{10}/{\alpha }^5={\alpha }^5$ y el corrimiento 2:

$${\sigma }^{(3)}(X)=(1+X+{\alpha }^5{X}^2)+{\alpha }^5{X}^2(1+X)=1+X+{\alpha }^5{X}^3$$

Son las filas pares de la tabla de la forma general, renumeradas. El grado final, 3, es el número de errores. Las raíces de $1+X+{\alpha }^5{X}^3$ son ${\alpha }^3$, ${\alpha }^{10}$ y ${\alpha }^{12}$; sus inversos, ${\alpha }^{12}$, ${\alpha }^5$ y ${\alpha }^3$, son los localizadores, y dan las posiciones 12, 5 y 3. Encontrarlas es la búsqueda de Chien. Peterson, con su fórmula cerrada para $t=3$, da el mismo polinomio; y comparados en los 575 patrones de hasta tres errores del (15,5), los dos coinciden siempre.

Antes de usar el widget, pensá la fila 1 de otra manera: ¿qué pasa si ahí se elige $\rho =-\frac{1}{2}$, el candidato que pierde? ¿Y si el que se elige mal es el de la fila 2?

Arranca con el patrón de la Tabla 6 y sólo las dos filas del arranque. Hacela entera eligiendo vos. Después volvé a empezar y, en la fila 2, elegí $\rho =0$; y otra vez con $\rho =-\frac{1}{2}$ en la fila 1. Probá con dos errores, en las posiciones 4 y 9, y con cuatro. Al final, prendé «sin inversiones» y mirá qué cambia y qué no.

05 · Lo que dice cada columna

$l_{\mu }$ es el largo del registro, y sólo crece. Al terminar es el grado de $\sigma (X)$, y el decodificador lo lee como el número de errores. La columna $2\mu -l_{\mu }$ está para elegir $\rho$ de un vistazo: entre las filas con $d_{\mu }\ne 0$, la del número más grande.

Con menos errores que $t$, las últimas discrepancias son cero. Con errores en 4 y 9, $d_0={\alpha }^{14}$, $d_1={\alpha }^{12}$ y $d_2=0$: la última iteración copia $1+{\alpha }^{14}X+{\alpha }^{13}{X}^2$, de grado 2. Un circuito de latencia fija la hace igual.

Elegir mal $\rho$ no siempre se nota, y por eso el criterio importa. Sobre el patrón de la Tabla 6:

El primer desvío llega bien por casualidad. Con $\rho =-\frac{1}{2}$ en la fila 1 el corrimiento es 3, y sale ${\sigma }^{(2)}=1+X+{\alpha }^5{X}^3$: genera los primeros cuatro síndromes, pero con $l_2=3$ donde alcanzaba 2. Ese registro de más resultó ser el localizador, y $d_2=0$. En la fila 2 ya no hay margen: un registro más largo que el necesario es un $\sigma (X)$ que no es el localizador, y el decodificador o declara falla o, con el de grado 5, encuentra una raíz que no es un error.

06 · La versión sin inversiones

La corrección lleva el factor $d_{\mu }{d}_{\rho }^{-1}$: una inversión en $GF(2^m)$ en cada iteración con discrepancia no nula, dentro del lazo. En hardware la inversión es la operación lenta, y queda en el camino crítico. La salida es que a $\sigma (X)$ se lo puede multiplicar por cualquier constante no nula sin que cambien sus raíces. Multiplicando la actualización por $d_{\rho }$:

$${\sigma }^{(\mu +1)}(X)=d_{\rho }{\sigma }^{(\mu )}(X)+d_{\mu }{X}^{2(\mu -\rho )}{\sigma }^{(\rho )}(X)$$

Sólo productos. Esta versión, de H. O. Burton en 1971, cambia la inversión y el producto del factor por un producto más en cada coeficiente de ${\sigma }^{(\mu )}$: todos los coeficientes hacen la misma cuenta, y por eso se arma como una fila de celdas iguales.

Que todo siga cerrando tiene un paso no evidente: las filas guardadas también quedan escaladas, cada una por su propio factor, y la discrepancia, que es lineal en $\sigma$, se escala con el mismo factor que su fila. Si la fila $\mu$ tiene factor $c_{\mu }$ y la $\rho$ tiene $c_{\rho }$, la cuenta sin inversiones da

$$c_{\rho }{d}_{\rho }·c_{\mu }{\sigma }^{(\mu )}+c_{\mu }{d}_{\mu }·X^{2(\mu -\rho )}{c}_{\rho }{\sigma }^{(\rho )}=c_{\mu }{c}_{\rho }{d}_{\rho }{\sigma }^{(\mu +1)}$$

el polinomio con inversiones por un factor no nulo. En la Tabla 6 las dos primeras correcciones salen con $d_{\rho }=1$ y no cambian nada; la tercera usa $d_1={\alpha }^5$:

$${\sigma }^{(3)}(X)={\alpha }^5+{\alpha }^{5}X+{\alpha }^{10}{X}^3={\alpha }^5(1+X+{\alpha }^5{X}^3)$$

En el ejemplo integral, el de la Tabla 8, el factor se acumula antes: la discrepancia de la fila 2 sale ${\alpha }^8$ en lugar de ${\alpha }^7$, porque la fila ya venía multiplicada por $\alpha$, y el polinomio final es ${\alpha }^2$ veces el de la tabla. Comparadas las dos versiones en los 575 patrones de uno, dos y tres errores del (15,5), el polinomio sin inversiones fue el otro por una constante no nula en todos, con las mismas raíces. Lo mismo en los 120 patrones del (15,7), los 496 del (31,21) y los 4991 del (31,16).

07 · Con más de t errores

El teorema de Massey promete el registro más corto, no el localizador. Con más de $t$ errores, el registro más corto que genera los $2t$ síndromes suele ser otro, y el algoritmo lo entrega igual, sin aviso. Lo que queda es mirar el resultado: el decodificador declara falla si el grado pasa de $t$, o si Chien no encuentra tantas raíces distintas como el grado. En los 1365 patrones de cuatro errores del (15,5):

Los 840 declarados y los 525 entregados en silencio son el mismo reparto que da cualquier decodificación incompleta de este código: los 525 son los patrones a distancia 3 de alguna de las 15 palabras de peso 7, y para ellos el registro más corto es un localizador legítimo, el de esa otra palabra. Entre los 90 de grado 4 hay 15 en que las cuatro raíces son justo las cuatro posiciones erradas, y el decodificador falla igual: con seis síndromes no hay garantía de que un registro de cuatro etapas sea el bueno, y no lo usa.

08 · La cuenta en Python

El algoritmo binario sobre $GF(2^4)$ con tablas de exponentes y logaritmos, en sus dos versiones: la única diferencia está en los dos coeficientes a y b de la actualización. Corre la Tabla 6 y después compara las dos versiones en los 575 patrones. Lo de los comentarios es la salida de correrlo.

from itertools import combinations

# GF(2^4) con 1 + X + X^4: EXP[i] es α^i como entero de 4 bits.
EXP, LOG = [1] * 30, {}
for i in range(1, 30):
    x = EXP[i - 1] << 1
    EXP[i] = x ^ 0b10011 if x & 0b10000 else x
for i in range(15):
    LOG[EXP[i]] = i

def mul(a, b):
    return 0 if a == 0 or b == 0 else EXP[LOG[a] + LOG[b]]

def inv(a):
    return EXP[(15 - LOG[a]) % 15]

def alfa(x):
    return '0' if x == 0 else '1' if x == 1 else 'α' if x == 2 else f'α^{LOG[x]}'

def sindromes(posiciones, t=3):
    S = [0] * (2 * t + 1)                  # S[i] es S_i; S[0] no se usa
    for i in range(1, 2 * t + 1):
        for j in posiciones:
            S[i] ^= EXP[(i * j) % 15]
    return S

def berlekamp(S, t=3, sin_inversiones=False):
    """La versión binaria: filas (μ, σ, d, l), con μ = -1/2 como -0.5."""
    filas = [(-0.5, [1], 1, 0)]
    sigma, l = [1], 0
    for mu in range(t + 1):
        d = None
        if mu < t:                         # d_μ = Σ σ_i S_{2μ+1-i}, con σ_0 incluido
            d = 0
            for i, c in enumerate(sigma[:2 * mu + 1]):
                d ^= mul(c, S[2 * mu + 1 - i])
        filas.append((mu, sigma, d, l))
        if not d:
            continue
        rho = max((f for f in filas[:-1] if f[2]), key=lambda f: 2 * f[0] - f[3])
        k = int(2 * (mu - rho[0]))         # el corrimiento 2(μ - ρ)
        if sin_inversiones:
            a, b = rho[2], d               # σ' = d_ρ σ + d_μ X^k σ^(ρ)
        else:
            a, b = 1, mul(d, inv(rho[2]))  # σ' = σ + d_μ d_ρ⁻¹ X^k σ^(ρ)
        nuevo = [mul(a, c) for c in sigma] + [0] * (len(rho[1]) + k)
        for i, c in enumerate(rho[1]):
            nuevo[i + k] ^= mul(b, c)
        while nuevo[-1] == 0:
            nuevo.pop()
        sigma, l = nuevo, max(l, rho[3] + k)
    return filas

S = sindromes([3, 5, 12])
print('S =', ', '.join(alfa(s) for s in S[1:]))
for nombre, si in (('con inversiones', False), ('sin inversiones', True)):
    print(nombre)
    for mu, sigma, d, l in berlekamp(S, sin_inversiones=si):
        print(f'  {mu:>4}  σ = {" ".join(alfa(c) for c in sigma):<14} d = {alfa(d) if d is not None else "—":<4} l = {l}')

# Todos los patrones de 1, 2 y 3 errores: ¿el σ sin inversiones es c·σ?
distintos = 0
total = 0
for w in (1, 2, 3):
    for pos in combinations(range(15), w):
        S = sindromes(pos)
        s1 = berlekamp(S)[-1][1]
        s2 = berlekamp(S, sin_inversiones=True)[-1][1]
        c = s2[0]
        total += 1
        distintos += [mul(x, inv(c)) for x in s2] != s1
print(total, 'patrones,', distintos, 'distintos')
# S = 1, 1, α^10, 1, α^10, α^5
# con inversiones
#   -0.5  σ = 1              d = 1    l = 0
#      0  σ = 1              d = 1    l = 0
#      1  σ = 1 1            d = α^5  l = 1
#      2  σ = 1 1 α^5        d = α^10 l = 2
#      3  σ = 1 1 0 α^5      d = —    l = 3
# sin inversiones
#   -0.5  σ = 1              d = 1    l = 0
#      0  σ = 1              d = 1    l = 0
#      1  σ = 1 1            d = α^5  l = 1
#      2  σ = 1 1 α^5        d = α^10 l = 2
#      3  σ = α^5 α^5 0 α^10 d = —    l = 3
# 575 patrones, 0 distintos

Los coeficientes van del grado 0 para arriba: 1 1 0 α^5 es $1+X+{\alpha }^5{X}^3$. La discrepancia se calcula con ${\sigma }_0$ adentro de la suma, que con inversiones vale 1 y sin inversiones no: así la misma línea sirve para las dos. Y la comparación del final divide por ${\sigma }_0$ sólo para comparar; el decodificador no lo necesita.

Para tener a mano

1. $\sigma (X)$ es el polinomio de conexión del registro más corto que genera $S_1,\dots ,S_{2t}$. Su grado es el número de errores.
2. Discrepancia: lo que falla el registro actual en la identidad siguiente. Si es cero, la fila se copia.
3. Si no, se corrige con la fila $\rho$ de $d_{\rho }\ne 0$ y mayor $\rho -l_{\rho }$: ${\sigma }^{(\mu +1)}={\sigma }^{(\mu )}+d_{\mu }{d}_{\rho }^{-1}{X}^{\mu -\rho }{\sigma }^{(\rho )}$, y $l_{\mu +1}=\max (l_{\mu },l_{\rho }+\mu -\rho )$.
4. El término agregado no toca las identidades viejas y cancela la nueva. El criterio de $\rho$ es lo que hace al registro el más corto.
5. Binario: como $S_{2i}=S_i^2$, los pasos impares tienen $d=0$. Quedan $t$ iteraciones, $\mu =-\frac{1}{2},0,\dots ,t$, con corrimiento $2(\mu -\rho )$ y criterio $2\rho -l_{\rho }$.
6. Tabla 6, errores en 3, 5 y 12: $d_0=1$, $d_1={\alpha }^5$, $d_2={\alpha }^{10}$; $\rho =-\frac{1}{2},0,1$; $\sigma =1+X+{\alpha }^5{X}^3$.
7. Sin inversiones: ${\sigma }^{(\mu +1)}=d_{\rho }{\sigma }^{(\mu )}+d_{\mu }{X}^{2(\mu -\rho )}{\sigma }^{(\rho )}$. El mismo $\sigma$ por una constante, las mismas raíces; ${\sigma }_0\ne 1$ y las discrepancias también escaladas.
8. Con más de $t$ errores entrega el registro más corto, que no es el localizador: con cuatro en el (15,5), 840 fallas declaradas y 525 palabras equivocadas.