Finding The Value Of Pi


Everyone must have heard of Pi, the irrational number that never ends. It goes by this :

3.14159265...

I have been learning about Pi since I was in 8th grade. All the time in school, I was curious about to which Pi extends.

Now, I know the extend – Infinity :P. I dedicated my precious Christmas vacation for finding the useless value of Pi.

Why ? Curiosity. Simply, the curiosity made me want to do this. Besides, it’s always fun to do something interesting.

The source code is on GitHub. Please Star it. Math Enthusiasts, fork it. 😉

Findings

1 Million digits and 5 Million digits were successfully found. You can see them here.

The current World Record is > 1 Trillion digits. (I should get a super computer soon).

Equation

Craig Wood has made a Python script to calculate N digits of Pi. It’s a very complicated script and it took me 2 days to understand it.

Pi equation

Pi equation

What I did is, create the client side and a WebSocket server. The Pi calculation is ran through Python and inserted to the database.

It’s a Live process. Meaning, a person on the client can see the step where the Pi is being calculated. If we look at the Python script :

import math
import os
import time
import MySQLdb
import base64
import zlib
import sys

digits_to_calculate = int(sys.argv[1])

db = MySQLdb.connect(host=os.environ["OPENSHIFT_MYSQL_DB_HOST"],    # your host, usually localhost
                     user=os.environ["OPENSHIFT_MYSQL_DB_USERNAME"],         # your username
                     passwd=os.environ["OPENSHIFT_MYSQL_DB_PASSWORD"],  # your password
                     db="ws")
cur = db.cursor()

# The time when program start
START_TIME = time.time()

cur.execute("""
UPDATE `pi` SET `value` = %s WHERE `key_name` = %s
""", (START_TIME, "start"))

def output(name, content):
  cur.execute("UPDATE `pi` SET `value` = %s WHERE `key_name` = %s", (content, name))
  db.commit()

sq_i = 1
def sqrt(n, one):
  global sq_i
  """
  Return the square root of n as a fixed point number with the one
  passed in.  It uses a second order Newton-Raphson convergence.  This
  doubles the number of significant figures on each iteration.
  """
  # Use floating point arithmetic to make an initial guess
  floating_point_precision = 10**16
  output("status", "sqrt," + str(sq_i))
  
  n_float = float((n * floating_point_precision) // one) / floating_point_precision
  x = (int(floating_point_precision * math.sqrt(n_float)) * one) // floating_point_precision
  n_one = n * one
  while 1:
    sq_i+=1
    output("status", "sqrt," + str(sq_i))
    x_old = x
    x = (x + n_one // x) // 2
    if x == x_old:
      break
  return x

i = 0
Qab = 0
def pi_chudnovsky_bs(digits):
    global i
    """
    Compute int(pi * 10**digits)

    This is done using Chudnovsky's series with binary splitting
    """
    output("status", "init,1")
    one = 10**digits
    output("status", "init,2")
    sqrtC = sqrt(10005*one, one)
    
    output("status", "init,3")
    C = 640320
    C3_OVER_24 = C**3 // 24
    
    def bs(a, b):
        global i, Qab
        """
        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)
        """
        
        """
        Subin Siby <subinsb.com>
        """
        
        if i != 0 and i != 1:
          calculated_digits = str(digits - (digits/i) + 7)
          output("status", "pi," + str(i) + "," +calculated_digits + "," + str(digits_to_calculate))
          
        """
        Calculate
        """
        if b - a == 1:
            # Directly compute P(a,a+1), Q(a,a+1) and T(a,a+1)
            if a == 0:
                Pab = Qab = 1
            else:
                Pab = (6*a-5)*(2*a-1)*(6*a-1)
                Qab = 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
        
        i+=1
        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, S = bs(0, N)
    
    return Q*426880*sqrtC // S
    
pi = pi_chudnovsky_bs(digits_to_calculate)
# Compress Pi
pi = zlib.compress(str(pi))

# Save to DB
output("pi", pi)

db.close()

we can see the call of function output. This function updates the respective values in DB. The “status” key holds the current status of Pi calculation. It may have values like :

Status

Stage

init,[I] The initialization stage. Creation of variable etc.
sqrt,[I] The loop inside sqrt() function. I indicates the current iteration value.
pi,[I],[C],[D] The loop inside bs() function.
I indicates the current iteration value.
C indicates the number of calculated digits (most likely inaccurate)
D indicates the number of digits to be calculated

The [I], [C], [D] above indicates an integer starting from 1.

Technical Spec

This was also an experiment of WebSocket. It was a test of Live result obtaining using WebSocket.

The server is “ws-subins.rhcloud.com” using the Francium DiffSocket library. You can see “ws-subins.rhcloud.com” shows a blank page as the WebSocket server is running. Server :

PHP 5.4.4
Python 2.6.6
MySQL 5.5.45

When someone requests for the full value of Pi, sending 5 MB file (5 Million digits) to client is quite a task. But, WebSocket handles it well.

Conclusion

The record is > 1 Trillion digits. It will surely require a super computer to calculate that much digits. I can see the processor load in my system when just calculating 1 Million digits.

If someone could get me a super computer, we could break the record 😉