16
\$\begingroup\$

The Fibonacci sequence is defined by the recurrence relation:

Fn = Fn−1 + Fn−2, where F1 = 1 and F2 = 1.

Hence the first 12 terms will be:

F1 = 1
F2 = 1
F3 = 2
F4 = 3
F5 = 5
F6 = 8
F7 = 13
F8 = 21
F9 = 34
F10 = 55
F11 = 89
F12 = 144

The 12th term, F12, is the first term to contain three digits.

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

from time import time


def fib(n):
    """returns index of the first fibonacci of length n digits."""
    fibs = [0, 1]
    while len(str(fibs[-1])) < n:
        fibs.append(fibs[-1] + fibs[-2])
    return len(fibs) - 1


if __name__ == '__main__':
    time1 = time()
    print(fib(1000))
    print(f'Time: {time() - time1} seconds.')
\$\endgroup\$
  • 4
    \$\begingroup\$ If you examine the Wikipedia article on Fibonacci numbers, in particular the section "Computation by Rounding, it should be possible to calculate the the first F_n > 10^1000 by taking base 10 logs of both sides. \$\endgroup\$ – James K Polk Jul 20 at 13:27
  • \$\begingroup\$ @JamesKPolk Exactly. \$\endgroup\$ – Eric Duminil Jul 21 at 13:18
18
\$\begingroup\$

It's possible to calculate the nth Fibonacci number directly thanks to Binet's formula.

Due to floating point errors, the formula doesn't give correct values if n is too large, but it works fine in order to calculate the number of digits.

from math import log10, floor, ceil

def fibonacci_digits(n):
    if n < 2:
        return 1
    ϕ = (1 + 5**0.5) / 2
    return floor(n * log10(ϕ) - log10(5) / 2) + 1

It's possible to calculate the inverse of this function. This way, you can find the first Fibonacci with k digits directly, without any loop:

def euler25(k):
    if k < 2:
        return 1
    ϕ = (1 + 5**0.5) / 2
    return ceil((k + log10(5) / 2 - 1) / log10(ϕ))

It seems to return the correct answer for any k in less than 1 μs :

>>> euler25(1000)
4782
>>> euler25(10_000)
47847
>>> euler25(100_000)
478495
>>> euler25(100_000_000_000_000_000)
478497196678166592

The digits are the same as in:

>>> 1 / log10(ϕ)
4.784971966781666
\$\endgroup\$
  • 1
    \$\begingroup\$ Oh god, φ as an identifier. Can't wait to start seeing bugs like φ = 1; assert ϕ == 1... \$\endgroup\$ – Peilonrayz Jul 21 at 4:06
  • \$\begingroup\$ @Peilonrayz: I agree it can be dangerous. I used two different symbols at first, without noticing it. :-/ Your example isn't a bug, it's a feature! \$\endgroup\$ – Eric Duminil Jul 21 at 7:28
  • 2
    \$\begingroup\$ You say feature, I say death trap. ;) Yeah I know about NFKC normalization, that's how I made my example legit ;P \$\endgroup\$ – Peilonrayz Jul 21 at 7:33
  • 2
    \$\begingroup\$ @Peilonrayz: i leads the pack with 21 collisions: "iᵢⁱℹⅈⅰⓘi𝐢𝑖𝒊𝒾𝓲𝔦𝕚𝖎𝗂𝗶𝘪𝙞𝚒". Next time I want to troll a colleague, I'll be sure to use every single synonym for i in a script. :D \$\endgroup\$ – Eric Duminil Jul 21 at 8:26
  • 2
    \$\begingroup\$ You've gained my eternal ire and awe. :O Also, 𝚄𝓃𝙞𝙘𝔬de 𝒊𝘴 𝚝ⓗ𝕖 ⓝ𝔢𝑤 |𝟯³𝟳  ˢ𝒑ⅇ𝖆𝙠.‎ - Yeah there's an impressive amount of collisions, which should make Unicode easier to support. \$\endgroup\$ – Peilonrayz Jul 21 at 8:37
17
\$\begingroup\$

Two comments:

  1. There isn't any reason to keep a list of all the fib numbers. You only need the last two numbers to compute the next one.

  2. A number with 1000 digits satisfies 1e1000 <= f < 1e1001. It might be faster to compare f to the numbers rather than to convert it to a string and get its length.

def fib_of_length(n):
    """returns index of the first fibonacci of length n digits.
       for n > 0.
    """

    bound = 10**(n-1)
    fib, previous_fib = 1,0
    index = 1
    while fib < bound:
        fib, previous_fib = fib + previous_fib, fib
        index += 1

    return index
\$\endgroup\$
  • \$\begingroup\$ It's reasonably fast, taking into account that you iterate over every Fibonacci number. \$\endgroup\$ – Eric Duminil Jul 21 at 8:11
15
\$\begingroup\$

This isn't a full review but something I found very interesting about this problem: I would have thought that keeping a list of every Fibonacci number would get expensive over larger numbers but that doesn't seem to be the case with my testing.

Instead it seems the bottleneck is the line:

