Berlekamp–Massey algorithm

Berlekamp–Massey algorithm
You are encouraged to solve this task according to the task description, using any language you may know.

The Berlekamp–Massey algorithm is an efficient algorithm for finding the shortest linear feedback shift register (LFSR) that generates a given sequence over a finite field. It determines the minimal polynomial of a linearly recurrent sequence, which is crucial in applications like error-correcting codes, cryptography, and sequence analysis.

Overview

Given a sequence of elements ${\displaystyle s_{0},s_{1},\dots ,s_{n-1}}$ from a finite field, the Berlekamp–Massey algorithm computes the shortest LFSR, characterized by its connection polynomial ${\displaystyle C(x)=1+c_{1}x+c_{2}x^{2}+\dots +c_{L}x^{L}}$, such that the sequence satisfies the recurrence relation:

${\displaystyle s_{i+L}=-\sum _{k=1}^{L}c_{k}s_{i+L-k}}$

for all valid indices. The degree ${\displaystyle L}$ of the polynomial is the linear complexity of the sequence, representing the length of the shortest LFSR.

History

The algorithm was independently developed by Elwyn Berlekamp in 1968 for decoding BCH codes and by James Massey in 1969 for cryptographic applications. It remains a cornerstone in coding theory and sequence analysis.

Algorithm Description

The Berlekamp–Massey algorithm iteratively constructs the connection polynomial by processing the input sequence element by element. It maintains a current polynomial ${\displaystyle C(x)}$ and updates it whenever a discrepancy is found between the predicted and actual sequence values. The key steps are:

1. Initialize ${\displaystyle C(x)=1}$, degree ${\displaystyle L=0}$, and a temporary polynomial ${\displaystyle B(x)=1}$. 2. For each sequence element ${\displaystyle s_{i}}$:

  * Compute the discrepancy ${\displaystyle \Delta =s_{i}+\sum _{k=1}^{L}c_{k}s_{i-k}}$.
  * If ${\displaystyle \Delta \neq 0}$:
    * Update ${\displaystyle C(x)\leftarrow C(x)-\Delta xB(x)}$.
    * If the discrepancy indicates a need for a longer LFSR (i.e., ${\displaystyle 2L\leq i}$), update ${\displaystyle L}$ and adjust ${\displaystyle B(x)}$.

3. Output ${\displaystyle C(x)}$ and ${\displaystyle L}$.

The algorithm runs in ${\displaystyle O(n^{2})}$ time, where ${\displaystyle n}$ is the sequence length, and uses ${\displaystyle O(n)}$ space.

Applications

• Error Correction: Used in decoding Reed–Solomon codes and BCH codes.
• Cryptography: Analyzes the linear complexity of sequences in stream ciphers to assess security.
• Sequence Analysis: Identifies patterns in sequences for signal processing and data compression.

Example

For the binary sequence ${\displaystyle 0,0,1,1,0,1,0}$ over ${\displaystyle GF(2)}$, the algorithm yields the connection polynomial ${\displaystyle C(x)=1+x+x^{3}}$, indicating a linear complexity of 3.

ALGOL 68

Translation of: FreeBASIC

