14 of 15 index j was already unsed

We are given the sequence $b_i$ for $i$ multiple of $m>1$ with $0\le i<2n\,m$.

  • If all the given $b_i$ are identical, we can't compute any missing $b_i$ (beyond a wild guess they are all $b_0$), stop.
  • Build the sequence of $2n$ known bits $c_j=b_{(j\,m)}$ with $0\le j<2n$.
  • Use Berlekamp–Massey to find a LFSR (polynomial and initial state) matching $c_j$.
  • Find the minimal period $p$ of the $c_j$ :
    • for a small polynomial, we can simulate operation of the LFSR until the state repeats (OK to degree say 40 on a PC)
    • or we can use baby-step/giant-step (OK to degree say 75 on a PC, see final section)
    • or better, call math to the rescue!
      • for irreducible polynomial, the minimal period $p$ is the same starting from any point other than zero, and is some divisor of $2^k-1$ (it is exactly $2^k-1$ if and only if the polynomial is primitive and the starting point non-zero). We can factor $2^k-1$ then pull out the factors that we can while still getting a period, to quickly find the minimal period.
      • if the polynomial is non-irreducible, we can factor it into a product of irreducible polynomials, and if their degrees are $k_r$, then the minimal period $p$ divides $\displaystyle\prod(2^{k_r}-1)$ for all starting point (but potentially varies according to starting point). We can compute and factor that period, and quickly find the minimal period as above.
  • If $\gcd(m,p)\ne1$, there are missing $b_i$ for which we have no clue, stop.
  • $b_i$ hopefully also has period $p$, thus $\forall j\in\Bbb N,\ c_j=b_{(j\,m\bmod p)}$ ; it follows that $\forall i\in\Bbb N, b_i=c_{(i\,m^{-1}\bmod p)}$

$m^{-1}\bmod p$ needs to be computed only once, e.g. with the half extended Euclidean algorithm. Notice that for large polynomial, we can still efficiently compute $c_j$ for large arbitrary $j$ (see this), once we have the polynomial and starting state.

If we want the polynomial yielding the $b_i$, an easy option is to use Berlekamp–Massey again on the first $2n$ bits $b_j$. Some polynomial arithmetic can also give the result directly from the polynomial for the LFSR generating the $c_j$. The initial state of the LFSR for the $b_j$ (in Fibonacci form) is of course the first $k$ bits $b_i$.


Addition: The baby-step/giant-step method works to find in $O(2^{k/2})$ steps the period of any generator of period $<2^k$ allowing fast-forward:

  • select some "giant step" $s$, e.g. $s=2^{\lceil (k+\log_2k)/2\rceil}$
  • compute and store in a table allowing fast search the state of the generator at step $i\;s$ for $i$ in $(0,2^k/s]$
  • search in this table the state of the generator at step $j$ in $[0,s)$ (with "baby step" 1).

When we find a match, that tells state $j$ matches state $i\;s$, thus a period is $i\;s-j$. The smallest period is a divisor of that, found by factoring the known period and removing those prime factors that can be while keeping what remains a period (which can be efficiently verified).

Because the values in the table are well distributed, using them modulo a hash table size works beautifully. We can use separate tables according to some bits of the state (like the upper half or nearly that), saving space (slightly less than what we need to store $i$). For $k=64$ and $s=2^{35}$ we need $2^{29}$ entries and 6GB RAM is ample.

It is essential that the baby step is fast, and that the giant step is reasonably fast. With an LFSR converted to Galois form where $S_j=S_0\,x^j\bmod P$, the baby step is a multiplication by $x\pmod P$ (of cost $O(k)$ bit operations, or just S = (P&-(S&1))^(S>>1) if $S$ fits a single computer word), whereas the giant step is a multiplication by $x^s\bmod P$; by precomputing $x^s\bmod P$ is can be down to about $k$ times more costly.

The fudge factor $\log_2k$ in the recommended $s$ assumes moreless that giant/baby step cost ratio, but $s$ can be further increased if memory is an issue. Also, if the period is much smaller than $2^k$, we might not need to make all the giant steps, saving both on time and memory.