Pi-calculating program

Question

I saw this question and answer about calculating pi on Stack Overflow, and I decided to write my own program for calculating pi. I used Python and only integers (I didn't want to use floating point numbers), and used the Gauss–Legendre algorithm because it was the simplest to implement (I considered using the Borwein's algorithm, but I didn't want to calculate third roots of numbers, and the Chudnovsky algorithm seemed a little complicated, although maybe I'll give it a try). I would like to know how efficient are the algorithms above, in \$\mathcal{O}(n)\$ (n is the number of digits to calculate), and how efficient is my program?

# Calculate next square root approximation of the number.
def calculate_next_square_root(number, square_root):
    next_square_root = ((number / square_root) + square_root) / 2
    return next_square_root

# Calculate the square root of the number.
def calculate_square_root(number, digits, add):
    # Start with 1, followed by (digits + add) zeros.
    square_root = 1 * (10 ** (digits + add))
    # Calculate next square root approximation of the number.
    next_square_root = calculate_next_square_root(number=number, square_root=square_root)
    while (next_square_root != square_root):
        # Replace square root with next square root.
        square_root = next_square_root
        # Calculate next square root approximation of the number.
        next_square_root = calculate_next_square_root(number=number, square_root=square_root)
    return square_root

# Calculate next pi approximation.
def calculate_next_pi(a, b, t, digits, add):
    next_pi = ((10 ** (digits + add)) * ((a + b) ** 2)) / (4 * t)
    return next_pi

# Calculate pi to 5,000 digits.
digits = 5000
add = 500
a = 10 ** (digits + add)
b = calculate_square_root(number=(10 ** ((digits + add) * 2)) / 2, digits=digits, add=add)
t = (10 ** ((digits + add) * 2)) / 4
p = 1
pi = -1 # pi must be different than next_pi
next_pi = calculate_next_pi(a=a, b=b, t=t, digits=digits, add=add)
n = 0
while (next_pi != pi):
    pi = next_pi
    next_a = (a + b) / 2
    next_b = calculate_square_root(number=a * b, digits=digits, add=add)
    next_t = t - (p * ((a - next_a) ** 2))
    next_p = 2 * p
    a = next_a
    b = next_b
    t = next_t
    p = next_p
    next_pi = calculate_next_pi(a=a, b=b, t=t, digits=digits, add=add)
    n += 1

# Remove the last 500 digits.
pi /= (10 ** add)

# Print the results.
print pi
print n

My program takes 52 lines of code (including comments), and I'm not interested in programs which take too many lines to implement, up to 100 or 200 lines are fine. I checked and it takes about one second to calculate 5,000 digits of pi (with add = 500), and a little less than two minutes to calculate 50,000 digits (with add = 5000). With 5,000 digits it took 14 iterations, and with 50,000 digits 17 iterations, but I didn't count how many iterations it took to calculate the square roots. Is my program efficiently using the Gauss–Legendre algorithm, or can it be made more efficient?

By the way, add doesn't have to be 500, even 4 or 5 is enough. But I wanted to make sure all the digits of pi are correct (and they are). With add = 0 the last 2 or 3 digits are not correct.

I took the square root program from a program I wrote in 2013, but I don't know the name of the algorithm I used and whether it's the most efficient algorithm or not.

Update: I remind you my question: I would like to know how efficient are the algorithms above, in \$\mathcal{O}(n)\$ (n is the number of digits to calculate), and how efficient is my program? I prefer not to import modules, but use pure Python (either 2 or 3) without any import. And of course if we start with a number closer to the square root (such as a or b) then it will take fewer iterations to calculate the square root.

Answer 1

This is a bit hard to answer, but let's try nevertheless:

