Berlekamp–Massey Algorithm in Haskell

Question

I was recently tasked to write an implementation of the Berlekamp–Massey algorithm (for GF2) from a given inductional proof and decided to fresh up on my Haskell, since converting a proof by induction into a recursive algorithm should be quite straight-forward.

The variable names have been taken from the given proof (and thus can't be changed, although they differ from other notations or are rather short); here's a map:

• ks: the input bit stream
• len: the length of ks
• n: degree of the current annihilating polynomial
• phi: said polynomial, where a "bit" set at position \$i\$ indicates that the the \$i\$-th power is "part" of the polynomial, ie.: $$ \varphi(X) = \sum_{i=0}^n \text{phi}[i] \cdot X^i $$
• s, psi: analogously the degree and the polynomial itself of the last iteration that had different values
• r: the maximum length of the prefix of ks that psi annihilates successfully
• dl: the discrepancy between the annihilation by phi and the \$l\$-prefix of ks
• eta: The new polynomial, calculated as: $$ \eta(X) = \varphi(X) + X^{l-r} \cdot \psi(X) $$

data State = State { n :: Int
                   , phi :: [Int]
                   , s :: Int
                   , psi :: [Int]
                   , r :: Int
                   }
    deriving Show

discr :: Int -> [Int] -> Int -> [Int] -> Int
discr l ks n phi = foldr (\ a k -> (a+k) `mod` 2) 0
    (zipWith (*) (take (n+1) phi) (reverse $ take (n+1) $ drop (l-n) ks))

buildEta :: [Int] -> [Int] -> Int -> Int -> Int -> [Int]
buildEta phi psi l r s
    = (take (l-r) phi) ++
      (zipWith (\ c b -> (c+b) `mod` 2) (take (s+1) $ drop (l-r) phi) (take (s+1) psi)) ++
      (drop (l-r+s+1) phi)

bk' :: Int -> [Int] -> State
bk' (-1) ks = State 0 (1:lenzeroes) 0 (1:lenzeroes) (-1)
    where lenzeroes = replicate (length ks) 0
bk' l ks
    | dl == 0   = State n phi s psi r
    | n*2 <= l  = State (l+1-n) eta n phi l
    | otherwise = State n eta s psi r
    where State n phi s psi r = bk' (l-1) (ks)
          dl = discr l ks n phi
          eta = buildEta phi psi l r s

bk :: [Int] -> (Int,[Int])
bk ks = (n,phi)
    where State n phi _ _ _ = bk' (length ks) ks

I'm not quite happy with the result though, especially the buildEta function and the use of length. I think the latter could be improved by reversing ks and then giving only the tail to the recursive call, but somehow I still need to build the phi and psi lists correctly, which I fail to do "nicely".

Style aside, the code works (as far as I know), eg.:

*Main> bk [1, 0, 1, 0, 0, 1, 1, 0, 1]
(5,[1,0,1,0,1,1,0,0,0])

Which means that the degree of the polynomial is \$n\$ and it has the following form: $$ \varphi = X^0 + X^2 + X^4 + X^5 $$ Or, as an LSFR, A=[1,1,0,1,0] which indeed generates ks with [1, 0, 1, 0, 0] as initialization sequence.

---

Edit: Had someone review my code offline and did some minor style fixes already:

• some point-free style where appropriate
• replacing some lambda with just (*)
• replacing sum ()mod2 with foldr and lambda
• found a way to rewrite the if-then-else into nice guards
• replaced the list-comprehension for zeroes with replicate
• line breaks :-)

Answer 1

A central question I ask myself on these is whether the explicit recursion follows an abstractable pattern. In this case, only bk' has explicit recursion, and its pattern is simple - the call tree is always linear and brings l incrementally down to -1.

buildEta applies a transformation to a particular slice of phi.

Why are you working with Ints? These are Bools.

bk :: [Bool] -> (Int,[Bool])
bk ks = (n,phi) where
  State n phi _ _ _ = foldr (bk' ks) (State 0 lenzeroes 0 lenzeroes (-1)) [l,l-1..0]
  lenzeroes = True : replicate l False
  l = length ks

bk' :: [Bool] -> Int -> State -> State
bk' ks l (State n phi s psi r)
  | not dl    = State n phi s psi r
  | n*2 <= l  = State (l+1-n) eta n phi l
  | otherwise = State n eta s psi r
  where dl = discr l ks n phi
        eta = overSlice (l-r) (s+1) (zipWith xor psi) phi

discr :: Int -> [Bool] -> Int -> [Bool] -> Bool
discr l ks n = foldr xor False . zipWith (&&) (reverse $ take (n+1) $ drop (l-n) ks)

overSlice :: Int -> Int -> ([a] -> [a]) -> ([a] -> [a])
overSlice from size f a = b ++ f d ++ e where
  (b, c) = splitAt from a
  (d, e) = splitAt size c