BEGIN # Berlekamp–Massey algorithm - translation of the FreeBASIC sample     #

    INT mdl = 2;

    PROC fp = ( INT a, k )INT:
         BEGIN
            INT result := 1, bbase := a MOD mdl, expo := k;
    
            WHILE expo > 0 DO
                IF ODD expo THEN result *:= bbase MODAB mdl FI;
                bbase *:= bbase MODAB mdl;
                expo OVERAB 2
            OD;
            result
        END # fp # ;

    PROC berlekamp massey = ( []INT a in )[]INT:
         BEGIN
            [ 0 : -1 ]INT empty row;
            REF[]INT ans coef := empty row, lst := empty row;
            []INT a = a in[ AT 0 ];
            INT   n = UPB a;
            INT   lst size := 0, w := 0, delta := 0;
    
            FOR i TO n DO
                INT tmp := 0;
                INT ans len = UPB ans coef;
                FOR j FROM 0 TO ans len DO
                    IF i - 1 - j >= 1 THEN
                        tmp +:= ( a[ i - 1 - j ] * ans coef[ j ] ) MODAB mdl
                    FI
                OD;
        
                INT  discrepancy = ( a[ i ] - tmp + mdl ) MOD mdl;
                IF   discrepancy = 0 
                THEN SKIP
                ELIF w = 0
                THEN
                    ans coef := HEAP[ 0 : i - 1 ]INT;
                    FOR j FROM 0 TO i - 1 DO ans coef[ j ] := 0 OD;
                    w     := i;
                    delta := discrepancy
                ELSE
                    [ 0 : UPB ans coef ]INT ahora := ans coef;

                    INT inv delta  = fp( delta, mdl - 2 );
                    INT mul        = ( discrepancy * inv delta ) MOD mdl;
                    INT needed len = lst size + i - w;
        
                    IF UPB ans coef + 1 < needed len THEN
                        ans coef := HEAP[ 0 : needed len - 1 ]INT;
                        ans coef[ 0 : UPB ahora AT 0 ] := ahora;
                        FOR j FROM UPB ahora + 1 TO UPB ans coef DO ans coef[ j ] := 0 OD
                    FI;

                    IF i - w - 1 >= 0 THEN
                        ans coef( i - w - 1 ) +:= mul MODAB mdl
                    FI;
        
                    FOR j FROM 0 TO lst size - 1 DO
                        INT idx = i - w + j;
                        IF idx < UPB ans coef  + 1 THEN
                            INT term to subtract = ( mul * lst[ j ] ) MOD mdl;
                            ans coef[ idx ] -:= term to subtract +:= mdl MODAB mdl
                        FI
                    OD;
        
                    IF UPB ans coef > UPB ahora THEN
                        lst := HEAP[ 0 : UPB ahora ]INT;
                        lst[ AT 0 ] := ahora;
                        lst size := UPB ahora + 1;
                        w := i;
                        delta := discrepancy
                    FI
                FI;
    
                FOR m FROM 0 TO UPB ans coef DO ans coef[ m ] +:= mdl MODAB mdl OD

            OD;
            ans coef
         END # berlekamp massey # ;

CO # uncomment if calculation of specific terms is required #
    PROC calculate term = ( INT m, []INT coef, h )INT:
         IF    m <= UPB h
         THEN ( h[ m ] + mdl ) MOD mdl
         ELIF INT k = UPB coef + 1;
              k = 0
         THEN 0
         ELSE [ 0 : k ]INT p coeffs;
              p coeffs[ 0     ] := ( mdl - 1 ) MOD mdl;
              p coeffs[ 1 : k ] := coef;
    
              INT degree k := k;
              [ 0 : degree k - 1 ]INT f, g; 
              f[ 0 ] := 1;
              FOR i TO UPB f DO f[ i ] := 0 OD;
              FOR i FROM 0 TO UPB g DO g[ i ] := 0 OD;
              IF k = 1 THEN
                  g[ 0 ] := p coeffs[ 1 ]
              ELSE
                  g[ 1 ] := 1
              FI;
    
              INT power := m;
              [ 0 : 2 * degree k - 1 ]INT tmp;
              INT term := 0;
    
              WHILE power > 0 DO
                 IF ODD power THEN
                    FOR i FROM 0 TO UPB tmp DO tmp[ i ] := 0 OD;
                    FOR i FROM 0 TO degree k - 1 DO
                        IF f[ i ] /= 0 THEN
                            FOR j FROM 0 TO UPB f DO
                                tmp[ i + j ] +:= ( f[ i ] * g[ j ] ) MODAB mdl
                            OD
                        FI
                    OD;
                    FOR i FROM UPB tmp BY - 1 TO UPB f + 1 DO
                        IF tmp[ i ] /= 0 THEN
                            term     := tmp[ i ];
                            tmp[ i ] := 0;
                            FOR j FROM 0 TO k DO
                                INT idx = i - j;
                                IF idx >= 0 THEN
                                    tmp[ idx ] +:= ( term * p coeffs[ j ] ) MODAB mdl
                                FI
                            OD
                        FI
                    OD;
                    f := tmp[ 0 : UPB f AT 0 ]
                 FI;
                 FOR i FROM 0 TO UPB tmp DO tmp[ i ] := 0 OD;
                 FOR i FROM 0 TO UPB g DO
                     IF g[ i ] /= 0 THEN
                         FOR j FROM 0 TO UPB g DO
                             tmp[ i + j ] +:= ( g[ i ] * g[ j ] ) MODAB mdl
                         OD
                     FI
                 OD;
                 FOR i FROM UPB tmp BY -1 TO degree k DO
                     IF tmp[ i ] /= 0 THEN
                         term := tmp[ i ];
                         tmp[ i ] := 0;
                         FOR j FROM 0 TO k DO
                             INT idx = i - j;
                             IF idx >= 0 THEN
                                 tmp[ idx ] +:= ( term * p coeffs[ j ] ) MODAB mdl
                             FI
                         OD
                     FI
                 OD;
                 g := tmp[ 0 : UPB g AT 0 ];
        
                 power OVERAB 2
              OD;
    
              INT final ans := 0;
              FOR i FROM 0 TO degree k - 1 DO
                  IF i + 1 <= UPB h THEN
                      final ans +:= ( h[ i + 1 ] * f[ i ] ) MODAB mdl
                  FI
              OD;
    
              ( final ans + mdl ) MOD mdl
         FI # calculate term # ;
