Berlekamp-Massey minimal LFSR issues

Question

I am having some issues getting the correct LFSR for my sequence(pattern), when I implement it as LFSR and the corresponding taps, it doesn't generate the sequence, any suggestions? The goal patt is {1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1};

My code following wikipedia's version for binary field (https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm):

#include <stdio.h>

int main()
{
    int patt[]={1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1};
    int n=sizeof(patt)/sizeof(int);

    int N=0, L=0, m=-1, b[n], c[n], d=0, t[n], j;
    b[0]=1;
    c[0]=1;
    float val;

    for(int i=1; i<n; i++){
        b[i]=0;
        c[i]=0;
        //printf("b[%d]=%d, c[%d]=%d; ",i,b[i],i,c[i]);
    }


    while (N < n){

        printf("N=%d, ",N);
        d=c[0]*patt[N];//initializing the value of d

        for(int i=1; i<=L; i++){
            //printf("d = %d + %d*%d, ",d,c[i],patt[N-L]);
            d=d ^ c[i]*patt[N-L];
            //printf("d=%d \n",d);
        }

        printf("d=%d\n", d);

        if (d==0){
            printf("c=c\n\n");
        }
        else{
            for(int i=0; i<n; i++){
                t[i]=c[i];
            }

            j=0;
            while(N-m+j<=n-1){
                printf("c[%d-%d+%d]=c[%d-%d+%d]^b[%d]; c[%d]=c[%d]^b[%d], %d=%d^%d; ", N, m, j, N, m, j, j, N-m+j, N-m+j, j, c[N-m+j], c[N-m+j], b[j]);
                c[N-m+j]=c[N-m+j]^b[j];//XOR operator: ^
                printf("c=%d\n",c[N-m+j]);
                j++;
            }
            printf("\n");
            val=N;
            val=val/2;
            printf("L=%d, N=%d, N/2=%f \n",L, N, val);

            if(L<= val){
                printf("updating L, m & b\n\n");
                L=N+1-L;
                m=N;

                for(int i=0; i<n; i++){
                    b[i]=t[i];
                }
            }
        }
        N++;
    } 

    int CiSi=c[L]*patt[0];;

    for(int i=1; i<L; i++){
        CiSi=CiSi ^ c[L-i]*patt[i];//XORing
    }

    printf("CiSi = %d;", CiSi);

    printf("c=");

    for(int i=0; i<n; i++){
        printf("%d ",c[i]);
    }

    return 0;
}

The answers per cycle:

N=0, d=1
c[0--1+0]=c[0--1+0]^b[0]; c[1]=c[1]^b[0], 0=0^1; c=1
c[0--1+1]=c[0--1+1]^b[1]; c[2]=c[2]^b[1], 0=0^0; c=0
c[0--1+2]=c[0--1+2]^b[2]; c[3]=c[3]^b[2], 0=0^0; c=0
c[0--1+3]=c[0--1+3]^b[3]; c[4]=c[4]^b[3], 0=0^0; c=0
c[0--1+4]=c[0--1+4]^b[4]; c[5]=c[5]^b[4], 0=0^0; c=0
c[0--1+5]=c[0--1+5]^b[5]; c[6]=c[6]^b[5], 0=0^0; c=0
c[0--1+6]=c[0--1+6]^b[6]; c[7]=c[7]^b[6], 0=0^0; c=0
c[0--1+7]=c[0--1+7]^b[7]; c[8]=c[8]^b[7], 0=0^0; c=0
c[0--1+8]=c[0--1+8]^b[8]; c[9]=c[9]^b[8], 0=0^0; c=0
c[0--1+9]=c[0--1+9]^b[9]; c[10]=c[10]^b[9], 0=0^0; c=0
c[0--1+10]=c[0--1+10]^b[10]; c[11]=c[11]^b[10], 0=0^0; c=0
c[0--1+11]=c[0--1+11]^b[11]; c[12]=c[12]^b[11], 0=0^0; c=0