len(str(fibs[-1])) < n:

I ran the function for n = 10,000 and it took roughly 32 seconds to run. If you instead do a comparison to a constant value such as:

threshold = 10 ** (n-1)
while fibs[-1] <= threshold:

The time drops to about a 10th of a second. If you change your list to only track the last three numbers that times is roughly halved again.

\$\endgroup\$
  • 1
    \$\begingroup\$ You should consider moving the actual computation out to a constant ... threshold = 10 ** (n - 1) ... then later you have fibs[-1] <= threshold .... so you only have to do that computation once. \$\endgroup\$ – rolfl Jul 22 at 11:26
  • 1
    \$\begingroup\$ Instead of just tacking on text at the end, fix your initial answer to show the good way first. Some Python implementations would CSE that loop-invariant 10 ** (n-1), but unlike C/C++ you apparently have to do that optimization manually. Anyway yes, converting a huge integer from binary to a base-10 string is very slow, requiring repeated division by 10 or equivalent. \$\endgroup\$ – Peter Cordes Jul 22 at 19:59
10
\$\begingroup\$

Project Euler has quite a few problems on the Fibonacci numbers. Other things also repeatedly crop up, like prime numbers, triangular numbers, and so on. In fact, quite a few problems build on earlier, easier, problems.

Because of this it is worth it to start designing your functions for reusability. Put these functions into modules, so you can reuse them later. This one should probably go to the sequences.py module or just a general utils.py.

Currently you have one function that does everything. Instead you will probably want a generator function that generates all Fibonacci numbers, which you can then use afterwards:

def fibonacci():
    """Yield the numbers from the Fibonacci sequence"""
    a, b = 1, 1
    while True:
        yield a
        a, b = b, a + b

The actual problem can then be solved in many ways, either like you did or by writing another generally applicable function, first:

def first(it, condition=bool):
    """Return the first truthy (default) or first value satisfying `condition` from `it`."""
    it = iter(it)
    return next(filter(condition, it))

The itertools module is well worth a read, although it is not needed here (I had used it in an earlier version). If you want to spend an entire afternoon, also have a look at the more_itertools module.

The actual solution is then:

if __name__ == "__main__":
    condition = lambda n: 10**999 <= n[1] < 10**1000
    print(first(enumerate(fibonacci(), 1), condition)[0])
\$\endgroup\$
1
\$\begingroup\$

As others have pointed out, if you are brute-forcing this:

  • there's no reason to keep the old values around in a list. You just need a += b; b += a; to go 2 steps, or in Python a,b = b,a+b to go one step.
  • binary integer -> base-10 digit string is very expensive, requiring repeated division by 10 or equivalent. (And that's division of the whole BigInteger by 10; it can't be done in chunks because the lowest decimal digit depends on all higher binary digits.) So compare against a threshold instead.

There's another trick nobody's mentioned:

Fib(n) grows fast enough that you can discard some of the the lowest digits occasionally without affecting the most-significant digits of your running total. Fib(n) grows close to exponentially; faster than any polynomial function.

For example, an Extreme Fibonacci code-golf question was to print the first 1000 digits of Fib(1 billion) with a program that could actually finish in a reasonable amount of time. One Python2 answer from @AndersKaseorg used this code:

## golfed for minimum source-code size in bytes, not readability or performance
a,b=0,1
i=1e9
while i:
 a,b=b,a+b;i-=1
 if a>>3360:a/=10;b/=10
print a
# prints first 1000 digits of Fib(1**9),
# plus some (possibly wrong) extra trailing digits, as allowed by that question

which discards the low decimal digit of a and b if a > 2 ** 3360. That's a 1012 digit number, so this keeps at least 1012 significant digits. (And slightly more in b which is the larger of the 2).

For your case, you only care about the magnitude, so basically 1 significant digit of precision. You probably only have to keep about 12 or 13 significant digits to make sure the leading significant digit is always correct and doesn't grow too soon or too late. Maybe fewer; you might get away with keeping your running totals small enough to fit into 32-bit integers. Although I think Python on 64-bit machines will use 64-bit integers.

For your case, you do care about how many total digits there are. So you need to count how many times you truncate.

For efficiency, you probably want to set a higher truncation threshold and divide by 100 or 1000 so you don't have to do as many divisions. You might also be able to bit-shift (divide by a power of 2) instead of dividing by 10, if handle your final threshold appropriately, but that's trickier.

(That Python version takes about 12.5 minutes on a 4.4GHz Skylake x86-64 with Arch Linux CPython 2.7. My hand-written asm answer on the same question takes about 70 seconds by using base-10^9 extended precision, making it efficient to discard groups of 9 decimal digits.)

You can also unroll the loop so you only check for truncation every few additions because compare isn't free either. Perhaps have an efficient main loop that stops one or two orders of magnitude below the target digit-count, then go one step at a time without further truncation.

Working out the implementation details is left as an exercise for the reader (because the most efficient solution for large n is the closed-form FP math version). Keeping the numbers small should make the run time scale linearly with n, instead of having each addition take longer as the numbers get larger.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.