Berlekamp-Massey
berlekamp_massey(A, m)
returns a vector C
of length \(n\) which satisfies
\[ \begin{aligned}
A_x &= \sum_{i=0}^{i=k-1} {C_i A_{x-k+i}} \\
&= C_0 A_{x-k} + C_1 A_{x-k+1} + \cdots + C_{k-1} A_{x-1}
\end{aligned} \]
with minimum \(n\) under prime modulo \(m\). It is safe to have the length of vals
as at least \(3n\).
Example
fn main() { // vals[x] = vals[x-3] + 2*vals[x-2] + 3*vals[x-1] let m: u64 = 1000000007; let mut vals: Vec<u64> = vec![1, 2, 3]; for x in 3..20 { vals.push((vals[x - 3] + 2 * vals[x - 2] + 3 * vals[x - 1]) % m); } let rec = berlekamp_massey(&vals, m); println!("{:?}", rec); // [1, 2, 3] } // Berlekamp-Massey // References // https://blog.naver.com/jinhan814/222140081932 // https://koosaga.com/231 fn rem_pow(mut base: i64, mut exp: i64, m: i64) -> i64 { let mut result = 1; while exp != 0 { if exp & 1 != 0 { result = (result * base) % m; } exp >>= 1; base = (base * base) % m; } result } /// Finds rec[n] which satisfies /// vals[d] = rec[0]vals[0] + rec[1]vals[1] + ... + rec[d-1]vals[d-1] /// with minimum n. fn berlekamp_massey(vals: &[u64], m: u64) -> Vec<u64> { let m = m as i64; let mut cur: Vec<i64> = Vec::new(); let (mut lf, mut ld) = (0, 0); let mut ls: Vec<i64> = Vec::new(); for i in 0..vals.len() { let mut t = 0; for (j, v) in cur.iter().enumerate() { t = (t + vals[i - j - 1] as i64 * v) % m; } if (t - vals[i] as i64) % m == 0 { continue; } if cur.len() == 0 { cur = vec![0; i + 1]; lf = i; ld = (t - vals[i] as i64) % m; continue; } let k = -(vals[i] as i64 - t) * rem_pow(ld, m - 2, m) % m; let mut c: Vec<i64> = vec![0; i - lf + ls.len()]; c[i - lf - 1] = k as i64; for (p, j) in ls.iter().enumerate() { c[i - lf + p] = -j * k % m; } if c.len() < cur.len() { c.extend((0..(cur.len() - c.len())).map(|_| 0)); } for j in 0..cur.len() { c[j] = (c[j] + cur[j]) % m; } if i - lf + ls.len() >= cur.len() { ls = cur; lf = i; ld = (t - vals[i] as i64) % m; } cur = c; } for i in 0..cur.len() { cur[i] = (cur[i] % m + m) % m; } cur.into_iter().rev().map(|x| x as u64).collect() }
Code
// Berlekamp-Massey
// References
// https://blog.naver.com/jinhan814/222140081932
// https://koosaga.com/231
fn rem_pow(mut base: i64, mut exp: i64, m: i64) -> i64 {
let mut result = 1;
while exp != 0 {
if exp & 1 != 0 {
result = (result * base) % m;
}
exp >>= 1;
base = (base * base) % m;
}
result
}
/// Finds rec[n] which satisfies
/// vals[d] = rec[0]vals[0] + rec[1]vals[1] + ... + rec[d-1]vals[d-1]
/// with minimum n.
fn berlekamp_massey(vals: &[u64], m: u64) -> Vec<u64> {
let m = m as i64;
let mut cur: Vec<i64> = Vec::new();
let (mut lf, mut ld) = (0, 0);
let mut ls: Vec<i64> = Vec::new();
for i in 0..vals.len() {
let mut t = 0;
for (j, v) in cur.iter().enumerate() {
t = (t + vals[i - j - 1] as i64 * v) % m;
}
if (t - vals[i] as i64) % m == 0 {
continue;
}
if cur.len() == 0 {
cur = vec![0; i + 1];
lf = i;
ld = (t - vals[i] as i64) % m;
continue;
}
let k = -(vals[i] as i64 - t) * rem_pow(ld, m - 2, m) % m;
let mut c: Vec<i64> = vec![0; i - lf + ls.len()];
c[i - lf - 1] = k as i64;
for (p, j) in ls.iter().enumerate() {
c[i - lf + p] = -j * k % m;
}
if c.len() < cur.len() {
c.extend((0..(cur.len() - c.len())).map(|_| 0));
}
for j in 0..cur.len() {
c[j] = (c[j] + cur[j]) % m;
}
if i - lf + ls.len() >= cur.len() {
ls = cur;
lf = i;
ld = (t - vals[i] as i64) % m;
}
cur = c;
}
for i in 0..cur.len() {
cur[i] = (cur[i] % m + m) % m;
}
cur.into_iter().rev().map(|x| x as u64).collect()
}