(*******************************************************)
(* The BERLEKAMP-MASSEY algorithm in Mathematica     *)
(*                                                    *)
(* Author: Klaus Pommerening                          *)
(*        Johannes-Gutenberg-Universität              *)
(*        Mainz, Germany                              *)
(* 3. Januar 2001                                     *)
(* last change: 9. August 2003                        *)(*******************************************************)
(* Constructs the shortest linear shift register that  *)
(*   generates a given bit sequence.                   *)
(*******************************************************)

(*******************************************************)
(* Global variables                                    *)
(*******************************************************)

(* Insert your input sequence here as u.               *)
(*   u should be a sequence of bits 0 or 1; in any     *)
(*   case inside the procedure u is reduced mod 2.     *)
(* Here is a typical example:                          *)
u = {0,0,1,1,0,1,1,1,0}

(* Output 1: The linear profile of u will be put into *)
(*   the global list linprof.                          *)
linprof = {0}

(* Output 2: The feedback polynomial of the            *)
(*   constructed linear shift register will be put     *)
(*   into the global variable phi, a polynomial in T.  *)
(*   For phi = 1 + a_1 T + ... + a_l T^l the LSR       *)
(*   recurrence is given by                             *)
(*      u_i = a_1 u_{i-1} + ... + a_l u_{i-l}.         *)
phi = 1

(*******************************************************)
(* The BM algorithm                                    *)
(* Preconditions:                                      *)
(*   The input list u is filled.                       *)
(*   The output list linprof is initialized as {0}.    *)
(*   The output polynomial phi is initialized as 1.    *)
(* The return value is the linear complexity of u      *)
(*   that is also output as the last entry of linprof. *)
(*******************************************************)

bma :=
Module[{        (* Define and initialize local variables: *)
psi = 1 ,                                  (* last feedback polynomial *)
eta = 0,                                     (* new feedback polynomial *)
l = 0,            (* linear complexity up to actual index *)
t = 0,                       (* linear complexity in last state *)
m = 0,                                              (* new linear complexity *)
r = -1 ,                                                                     (* last index *)
d = 0,  (* discrepancy between predicted bit and    *)
                                   (* true bit -- always 0 or 1 in F_2 *)
nn, n, b, feedbacklist, i        (* auxiliary variables *)
},
nn = Length[u];
For [n = 0,  n <  nn, n++,
b = Take[u,{n-l+1,n}];        (* coeff. for feedback *)
feedbacklist = {};                                       (* feenback taps *)
For[i=1, i ≤ l, i++,
feedbacklist = Prepend[feedbacklist, Coefficient[phi,T,i]]];
d = Mod[u[[n+1]] - feedbacklist . b, 2];
If[d == 1,
eta = PolynomialMod[phi - T^(n-r) psi, 2];
If[2*l ≤ n,
m = n+1-l;
t = l;
l = m;
psi = phi;
r = n
]; (* END if 2*l ≤ n -----------------------------*)
phi = eta
]; (* END if d == 1 --------------------------------*)
linprof = Append[linprof, l]
]; (* END for n ... ---------------------------------*)
Return[l]
] (* END Module --------------------------------------*)

(*---------- Example application ----------------------*)
ls = bma

linprof

phi

ListPlot[Rest[linprof], AxesOrigin -> {0,0}]