N := 4; s := Array(0..3,[1,1,0,0]); C := 1; B := 1; L := 0; m := 1; b := 1; p := 2; for n from 0 to N-1 do d := s[n] + add( coeff(C,x,i)*s[n-i], i=1..L ) mod p; if d=0 then m := m+1; elif 2*L<=n then T := C; C := C-modp(d/b,p)*expand( x^m*B ) mod p; L := n+1-L; B := T; b := d; m := 1; else C := C-modp(d/b,p)*expand( x^m*B ) mod p; m := m+1; fi; od; L; C; B; # Check for n from L to N-1 do s[n]+add(coeff(C,x,i)*s[n-i],i=1..L) mod p; od; # Maple code for the Berlekamp-Massey algorithm # Adapted from www.cs.wisc.edu/~cs435-1/bermas.m # Transliteration of # Massey, "Shift-Register Synthesis and BCH Decoding," # IEEE Trans. Inform. Theory, 15(1):122-127, 1969. # Input: P, either 0 or a prime # If P>0 then we work over the field K = Z/Z[P] (mod P) # else we work over the field K = Q (rationals) # N, a positive integer # s, a list of >= 2*N terms in K # x, a formal variable # Returns: Unique monic annihilator of minimum degree, over K[x]. BM := proc(s, N, P, x) local C,B,T,L,k,i,n,d,b,safemod; ASSERT(nops(s) = 2*N); safemod := (exp, P) -> `if`(P=0, exp, exp mod P); B := 1; C := 1; L := 0; k := 1; b := 1; for n from 0 to 2*N-1 do d := s[n+1]; for i from 1 to L do d := safemod(d + coeff(C,x^i)*s[n-i+1], P); od; if d=0 then k := k+1 fi; if (d <> 0 and 2*L > n) then C := safemod(expand(C - d*x^k*B/b), P); k := k+1; fi; if (d <> 0 and 2*L <= n) then T := C; C := safemod(expand(C - d*x^k*B/b), P); B := T; L := n+1-L; k := 1; b := d; fi; od; return C; end: # Compute lambda(z) usign the Euclidean algorithm BMEA := proc(a,p,x) local m,R0,R1,V0,V1,i; n := iquo( nops(a), 2 ); m := 2*n-1; R0 := x^(2*n); R1 := add( a[m+1-i]*x^i, i=0..m ); #lprint("R0=",R0); #lprint("R1=",R1); V0 := 0; V1 := 1; while n <= degree(R1,x) do R0,R1 := R1,Rem(R0,R1,x,'Q') mod p; #lprint("R0=",R0); #lprint("R1=",R1); V0,V1 := V1,Expand(V0-Q*V1) mod p; #lprint("V0=",V0); #lprint("V1=",V1); od; i := 1/lcoeff(V1,x) mod p; i*V1 mod p; end: P := 103: d := 4: num := 21+83*x+90*x^2+4*x^3: # degree < d den := 1+11*x+23*x^2+58*x^3+69*x^4: # monic, degree <= d f := series(num/den, x=0, 4*d) mod P: s := [seq(coeff(f, x, i), i=0..13)]; g := BM(s, d, P, x); h := BMEA(s,P,x); Factor(h) mod P; with(numtheory): p := 997; f := randpoly(x,degree=100,terms=5) mod p; f := unapply(f,x); alpha := primroot(p); a := [seq( alpha^i mod p, i=0..19 )]; s := [seq( f(alpha^i mod p) mod p, i=0..19 )]; g := BM(s, nops(s)/2, p, x); Lambda := BMEA(s,p,x); R := Roots(Lambda) mod p; R := [seq( r[1], r=R )]; L := [seq( mlog(r,alpha,p), r=R )]; F := unapply( add( c[i]*x^L[i], i=1..nops(R) ), x ); n := nops(R); eqns := {seq( F(alpha^i mod p) mod p = s[i+1], i=0..n-1)}; sols := msolve(eqns,p); eval( F(x), sols ); f(x) mod p; roots=R; k := [seq( alpha^L[i] mod p, i=1..4 )]: 'k'=k; M := mul(x-k[i],i=1..n): M := Expand(M) mod p: 'M'=M; 'L'=Lambda; P := Array(1..n): X := Vector(n): U := Vector(n); for j to n do q := Quo(M,x-k[j],x) mod p; print(degree(q)); U[j] := 1/eval(q,x=k[j]) mod p; P[j] := q; X[j] := U[j]*add( s[i]*coeff(P[j],x,i-1), i=1..n ) mod p; od: 'X' = X; F(x); f(x); C := proc(i,j) coeff(P[i],x,j-1) end: M := Matrix(n,n,C); w := Vector(s[1..4]); c := M.w mod p: for i to n do c[i] := c[i]*U[i] mod p od: c=X; #alpha := 7; #p := 997; #seq( alpha^i mod p, i=0..9 ); #f := (x-alpha^3)*(x-alpha^4) mod p; #a := [seq( eval(f,x=alpha^i mod p) mod p, i=0..19 )]; #for i from 2 to 10 by 2 do Lambda := BMEA(s[1..i],p,x); Roots(Lambda) mod p; od; with(numtheory): f := x+x^10; p := 11; alpha := 2; a := [seq( eval(f,x=alpha^i mod p), i=0..5 )]; Lambda := BMEA(a,p,x); R := Roots(Lambda) mod p; seq( mlog(r[1],alpha,p), r=R );