LFSR Jump Ahead Algorithm
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!