Introduction

Fibonacci is so simple, but so complex as to seem organic. It is a great homework exercise for Divacon.

Below are three derivations of a Fibonacci(N) function, (1) for \(N=2^a\), (2) for arbitrary \(N\), and (3) using Divacon to tighten things up.

The n'th (zero based) Fibonacci number F(n) is defined as:

F(0) = F(1) = 1
F(n) = F(n-1) + F(n-2)

Since each requires its immediate predecessors to calculate it, F(n) for large \(n\) seems to require O(n) operations. F(n) for \(n>2^{1000}\) could take a very long time indeed.

George Mou says he can do it in O(log(n)) operations.

That would be O(1000) time to calculate F(\(2^{1000}\)): perhaps a second or two on this web server (since the results below seem to be sub-millisecond), instead of

O(\(2^{1000}\)) ~ O(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000).

On a system achieving one addition per clock cycle, and terahertz speeds (faster than may ever be possible given the fade-out of Moore's Law), and with the age of the universe at \(4x10^{17}\) seconds, or roughly on the order of \(2^{60}\) seconds ago, it would take \(~2^{1000-40-60}=2^{900}\) times the age of the universe to solve such a problem directly. Just sayin'.

George suggested I use Divacon to do the same as an exercise. Here, then, is my homework.

Domain, Co-Domain, Operation, Co-Operation

The domain of the Fibonacci sequence is the set of Fibonacci numbers themselves, \(\forall\ n\in I\ \ F(n)\).

A convenient codomain would seem to be the set of pairs of Fibonacci summands: F(\(n\)-1), F(\(n\)-2), which sum to F(\(n\)). Can we please have a notation which hides the index arithmetic?

$$ \begin{align} for\ n\ \in \ I: & \\ \text{let}\ N'\ \ &\text{≜}\ \text{F(}n-1)\\ \text{let}\ N''\ &\text{≜}\ \text{F(}n-2) \\ \text{let}\ CAP(n) = N\ \ \ &\text{≜} \begin{bmatrix} \text{F(}n-1)\\ \text{F(}n-2) \end{bmatrix} = \begin{bmatrix} N' \\ N'' \end{bmatrix} \end{align}$$

What operation shall we associate with these vectors? Mapping from codomain to domain is accomplished by adding the summands. But how to combine them? Is there an \(N\# K\) in the codomain corresponding to f(n+k) in the domain?

This Fibonacci recurrence shifts from the n'th to the n+k'th (via a two-sample auto-convolution, meaning that the \(k\)'s count up while the \(n\)'s count down, and they are both taken from the same series):

 F(n+k) = F(k)*F(n-2)+F(k+1)*F(n-1)

This enables a doubling of the index from F(\(n\)) to F(\(2n\)) by setting \(k=n\). A doubling: a bit shift.

Arrays of Fibonacci results can be arranged by doublings to any desired rank. Allow me to overload the square bracket notation, to specify the rank of such a doubling. \(F[a]\), then, is an array of length \(2^{(a+1)}\):

 index: 0 1 2 3 4 5  6  7  8  9  10 11  12  13  14  15
F[0] = [1 1]
F[1] = [1 1 2 3]
F[2] = [1 1 2 3 5 8 13 21]
F[3] = [1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987]

F[\(a\)](n) (with value F(n)) will refer to the n'th element of the \(a\)'th of these doubled arrays, so for example
F[3](8)=34. Arrays F[a] are half redundant since F[a] is always the first half of F[a+1] and F(i)=F[a](i) wherever F[a](i) is defined. But it may help see and think clearly about the indexes involved in going from the end of one level to the end of the next. It is remarkable how confused I can become with such a simple function as Fibonacci!

Do you believe me when I say:

F[\(a\)] contains elements F[\(a\)](0) .. F[\(a\)]\((2^{(a+1)}-1)\).

To develop an intuition for the dependencies here, you might graphically draw out the Fibonacci numbers with their summands and daughters in a tree structure of dependencies. The subtrees responsible for a given sum recur, once then twice, then three, five, and eight times, in succeeding-Fibonacci-number dependency trees. If you skip a single index down in the array, you will find there are one of the first-smaller, two of the second-smaller, three of the third-smaller subtrees, \(F(n)\) of the \(n\)'th smaller subtrees, no matter what the value of \(n\). It's really a lovely self-inclusive mosaic pattern however you look at it. (I thought about drawing a fractal with it, but here's what I came up with, instead.)

Fortunately the F(\(n+k\)) recurrence above gives us a handle on this chaotic fractal. With it we may calculate, given the 2-vector \(N\) characterizing the \(n\)'th Fibonacci number, a second 2-vector, call it \((N+N)\) or \((2N)\) which characterizes the \(2n\)'th Fibonacci number, or shall we say, is that Fibonacci number in the codomaina. We will need an operation on 2-vectors yielding a 2-vector. The double-cross symbol, \(\#\), evokes two pairwise combinations: perfect.

Although a pair of products of 1x2 and 2x1 matrices amounts to the same thing I'll use two dot products in the following definition and derivation:

$$ \begin{align} N\ \#\ K\ \ &\text{≜} \begin{bmatrix} N \cdot (K+1)\\ N \cdot K \end{bmatrix} = \begin{bmatrix} N'(K+1)' + N''(K+1)'' \\ N'K' + N''K'' \end{bmatrix}\\ & =\begin{bmatrix} f(n-1)f(k+1-1) + f(n-2)f(k+1-2) \\ f(n-1)f(k-1) + f(n-2)f(k-2) \end{bmatrix} =\begin{bmatrix} f(n-1)f(k) + f(n-2)f(k-1) \\ f(n+(k-2)) \end{bmatrix} \\ & =\begin{bmatrix} f(n+(k-1)\\ f(n+(k-2)) \end{bmatrix} =\begin{bmatrix} (N+K)'\\ (N+K)'' \end{bmatrix} = CAP(n+k) \end{align}$$

That is, the codomain combination (\(\#\)) calculated according to the definition-given formula, becomes the domain combination (\(+\)) of indexes in the Fibonacci sequence.

The nervous reader, seeking to fit these calculations into a well-founded and incremental process (say, doubling the indexes) may wish to check that it doesn't assume its output before calculating it. Observe that although \((K+1)''=f(k+1-2)=f(k-1)\) is among the available components of \(N\) and \(K\), namely \(K'=f(k-1)\), we do not immediately find \((K+1)'=f(k)\) among those available components, which might suggest we are not ready to make this calculation. However the Fibonacci definition itself provides a substitution: \(f(k)=f(k-1)+f(k-2) = K' + K''\). This lets us proceed to double or shift an index while using only available components in \(N\) and \(K\) to calculate the shifted vector.

In a doubling process which at each step arrives at \(F[a]\) with \(n=2^a+1\) elements, although we may have calculated values \(F(k)\) up to \(k \le n\), we cannot assume availability of \(F(k)\) for \(k>n\). Since we mustn't get ahead of ourselves, our maximum shift will be achieved with \(k=n\), which gives a doubling. Great.

$$ \begin{align} N\ \#\ N\ \ &=\begin{bmatrix} N \cdot (N+1)\\ N \cdot N \end{bmatrix} &= \begin{bmatrix} N'(N+1)' + N''(N+1)''\\ N'N' + N''N'' \end{bmatrix} &= \begin{bmatrix} N' (N' + N'') + N'' N'\\ N'N' + N''N'' \end{bmatrix} &= \begin{bmatrix} (2N)'\\ (2N)'' \end{bmatrix} \end{align}$$

Note the doubled array F[\(a\)+1] has length \(m=2n=2^{(a+2)}\). And its final-element 2-vector (2N) depends (only!) on the final-element 2-vector (N) of the F[\(a\)] array. Who knew that such potentially remote entries in the Fibonacci series had such close relationships?!

Shall we look at an example? Consider the case of \(a=3\) so \(n=8, a+1=4, m=2n=16\). In the pattern we seek F(m-1) and F(m-2), so let's calculate F[4][14] and F[4][15] from F[3][6] and F[3][7].

Fib[4][14] = Fib[3][6]*Fib[3][6] + Fib[3][7]*Fib[3][7]
	   = 13*13 + 21*21
	   = 169 + 441
	   = 610

and

Fib[4][15] = Fib[3][7]*(2*Fib[3][6] + Fib[3][7])
	   = 21 * (2*13 + 21)
	   = 21 * (26+21)
	   = 21 * 47
	   = 987

To calculate each new 2-vector, (e.g., F(\(2^a\)-1) and F(\(2^a\)-2)) requires 2 multiplies and one add for each element.

For F(2^1000-1) requires 4000 multiplies and 2000 adds or in general O(log(N)) operations.

The following table is re-generated every time this page is loaded:

n	F(n-2)	F(n-1)	Δ hrtime
4	F(2) = 2	F(3) = 3	0.053311ms
8	F(6) = 13	F(7) = 21	0.00241ms
16	F(14) = 610	F(15) = 987	0.00137ms
32	F(30) = 1346269	F(31) = 2178309	0.00097ms
64	F(62) = 6557470319842	F(63) = 10610209857723	0.00097ms
128	F(126) = 1.5557697022053E+26	F(127) = 2.5172882568355E+26	0.0019ms
256	F(254) = 8.7571595343019E+52	F(255) = 1.4169381771406E+53	0.00152ms
512	F(510) = 2.7745922289306E+106	F(511) = 4.489384531331E+106	0.00167ms
1024	F(1022) = 2.7852935506996E+213	F(1023) = 4.5066996336778E+213	0.00165ms
2048	F(2046) = INF	F(2047) = INF	0.00122ms

These results are limited by PHP LONG INT and DOUBLE FLOAT computation, but it seems the algorithm and results are correct as far as the eye can see, comparing to some reference results:

F[14]=610, F[15]=987
F[30]=1346269, F[31]=2178309
F[62]=6557470319842, F[63]=10610209857723
F[126]=155576970220531065681649693, F[127]=251728825683549488150424261
...
F[1022] = 278529355069959292393881241266809350935330735212370380691318266898736\\
	  950320346518362561675961332445274995854966996688219111789542501520845\\
	  546940373127265215824082562848481813148554423082730494051913219529946\\
	  6733282
F[1023] = 450669963367781981310438323572888604936786059621860483080302314960003\\
	  064570872139624879260914103039624487326658034501121953020936742558101\\
	  987106764609420026228520234665586889971108924677841335400410363155392\\
	  5405243

---

This PHP code generates the above table (live! those durations were re-measured just now).

function Fibonacci($n) {
  $A = log($n,2)-1;
  $Fn_1 = $Fn_2 = 1;
  for ($a=1; $a<=$A; $a++) {
    $Fm_2 = ($Fn_2 * $Fn_2) + ($Fn_1 * $Fn_1);
    $Fm_1 = $Fn_1 * (2 * $Fn_2 + $Fn_1);
    $Fn_1 = $Fm_1;
    $Fn_2 = $Fm_2;
  }
  return "F(" . ($n-2) . ") = $Fm_2, F(" . ($n-1) . ") = $Fm_1\n";
}

Containing an actual log(n) function call and an iteration counting up to that many, this is clearly O(log(n)) in complexity.

So far so good, but so far this only works on powers of 2.

Fib(\(p\)) for \(p\) not a power of 2

Every counting number \(p\) is a sum of powers of 2, like 5 is a sum of \(2^2 + 2^0\), which are the binary digits in \(p\): decimal 5 = binary 101.

Considering \(p\) as a data structure comprising various properties, we might give names to

• the list of indexes of all the "1" bits in \(p\) (pbits),
• the count of how many there are (pbits.length), and
• the highest index of all the "on" bits in p (pbits.max).

For example, for \(p=17=2^0+2^4\), pbits=[0,4], pbits.length=2, pbits.max=4.

Each bit being a power of 2 also corresponds to a doubling of the Fibonacci array, from F[\(a\)] to F[\(a\)+1], etc. And each \(N\)-sized array has its 2-vector, \([\ N'\ N''\ ]^T\).

Now we approach the nub of the question, and the answer is Yes:

Let \(n\) be a power of 2, \(n=2^{(a+1)}\), so F(n-1) is the last element in F[\(a\)].
Let \(0\le k\lt n\). Then F(n+k) is in the 2nd half of F[\(a+1\)]. Sorry for the index jiu jitsu! Please compare definitions above to satisfy yourself.

Using the recurrence equation, F(n+k) depends on values \(N', N'', K', K''\), all within F[\(a\)] or (redundantly) in the first half of F[\(a\)+1].

In a divide-and-conquer approach we might knock off the highest bit of \(p\) by collecting for its highest bit \(2^a=n\) the values \(N', N''\), and then, for one on bit at a time of rank, say, b, such that \(k=2^b\) and \(b We could also proceed in the other direction, accumulating K',K'' for \(k\) comprised of lower bits and calculating (N+K)' and (N+K)'' for the next encountered bit \(n=2^a\) for some larger \(a\).

Either direction requires that we separately solve the trivial cases for n=0 itself as well as the 0th bit n=1, and then with those we can repeatedly build up the Fibonacci number 2-vectors for powers of two.

The modified recurrence allows us to combine the parameters for any two \(n\) and \(k\) elements to create the parameters for the \(n+k\)'th element. That, combined with parameters of the next on bit, lets us repeatedly combine pairs up to the last one, at which the summary is to say, \(P = [ P' P'']^T\).

To summarize the process, we will scroll up from the 0'th to the highest bit in \(p\) collecting intermediate results which are our parameters, N', N'' for each level.

Then we will scroll across again, incrementally incorporating previous \(k\)'th results \(K', K''\), with on-bit-triggered single-bit results N', N'' for each of the on bits, using the recurrence for f(n+k).

Here is an implementation in PHP:

function as_bits($p) {
  $res = [];
  for ($i=$j=1;$i<$p;$i<<=1,$j++)
    ;
  $j--; 
  for (  ; $j>=0; $j--,$i>>=1  ) {
    if ($p&$i) {
      $p-=$i;
      array_push($res,$j); // smaller j's go to the end of the results array.
    }
  }
  return $res;
}

function CombinePairs($pk, $pn) {
  $pth = ($pk["pen"] + $pk["ult"]) * $pn["pen"] + ($pk["pen"] + 2*$pk["ult"]) * $pn["ult"];
  $res = [];
  $res["ult"] = ($pk["ult"] * $pn["pen"]) + (($pk["pen"]+$pk["ult"])*$pn["ult"]);
  $res["pen"] = $pth - $res["ult"];
  return ($res);
}

function AnyFibonacci($p) {
  $pbits = as_bits($p);

  // create per-bit intermediate results
  $Fpair = [];
  $Fpair[0] = array("ult" => 1,"pen" => 0); 
  $Fpair[1] = array("ult" => 1,"pen" => 1);
  if ($pbits[0]>0) {
    for ($a=2; $a<=$pbits[0]; $a++) {
      $Fpair[$a] = CombinePairs($Fpair[$a-1],$Fpair[$a-1]);
    }
  }

  // Combine from a=1 up to the highest in pbits.
  $len = count($pbits);
  if ($len==0) // p has all bits zero
    return 1;
  else if ($len==1) { // p has only one bit set to one
    return $Fpair[$pbits[0]]["ult"]+$Fpair[$pbits[0]]["pen"];
  } else if ($len>1) { // p has more than one bit set to one.
    $intres = CombinePairs($Fpair[$pbits[$len-1]],$Fpair[$pbits[$len-2]]);
    for ($i = $len-3; $i >= 0; $i--) {
      $intres = CombinePairs($Fpair[$pbits[$i]],$intres);
    }
    $res = $intres["ult"] + $intres["pen"];
    return $res;
  }
}