CO

    PROC solve = VOID:
         BEGIN

            []INT h input = []INT( 0, 0, 1, 1, 0, 1, 0 )[ AT 0 ];
            INT n = UPB h input + 1;
            [ 0 : n ]INT h;
 
            h[ 0 ] := 0; FOR i FROM LWB h TO UPB h - 1 DO h[ i + 1 ] := h input[ i ] OD;
    
            []INT ans coef = berlekamp massey( h );
    
            print( ( "Coefficients:" ) );
            FOR i FROM 0 TO UPB ans coef DO
                print( ( " ", whole( ( ans coef[ i ] + mdl ) MOD mdl, 0 ) ) )
            OD;
            print( ( newline ) )
    
            CO # Uncomment to calculate a specific term #
               ;
               INT m      = 10;
               INT result = calculate term( m, ans coef(), h() );
               print( ( "Term ", whole( m, 0 ) ) );
               print( ( ": ", whole( ( result + mdl ) MOD mdl, 0 ), newline ) )
            CO

         END # solve # ;

    solve

END

Output:

Coefficients: 1 1 0 1

FreeBASIC

Translation of: Python

Const As Integer MDL = 2

Function fp(a As Integer, k As Integer) As Integer
    Dim As Integer result = 1
    Dim As Integer bbase = a Mod MDL
    Dim As Integer expo = k
    
    While expo > 0
        If (expo And 1) Then result = (result * bbase) Mod MDL
        bbase = (bbase * bbase) Mod MDL
        expo \= 2
    Wend
    
    Return result
End Function

Sub Berlekamp_Massey(a() As Integer, ans_coef() As Integer)
    Dim As Integer n = Ubound(a)
    Dim As Integer lst_size = 0
    Dim As Integer w = 0
    Dim As Integer delta = 0
    Dim As Integer i, j, tmp, idx
    
    Redim ans_coef(-1)
    Dim As Integer lst()
    Redim lst(-1)
    
    For i = 1 To n
        tmp = 0
        Dim As Integer ans_len = Ubound(ans_coef)
        For j = 0 To ans_len
            If i - 1 - j >= 1 Then
                tmp = (tmp + a(i - 1 - j) * ans_coef(j)) Mod MDL
            End If
        Next j
        
        Dim As Integer discrepancy = (a(i) - tmp + MDL) Mod MDL
        
        If discrepancy = 0 Then Continue For
        
        If w = 0 Then
            Redim ans_coef(i - 1)
            For j = 0 To i - 1
                ans_coef(j) = 0
            Next j
            
            w = i
            delta = discrepancy
            Continue For
        End If
        
        Dim As Integer Ahora(Ubound(ans_coef))
        For j = 0 To Ubound(ans_coef)
            Ahora(j) = ans_coef(j)
        Next j
        
        Dim As Integer inv_delta = fp(delta, MDL - 2)
        Dim As Integer mul = (discrepancy * inv_delta) Mod MDL
        Dim As Integer needed_len = lst_size + i - w
        
        If Ubound(ans_coef) + 1 < needed_len Then
            Redim Preserve ans_coef(needed_len - 1)
            For j = Ubound(Ahora) + 1 To Ubound(ans_coef)
                ans_coef(j) = 0
            Next j
        End If
        
        If i - w - 1 >= 0 Then
            ans_coef(i - w - 1) = (ans_coef(i - w - 1) + mul) Mod MDL
        End If
        
        For j = 0 To lst_size - 1
            idx = i - w + j
            If idx < Ubound(ans_coef) + 1 Then
                Dim As Integer term_to_subtract = (mul * lst(j)) Mod MDL
                ans_coef(idx) = (ans_coef(idx) - term_to_subtract + MDL) Mod MDL
            End If
        Next j
        
        If Ubound(ans_coef) + 1 > Ubound(Ahora) + 1 Then
            Redim lst(Ubound(Ahora))
            For j = 0 To Ubound(Ahora)
                lst(j) = Ahora(j)
            Next j
            lst_size = Ubound(Ahora) + 1
            
            w = i
            delta = discrepancy
        End If
    Next i
    
    For i = 0 To Ubound(ans_coef)
        ans_coef(i) = (ans_coef(i) + MDL) Mod MDL
    Next i