L=0, N=0, N/2=0.000000 
updating L, m & b

N=1, d=0
c=c

N=2, d=1
c[2-0+0]=c[2-0+0]^b[0]; c[2]=c[2]^b[0], 0=0^1; c=1
c[2-0+1]=c[2-0+1]^b[1]; c[3]=c[3]^b[1], 0=0^0; c=0
c[2-0+2]=c[2-0+2]^b[2]; c[4]=c[4]^b[2], 0=0^0; c=0
c[2-0+3]=c[2-0+3]^b[3]; c[5]=c[5]^b[3], 0=0^0; c=0
c[2-0+4]=c[2-0+4]^b[4]; c[6]=c[6]^b[4], 0=0^0; c=0
c[2-0+5]=c[2-0+5]^b[5]; c[7]=c[7]^b[5], 0=0^0; c=0
c[2-0+6]=c[2-0+6]^b[6]; c[8]=c[8]^b[6], 0=0^0; c=0
c[2-0+7]=c[2-0+7]^b[7]; c[9]=c[9]^b[7], 0=0^0; c=0
c[2-0+8]=c[2-0+8]^b[8]; c[10]=c[10]^b[8], 0=0^0; c=0
c[2-0+9]=c[2-0+9]^b[9]; c[11]=c[11]^b[9], 0=0^0; c=0
c[2-0+10]=c[2-0+10]^b[10]; c[12]=c[12]^b[10], 0=0^0; c=0

L=1, N=2, N/2=1.000000 
updating L, m & b

N=3, d=0
c=c

N=4, d=0
c=c

N=5, d=0
c=c

N=6, d=1
c[6-2+0]=c[6-2+0]^b[0]; c[4]=c[4]^b[0], 0=0^1; c=1
c[6-2+1]=c[6-2+1]^b[1]; c[5]=c[5]^b[1], 0=0^1; c=1
c[6-2+2]=c[6-2+2]^b[2]; c[6]=c[6]^b[2], 0=0^0; c=0
c[6-2+3]=c[6-2+3]^b[3]; c[7]=c[7]^b[3], 0=0^0; c=0
c[6-2+4]=c[6-2+4]^b[4]; c[8]=c[8]^b[4], 0=0^0; c=0
c[6-2+5]=c[6-2+5]^b[5]; c[9]=c[9]^b[5], 0=0^0; c=0
c[6-2+6]=c[6-2+6]^b[6]; c[10]=c[10]^b[6], 0=0^0; c=0
c[6-2+7]=c[6-2+7]^b[7]; c[11]=c[11]^b[7], 0=0^0; c=0
c[6-2+8]=c[6-2+8]^b[8]; c[12]=c[12]^b[8], 0=0^0; c=0

L=2, N=6, N/2=3.000000 
updating L, m & b

N=7, d=0
c=c

N=8, d=1
c[8-6+0]=c[8-6+0]^b[0]; c[2]=c[2]^b[0], 1=1^1; c=0
c[8-6+1]=c[8-6+1]^b[1]; c[3]=c[3]^b[1], 0=0^1; c=1
c[8-6+2]=c[8-6+2]^b[2]; c[4]=c[4]^b[2], 1=1^1; c=0
c[8-6+3]=c[8-6+3]^b[3]; c[5]=c[5]^b[3], 1=1^0; c=1
c[8-6+4]=c[8-6+4]^b[4]; c[6]=c[6]^b[4], 0=0^0; c=0
c[8-6+5]=c[8-6+5]^b[5]; c[7]=c[7]^b[5], 0=0^0; c=0
c[8-6+6]=c[8-6+6]^b[6]; c[8]=c[8]^b[6], 0=0^0; c=0
c[8-6+7]=c[8-6+7]^b[7]; c[9]=c[9]^b[7], 0=0^0; c=0
c[8-6+8]=c[8-6+8]^b[8]; c[10]=c[10]^b[8], 0=0^0; c=0
c[8-6+9]=c[8-6+9]^b[9]; c[11]=c[11]^b[9], 0=0^0; c=0
c[8-6+10]=c[8-6+10]^b[10]; c[12]=c[12]^b[10], 0=0^0; c=0

