5
\$\begingroup\$

I have implemented a calculation of Pi in Python, just for the lulz. But now I wonder if I could improve this more. Could you help me out? I do not program a lot in Python.

from time import time
from decimal import *

def faculty(m):
    k=m
    n=1
    while(n<m):
        k=k*n
        n+=1
    if(k>0):
        return k
    else:
        return 1

def ramanujan(b, n):
    return b+Decimal(faculty(4*n)*(1103+26390*n))/Decimal((faculty(n)**4)*396**(4*n))

def chudnovsky(b, n):
    return b+(((-1)**n)*(faculty(6*n))*(13591409+(545140134*n)))/(faculty(Decimal(3)*n)*(faculty(n)**3)*640320**(3*n+(Decimal(3)/Decimal(2))))

def calc(x, a, func):
    b=Decimal(0)
    n=Decimal(0)
    while(n<x):
        b = func(b, n)
        n=n+1
    return ((a*b)**(-1))

def calcrama():
    return calc(20, Decimal((2*Decimal(2).sqrt())/9801), ramanujan)

def calcchud():
    return calc(20, 12, chudnovsky)

def save(name, func):
    fobj = open(name, "w")
    t = time()
    pi = func()
    t = time() - t
    fobj.write("Time: "+str(t)+"\nPi: "+str(pi))
    fobj.close()

getcontext().prec = 1000

save("rama.txt", calcrama)
save("chud.txt", calcchud)
\$\endgroup\$
7
\$\begingroup\$

I find that your program is obfuscated and therefore incomprehensible to anyone except for perhaps an expert in those calculations methods. You should pick more informative variable names than a, b, k, m, n, and x.

faculty(m) should be renamed factorial(m). Ideally, you wouldn't need it — see below.

If you could verify that the denominator is always a divisor of the numerator, you should use // integer division.

Use more whitespace for readability. Here, I would break with Python guidelines and be more generous with whitespace than usual.

Verbose comments, especially docstrings, would also be appreciated for complicated mathematics:

def chudnovsky(b, n):
    """
    Based on the identity

                    n
      1     inf (-1)  (6n)! (13591409 + 545140134n)
    ----- = SUM -----------------------------------
    12 PI   n=0              3       3n + 1.5
                   (3n)! (n!)  640320

    Returns the nth term of the sum.
    """
    return b + (
        ((-1)**n) * (factorial(6 * n)) * 
        (13591409 + (545140134 * n))
    ) // (
        factorial(Decimal(3) * n) *
        (factorial(n)**3) *
        640320**(3 * n + (Decimal(3) / Decimal(2)))
    )

Many parts of that expression would be more efficiently computed based on their counterparts in the previous term, instead of starting from scratch. You might want to build an generator that caches those intermediate results.

Furthermore, you have very large numbers in both your numerators and denominators, especially factorials. You can help keep those small by cancelling out the factorials in the numerator and denominator as soon as possible.

def ChudnovskyTerms():
    """
    Based on the identity

                      n
      1    inf 12 (-1)  (6n)! (13591409 + 545140134n)
    ---- = SUM --------------------------------------
     PI    n=0                3       3n + 1.5
                    (3n)! (n!)  640320

    yields successive terms of the sum.
    """
    term = 0
    sign = 1
    factorial_6_3 = 1         # = factorial(6 * n) / factorial(3 * n)
    numerator_sum = 13591409
    n_factorial_cubed = 1
    denominator = 640320**1.5

    while True:
        yield 12 * sign * factorial_6_3 * …

        term += 1
        sign *= -1
        for i in range(6):
            factorial_6_3 *= (6 * term - i)
        for i in range(3):
            factorial_6_3 //= (3 * term - i)
        numerator_sum += 545140134
        # etc.
\$\endgroup\$
4
  • \$\begingroup\$ So wrapping it in a class increases the speed significantly? \$\endgroup\$
    – NaCl
    Jan 20 '14 at 23:37
  • 2
    \$\begingroup\$ Keeping intermediate results rather than discarding them increases the speed. Packaging the code in a generator gives it a reasonable interface, and the overhead wouldn't be worse than what you had before. \$\endgroup\$ Jan 20 '14 at 23:45
  • \$\begingroup\$ self.term looks wrong to me: what's self? \$\endgroup\$ Jan 22 '14 at 16:02
  • \$\begingroup\$ @GarethRees Thanks for spotting that. It was left over from an old version before I rewrote it as a generator. \$\endgroup\$ Jan 22 '14 at 16:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.