End Sub

Function calculate_term(m As Integer, coef() As Integer, h() As Integer) As Integer
    Dim As Integer i, j, k
    
    If m <= Ubound(h) Then Return (h(m) + MDL) Mod MDL
    k = Ubound(coef) + 1
    If k = 0 Then Return 0
    
    Dim As Integer p_coeffs(k)
    p_coeffs(0) = (MDL - 1) Mod MDL
    For i = 0 To k - 1
        p_coeffs(i + 1) = coef(i)
    Next i
    
    Dim As Integer degree_k = k
    Dim As Integer f(degree_k - 1), g(degree_k - 1)
    
    f(0) = 1
    For i = 1 To degree_k - 1
        f(i) = 0
    Next i
    
    If k = 1 Then
        g(0) = p_coeffs(1)
        For i = 1 To degree_k - 1
            g(i) = 0
        Next i
    Else
        g(0) = 0
        g(1) = 1
        For i = 2 To degree_k - 1
            g(i) = 0
        Next i
    End If
    
    Dim As Integer power = m
    Dim As Integer tmp(2 * degree_k - 1)
    Dim As Integer term, idx
    
    While power > 0
        If (power And 1) Then
            For i = 0 To 2 * degree_k - 1
                tmp(i) = 0
            Next i
            
            For i = 0 To degree_k - 1
                If f(i) = 0 Then Continue For
                For j = 0 To degree_k - 1
                    tmp(i + j) = (tmp(i + j) + f(i) * g(j)) Mod MDL
                Next j
            Next i
            
            For i = 2 * degree_k - 1 To degree_k Step -1
                If tmp(i) = 0 Then Continue For
                term = tmp(i)
                tmp(i) = 0
                For j = 0 To k
                    idx = i - j
                    If idx >= 0 Then
                        tmp(idx) = (tmp(idx) + term * p_coeffs(j)) Mod MDL
                    End If
                Next j
            Next i
            
            For i = 0 To degree_k - 1
                f(i) = tmp(i)
            Next i
        End If
        
        For i = 0 To 2 * degree_k - 1
            tmp(i) = 0
        Next i
        
        For i = 0 To degree_k - 1
            If g(i) = 0 Then Continue For
            For j = 0 To degree_k - 1
                tmp(i + j) = (tmp(i + j) + g(i) * g(j)) Mod MDL
            Next j
        Next i
        
        For i = 2 * degree_k - 1 To degree_k Step -1
            If tmp(i) = 0 Then Continue For
            term = tmp(i)
            tmp(i) = 0
            For j = 0 To k
                idx = i - j
                If idx >= 0 Then
                    tmp(idx) = (tmp(idx) + term * p_coeffs(j)) Mod MDL
                End If
            Next j
        Next i
        
        For i = 0 To degree_k - 1
            g(i) = tmp(i)
        Next i
        
        power \= 2
    Wend
    
    Dim As Integer final_ans = 0
    For i = 0 To degree_k - 1
        If i + 1 <= Ubound(h) Then
            final_ans = (final_ans + h(i + 1) * f(i)) Mod MDL
        End If
    Next i
    
    Return (final_ans + MDL) Mod MDL
End Function