L=5, N=8, N/2=4.000000 
N=9, d=0
c=c

N=10, d=0
c=c

N=11, d=0
c=c

N=12, d=1
c[12-6+0]=c[12-6+0]^b[0]; c[6]=c[6]^b[0], 0=0^1; c=1
c[12-6+1]=c[12-6+1]^b[1]; c[7]=c[7]^b[1], 0=0^1; c=1
c[12-6+2]=c[12-6+2]^b[2]; c[8]=c[8]^b[2], 0=0^1; c=1
c[12-6+3]=c[12-6+3]^b[3]; c[9]=c[9]^b[3], 0=0^0; c=0
c[12-6+4]=c[12-6+4]^b[4]; c[10]=c[10]^b[4], 0=0^0; c=0
c[12-6+5]=c[12-6+5]^b[5]; c[11]=c[11]^b[5], 0=0^0; c=0
c[12-6+6]=c[12-6+6]^b[6]; c[12]=c[12]^b[6], 0=0^0; c=0

L=5, N=12, N/2=6.000000 
updating L, m & b

CiSi = 0;

CiSi = 0; Testing the values mentioned as a result of the algorithm looks correct because is equal to zero, but

c=1 1 0 1 0 1 1 1 1; excluding the last 4 zeros due to their values as zeros

c=1 1 0 1 0 1 1 1 1, this are the coefficients of the polynomial starting from left to right: C0,..., Ck: 1 + x^2 + x^4 + x^5 + x^6 + x^7

When I implement this values the results are NOT RIGHT

Implementation of the LFSR with the taps in the corresponding positions, x0 is excluded according to https://en.wikipedia.org/wiki/Linear-feedback_shift_register and also Wang Laung-Terng, Wu Cheng-Wen and W. Xiaoqing, "VLSI Test Principles and Architectures - Design for Testability", 2006:

%Matlab source code
clear all;
seed=[1 1 0 0 0 0 1 0];
seed_sz=size(seed);
%Loop to initialize a array
for i=1:50
    A{i}=1:seed_sz(1,2);
    A{i}(1,1:end)=0;
end

filename='LFSR rightshift no x0 c program.xlsx';

for i=1:50
    A{i}=seed;
    xlswrite(filename,A{i},'1',['A',int2str(i)]);
    XOR_output=xor(seed(1,8),seed(1,7));
    XOR_output=xor(XOR_output,seed(1,6));
    XOR_output=xor(XOR_output,seed(1,5));
    XOR_output=xor(XOR_output,seed(1,3));
    XOR_output=xor(XOR_output,seed(1,1));

    %Right shift the seed
    seed=circshift(seed,1);

    seed(1,1)=XOR_output;
end

Flow chart of the algorithm adapted from wikipedia

I added a supplemental answer to fgrieu's answer to show the code for Fibonacci and Galois LFSRs, both right shifting and left shifting examples. All 4 output the string pattern, plus two more bits, 0 1, which completes the repeating 15 bit pattern. — rcgldr, May 26 '18 at 01:42

fgrieu · Accepted Answer · 2018-05-25T11:41:04.227

Compared to the Wikipedia pseudocode that it aims at implementing, the question's code has two clear discrepancies (at least the second is a fatal bug):

d=c[0]*patt[N] should be d=patt[N] to match d ← s_N…
d=d ^ c[i]*patt[N-L] should be d=d ^ c[i]*patt[N-i] to match the other terms of d ← s_N ⊕ c₁s_N−1 ⊕ c₂s_N−2 ⊕ … ⊕ c_Ls_N−L

As an aside the code does not show the final L, that is the length of the minimal LFSR for the stream. But that L is also the index of the last 1 in the output, so we can get away with that omission.

