Improving Implementation of Binet's Formula for finding nth Fibonacci Number in Python 3

Question

I tried to Implement Binet's formula for finding nth Fibonacci Number in Python 3.

def nth_fib(n):
    # this function returns fibonacci number of
    # the given term by using Binet's Formula
    sq5 = 5 ** 0.5
    phi = (sq5 + 1) / 2
    fib = (phi ** n) - (-phi ** -n)
    fib //= sq5
    return int(fib)

The problem with this implementation:

The biggest value it can handle is 1474. Passing values greater than 1474 throws the following Exception.

OverflowError: (34, 'Numerical result out of range')

How can this solution be improved for handling inputs ranging from 0 to 10^14?

Reference:

Binet's Formula for the nth Fibonacci number

Is it related to this? https://stackoverflow.com/questions/61874703/finding-last-digit-of-sum-from-m-to-n-fibonacci-numbers-0-%e2%89%a4-%e2%89%a4-%e2%89%a4-1014 — , May 18 '20 at 17:28
The Binet formula won't help if you need the exact Fibonacci number. However, there are duplication formulas. You can find them by applying [binary exponentiation](https://en.wikipedia.org/wiki/Exponentiation_by_squaring) to the [matrix form](https://en.wikipedia.org/wiki/Fibonacci_number#Matrix_form) of Fibonacci numbers. — , May 18 '20 at 17:30
Then I will ask the same question: where does this question come from? As to why Binet won't help: because the square root is an irrational number. If you want the exact result you will need multiprecision floating point. It will probably be faster by other means. Anyway, fib(10^14) is so large that you will have to focus on only part of the information: first digits (then floating point may help), last digits (then it's a form of modular arithmetic), sum of digits, etc. — , May 18 '20 at 17:45

hobbes · Answer 1 · 2020-05-19T17:00:18.677

This is happening because arbitrarily large exponentiation only works with integers, which are bignums in Python. Floating point numbers will fail. For example:

In [41]: phi
Out[41]: 1.618033988749895

In [42]: 5 ** 1500
Out[42]: 28510609648967058593679017274152865451280965073663823693385003532994042703726535281750010939152320351504192537189883337948877940498568886988842742507258196646578577135043859507339978111500571726845535306970880115202339030933389586900213992268035185770649319797269196725831118636035211367342502592161612681404558896878205505259742673921998666848316296574456143285153407461693074529608060405705703190247031916733545429301523565202628619442784043773875799299799772062596279270685668750358350581239751392647377917727924073955752619811973924353072146897222054396284190793435454619462166959138549077025548151961129557730113226497053327025918024691450322204632795881761117317264715060152457060422911440809597657134113164654343933125576083446389585308532864118204843115878436344284086952443434298108182889069338971572783051504615283483170635029160778619107133456847839866260715887917144004772675646444499010890878045793828781976559446412621993167117009741097351499347086624666372905178820086046962818676294533224769602031134496655795373953878879547119140625

In [43]: phi ** 1500
---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
<ipython-input-43-38afd4fed496> in <module>()
----> 1 phi ** 1500

OverflowError: (34, 'Numerical result out of range')

The solution is to use the Decimal class, which can handle arbitrary precision floating point operations, including exponentiation:

In [47]: from decimal import *

In [48]: getcontext().power(Decimal(phi), Decimal(1500))
Out[48]: Decimal('3.030123816655090678595267922E+313')

With that in mind, a rewritten nth_fib function might look like this (I've consolidated some of the math and removed the integer division to avoid a type error):

from decimal import *

def nth_fib(n):
    # Tweak your Decimal context depending on your needs.
    ctx = Context(prec=60, rounding=ROUND_HALF_EVEN)
    sq5 = Decimal(5 ** 0.5)
    phi = Decimal(sq5 + 1) / 2
    fib = (ctx.power(phi, Decimal(n)) - ctx.power(-phi, -n)) / sq5
    return int(fib)

print(nth_fib(5))
print(nth_fib(1500))

Which should give output like:

5
13551125668563783707293396130000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

As noted in the comments, any operations on floats accumulate errors, the scale of which depends on the operations performed and the frequency thereof.

I hope this helps. Cheers!

*...Python's default exponentiation operator \*\* works on integers...*. I'm not sure what you meant to say there but the \*\* works with any numeric operands — President James K. Polk, May 18 '20 at 14:53
It was late when I wrote that. I meant that arbitrarily large exponentiation will only work with integers. I'll update my answer. — hobbes, May 18 '20 at 14:59
@GaganDeepSingh Sorry, I forgot the original `int()` cast. I've updated my answer. — hobbes, May 18 '20 at 16:01
@hobbes Numbers being returned by this function are incorrect. (Beginning from nth Term 77). :( — Gagan Deep Singh, May 18 '20 at 17:19
@GaganDeepSingh It's possible there are rounding errors or floating point precision errors. The general solution stands, though: if you want to perform arbitrarily large exponentiation with floats, you'll need to use the `Decimal` module. — hobbes, May 18 '20 at 17:58
Good answer. @GaganDeepSingh: Obviously as `n` increases you will need increasing precision. Figuring out what precision is needed as a function of n is probably not too hard as just a rough over-estimate will to. — President James K. Polk, May 19 '20 at 01:32

Improving Implementation of Binet's Formula for finding nth Fibonacci Number in Python 3

1 Answers1