Sub Solve()
    Dim As Integer h_input(6) = {0, 0, 1, 1, 0, 1 ,0}    
    Dim As Integer i, n
    n = Ubound(h_input) + 1
    
    Dim As Integer h(n)
    h(0) = 0
    For i = 0 To n - 1
        h(i + 1) = h_input(i)
    Next i
    
    Dim As Integer ans_coef()
    Berlekamp_Massey(h(), ans_coef())
    
    Print "Coefficients:";
    For i = 0 To Ubound(ans_coef)
        Print " " & ((ans_coef(i) + MDL) Mod MDL);
    Next i
    Print
    
    ' Uncomment to calculate a specific term
    ' Dim As Integer m = 10
    ' Dim As Integer result = calculate_term(m, ans_coef(), h())
    ' Print "Term " & m & ": " & ((result + MDL) Mod MDL)
End Sub

Solve()

Sleep

Output:

Coefficients: 1 1 0 1

Python

MOD = 2
def fp(a, k):
    return pow(a, k, MOD)
def berlekamp_massey(a):
    n = len(a) - 1  
    ans_coef = []  
    lst = []  
    w = 0  
    delta = 0  
    for i in range(1, n + 1):
        tmp = 0
        for j in range(len(ans_coef)):
            if i - 1 - j >= 1:
                tmp = (tmp + a[i - 1 - j] * ans_coef[j]) % MOD
        discrepancy = (a[i] - tmp + MOD) % MOD
        if discrepancy == 0:
            continue  
        if w == 0:
            ans_coef = [0] * i  
            w = i
            delta = discrepancy
            continue
        now = list(ans_coef)
        mul = discrepancy * fp(delta, MOD - 2) % MOD
        needed_len = len(lst) + i - w
        if len(ans_coef) < needed_len:
            ans_coef.extend([0] * (needed_len - len(ans_coef)))
        if i - w - 1 >= 0:
            ans_coef[i - w - 1] = (ans_coef[i - w - 1] + mul) % MOD
        for j in range(len(lst)):
            idx = i - w + j
            if idx < len(ans_coef):  
                term_to_subtract = (mul * lst[j]) % MOD
                ans_coef[idx] = (ans_coef[idx] - term_to_subtract + MOD) % MOD
        if len(ans_coef) > len(now):  
            lst = now
            w = i
            delta = discrepancy
    return [(x + MOD) % MOD for x in ans_coef]
def calculate_term(m, coef, h):
    k = len(coef)
    if m < len(h):
        return (h[m] + MOD) % MOD
    if k == 0:  
        return 0
    p_coeffs = [0] * (k + 1)
    p_coeffs[0] = (MOD - 1) % MOD  
    for i in range(k):
        p_coeffs[i + 1] = coef[i]
    def poly_mul(a, b, degree_k, p_poly):
        res = [0] * (2 * degree_k)
        for i in range(degree_k):
            if a[i] == 0: continue
            for j in range(degree_k):
                res[i + j] = (res[i + j] + a[i] * b[j]) % MOD
        for i in range(2 * degree_k - 1, degree_k - 1, -1):
            if res[i] == 0: continue
            term = res[i]
            res[i] = 0  
            for j in range(degree_k + 1):
                idx = i - j
                if idx >= 0:
                    res[idx] = (res[idx] + term * p_poly[j]) % MOD
        return res[:degree_k]
    f = [0] * k
    g = [0] * k
    f[0] = 1  
    if k == 1:
        if k == 1:
            g[0] = p_coeffs[1]  
        else:
            g[1] = 1  
    else:  
        g[1] = 1  
    power = m
    while power > 0:
        if power & 1:
            f = poly_mul(f, g, k, p_coeffs)
        g = poly_mul(g, g, k, p_coeffs)
        power >>= 1
    final_ans = 0
    for i in range(k):
        if i + 1 < len(h):
            final_ans = (final_ans + h[i + 1] * f[i]) % MOD
    return (final_ans + MOD) % MOD
def solve():
    h_input = [0,0,1,1,0,1,0]
    n=len(h_input)
    h = [0] + h_input
    ans_coef = berlekamp_massey(h)
    print(*((x + MOD) % MOD for x in ans_coef))
    # m=10
    # result = calculate_term(m, ans_coef, h)
    # print((result + MOD) % MOD)
if __name__ == "__main__":
    solve()

Output:

1 1 0 1

Wren

Translation of: Python

Library: Wren-math

