Python: "long int too large to convert to float" when calculating pi

Question

I get this error when using a python script that calculates pi using the Gauss-Legendre algorithm. You can only use up to 1024 iterations before getting this:

    C:\Users\myUsernameHere>python Desktop/piWriter.py
    End iteration: 1025
    Traceback (most recent call last):
      File "Desktop/piWriter.py", line 15, in <module>
        vars()['t' + str(sub)] = vars()['t' + str(i)] - vars()['p' + str(i)] * math.
    pow((vars()['a' + str(i)] - vars()['a' + str(sub)]), 2)
    OverflowError: long int too large to convert to float

Here is my code:

import math

a0 = 1
b0 = 1/math.sqrt(2)
t0 = .25
p0 = 1

finalIter = input('End iteration: ')
finalIter = int(finalIter)

for i in range(0, finalIter):
        sub = i + 1
        vars()['a' + str(sub)] = (vars()['a' + str(i)] + vars()['b' + str(i)])/ 2
        vars()['b' + str(sub)] = math.sqrt((vars()['a' + str(i)] * vars()['b' + str(i)]))
        vars()['t' + str(sub)] = vars()['t' + str(i)] - vars()['p' + str(i)] * math.pow((vars()['a' + str(i)] - vars()['a' + str(sub)]), 2)
        vars()['p' + str(sub)] = 2 * vars()['p' + str(i)]
        n = i

pi = math.pow((vars()['a' + str(n)] + vars()['b' + str(n)]), 2) / (4 * vars()['t' + str(n)])
print(pi)

Ideally, I want to be able to plug in a very large number as the iteration value and come back a while later to see the result.

Any help appreciated! Thanks!

Solution 1

Floats can only represent numbers up to sys.float_info.max, or 1.7976931348623157e+308. Once you have an int with more than 308 digits (or so), you are stuck. Your iteration fails when p1024 has 309 digits:

179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137216L

You'll have to find a different algorithm for pi, one that doesn't require such large values.

Actually, you'll have to be careful with floats all around, since they are only approximations. If you modify your program to print the successive approximations of pi, it looks like this:

2.914213562373094923430016933707520365715026855468750000000000
3.140579250522168575088244324433617293834686279296875000000000
3.141592646213542838751209274050779640674591064453125000000000
3.141592653589794004176383168669417500495910644531250000000000
3.141592653589794004176383168669417500495910644531250000000000
3.141592653589794004176383168669417500495910644531250000000000
3.141592653589794004176383168669417500495910644531250000000000

In other words, after only 4 iterations, your approximation has stopped getting better. This is due to inaccuracies in the floats you are using, perhaps starting with 1/math.sqrt(2). Computing many digits of pi requires a very careful understanding of the numeric representation.

Solution 2

As noted in previous answer, the float type has an upper bound on number size. In typical implementations, sys.float_info.max is 1.7976931348623157e+308, which reflects the use of 10 bits plus sign for the exponent field in a 64-bit floating point number. (Note that 1024*math.log(2)/math.log(10) is about 308.2547155599.)

You can add another half dozen decades to the exponent size by using the Decimal number type. Here is an example (snipped from an ipython interpreter session):

In [48]: import decimal, math    
In [49]: g=decimal.Decimal('1e12345')    
In [50]: g.sqrt()
Out[50]: Decimal('3.162277660168379331998893544E+6172')
In [51]: math.sqrt(g)
Out[51]: inf

This illustrates that decimal's sqrt() function performs correctly with larger numbers than does math.sqrt().