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.

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
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 :

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 😉