LFSR Jump Ahead Algorithm A related question has been asked and answered here.  My question is specifically about figuring out C for very large J.  I have searched the web and found many scholarly articles on "jump ahead" and how to efficiently do so, but I am looking for a more "layman's" algorithm.

Here is the background information.  I am working through understanding the areas of weakness in using LFSR in cryptography.  To proof the concept, I have an encrypted message from a stream cipher using an N-bit LFSR.  The plain-text is ASCII, so I know the plain-text bits recurring at a specific interval (8), call it m.  From that, I get a set of non-consecutive key bits.  I have used the Berlekamp-Massey on the first 128 such bits at index multiples of m, then used the resulting coefficients and the first 64 bits passed to BM to initialize an LFSR.  I've run through the rest of the message proving to myself that the bits at multiple m are all predictable.

Now I would like to prove the remaining piece to myself, namely to get back the LFSR coefficients for the message, with only having a subset of the key bits and an associated LFSR for those bits.  In the question I made reference to, I am stuck trying to figure out $b_i=c_{(i\,m^{-1}\bmod p)}$.  In other words, $c_{(i\,m^{-1}\bmod p)}$ is a hugh jump forward with the LFSR, too large to just step to that point in the state.  As, I said, I found articles that discuss LFSR "jump ahead", but I think I need something more explicit on the approach.

For instance, a simple method (and slow), is to precompute the matrix power $A := F^J$(A Fast Jump Ahead Algorithm for Linear Recurrences in a Polynomial Space).  But I find no information on how to create the matrix $A$ from the coefficients.  To prove the concept, I don't need a fast algorithm, just something I can apply to close the loop.

An algorithm (Python, C, Java, et al), or implementation would be most helpful.  But even an article that give more complete details on the approach could serve the purpose.  Can you recommend something?

Here is a specific example using Fibonacci mode to establish conventions:

Coeffients: 1111100000010100101010001111110101111100110010110011101111001001
State:      0000111101011011011000011101100111001010000011001001100100011110
Output:     000011110101101101100001110110011100101000001100100110010001111010110011001100101111011110110001010001001110010011101010110101001001110000000101111110001001110001000010100111110111001001100001111001010011110110110010010110010001101000100001111101110011101100011001110100010111101111010101101001001000011100011001010110100010001111011110000110010100000000011100011111011110011111100110010110001100011101001111010101101101100010100110110110100000100011111001011111001111000010011011101000111110111101001111001101011101101001111110001111011111100100101010011000101011001010111001001001..


Shortened; see here for the origin.
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&gt;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\&amp;(0-b) when a\&amp;-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&lt;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 &lt;stdint.h&gt;
#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&lt;&lt;63)&gt;&gt;(j&amp;63);    // takes care of the low-order 6 bits of j
    if((j&gt;&gt;=6)!=0)
        for(v=p;;) {                // for remaining bits of j, if any
            if ((j&amp;1)!=0)           // perform x = x*v mod p
                for(u=x, x&amp;=-((w=v)&gt;&gt;63); w+=w; x^=(u=(p&amp;-(u&amp;1))^(u&gt;&gt;1))&amp;-(w&gt;&gt;63));
            if ((j&gt;&gt;=1)==0) break;  // we have handled all bits of j
            // perform v = v*v mod p
            for(v&amp;=-((u=w=v)&gt;&gt;63); w+=w; v^=(u=(p&amp;-(u&amp;1))^(u&gt;&gt;1))&amp;-(w&gt;&gt;63));
            }
    // perform p*reflected(s) and keep coefficients 127 to 64 of that into v
    for(w=(p+p)|1, v=0; v^=w&amp;-(s&gt;&gt;63), w+=w; s+=s);
    // perform u = v*x  mod p
    for(u=v&amp;-(x&gt;&gt;63); x+=x; u^=(v=(p&amp;-(v&amp;1))^(v&gt;&gt;1))&amp;-(x&gt;&gt;63));
    // here u is the Galois state after j steps; step it into x, reflecting
    for(w=64, x=0; x+=x+(v=u&amp;1), u=(p&amp;-v)^(u&gt;&gt;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 &lt;stdio.h&gt;
int main(void) {
    uint64_t j,
        p=0xf814a8fd7ccb3bc9,   // 1111100000010100101010001111110101111100110010110011101111001001
        s=0x0f5b61d9ca0c991e,   // 0000111101011011011000011101100111001010000011001001100100011110
        t=s,x;
    for(j=0; j&lt;800; ++j) {
        printf("%d", (int)(t&gt;&gt;63));
        // evolve the state t step by step
        x=p&amp;t; x^=x&gt;&gt;32; x^=x&gt;&gt;16; x^=x&gt;&gt;8; x^=x&gt;&gt;4; x^=x&gt;&gt;2; x^=x&gt;&gt;1; t=(t+t)|(x&amp;1);
    }
    printf("\n");
    // same sequence with direct access
    for(j=0; j&lt;800; ++j) printf("%d", (int)(fffibo64(p,s,j)&gt;&gt;63));
    printf("\n");
    return 0;
}


Try it online!