Library: Wren-fmt

import "./math" for Int
import "./fmt" for Fmt

var MOD = 2

var fp = Fn.new { |a, k| Int.modPow(a, k, MOD) }

var berlekampMassey = Fn.new { |a|
    var n = a.count - 1
    var ansCoef = []
    var lst = []
    var w = 0
    var delta = 0
    for (i in 1...n + 1) {
        var tmp = 0
        for (j in 0...ansCoef.count) {
            if (i - 1 - j >= 1) tmp = (tmp + a[i - 1 - j] * ansCoef[j]) % MOD
        }
        var discrepancy = (a[i] - tmp + MOD) % MOD
        if (discrepancy == 0) continue 
        if (w == 0) {
            ansCoef = List.filled(i, 0)
            w = i
            delta = discrepancy
            continue
        }
        var now = ansCoef.toList
        var mul = discrepancy * fp.call(delta, MOD - 2) % MOD
        var neededLen = lst.count + i - w
        if (ansCoef.count < neededLen) {
            ansCoef.addAll(List.filled(neededLen - ansCoef.count, 0))
        }
        if (i - w - 1 >= 0) {
            ansCoef[i - w - 1] = (ansCoef[i - w - 1] + mul) % MOD
        }
        for (j in 0...lst.count) {
            var idx = i - w + j
            if (idx < ansCoef.count) {
                var termToSubtract = (mul * lst[j]) % MOD
                ansCoef[idx] = (ansCoef[idx] - termToSubtract + MOD) % MOD
            }
        }
        if (ansCoef.count > now.count) {
            lst = now
            w = i
            delta = discrepancy
        }
    }
    return ansCoef.map { |x| (x + MOD) % MOD }.toList
}

var calculateTerm = Fn.new { |m, coef, h|
    var k = coef.count
    if (m < h.count) return (h[m] + MOD) % MOD
    if (k == 0) return 0
    var pCoeffs = List.filled(k + 1, 0)
    pCoeffs[0] = (MOD - 1) % MOD
    for (i in 0...k) pCoeffs[i + 1] = coef[i]

    var polyMul = Fn.new { |a, b, degreeK, pPoly|
        var res = List.filled(2 * degreeK, 0)
        for (i in 0...degreeK) {
            if (a[i] == 0) continue
            for (j in 0...degreeK) res[i + j] = (res[i + j] + a[i] * b[j]) % MOD
        }
        for (i in (2 * degreeK - 1)...(degreeK - 1)) {
            if (res[i] == 0) continue
            var term = res[i]
            res[i] = 0
            for (j in 0..degreeK) {
                var idx = i - j
                if (idx >= 0) res[idx] = (res[idx] + term * pPoly[j]) % MOD
            }
        }
        return res[0...degreeK]
    }

    var f = List.filled(k, 0)
    var g = List.filled(k, 0)
    f[0] = 1
    if (k == 1) g[0] = pCoeffs[1] else g[1] = 1
    var power = m
    while (power > 0) {
        if (power % 2 == 1) f = polyMul.call(f, g, k, pCoeffs)
        g = polyMul.call(g, g, k, pCoeffs)
        power = power >> 1
    }
    var finalAns = 0
    for (i in 0...k) if (i + 1 < h.count) finalAns = (finalAns + h[i + 1] * f[i]) % MOD
    return (finalAns + MOD) % MOD
}

var solve = Fn.new {
    var hInput = [0, 0, 1, 1, 0, 1, 0]
    var n = hInput.count
    var h = [0] + hInput
    var ansCoef = berlekampMassey.call(h)
    var coeffs = ansCoef.map { |x| (x + MOD) % MOD }.toList
    System.print("Coefficients: %(coeffs) (lowest to highest degree)")
    System.write("Hence connection polynomial: ")
    Fmt.pprint("$d", coeffs[-1..0], "", "x")
    System.print("and linear complexity: %(coeffs.count - 1)")
    /* Uncomment the following 3 lines to calculate and show a term. */ 
    // var m = 10
    // var result = calculateTerm.call(m, ansCoef, h)
    // System.print((result + MOD) % MOD)
}

solve.call()

Output:

Coefficients: [1, 1, 0, 1] (lowest to highest degree)
Hence connection polynomial: x³ + x + 1
and linear complexity: 3