DEV Community

Parambir Singh
Parambir Singh

Posted on • Edited on

TIL: Calculating n digits of pi using Chudnovsky Algorithm

Today I stumbled upon Chudnovsky Algorithm to calculate the value of π to N digits of precision. Python code for this algorithm looks like the following:

import decimal

def compute_pi(n):
    decimal.getcontext().prec = n + 1
    C = 426880 * decimal.Decimal(10005).sqrt()
    K = 6.
    M = 1.
    X = 1
    L = 13591409
    S = L

    for i in range(1, n):
        M = M * (K ** 3 - 16 * K) / ((i + 1) ** 3)
        L += 545140134
        X *= -262537412640768000
        S += decimal.Decimal(M * L) / X

    pi = C / S
    return pi
Enter fullscreen mode Exit fullscreen mode

I wrote a simple CLI to test this out:

#!/usr/bin/env python3

import argparse
from lib import compute_pi

parser = argparse.ArgumentParser(description='Print n digits of pi.')
parser.add_argument('num_digits', metavar='N', type=int, help='no. of digits of pi to print')
args = parser.parse_args()

print(compute_pi(args.num_digits))
Enter fullscreen mode Exit fullscreen mode

Here's the CLI in action:

> ./pi.py 30
3.141592653589741586517727315459

> ./pi.py 300
3.141592653589741586517727315457865035780252691261563179943288214795808630531389642185274931230804430454419117074147967105366083976712333542218321180274249883145873143454428446008580088034341219473373000151443532721504141865178673966393142941520166862874509797611548477147655085787688540025728361617601
Enter fullscreen mode Exit fullscreen mode

The complete source code is available in this Github repository.

Top comments (4)

Collapse
 
lysergikproductions profile image
LysergikProductions • Edited

I rewrote the code for use directly in a Terminal by simply running $ python3 [scriptTitle].py [number of iterations to perform]

This code also prevents unnecessary iterations by checking when the last 10 iterations have not updated the accuracy of the requested number of digits.

"""
pi.py
# 2020 (C) Nikolas A. Wagner
# License: GNU GPLv3

# Build_010

                   -- Purpose --
 Stress test CPU for its single-threaded performance

 """
 # setting global variables
import time; import decimal
import os; import sys

txtPi = ""; final_k = 1; its = 1; n = 0
decimalArgs = decimal.Decimal(0)
validArgs = True; startTime = None

# evaluate the command line arguments and set flags for MAIN
def INIT():
    os.system('clear && printf "Initializing..\n"')
    global decimalArgs; global validArgs; global startTime

    try:
        args = int(sys.argv[1]); print(args)
        decimalArgs = decimal.Decimal(args)
    except:
        print("""I can't do that! I need a whole number with no extra symbols!\n""")
        validArgs = False; return
    finally:
        if decimalArgs < 34:
            decimalArgs = 34
            print('Sorry, but I can only be accurate from 34 digits onward.')
            os.system('sleep 2')

        elif decimalArgs > 10000000:
            decimalArgs = 10000000
            print('For the sake of your Terminal, I will limit you to 10 MIL. Edit code to change this if you insist.')
            os.system('sleep 2')

    startTime = time.time()

INIT()

# define functions; finally calling MAIN()
def getUptime():
    return time.time() - startTime

def compute_pi(n):
    global final_k; global txtPi
    decimal.getcontext().prec = n + 1
    its = n; n += 1

    x = 0
    C = 426880 * decimal.Decimal(10005).sqrt()

    K = 6.0; M = 1.0; X = 1
    L = 13591409; S = L

    for k in range(1, n):
        if x == 0:
            print('calculating k= {0}..\n'.format(k - 1))

        lastPi = decimal.Decimal(C / S)
        M = M * (K ** 3 - 16 * K) / ((k + 1) ** 3)
        L += 545140134
        X *= -262537412640768000
        S += decimal.Decimal(M * L) / X
        pi = decimal.Decimal(C / S)

        os.system('clear')
        print('{0}\n'.format(pi))

        if lastPi != pi:
            print('Accuracy improved from iteration {0}'.format(k))
            final_k = k
        else:
            x += 1; n+= 1
            print("Pi not changed {0} time/s".format(x))

            if x >= 10:
                os.system('history -c && clear')
                print(pi); txtPi = str(pi)
                print("\n{0} digits of precision achieved in {1} iterations!".format(its, final_k))
                return ""

def MAIN():
    global validArgs

    if validArgs == True:
        result = compute_pi(int(decimalArgs)); print(result)
        f = open('pi.txt', 'w+'); f.write(txtPi)

        uptime = str(getUptime())
        print("This took {0} seconds".format(uptime))

        print("""
    The last one to three digits are tentative
    The decimal module rounds digits of precision rather than trimming them
""")

    else:
        print('INIT failed; please try again!')

MAIN()
Collapse
 
tigersharks profile image
TigerSharks • Edited

This result does not seem correct. I think the first 39 digits of pi are:

3.14159265358979323846264338327950288420

This does not agree with your result,

3.14159265358974158651772731545786503578
000000000000000^
Difference starts here. Sorry, I did not debug your code. I'm new to Python and wrote this program as an exercise. My results are different and agree with what's published for pi online.

Also, you don't need to compute n-1 loops of the Chudnovsky series to get n digits of precision. To compute 30 digits of precision (31 total) you only need to compute Sigma 4 times -- that includes the freebie you get to start at k=0.

It's clear that something is wrong with your code because your 30th position is a 9, but the 30th position in the 301 digit calculation is a 7. It's not possible to round a 7 to a 9.

I'm curious what went wrong.

Collapse
 
vsierrajc profile image
Juan Carlos Sierra

I check the pi value using the benchmark software SUPER_PI_MOD-1.5 and i get pi with four millions digits
PI=3.1415926535 8979323846 2643383279 5028841971 6939937510
5820974944 5923078164 0628620899 8628034825 3421170679
8214808651 3282306647 0938446095 5058223172 5359408128
4811174502 8410270193 8521105559 6446229489 5493038196
4428810975 6659334461 2847564823 3786783165 2712019091
the first result is correct

Collapse
 
tigersharks profile image
TigerSharks • Edited

In the original post, Parambir shows the output for 30 and 300 digits of pi. It's in the code block right after the "Here's the CLI in action:"

./pi.py 30
3.141592653589741586517727315459

./pi.py 300
3.1415926535897415865177273154578650357802526912....

The results as shown in the author's work are wrong, and don't agree with your results.