4

I was doing a project Euler problem where I needed to find several Fibonacci numbers, but their index was so large that I could not use the typical recursive method. Instead, I used Binet's rule:

$$ F_n = \frac{\phi^n - \psi^n}{\sqrt{5}}$$

But since phi and psi and root 5 are all irrational, I didn't know how many digits past the decimal place I needed to drag them out, especially since the Fibonacci numbers were hundreds or thousands of digits long. I ended up making the irrational numbers obscenely accurate to 1000 or 10000 decimal places just to be safe.

Is there some kind of rule like "x amount of precision before will result in y amount of precision in the answer"?

I have already solved the problem just fine, I was just wondering about the precision part.

Ryan
  • 1,200
  • 1
    Note that $F_n$ can be calculated by rounding $\frac{\phi^n}{\sqrt{5}}$ correctly to the next integer, where $\phi=\frac{\sqrt{5}+1}{2}$, so you only need the power of one irrational number. Additionally, there are programs that can easily handle natural numbers of the magnitude you want to calculate, so my suggestion would be to use the recursion – Peter Jan 03 '17 at 00:01
  • @peter Yeah I know that, but as n increases, there comes a point where $F_n$ is not even accurate to +-.5 with just floating point decimal precision – Ryan Jan 03 '17 at 00:10
  • A program that can handle numbers with many decimal digits usually can handle very large natural numbers as well. At some point, even the most powerful programs will fail to calculate $F_n$ exactly, but it is still easy to get very good approximations using the powers of $\phi$ – Peter Jan 03 '17 at 00:14
  • 1
    Have you tried keeping a table of sparse values and then applying the Fast Doubling method on this page? I have found this algorithm extremely fast in practice. (Note that using $F(3k)$ to get a tripling method wasn't as fast actually... can't remember why). Given that we are on computers and working with the number $2$, significant speed-ups can occur. A naive application alone will halve the index of your Fibonacci numbers, which is exponentially (correct word for once!) faster. – Brevan Ellefsen Jan 03 '17 at 00:15
  • Looking at the notes, calculating $F{1000000}$ takes $1$ second and calculating $F_{4000000}$ takes only 7 seconds using this approach. If you are in a time crunch and your indices are of size 10000000 (10 million in the USA) or higher you are going to encounter some problems getting an exact approach anyway. – Brevan Ellefsen Jan 03 '17 at 00:23
  • @BrevanEllefsen: looking at what notes? – Rob Arthan Jan 03 '17 at 01:02
  • To all, I already solved the problem. I was just interested in the precision aspect of the problem that I mentioned in the question – Ryan Jan 03 '17 at 01:07
  • @BrevanEllefsen: we may even improve such fast doubling procedure in order that the computation of $F_n$ only involves a few squarings. Code disclosed in my answer. – Jack D'Aurizio Jan 03 '17 at 01:14

2 Answers2

5

You do not really need to work in floating point: Fibonacci numbers fulfill interesting duplication identities, so there are fast algorithms in integer arithmetic, as Brevan Ellefsen mentioned.

Here it is some $\text{Maxima}$ code for the fast computation of the $n$-th Fibonacci number.

 /* Fast Fibonacci Numbers Computation */
 fastfib(n):=block([O:2,L:?integer\-length(n),i,f0:2,f1:1,u,v],
   for i:1 thru L do (
    u:f0^2, v:f1^2,
    if ?logbitp(L-i,n) then (f1:v+O, f0:f1-u+O, O:-2  )
                       else (f0:u-O, f1:v+O-f0, O:2   )),
  return((?ash(f1,1)-f0)/5)
  );

This is included in the $\text{gf}$ package. In order to compute $F_n$, we just need to compute $2\log_2(n)$ squares. With Karatsuba algorithm (invoked in $\text{Maxima}$ through the $\text{fasttimes}$ command) or similar approaches the computation of $m^2$ requires just $\log_2(m)^{1.585}$ elementary operations.

This is pretty efficient. Here it is some test on my machine:

Maxima 5.32.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.10 (a.k.a. GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) load("gfpatch.mac");
(%o1)                             gfpatch.mac
(%i2) showtime:true;
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%o2)                                true
(%i3) fastfib(10000000)$
Evaluation took 0.0800 seconds (0.0900 elapsed)

Back to the actual question: since $\phi^n=\exp(n\log\phi)$ and $\frac{d}{dx}x^n = n x^{n-1}$, we roughly have $(\phi+\varepsilon)^n \approx \phi^n+n\varepsilon\phi^{n-1}$, so a good accuracy is achieved as soon as $\varepsilon\ll\frac{1}{n\phi^{n-1}}$, and we need to store $\approx 0.7n$ binary digits of $\varphi$ to compute $\varphi^n$ with decent accuracy.

Jack D'Aurizio
  • 353,855
  • 2
    Nice, haven't seen this before. Will keep in mind for future applications! Worth noting though that this doesn't answer the question the OP seemingly intended to ask, namely how many decimal digits of $\phi$ must be stored to calculate $\phi^n$ with an error of magnitude less than or equal to $\frac{1}{2}$?. Upvoted regardless - this is a far better solution for anything the OP likely needs to calculate – Brevan Ellefsen Jan 03 '17 at 01:21
  • @BrevanEllefsen: you are right, answer improved. – Jack D'Aurizio Jan 03 '17 at 01:34
2

You could definitely improve the calculations by considering only $\phi$ $$F_n=\left[\frac{\phi^n}{\sqrt{5}}\right]$$ Re the error, taking a $\phi^{*} \in \mathbb{Q}$ such that $|\phi - \phi^{*}|<\frac{1}{n2^{n-1}10^m}$, leads to: $$|\phi^n - (\phi^{*})^n|=|\phi - \phi^{*}|\cdot |\phi^{n-1}+\phi^{n-2}\phi^{*}+...+(\phi^{*})^{n-1}|\leq |\phi - \phi^{*}|\cdot n2^{n-1}<\frac{1}{10^m}$$

On another note, the very same wikipedia article suggests a procedure that:

These last two identities provide a way to compute Fibonacci numbers recursively in $O(\log{n})$ arithmetic operations ...

rtybase
  • 16,907