Monday, March 14, 2016

Using Your Computer to Calculate Pi

On March 21, 2015, Rajveer Meena recited 70,000 digits of pi from memory in just under 10 hours in Vellore, India, for the Guinness World Record of most digits of pi memorized.  Apart from utter amazement by this feat of memory, have you ever wondered how mathematicians can even calculate such a long list numbers for pi?


There are actually several methods for calculating pi.  You might recall from your math classes that pi is the ratio between the circumference and diameter of any circle, and so one way to calculate pi is to construct a circle, directly measure its circumference and diameter, and divide the two.  Unfortunately, this is only precise enough to obtain a few decimal places of pi.  Another way to calculate pi is to find the perimeter or areas of inscribed or circumscribed polygons of circles.  This was the technique favored by mathematicians before the invention of calculators and computers, but its precision is still limited to “just” a few hundred digits of pi.

inscribed octagon

Today, computers can be used to calculate millions of digits of pi by using infinite series formulas.  There are several different ways to do this, but one of the most efficient methods that is relatively straightforward is to combine Machin’s formula:


with the arc-cotangent power series formula:


The computer program for this is as follows:

01
# Pi Calculator
02
# Python 2.7.3
03
# After running, type "pi(n)" where n is the number of decimals for pi.  For
04
#  example, if you would like to calculate 100 decimals for pi, type "pi(100)".
05

06
# import python libraries
07
from decimal import Decimal, getcontext
08
from time import time, strftime
09
import datetime
10

11
# arccot function using power formula arccot = 1/x - 1/(3x^3) + 1/(5x^5) ...
12
def arccot(x, digits):
13
    # set precision and starting values
14
    getcontext().prec = digits
15
    total = 0
16
    n = 1
17
    # loop while new term is large enough to actually change the total
18
    while Decimal((2 * n - 1) * x ** (2 * n - 1)) < Decimal(10 ** digits):
19
        # find value of new term
20
        term = ((-1)**(n - 1)) * 1 / Decimal((2 * n - 1) * x ** (2 * n - 1))
21
        # add the new term to the total
22
        total += term
23
        # next n
24
        n += 1
25
    # return the sum
26
    return total
27

28
# pi function
29
def pi(decimals):
30
    # start timer
31
    timestart = time()
32

33
    # find pi using Machin's Formula pi = 4 * (4 * arccot(5) - arccot(239))
34
    #  and the power formula for arccot (see arccot function above)
35
    print "pi = " + str(Decimal(4 * (4 * arccot(5, decimals + 3) - arccot(239,
36
        decimals + 3))).quantize(Decimal(10) ** (-decimals)))
37

38
    # display elapsed time
39
    timeelapsedint = round(time() - timestart, 2)
40
    timeelapsedstr = str(datetime.timedelta(seconds = round(
41
        timeelapsedint, 0)))
42
    print "runtime: " + timeelapsedstr + " or " + str(
43
        timeelapsedint) + " seconds."

Here’s what happens when you use the program to find the first 1,000 decimal places for pi on a dual-core 1.67 GHz processor (a mediocre-speed laptop for 2016):

>> pi(1000)
pi = 3.14159265358979323846264338327950288419716939937510
582097494459230781640628620899862803482534211706798214808
651328230664709384460955058223172535940812848111745028410
270193852110555964462294895493038196442881097566593344612
847564823378678316527120190914564856692346034861045432664
821339360726024914127372458700660631558817488152092096282
925409171536436789259036001133053054882046652138414695194
151160943305727036575959195309218611738193261179310511854
807446237996274956735188575272489122793818301194912983367
336244065664308602139494639522473719070217986094370277053
921717629317675238467481846766940513200056812714526356082
778577134275778960917363717872146844090122495343014654958
537105079227968925892354201995611212902196086403441815981
362977477130996051870721134999999837297804995105973173281
609631859502445945534690830264252230825334468503526193118
817101000313783875288658753320838142061717766914730359825
349042875546873115956286388235378759375195778185778053217
12268066130019278766111959092164201989
runtime: 0:00:04 or 3.73 seconds.
  
So with less than 50 lines of code, you can have your computer calculate 1,000 decimal places of pi in less than 4 seconds!  At this rate, you will have 1,000,000 decimal places of pi in just over an hour!  Not bad, considering the very first computer, the ENIAC, took 70 hours to calculate pi to 2,037 decimal places in 1949.  Computers have come a long way since then!


With better computers and more efficient (but more complicated) infinite series formulas, pi has been calculated to over 13.3 trillion digits!  This number is unfathomable to most of us.  If you were to recite this many digits of pi at the same speed as the 2015 world record holder Meena (about 2 digits per second), it would take you over 200,000 years to finish!  So why bother to find so many digits of pi?   Most mathematicians will give the same reason that mountain climbers use for climbing Mount Everest: because it’s there.


9 comments:

  1. I am slightly confused with your power series formula above the code. The middle section appears to be the Laurent series and the right side seems to be either missing a minus sign in the exponent or to be the power series of arctangent. Could you enlighten me?

    ReplyDelete
  2. This comment has been removed by the author.

    ReplyDelete
  3. This commended has been removed by the author.

    ReplyDelete
  4. 4 * arctangent(1) =
    3.1415926535897932384626433832795
    This is the method I use to generate huge pi approximations.

    ReplyDelete