CODEWITHSUNDEEP
pythonalgorithmmathfloating-point
user2998454
user2998454

Reputation: 155

Python miscalculating for large floats

My algorithm returns approximations for pi. It works by using the equation 

a_2n= math.sqrt(-2*(math.sqrt((n**2)-4*(a_n**2))-n)*n)/2

where a_n is the approx when n is a certain number and a_2n is the approx when n is twice that. It starts with n = 4 and a_n is 2, and applies the formula until n is sufficiently high. The higher n is the more accurate the calculation, only it suddenly stops converging to pi when n>=2^22

Here's the full code:

import math
def pi(n):
    value = 2
    round = 4
    loop = 2
    while round != n:
        value= a2n(value,round)
        round = round*2
        loop +=1
        print("n=2^", loop, " pi = ", value)
    return value


def a2n(an,n):
    return  math.sqrt(-2*(math.sqrt((n**2)-4*(an**2))-n)*n)/2


print(pi(2**25))

I'm confident the math is fine, so I'm thinking python is having trouble with the larger numbers. It goes from '3.141596553704' to '3.14167426502175' and gets worse from there.

Upvotes: 0

Views: 261

Answers (2)

Mark Dickinson
Mark Dickinson

Reputation: 30561

The iterative formula you're using is somewhat ill-conditioned: as the iteration progresses, the quantity n becomes very much larger than an. So in the expression

math.sqrt(n**2-4*an**2)-n

the result of the square root will be close to n, so the outer subtraction is a subtraction of two quantities that are nearly equal (in a relative sense). Now if you're computing with regular Python floats, which have 16-ish decimal digits of precision, that subtraction is going to be giving you a result that's accurate to only a handful of digits as the iteration progresses. See the Wikipedia page on loss of significance for a more general discussion of this problem.

Short story: to get d digits of pi using your formula as originally written, you're going to need to work with many more than d digits in the intermediate computations. With a little bit of work, you can show that you're going to need to work with a little more than 2d digits of precision internally in order to get d accurate digits of pi. And even then you have to be careful only to use exactly as many iterations as you need: too many iterations and the accuracy will be lost again, no matter how much intermediate precision you use.

But there's a much better alternative here than doubling the intermediate precision, and that's to rewrite your formula so that it avoids the loss of significance in the first place. If you multiply sqrt(n**2-4*an**2)-n by the conjugate expression sqrt(n**2-4*an**2)+n, you get -4*an**2. So the original difference can be rewritten as -4*an**2/(sqrt(n**2-4*an**2)+n). Plugging that into your original formula and simplifying a bit leads to an iteration step that looks like this:

def a2n(an, n):
    return an*math.sqrt(2*n/(math.sqrt(n*n-4*an*an)+n))

Mathematically, that's doing exactly the same calculation as your a2n function, but it's a whole lot better behaved from a numerical point of view.

If you use this iteration step in place of the original, you'll see much smaller rounding error, and you should be able to get up to 15 digits of accuracy just with Python floats. Indeed, running your code with this new iteration, I get a value of 3.1415926535897927 after 30 iterations, which is only off the best possible double-precision approximation to pi by a single ulp (unit in the last place).

To get more digits than this, we'll need to use the decimal module. Here's a snippet, based on your code, but using my suggested modified computation for the iterations, using the decimal module to get a value for pi that's accurate to 51 significant digits. It's using an internal precision of 55 significant digits, to allow for accumulation of rounding errors.

from decimal import getcontext

context = getcontext()
context.prec = 55    # Use 55 significant digits for all operations.
sqrt = context.sqrt  # Will work for both ints and Decimal objects,
                     # returning a Decimal result.

def step(an, n):
    return an*sqrt(2*n/(sqrt(n*n-4*an*an)+n)), n*2

def compute_pi(iterations):
    value, round = 2, 4
    for k in range(iterations):
        value, round = step(value, round)
    return value

pi_approx = compute_pi(100)
print("pi = {:.50f}".format(pi_approx))

And here's the result:

pi = 3.14159265358979323846264338327950288419716939937511

With the original formula, this would have needed a precision exceeding 100 for the intermediate calculations.

Upvotes: 4

Marcus Müller
Marcus Müller

Reputation: 36327

It's probably not python itself that has this problems, but at 2^22 you run into numerical problems with the types python uses internally.

For the 64 bit float that python uses, you get 52bit of mantissa; that limits your accuracy to around 16 (decimal) digits behind the comma.

EDIT: I don't really see how going back to Matlab helps -- Matlab internally usually also uses the same 64 bit float type by default.

Upvotes: 2

Related Questions