1. From a computational point of view, using a scripting language like python for iterative algorithms doesn't make the least sense -- your PC is going to be far more occupied by translating your python code to executable structures than it will be with executing the core algorithm. For example, every number in python (assuming a CPython implementation of the runtime environment) is an object -- which makes number-number multiplication far more expensive than the actual multiplication. Now, I assume your question is "if I did this in a compiled language, with fixed length numerical types, would this algorithm perform optimally"?
2. approximating a square root in software is bound to be horrendously expensive compared to just calculating the floating point square root; modern general purpose CPUs even have instructions that can take the square root of several numbers at once, in just a few cycles!
3. I just have to stress this once more: "efficient" doesn't really make sense in this context, without actually specifying your measurement for efficiency. If your measurement is "how long does it take to have N significant digits", most probably even the most naive C implementation would beat the hell out of your Python program -- just because python is such a capable scripting language, which requires a lot of indirection. You could specify things like "minimum number of multiplications plus additions", but you'd totally neglect possibly different cost for these operations, and, much worse, memory movement complexity etc. Now, if you're actually into this kind of discussion, think about what kind of architecture you actually try to optimize for: A single threaded computer with uniform multiply, add, memory address loading and storage costs has never really existed; the norm today are hyperthreaded multicore processors with a limited number of SIMD registers and multilevel caches, which can make algorithms which "look good on paper" perform terribly because they just misuse the capabilities of the architecture, whilst stupid implementation might work exceptionally well because of possible compile- and execution-time optimization.

Answer 2

I think the most efficient method is the Chudnovsky algoritm (100 million digits of Pi, in under 10 minutes!) To know the math behind it:

"""
Python3 program to calculate Pi using python long integers, BINARY
splitting and the Chudnovsky algorithm

"""

import math
from gmpy2 import mpz
from time import time

def pi_chudnovsky_bs(digits):
    """
    Compute int(pi * 10**digits)

    This is done using Chudnovsky's series with BINARY splitting
    """
    C = 640320
    C3_OVER_24 = C**3 // 24
    def bs(a, b):
        """
        Computes the terms for binary splitting the Chudnovsky infinite series

        a(a) = +/- (13591409 + 545140134*a)
        p(a) = (6*a-5)*(2*a-1)*(6*a-1)
        b(a) = 1
        q(a) = a*a*a*C3_OVER_24

        returns P(a,b), Q(a,b) and T(a,b)
        """
        if b - a == 1:
            # Directly compute P(a,a+1), Q(a,a+1) and T(a,a+1)
            if a == 0:
                Pab = Qab = mpz(1)
            else:
                Pab = mpz((6*a-5)*(2*a-1)*(6*a-1))
                Qab = mpz(a*a*a*C3_OVER_24)
            Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a)
            if a & 1:
                Tab = -Tab
        else:
            # Recursively compute P(a,b), Q(a,b) and T(a,b)
            # m is the midpoint of a and b
            m = (a + b) // 2
            # Recursively calculate P(a,m), Q(a,m) and T(a,m)
            Pam, Qam, Tam = bs(a, m)
            # Recursively calculate P(m,b), Q(m,b) and T(m,b)
            Pmb, Qmb, Tmb = bs(m, b)
            # Now combine
            Pab = Pam * Pmb
            Qab = Qam * Qmb
            Tab = Qmb * Tam + Pam * Tmb
        return Pab, Qab, Tab
    # how many terms to compute
    DIGITS_PER_TERM = math.log10(C3_OVER_24/6/2/6)
    N = int(digits/DIGITS_PER_TERM + 1)
    # Calclate P(0,N) and Q(0,N)
    P, Q, T = bs(0, N)
    one_squared = mpz(10)**(2*digits)
    sqrtC = (10005*one_squared).sqrt()
    return (Q*426880*sqrtC) // T

# The last 5 digits or pi for various numbers of digits
check_digits = {
        100 : 70679,
       1000 :  1989,
      10000 : 75678,
     100000 : 24646,
    1000000 : 58151,
   10000000 : 55897,
}

