Project Euler Solution 25: 1000-digit Fibonacci number

Code & Zahlen

In today's installment of the Project Euler series we have Problem 25: 1000-digit Fibonacci number, where we have to find the Fibonacci number with 1000 digits.

The Fibonacci sequence is defined by the recurrence relation:

$F_n = F_{n−1} + F_{n−2}$, where $F_1 = 1$ and $F_2 = 1$.

Hence the first 12 terms will be:

  • $F_1 = 1$
  • $F_2 = 1$
  • $F_3 = 2$
  • $F_4 = 3$
  • $F_5 = 5$
  • $F_6 = 8$
  • $F_7 = 13$
  • $F_8 = 21$
  • $F_9 = 34$
  • $F_{10} = 55$
  • $F_{11} = 89$
  • $F_{12} = 144$

The 12th term, $F_{12}, is the first term to contain three digits.

What is the index of the first term in the Fibonacci sequence to contain 1000 digits?

I am not sure about this problem. In Solution 2: Even Fibonacci numbers we constructed a Fibonacci sequence generator which has O(1) time per number. We can just use this one to enumerate the whole sequence and stop whenever we have found one with 1000 digits:

def solution_iterative() -> int:
    for i, f in enumerate(fibonacci_generator(), 2):
        if len(str(f)) >= 1000:
            return i

This runs in 26 ms and is quite sufficient for this problem.

Perhaps we can use the diagonalization of the the Fibonacci sequence to get it even faster? Taking the derivation from the website, we can write the number $F_n$ as such: $$ F_n = \frac{1}{\sqrt 5}(\phi^n - \psi^n) \,, $$ where $$ \phi = \frac{1 + \sqrt 5}{2} \,, \qquad \psi = \frac{1 - \sqrt 5}{2} \,. $$

We could just put in $F_n = 10^{1000}$ and solve for $n$. But how would one do that? The problem is that $\psi$ is negative and flips it sign with every iteration. One could just take a look at even or off $n$ and then analytically continue that function.

Then one would use a numerical equation solver to solve this equation for $n$ and just take a look at the Fibonacci numbers in that region. This could be done with an open method like Newton or secant. Of one could establish an upper bound and then use bisection.

The problem here are the IEEE floating point numbers with 64 bit. They can only go up to around $10^{307}$, which is not sufficient for the $10^{1000}$ that we want. Therefore we cannot compute these high numbers directly with the precision that we have. We would need to have around 300 bits of precision to pull that off. And we don't have this, the largest that NumPy offers are 256 bit. One could then use an arbitrary precision library to get this, but that will hardly be faster than just enumerating about 5000 Fibonacci numbers (the solution is 4,782).

So in the end, the straightforward solution is the best one.