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.

• I'm using Fibonnaci mode. I will update question with a specific example. – Les Jul 9 '18 at 15:07

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: 1. 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$$); 2. Compute $$P\,Q_0$$ and keep bits of order $$2n-1$$ to $$n$$, yielding $$S_0$$ from order $$n-1$$ to $$0$$; 3. Compute $$S_j=S_0\,x^j\bmod P$$; 4. Compute the quotient $$Q_j$$ of the polynomial division of $$S_j\,x^n$$ by $$P$$; 5. 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! • Some questions: 1) Are$P$in Galois and$P$in Fibonacci the same? The coefficients in question were obtained by BMA using the State in question, so if I add one more bit (equal to 1) to the right of the coefficients, then I would have the$P$polynomial, is that correct? 3) Can I summarize your solution as "transform Fibonacci to Galois and solve for Galois"? And 4), if so, would I do better to just convert my whole problem to Galois mode? Thanks for your time and patience. – Les Jul 10 '18 at 0:18 • @Les: 1) yes, with the convention I use,$P$is the same in Galois and Fibonacci; BUT that makes$P$in Fibonacci contrary to your/usual convention. Your "one more bit" is term 64 in my convention, 0 in your/usual convention, and yes that makes the full poly. 3)[sic] Yes. I'm not saying that's the only way. 4) Fibonacci is fine, conversion is a small part of the cost, and reflexion can actually require no code. I hope to post a short C routine, handling$n=64$– fgrieu Jul 10 '18 at 0:49 • I found a major hurldle with changing my approach to Galois, I cannot find a way to translate the Fibonacci "seed" (which is just the first N bits of the sequence) to a Galois seed. BM doesn't help, all implementations of BM use Fibonacci, and all scholarly articles describe BM in terms of Fibonacci. Without a translated seed for Galois, I have no$S_0$G-state and cannot do the$S_j:=S_0x^j mod P$transformation. Unless, I can find a way to translate F-seed to G-seed, I will turn back to my original approach of finding the$A:=F^J$matrix. – Les Jul 10 '18 at 20:59 • -i have reviewed the code and it's output. At the top, comments say you are fast forwarding and outputting a Fibonacci state at step J, but at the bottom, you say you have computed the Galois state at J. It looks like the Galois state, which i can't use since I don't know and can't find the Galois starting state for my sequence. – Les Jul 10 '18 at 22:47 • @Lee: the answer does the transformation from Fibonacci to Galois at steps 1 and 2; and back from Galois to Fibonacci at steps 4 and 5. The function fffibo64 has a Fibonacci-only interface and is directly usable to predict the sequence (as demoed in main), and performs Fibonacci to Galois incompute p*reflected(s) and keep coefficients 127 to 64 of that into v (where the expression w=(p+p)|1 re-introduces the coefficient of$P$missing in p), and back in code on the line under here u is the Galois state after j steps; step it into x, reflecting. – fgrieu Jul 11 '18 at 6:46 The matrix representation seems preferable in your case because • You indicate that you already know how to solve your problem in this form (by matrix exponentiation) • You said you are not concerned about performance The matrix representation is preferable in general becase • You can apply the tools of linear algebra to understand it • LFSR representation (Galois or Fibonacci) becomes irrelevant • It generalises beyond$GF(2)$Obtaining the matrix representation The matrix representation of a linear function can be obtained by evaluating it at points corresponding to rows of the identity matrix: $$A_0 = f(1,0,0,\ldots,0)$$ $$A_1 = f(0,1,0,\ldots,0)$$ $$A_2 = f(0,0,1,\ldots,0)$$ $$A_n = f(0,0,0,\ldots,1)$$ In other words, you just initialize the LFSR state to$0x0001$,$0x002$,$0x0004$, .etc, and each output becomes one row of the matrix. Then, given an arbitrary input vector$x$, you can obtain the same result as$f(x)$by matrix multiplication with$A$. Fast-forwarding the state Since LFSR iteration is just repeated application of$f$we have $$f(x) = Ax$$ $$f(f(x)) = AAx$$ $$f(f(f(x))) = AAAx$$ And in general $$f^n(x) = {A^n}x$$ For large$n$we can carry out the exponentiation in$O(\log(n))$operations, using the same exponentiation by squaring algorithm that is used in RSA encryption. For example$A^8 = A^{2\cdot2\cdot2} = ((A^2)^2)^2$. • Yes. That method has the advantage of working with any linear (extensible to affine) transformation without having to analyses its details. However, it is$O(m^3\,\log(n))$time,$O(m^2)$memory (with$m$the LFSR size and that answer's use of$n$), due to the cost/space of the matrix multiplication. It was fun to find an$O(m^2\,\log(n))$time,$O(m)\$ memory algorithm for the Fibonacci case; I only knew this was possible for the Galois case. – fgrieu Jul 10 '18 at 11:05
• @conchild - I think this may be the approach. I have tried switching to Galois to begin with and have hit a road block at determining the correct Seed. – Les Jul 10 '18 at 21:00
• @conchild - But with a Fibonacci LFSR, the bits of the initial state are the 1st N bits of the output. That would produce an A = I, which gets me nowhere (AAA=A=I). Or am I missing something? However, I could use those values (1,2,4,8,...) to initialize state of a Galois LFSR and get very different rows for A. What do you think? – Les Jul 11 '18 at 0:13
• @Les I think you should try it with Fibonacci form, observe that it works, then reason about it afterwards. You won't get A=I, but it might be close to the identity matrix. – conchild Jul 11 '18 at 7:05