(*******************************************************)
(* 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 q 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}]