Galois LFSR
In a Galois LFSR with polynomial $P$ of degree $n$, the state is a vector of $n$ bits assimilated to the binary coefficients of terms of degree $n-1$ to $0$ of a polynomial; we note both the state and its polynomial $S_j$ when at step $j\ge0$ .
The next state is computed as $S_{j+1}=S_j\,x\bmod P$. Equivalently, the state is shifted towards the high-order coefficient, which is output and dropped from the state; a 0 is introduced as the low-order coefficient; then if the output was 1, state bits at indexes corresponding to non-zero terms of $P$ (ignoring the high-order term $x^n$) are complemented.
It follows that in a Galois LFSR,
$$S_j=S_0\,x^j\bmod P$$
which allows efficient direct access. Polynomial modular exponentiation can uses the same methods as integer one except using carry-less arithmetic, with the simple algorithm having cost $O(\log(j)\,n^2)$.
Also, the output of the Galois LFSR in $j$ steps is the coefficients of the quotient $Q$ of the polynomial division of $S_0\,x^j$ by $P$, starting from the coefficient of degree $j-1$ for the first output bit, down to zero for the $j^\text{th}$; and it holds $P\,Q+S_j=S_0\,x^j$.
Fibonacci LFSR
The question uses a Fibonacci LFSR, as apparent by having the state at the beginning of the output. Afterwards, the next output bit is the parity of the bitwise-AND of the $n$ previous bits and the coefficients defining the polynomial. I'll use the convention that the previously produced bit is ANDed with the coefficient of degree $n-1$ of the polynomial. With that convention, the low-order (constant) coefficient of the polynomial is the leftmost bit of the question's Coeffients, and the polynomial $P$ is $1+x+x^2+x^3+x^4+x^{11}+x^{13}+x^{16}+\dots+x^{55}+x^{56}+x^{57}+x^{60}+x^{63}+x^{64}$
with the term $x^{64}$ not shown in the question (and never combined with anything by bitwise-AND when stepping the LFSR).
Caution: standard literature on Fibonacci LFSRs and Berlekamp-Massey often use the so-called reflected or reciprocal polynomial, where coefficients of order $i$ and $n-i$ are exchanged (which might change the order of the polynomial). The present answer does not use this form for the polynomial.
Without proof: with our convention, when the low-order coefficient of the polynomial is $1$ (as customary in Berlekamp-Massey), Galois and Fibonacci LFSRs using the same polynomial produce identical output sequences, only from different starting point and with different states (when the lowest-order non-zero coefficient of the polynomial has order $l>0$, the Fibonacci LFSRs reaches the same periodic sequence as the Galois LFSR after $l$ steps at most, when the Galois LFSR always enters that periodic sequence right from the start; we have the option of removing the first $l$ bits and reducing the polynomial's degree to $n-l$, so that the low-order coefficient of the polynomial becomes 1$).
We can always convert a Galois state to Finnonacci state yielding the same output by stepping the Galois LFSR $n$ times; its output (in chronological order) is the desired Fibonacci state (produced by increasing order of the corresponding bits in the polynomial, that is reflected compared to $Q$ obtained by polynomial division of $S_0\,x^n$ by $P$, by exchanging coefficients of order $i$ and $n-1-i$); that's because the Fibonacci LFSR first outputs its state (starting frow low-order, in our convention).
By the same reasoning, when the low-order coefficient of the polynomial is $1$, the Galois state $S_0$ corresponding to a certain Fibonacci state must be such that the corresponding output $Q_0$ of the Galois LFSR in the next $n$ bits is the Fibonacci state (produced by increasing order of the corresponding bits in the polynomial). It thus holds that $P\,Q_0+S_n=S_0\,x^n$. It follows that the coefficients of the Galois state $S_0$ from $n-1$ to $0$ are the coefficients of $P\,Q_0$ from $2n-1$ to $n$, where $Q_0$ is the Fibonacci state reflected by exchanging coefficients of order $i$ and $n-1-i$.
When the low-order coefficient of the polynomial is $1$, to compute the state of a Fibonacci LFSR after $j$ steps, we can thus convert its state to Galois form, fast-forward that, and convert back.
With $\widetilde{F_j}$ the Fibonacci state at step $j$ with reflection by exchanging terms $i$ and $n-1-i$, and $\lfloor U/V\rfloor$ the quotient of the polynomial division of $U$ by $V$
$$\widetilde{F_j}=\left\lfloor{\left(\left\lfloor{\widetilde{F_0}\;P}/{x^n}\right\rfloor\,x^j\bmod P\right)x^n}/P\right\rfloor$$
In more details:
- Reflect the Fibonacci initial state, yielding $Q_0=\widetilde{F_0}$ (where the first outputs $0$ to $n-1$ are the terms of order $n-1$ to $0$);
- Compute $P\,Q_0$ and keep bits of order $2n-1$ to $n$, yielding $S_0$ from order $n-1$ to $0$;
- Compute $S_j=S_0\,x^j\bmod P$;
- Compute the quotient $Q_j$ of the polynomial division of $S_j\,x^n$ by $P$;
- Reflect $Q_j$ to compute the Fibonacci state $F_j$ after $j$ steps (which from order $0$ to $n-1$ is also outputs $j$ to $j+n-1$).
Note: Steps 4. and 5. are equivalent to stepping the Galois LFSR $n$ times starting from $S_j$, forming the Fibonacci state from order $0$ to $n-1$ (left to right in the question).
An alternative matrix method considers the transformation resulting from stepping the LFSR (Fibonacci or Galois) as multiplication of a (sparse) matrix. That matrix can be raised to the $j^\text{th}$ to obtain the matrix corresponding to stepping $j$ times, but the matrix is no longer sparse, and (using schoolbook matrix multiplication algorithms) the time grows $O(\log(j)\,n^3)$ rather than $O(\log(j)\,n^2)$, with space $O(n^2)$ rather than $O(n)$.
Here is hopefully C99-conforming code specialized for $n=64$ and Fibonacci form, heavily optimized for speed. Notably, the code:
- Processes in parallel over 64-bit and uses no table, which I guesstimate saves like a factor of 100 in speed compared to bit-by-bit on a 64-bit CPU, on top of having a cost growing with $j$ as $O(\log(j))$.
- Is branch-predictor friendly (the only test in inner loops is their end condition); this is obtained by masking rather than testing, often extending a low-order bit to 64-bit using negation:
-(uint64_t)1 is 0xFFFFFFFFFFFFFFFF while -(uint64_t)0 is 0x0000000000000000, and that's standard-conformant.
Off-topic rant: If you are a compiler writer in certain large company, please make it unnecessary to use a workaround; no I won't bow to a\&(0-b) when a\&-b is more recognizable.
- Stores the low-order coefficient of the polynomial $P$ in the high-order bit of variable $p$: that's a common implementation technique for Galois LFSRs, and eases conversion from bits in the order of the question to integer.
The code first performs the computation of $x^j$ in step 3, by binary exponentiation scanning $j$ from right to left (the low-order 6 bits are optimized out because $x^j$ requires a single shift for $j<2^6$); then steps 1 and 2 (the initial reflection reduces to scanning bits of $s$ in the appropriate order); then the rest of 3; then steps 4 and 5 are handled by stepping the Galois LFSR (the final reflection reduces to accumulating the result in the appropriate order).
I wanted to fit the function's code including comments in the ≈25 lines requiring no scrolling, so I (ab)used the C syntax, like using multiple operators , and separators ; on the same line.
// fffibo64: Fast-Forward a 64-bit Fibonacci LFSR
// Inputs:
// p The polynomial, of degree 64; high-order bit of p is the poly's
// constant term, and must be 1; coefficient x^^64 is implicit
// s Initial Fibonacci state (that is outputs 0 to 63 starting at high-order)
// j Number of steps
// Output: Fibonacci state after j steps (that is outputs j to j+63 starting high-order)
#include <stdint.h>
#ifdef _MSC_VER // this vendor has its own views on C, workaround that
#pragma warning (push)
#pragma warning (disable : 4146) // avoids "unary minus operator applied to unsigned type"
#endif
uint64_t fffibo64(uint64_t p, uint64_t s, uint64_t j) {
uint64_t u,v,w,x; // workhorse temporaries
// compute x^^j mod p into x, with optimization for low-order 6 bits
x=((uint64_t)1<<63)>>(j&63); // takes care of the low-order 6 bits of j
if((j>>=6)!=0)
for(v=p;;) { // for remaining bits of j, if any
if ((j&1)!=0) // perform x = x*v mod p
for(u=x, x&=-((w=v)>>63); w+=w; x^=(u=(p&-(u&1))^(u>>1))&-(w>>63));
if ((j>>=1)==0) break; // we have handled all bits of j
// perform v = v*v mod p
for(v&=-((u=w=v)>>63); w+=w; v^=(u=(p&-(u&1))^(u>>1))&-(w>>63));
}
// perform p*reflected(s) and keep coefficients 127 to 64 of that into v
for(w=(p+p)|1, v=0; v^=w&-(s>>63), w+=w; s+=s);
// perform u = v*x mod p
for(u=v&-(x>>63); x+=x; u^=(v=(p&-(v&1))^(v>>1))&-(x>>63));
// here u is the Galois state after j steps; step it into x, reflecting
for(w=64, x=0; x+=x+(v=u&1), u=(p&-v)^(u>>1), --w;);
return x;
}
#ifdef _MSC_VER
#pragma warning (pop)
#endif
// demo of the above, producing the output stream bit by bit both sequentialy, and with
// direct access, with (hopefully) identical result.
#include <stdio.h>
int main(void) {
uint64_t j,
p=0xf814a8fd7ccb3bc9, // 1111100000010100101010001111110101111100110010110011101111001001
s=0x0f5b61d9ca0c991e, // 0000111101011011011000011101100111001010000011001001100100011110
t=s,x;
for(j=0; j<800; ++j) {
printf("%d", (int)(t>>63));
// evolve the state t step by step
x=p&t; x^=x>>32; x^=x>>16; x^=x>>8; x^=x>>4; x^=x>>2; x^=x>>1; t=(t+t)|(x&1);
}
printf("\n");
// same sequence with direct access
for(j=0; j<800; ++j) printf("%d", (int)(fffibo64(p,s,j)>>63));
printf("\n");
return 0;
}
Try it online!
