Prevent Rounding to Zero in Python

Question

I have a program meant to approximate pi using the Chudnovsky Algorithm, but a term in my equation that is very small keeps being rounded to zero.

Here is the algorithm:

import math
from decimal import *
getcontext().prec = 100

pi = Decimal(0.0)
C = Decimal(12/(math.sqrt(640320**3)))
k = 0
x = Decimal(0.0)
result = Decimal(0.0)
sign = 1
while k<10:
    r = Decimal(math.factorial(6*k)/((math.factorial(k)**3)*math.factorial(3*k)))
    s = Decimal((13591409+545140134*k)/((640320**3)**k))
    x += Decimal(sign*r*s)
    sign = sign*(-1)
    k += 1
result = Decimal(C*x)
pi = Decimal(1/result)


print Decimal(pi)

The equations may be clearer without the "decimal" terms.

import math

pi = 0.0
C = 12/(math.sqrt(640320**3))
k = 0
x = 0.0
result = 0.0
sign = 1
while k<10:
    r = math.factorial(6*k)/((math.factorial(k)**3)*math.factorial(3*k))
    s = (13591409+545140134*k)/((640320**3)**k)
    x += sign*r*s
    sign = sign*(-1)
    k += 1
result = C*x
pi = 1/result

print pi

The issue is with the "s" variable. For k>0, it always comes to zero. e.g. at k=1, s should equal about 2.1e-9, but instead it is just zero. Because of this all of my terms after the first =0. How do I get python to calculate the exact value of s instead of rounding it down to 0?

Answer 1

Try:

s = Decimal((13591409+545140134*k)) / Decimal(((640320**3)**k))

The arithmetic you're doing is native python - by allowing the Decimal object to perform your division, you should eliminate your error.

You can do the same, then, when computing r.

Answer 2

A couple of comments.

If you are using Python 2.x, the / returns an integer result. If you want a Decimal result, you convert at least one side to Decimal first.

math.sqrt() only return ~16 digits of precision. Since your value for C will only be accurate to ~16 digits, your final result will only be accurate to 16 digits.

Answer 3

If you're doing maths in Python 2.x, you should probably be putting this line into every module:

from __future__ import division

This changes the meaning of the division operator so that it will return a floating point number if needed to give a (closer to) precise answer. The historical behaviour is for x / y to return an int if both x and y are ints, which usually forces the answer to be rounded down.

Returning a float if necessary is generally regarded as a better way to handle division in a language like Python where duck typing is encouraged, since you can just worry about the value of your numbers rather than getting different behaviour for different types.

In Python 3 this is in fact the default, but since old programs relied on the historical behaviour of the division operator it was felt the change was too backwards-incompatible to be made in Python 2. This is why you have to explicitly turn it on with the __future__ import. I would recommend always adding that import in any module that might be doing any mathematics (or just any module at all, if you can be bothered). You'll almost never be upset that it's there, but not having it there has been the cause of a number of obscure bugs I've had to chase.

Answer 4

I feel that the problem with 's' is that all terms are integers, thus you are doing integer maths. A very simple workaround, would be to use 3.0 in the denominator. It only takes one float in the calculation to get a float returned.