6
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Wow, does python not have a square root function? \$\endgroup\$ – ChiefTwoPencils Feb 8 '15 at 17:31
  • \$\begingroup\$ echo 'scale=1000; 4*a(1)' | bc -l \$\endgroup\$ – Jan-Philip Gehrcke Feb 8 '15 at 17:33
  • 2
    \$\begingroup\$ You should call this program pypi or something. =) \$\endgroup\$ – 200_success Feb 8 '15 at 18:06
  • \$\begingroup\$ @ChiefTwoPencils I think it does, but it doesn't work for long integers. If I replace calculate_square_root with **0.5 I get this error message: "OverflowError: long int too large to convert to float". \$\endgroup\$ – Uri Feb 8 '15 at 18:56
  • 1
    \$\begingroup\$ Python 2 or 3? It is important because division is different. \$\endgroup\$ – Caridorc Feb 8 '15 at 21:47
4
\$\begingroup\$

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.
\$\endgroup\$
  • \$\begingroup\$ No. Most of this program's run time is spent inside bigint arithmetic functions which are implemented in C anyway. \$\endgroup\$ – Janne Karila Feb 9 '15 at 13:05
  • \$\begingroup\$ @JanneKarila: He does a sqrt approximation iteratively in python, definitely not in C \$\endgroup\$ – Marcus Müller Feb 9 '15 at 13:38
  • \$\begingroup\$ Yes but the number of iterations is rather small: 210 calls to calculate_next_square_root. \$\endgroup\$ – Janne Karila Feb 9 '15 at 14:21
  • \$\begingroup\$ With small integers a call to calculate_next_square_root takes 300 ns, but with OP's numbers 3 ms. From that I'd estimate Python's overhead to be 0.01 %. \$\endgroup\$ – Janne Karila Feb 9 '15 at 14:33
  • \$\begingroup\$ ((10 ** (digits + add)) * ((a + b) ** 2)) / (4 * t) has ten CPython PyObjects. For each one of these, at least a type lookup, a lookup for the correct operator and the value lookup. That's three times the indirection operations you'd need to do if this was a statically typed compiled language, not even counting the fact that Python somehow has to figure out what these strings mean. \$\endgroup\$ – Marcus Müller Feb 9 '15 at 15:08
3
\$\begingroup\$

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/

\$\endgroup\$
  • \$\begingroup\$ Which version of Python is this program for? I tried with Python 2.7.9 and I got this error message: "ImportError: No module named gmpy2". \$\endgroup\$ – Uri Feb 8 '15 at 18:21
  • \$\begingroup\$ @Uri You will have to download it from here:pypi.python.org/pypi/gmpy2 \$\endgroup\$ – Mohit Bhasi Feb 8 '15 at 18:24
  • \$\begingroup\$ I tried to install gmpy2 but it failed. Your answer is nice, but I'm looking for a pure integer operating program which doesn't have to import modules such as math and gmpy2. \$\endgroup\$ – Uri Feb 8 '15 at 18:50
  • \$\begingroup\$ @Uri: that's important information! I still don't get your use case, though. Purely integer sounds like you want to adapt this algorithm to something without floating point instructions (which, honestly, makes little sense when calculating pi), but you're using python. \$\endgroup\$ – Marcus Müller Feb 8 '15 at 18:53
  • 1
    \$\begingroup\$ @MohitBhasi I used your program without importing math and mpz, I replaced math.log10(C3_OVER_24/6/2/6) with the constant 14.1816474627 and (10005*one_squared).sqrt() with calculate_square_root(number=(10005*one_squared), digits=digits + 2, add=0), added my square root functions, removed mpz from the code (replaced it with "") and it took 20.65 seconds to calculate the first 100,000 digits of pi (Last 5 digits 24646 OK), which is much faster than my own program. \$\endgroup\$ – Uri Feb 8 '15 at 21:18
2
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Why do you prefer calling calculate_next_pi(a, b, t, digits, add) instead of calculate_next_pi(a=a, b=b, t=t, digits=digits, add=add)? I always prefer to call functions the second way (with parameter names), to avoid calling the wrong function or if the function parameters change I will get an exception. It's possible in Python but not in other languages like JavaScript (I also program in JavaScript). \$\endgroup\$ – Uri Feb 10 '15 at 5:39
  • \$\begingroup\$ My way is less verbose. \$\endgroup\$ – Caridorc Feb 10 '15 at 12:38
2
\$\begingroup\$

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.

\$\endgroup\$
  • 1
    \$\begingroup\$ Using b instead of 1 as an initial value is a big win; this is the 0th-order estimate. I wonder if the first-order estimate (a+b)/2 would be better, especially since it's already computed. \$\endgroup\$ – Charles Mar 12 '15 at 18:41

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.