Berlekamp Massey

src/lfsr.rs

@@ -86,6 +86,107 @@ impl Lfsr {
         assert_eq!(new_state.len(), self.state.len(), "State length mismatch");
         self.state = new_state;
     }
+
+    /// Berlekamp–Massey over GF(2^8).
+    ///
+    /// Returns the retroaction polynomial coefficients `c` such that:
+    ///   s[n] = c[1]*s[n-1] + c[2]*s[n-2] + ... + c[L]*s[n-L]
+    ///
+    /// The coefficients live in GF(2^8)
+    ///
+    /// The returned polynomial always has `c[0] = 1`.
+    pub fn berlekamp_massey(sequence: &[u8]) -> Vec<u8> {
+        assert!(!sequence.is_empty(), "sequence must not be empty");
+
+        // C(x): current retroaction polynomial
+        let mut c = vec![1u8];
+        // B(x): copy of the last "best" polynomial
+        let mut b = vec![1u8];
+
+        let mut l: usize = 0; // current linear complexity
+        let mut m: usize = 1; // number of steps since last update
+        let mut bb: u8 = 1; // last non-zero discrepancy
+
+        for n in 0..sequence.len() {
+            // discrepancy d = s[n] + sum_{i=1..L} c[i] * s[n-i]
+            let mut d = sequence[n];
+            for i in 1..=l {
+                d ^= gf_mul(c[i], sequence[n - i]);
+            }
+
+            if d == 0 {
+                m += 1;
+                continue;
+            }
+
+            let coef = gf_mul(d, gf_inv(bb));
+            let t = c.clone();
+
+            if c.len() < b.len() + m {
+                c.resize(b.len() + m, 0);
+            }
+
+            for i in 0..b.len() {
+                c[i + m] ^= gf_mul(coef, b[i]);
+            }
+
+            if 2 * l <= n {
+                l = n + 1 - l;
+                b = t;
+                bb = d;
+                m = 1;
+            } else {
+                m += 1;
+            }
+        }
+
+        c.truncate(l + 1);
+        c
+    }
+
+    pub fn predict_next_from_poly(history: &[u8], poly: &[u8]) -> u8 {
+        assert!(poly.len() >= 2, "polynomial must have degree at least 1");
+        assert_eq!(
+            history.len(),
+            poly.len() - 1,
+            "history must contain exactly degree-many previous samples"
+        );
+
+        let mut next = 0u8;
+        let degree = poly.len() - 1;
+
+        for i in 1..=degree {
+            next ^= gf_mul(poly[i], history[degree - i]);
+        }
+
+        next
+    }
+}
+
+fn gf_mul(mut a: u8, mut b: u8) -> u8 {
+    let mut p = 0u8;
+    for _ in 0..8 {
+        if (b & 1) != 0 {
+            p ^= a;
+        }
+        let hi = a & 0x80; // hi for high bit
+        a <<= 1;
+        if hi != 0 {
+            a ^= 0x1B; // reduction by 0x11B
+        }
+        b >>= 1;
+    }
+    p
+}
+
+fn gf_inv(x: u8) -> u8 {
+    assert!(x != 0, "cannot invert zero in GF(2^8)");
+    for y in 1u16..=255 {
+        if gf_mul(x, y as u8) == 1 {
+            return y as u8;
+        }
+    }
+    unreachable!("every non-zero element in GF(2^8) has an inverse");
 }
 
 #[cfg(test)]
@@ -170,4 +271,21 @@ mod tests {
         let mut lfsr = Lfsr::new(3, vec![0], vec![1, 2, 3]);
         lfsr.reset(vec![1, 2]);
     }
+
+    #[test]
+    fn test_berlekamp_massey() {
+        // polynomial: 1 + x + x^3 + X^4 + X^7 + x^10
+        let mut lfsr = Lfsr::new(10, vec![0, 1, 3, 4, 7], vec![1, 0, 0, 1, 0, 0, 1, 0, 0, 1]);
+
+        let sequence: Vec<u8> = (0..8).map(|_| lfsr.next()).collect();
+        let poly = Lfsr::berlekamp_massey(&sequence);
+
+        // polynomial 1 + X^3
+        assert_eq!(poly, vec![1, 0, 0, 1]);
+
+        let history = &sequence[sequence.len() - (poly.len() - 1)..];
+        let predicted = Lfsr::predict_next_from_poly(history, &poly);
+
+        assert_eq!(predicted, lfsr.next());
+    }
 }