With these changes the code outputs c=1 0 1 0 1 1 0 0 0 0 0 0 0, that is an LFSR with L=5 bits and the recurrence s_i ← s_i-2 ⊕ s_i-4 ⊕ s_i-5, or equivalently s_i+5 ← s_i+3 ⊕ s_i+1 ⊕ s_i. That indeed matches the sequence! Applied to the 5 first given terms it computes the next 8 ones:

s[ 0] :=                                     1
s[ 1] :=                                     1
s[ 2] :=                                     0
s[ 3] :=                                     0
s[ 4] :=                                     0
s[ 5] := s[ 3] ^ s[ 1] ^ s[ 0] = 0 ^ 1 ^ 1 = 0
s[ 6] := s[ 4] ^ s[ 2] ^ s[ 1] = 0 ^ 0 ^ 1 = 1
s[ 7] := s[ 5] ^ s[ 3] ^ s[ 2] = 0 ^ 0 ^ 0 = 0
s[ 8] := s[ 6] ^ s[ 4] ^ s[ 3] = 1 ^ 0 ^ 0 = 1
s[ 9] := s[ 7] ^ s[ 5] ^ s[ 4] = 0 ^ 0 ^ 0 = 0
s[10] := s[ 8] ^ s[ 6] ^ s[ 5] = 1 ^ 1 ^ 0 = 0
s[11] := s[ 9] ^ s[ 7] ^ s[ 6] = 0 ^ 0 ^ 1 = 1
s[12] := s[10] ^ s[ 8] ^ s[ 7] = 0 ^ 1 ^ 0 = 1

Interpretation of the program's output:

the rightmost 1 in the output tells the width of the LFSR, and the rest should be ignored;
the 1 digits in the output, except the leftmost, correspond to the terms to XOR in a Fibonnaci LFSR, from most recent to oldest;
the leftmost output is 1 and corresponds to the next term computed;
the 1 digits in reading order correspond to the terms of the Fibonacci polynomial from 1 to x^L (here, 1+x²+x⁴+x⁵), or equivalently to the terms of the Galois polynomial from x^L to 1 (here, x⁵+x³+x¹+1)

The initial state in Fibonacci convention is simply the first L terms of the given s_i, that is patt[i] in the question's code.

Here is a streamlined version of the code with less and clearer output, perusing for to clarify the intent, sticking to variable names in the original pseudocode, compatible with more C compilers, _{^{using Boolean operators when possible, staying away from floating point, and making loops with minimalist end conditions}}. It seems to work just fine in the few tests I made.

// Berlekamp-Massey algorithm per https://en.wikipedia.org/w/index.php?title=Berlekamp%E2%80%93Massey_algorithm&oldid=808089047#The_algorithm_for_the_binary_field
#include <stdio.h>
int main(void) {
    int s[]={1,1,0,0,0,0,1,0,1,0,0,1,1}; // bits of the stream to analyse
    #define n (sizeof(s)/sizeof(*s))     // how many bits there are
    int b[n], c[n], t[n], d, j, N, L=0, m=-1;
    for(j=n; --j>0;)
        b[j]=c[j]=0;
    b[0]=c[0]=1;
    for(N=0; N<n; ++N) {                 // For N=0 step 1 while N<n
        d=s[N];                          // first term  of discrepancy
        for(j=L; j>0; --j)               // other terms of discrepancy
            d ^= c[j]&s[N-j];
        if (d!=0) {                      // non-zero discrepancy
            for(j=n; --j>=0;)            // copy c to t
                t[j]=c[j];
            for(j=n-N+m; --j>=0;)        // XOR b (reversed) into c
                c[N-m+j] ^= b[j];
            if(L+L<=N) {                 // if L<=N/2
                L=N+1-L;
                m=N;
                for(j=n; --j>=0;)        // copy t to b
                    b[j]=t[j];
            }
        }
    }
    printf("s =");                       // show input
    for(j=0; j<n; j++)
        printf(" %d",s[j]);
    printf("\nc ="); 
    for(j=0; j<=L; j++)                  // show result
        printf(" %d",c[j]);
    printf("\nL = %d\n",L);              // show degree of polynomial
    return 0;
}

