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!

Answer 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.

Answer 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().

Answer 3

As noted above, getting lots of digits is going to be tricky, but looking at all those vars hurts my eyes. So here's a version of your code after (1) replacing your use of vars with dictionaries, and (2) using ** instead of the math functions:

a, b, t, p = {}, {}, {}, {}
a[0] = 1
b[0] = 2**-0.5
t[0] = 0.25
p[0] = 1

finalIter = 4

for i in range(finalIter):
    sub = i + 1
    a[sub] = (a[i] + b[i]) / 2
    b[sub] = (a[i] * b[i])**0.5
    t[sub] = t[i] - p[i] * (a[i] - a[sub])**2
    p[sub] = 2 * p[i]
    n = i

pi_approx = (a[n] + b[n])**2 / (4 * t[n])

Instead of playing games with vars, I've used dictionaries to store the values (the link there is to the official Python tutorial) which makes your code much more readable. You can probably even see an optimization or two now.

As noted in the comments, you really don't need to store all the values, only the last, but I think it's more important that you see how to do things without dynamically creating variables. Instead of a dict, you could also have simply appended the values to a list, but lists are always zero-indexed and you can't easily "skip ahead" and set values at arbitrary indices. That can occasionally be confusing when working with algorithms, so let's start simple.

Anyway, the above gives me

>>> print(pi_approx)
3.141592653589794
>>> print(pi_approx-math.pi)
8.881784197001252e-16

Answer 4

A simple solution is to install and use the arbitrary-precisionmpmath module which now supports Python 3. However, since I completely agree with DSM that your use ofvars()to create variables on the fly is an undesirable way to implement the algorithm, I've based my answer on his rewrite of your code and [trivially] modified it to make use ofmpmath to do the calculations.

If you insist on usingvars(), you could probably do something similar -- although I suspect it might be more difficult and the result would definitely harder to read, understand, and modify.

from mpmath import mpf  # arbitrary-precision float type

a, b, t, p = {}, {}, {}, {}
a[0] = mpf(1)
b[0] = mpf(2**-0.5)
t[0] = mpf(0.25)
p[0] = mpf(1)

finalIter = 10000

for i in range(finalIter):
    sub = i + 1
    a[sub] = (a[i] + b[i]) / 2
    b[sub] = (a[i] * b[i])**0.5
    t[sub] = t[i] - p[i] * (a[i] - a[sub])**2
    p[sub] = 2 * p[i]
    n = i

pi_approx = (a[n] + b[n])**2 / (4 * t[n])
print(pi_approx)  # 3.14159265358979