Fast and succinct random access into the Fibonacci Sequence.

Suppose that we wanted to query for values of Fibonacci numbers up to $F_{n = 2^{m+1}-1}$ for some fixed upfront $m$ . Clearly the fastest solution -- $\mathcal O(1)$ query and $\mathcal O(n)$ space -- is via direct tabulation.

Conversely, one very fast and popular method (binary; via identities for $F_{2k}$ and $F_{2k+1}$ requires only $\mathcal O(1)$ space:

typedef unsigned _BitInt(1024) U;
static U fib_fd(unsigned n){
  U a=0,b=1; /*F(0), F(1)*/
  int msb = 31 - __builtin_clz(n | 1u);
  for(int i=msb;i>=0;i--){
    U c = a * (2*b - a); /*F(2k)*/
    U d = a*a + b*b;     /*F(2k+1)*/
    if((n>>i)&1u){a=d; b=c+d;}else{a=c; b=d;}}
  return a;}

We can meet in the middle -- $\mathcal O(\text{popcount }n) \le \mathcal O(\log n)$ query and $\mathcal O(\log n)$ space -- as follows. We first pre-compute a table of $F_{2^k}\pm\{-1,0,1\}$ for $k=1..n$ via some clever algorithm. This requires $\log(n)$ auxiliary space space under a naive word model. Then notice $F_{a+b}=F_{b+1}\cdot F_{a}+F_{b}\cdot F_{a-1}$ .

Then we can write $F_n$ as $F_{\sum c_i 2^i}$ where $c_i$ are binary digits of n.

For $n < 1024$ we thus have:

typedef unsigned _BitInt(1024) U;
U a[]={0x0,0x1,0x2,0xd,0x262,0x148addwb,
0x5f6c7b064e2wb,0x80b0b47480d1a39bcca01dwb,
0xea0eeac2531f02ba57d6b8d2eef01bed585606b77ee2wb,
0x3063fbef018d3b2e449653f439a7c8916b8457f2e3f013074959cc9a5aacee984cdd0dc5cebd7e82fb2a9951dwb};
U b[]={0x1,0x1,0x3,0x15,0x3db,0x213d05wb,
0x9a661ca20bbwb,0xd039a6fbf25547aedf2ac5wb,
0x17ab6d825584b020438dc6f8aab80eca0ab1ea2592c3bwb,
0x4e4c2de9f74070d76307f6e4cdce1c3dab16355e1748e0dae9b5862a653ba81c62bbff81e69cb05f8e20081c5wb};
U c[]={0x1,0x2,0x5,0x22,0x63d,0x35c7e2wb,
0xf9d297a859dwb,0x150ea5b707326eb4aabcae2wb,
0x264c5c2e7ab6a04be90b3285d9a71088e0374a910ab1dwb,
0x7eb029d8f8cdac05a79e4ad90775e4cf169a8d50fb38f3e2330f52c4bfe896b4af990d47b55a2ee2894aa16e2wb};
U fib(U n){
  U p=0,q=1;
  for(int i=0;i<10;i++) if((n>>i)&1){
    U r=c[i]*p+b[i]*q, s=b[i]*p+a[i]*q; p=r; q=s;}
  return p;}

Testing on random values in the interval $[0, 1024)$ shows:

% ./fibs
correct for n=0..1023
table method:   5.021384 s, 502.14 ns/query
fast doubling:  8.426360 s, 842.64 ns/query
sink = 0x0