Possible Duplicate:
Handling very large numbers in Python
I have a python function to generate fibonacci numbers:
def fib(n):
return ((1+math.sqrt(5))**n - (1-math.sqrt(5))**n)/(2**n*math.sqrt(5))
I can feed the fib function numbers up to 700, where it starts to
OverflowError: (34, 'Numerical result out of range')
Do I need to use a special type like long to get around this?
The problem is that you're using doubles to calculate the value and the doubles are overflowing. Doubles give exact solutions only to about the 85th Fibonacci number.
If you want a fast and accurate calculation you're better off using an algorithm based on a better recurrence relationship, and using python bignum integers.
In particular you can use:
fib(2*n) = fib(n)^2 + fib(n-1)^2
fib(2*n-1) = fib(n)*(2*fib(n-1)+fib(n))
Or the equivalent matrix exponentiation formula (excuse the ugly formatting)
[ F_n F_{n-1} ] [ 1 1 ] ^N
[ ] = [ ]
[ F_{n-1} F_{n-2} ] [ 1 0 ]
Both of these result in algorithms that require O(log(N)) calculations rather than O(N).
Here's a complete solution in pseudo-code
If you do want to perform your calculations using doubles and the explicit formulae, then the formulae can be tweaked to give something faster that doesn't overflow until about the 1500th fibonacci number, and remains the same accuracy as your version. IIRC it is:
def fib(n):
return round( ((1+math.sqrt(5))/2)**n / math.sqrt(5) )
It's easy to isolate the error
>>> (1+math.sqrt(5))**700
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
OverflowError: (34, 'Numerical result out of range')
This method won't work well as floating point numbers don't have enough precision
for example, here
>>> (1+math.sqrt(5))**600
1.024664165563927e+306
you are working with only the first 15 or so digits. The remaining 291 will be treated as zeros when you do any arithmetic
See wikipedia for more information about accuracy problems with floating point numbers
You can always try this approach:
def fib(n, memo={0:0, 1:1}):
if n not in memo:
memo[n] = fib(n-1) + fib(n-2)
return memo[n]
print fib(800)
Output:
69283081864224717136290077681328518273399124385204820718966040597691435587278383112277161967532530675374170857404743017623467220361778016172106855838975759985190398725
If you actually want to use that algorithm, and you want to work beyond the limits of the built-in float, then yes, you need a different type.
If all you want is to get an approximate answer instead of an exception, that's easy; you can get infinite range out of the box. But if you also want to eliminate the rounding errors, you can't have infinite precision (that would take infinite time/space), so you have to know how to work out the precision needed for your range of inputs. (I'll leave that as an exercise for the reader.)
The standard library type decimal.Decimal may be all you need. It provides arbitrary-precision fixed- or floating-point decimal arithmetic according to the IEEE-854 standard. There are many cases for which it's unusable because it doesn't provide enough mathematical functions, but you only need basic arithmetic and sqrt, which are just fine. It can also be slow for huge numbers, but if you just want to calculate fib on a few three-digit numbers it's more than sufficient.
When Decimal is insufficient, there are a number of third-party modules, usually wrapping industry-standard C libraries like gmp/mpfr, such as bigfloat.
Here's how to get infinite range, but with rounding errors on roughly the same scale as the built-in float:
>>> s5 = decimal.Decimal(5).sqrt()
>>> def fib(n):
... return ((1+s5)**n - (1-s5)**n)/(2**n*s5)
>>> fib(800)
Decimal('6.928308186422471713629008226E+166')
>>> int(fib(800))
69283081864224717136290082260000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000L
>>> s5 = bigfloat.sqrt(5)
>>> def fib(n):
... return ((1+s5)**n - (1-s5)**n)/(2**n*s5)
>>> fib(800)
BigFloat.exact('6.9283081864226567e+166', precision=53)
>>> int(fib(800))
69283081864226566841137772774650010139572747244991592044952506898599601083170460360533811597710072779197410943266632999194601974766803264653830633103719677469311107072L
But notice that neither of these are actually the answer you'd get if you did the math perfectly; you've lost 24 digits to rounding errors. (The reason the values are different is that bigfloat is rounding in base 2, decimal in base 10.)
To fix that, you need more precision. All libraries provide some way to change the precision; bigfloat has more convenient options than most, but none are too onerous:
>>> decimal.getcontext().prec = 300
>>> s5 = decimal.Decimal(5).sqrt()
>>> def fib(n):
... return ((1+s5)**n - (1-s5)**n)/(2**n*s5)
>>> fib(800)
69283081864224717136290077681328518273399124385204820718966040597691435587278383112277161967532530675374170857404743017623467220361778016172106855838975759985190398725.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000048
>>> def fibp(n, p):
... with bigfloat.precision(p):
... s5 = bigfloat.sqrt(5)
... return ((1+s5)**n - (1-s5)**n)/(2**n*s5)
>>> fibp(800, 125)
BigFloat.exact('6.92830818642247171362900776814484912138e+166', precision=125)
>>> int(fibp(800, 125))
69283081864224717136290077681448491213794574774712670552070914552025662674717073354503451578576268674564384721027806323979200718479461097490537109958812524476157132800L