//  The above code outputs the following:
//  s = 1 1 0 0 0 0 1 0 1 0 0 1 1
//  c = 1 0 1 0 1 1
//  L = 5

I have a critic on the Wikipedia pseudocode and its transcription to code: indexes are going in reverse directions in the sequence under analysis s_i and the polynomial being constructed c_i; that seems to only create complications.

In the case of the implementation of the output polynomial, I am doing test as Fibonacci LFSR, but each of the cycles are not producing the expected terms, is there a trick for this implementation or what's the correct way? The LFSR diagram is:[diagram](https://drive.google.com/open?id=1K72mb6RhCOr5F4o5ZQzPDrcOcFt-IxRZ). The results per cycle: [cycles](https://drive.google.com/open?id=1xNNLV4mmbZH2w8e7w2e8pJsWhf9CqLyf). — Leo, May 25 '18 at 10:37
@Leonel Hernández: the diagram of the LFSR is OK, but you initialized it wrong. The most recent output is on the left (as stated in the answer), therefore the initial state you want is 00011 in reading order, obtained by reversing the first 5 terms in reading/chronological order. I added how to apply the recurrence. — fgrieu, May 25 '18 at 11:46
Hi, I am calculating the big O space complexity of the Berley-kamp, so far I've read couple of papers demonstrating that the time complexity is O(n^2). Would you say is the same for the space? — Leo, Jun 26 '19 at 16:19
@Leonel Hernández: hint: the above code works for large `n` (millions on a modern PC). The memory that it uses is declared in lines 4 and 6, starting with `int`. — fgrieu, Jun 26 '19 at 17:27
well, in that case would be 4*n right? In the big O complexity usually the constants (in this case 4) are diminished, so at the end is only O(n), right? — Leo, Jun 27 '19 at 09:16

rcgldr · Answer 2 · 2018-05-26T13:42:53.597

Supplement to fgrieu's answer. Code for Fibonacci and Galois LFSR, showing both left and right shift algorithms. These will produce the original string plus 2 more bits (0 1), which is a 15 bit repeating pattern. The starting LFSR values are different in order for the output pattern to match the original string pattern.

    // Fibonacci LFSR right shift
    printf("fr=");
    d = 0x03;
    do {
        printf(" %1x", d & 1);
        m = ((d>>3)^(d>>1)^(d>>0))&1;
        d  = (d >> 1)|(m << 4);
    } while (d != 0x3);
    printf("\n");

    // Fibonacci LFSR left shift
    printf("fl=");
    d = 0x18;
    do {
        printf(" %1x", d >> 4);
        m = ((d>>4)^(d>>3)^(d>>1)) & 1;
        d = ((d << 1)&0x1e) | m;
    } while (d != 0x18);
    printf("\n");

    //  Galios LFSR right shift
    printf("gr=");
    d = 0x1f;
    do {
        printf(" %1x", d & 1);
        m = d & 1;
        d >>= 1;
        if (m)
            d ^= 0x1a;
    } while (d != 0x1f);
    printf("\n");

    //  Galios LFSR left shift
    printf("gl=");
    d = 0x1f;
    do {
        printf(" %1x", d >> 4);
        d <<= 1;
        if (d & 0x20)
            d ^= 0x2b;
    } while (d != 0x1f);
    printf("\n");

output:

fr= 1 1 0 0 0 0 1 0 1 0 0 1 1 0 1
fl= 1 1 0 0 0 0 1 0 1 0 0 1 1 0 1
gr= 1 1 0 0 0 0 1 0 1 0 0 1 1 0 1
gl= 1 1 0 0 0 0 1 0 1 0 0 1 1 0 1

Berlekamp-Massey minimal LFSR issues

2 Answers2