if __name__ == "__main__":
    digits = 100
    pi = pi_chudnovsky_bs(digits)
    print(pi)
    #raise SystemExit
    for log10_digits in range(1,9):
        digits = 10**log10_digits
        start =time()
        pi = pi_chudnovsky_bs(digits)
        print("chudnovsky_gmpy_mpz_bs: digits",digits,"time",time()-start)
        if digits in check_digits:
            last_five_digits = pi % 100000
            if check_digits[digits] == last_five_digits:
                print("Last 5 digits %05d OK" % last_five_digits)
            else:
                print("Last 5 digits %05d wrong should be %05d" % (last_five_digits, check_digits[digits]))

For further details:http://www.craig-wood.com/nick/articles/pi-chudnovsky/

Answer 3

Others have talked about performance, now some style/readibility improvements:

---

a = 10 ** (digits + add)
b = calculate_square_root(number=(10 ** ((digits + add) * 2)) / 2, digits=digits, add=add)
t = (10 ** ((digits + add) * 2)) / 4
p = 1

You should avoid one letter variable names, prefer longer more descriptive names.

---

next_pi = calculate_next_pi(a=a, b=b, t=t, digits=digits, add=add)

You call functions in a very verbose way, you should prefer:

next_pi = calculate_next_pi(a, b, t, digits, add)

---

# Calculate next square root approximation of the number.
def calculate_next_square_root(number, square_root):

Your comment just repeats the obvious instead say which algorithm you chose and maybe also why you chose that algorithm over another.

---

def calculate_next_square_root(number, square_root):
    next_square_root = ((number / square_root) + square_root) / 2
    return next_square_root

def calculate_next_pi(a, b, t, digits, add):
    next_pi = ((10 ** (digits + add)) * ((a + b) ** 2)) / (4 * t)
    return next_pi

Omit needless words:

def calculate_next_square_root(number, square_root):
    return ((number / square_root) + square_root) / 2

def calculate_next_pi(a, b, t, digits, add):
    return ((10 ** (digits + add)) * ((a + b) ** 2)) / (4 * t)

---

digits = 5000
add = 500
a = 10 ** (digits + add)
b = calculate_square_root(number=(10 ** ((digits + add) * 2)) / 2, digits=digits, add=add)
t = (10 ** ((digits + add) * 2)) / 4
p = 1
pi = -1 # pi must be different than next_pi
next_pi = calculate_next_pi(a=a, b=b, t=t, digits=digits, add=add)
n = 0
while (next_pi != pi):
    pi = next_pi
    next_a = (a + b) / 2
    next_b = calculate_square_root(number=a * b, digits=digits, add=add)
    next_t = t - (p * ((a - next_a) ** 2))
    next_p = 2 * p
    a = next_a
    b = next_b
    t = next_t
    p = next_p
    next_pi = calculate_next_pi(a=a, b=b, t=t, digits=digits, add=add)
    n += 1

This big fat 'top-level' code should be a function, also because Python code runs faster in functions.

---

Mathematical code of this kind really benefits from automated testing, a simple example of automatic testing is:

import doctest

def add(a,b):
    """
    >>> add(3,9)
    12
    """

doctest.testmod()

This allows you to optimize without fear of breaking things.

Answer 4

I didn't count how many iterations it took to calculate the square roots.

That's easy to do because each iteration involves a function call. Simply profile the program by issuing this command from shell:

python -m cProfile filename.py

A report similar to this should appear. From the ncalls column you see calculate_next_square_root is called 210 times.

         242 function calls in 0.858 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.034    0.034    0.858    0.858 cr1a.py:2(<module>)
      210    0.682    0.003    0.682    0.003 cr1a.py:2(calculate_next_square_root)
       15    0.137    0.009    0.137    0.009 cr1a.py:20(calculate_next_pi)
       15    0.004    0.000    0.687    0.046 cr1a.py:7(calculate_square_root)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

I took the square root program from a program I wrote in 2013, but I don't know the name of the algorithm I used and whether it's the most efficient algorithm or not.

It is Newton's method. You could reduce the number of iterations with a better initial guess for the square root (instead of 1). You are computing the square root of a*b where a and b converge towards the same value. Thus either a or b would be a good initial guess for the square root. Using b I see the number of calls go down to 111.