4 Ways to Calculate Pi

pi_day.py

from functools import reduce
from math import factorial, comb, ceil, log2
from decimal import getcontext, Decimal as Dec


def bernoulli(n):
    bs = [Dec(1)]
    for m in range(1, n+1):
        bs.append(1 - sum(comb(m, k)*b / (m - k + 1) for k, b in zip(range(m), bs)))
    return abs(bs[-1])


# Formula taken from Plouffe (2022): https://arxiv.org/abs/2201.12601
def pi(n):
    getcontext().prec = n
    return (2*factorial(n)/((bernoulli(n))*2**n*reduce(Dec.__mul__, [1-(1/Dec(x)**n) for x in [2, 3, 5, 7]])))**(1/Dec(n))


# Bailey–Borwein–Plouffe formula for reference.
def bbp(n):
    getcontext().prec = n
    return sum(1/Dec(16**k) * (4/Dec(8*k+1)-2/Dec(8*k+4)-1/Dec(8*k+5)-1/Dec(8*k+6)) for k in range(n))


# Gauss-Legendre formula
def gauss_legendre(prec):
    getcontext().prec = prec
    half = Dec(1/2)
    a = Dec(1)
    b = Dec(a / Dec(2)**Dec(1/2))
    t = Dec(a / 4)
    for x in map(lambda x: Dec(2**x), range(ceil(log2(prec)))):
        y=a
        a+=b
        a*=half
        t-=x*(y-a)**2
        b*=y
        b **= half
    return (a + b)**2 / (4 * t)


if __name__ == "__main__":
    precision = 1000
    print(b := bbp(precision))
    print(c := gauss_legendre(precision))
    print(c := pi(precision))


print(''.join(b:=open(__file__).readlines()[49:]).count("#")/len(b)/4)
            #######
       #################
    #######--------########
  ######-------------#######
 ######-----------------######
 #####------------------######
######-------------------######
#####--------------------######
#####---------!\---------######
######--------!$#\-------######
 ######-------!#$#$#\---######
 ######-------!#$#$##$#\#####
  #######-----!$#$#$#$#$#\###
    ##########!$#$#$#$#$##\
        ######!$#$#$###
            ##!####

In celebration of Pi Day I decided to compile together some of my previous algorithms in Python for calculating Pi, in addition to two new ones. The first algorithm is the one published by Plouffe in 2022, based on Bernoulli numbers. After that we have the classic BBP, then Gauss-Legendre, and then my own bastardized version of the classic 1988 IOCCC entry by Westley.

Example output:

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420193
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420201
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199
3.140625

Hi Julius,

This is extremely cool, I was looking into the Plouffe algorithm from 2022 and was wondering if you had an implementation which returned only the nth digit?

I.e. func(0) = 3, func(2) = 4, etc?

You seem extremely knowledgeable about the calculation of Pi so I was just wondering if you knew whether the above would infact be the fastest way to calculate the nth decimal digit of Pi or at what position n it becomes quicker than existing algorithms which calculate all preceding digits to determine the nth value.

@LukePrior I was actually looking into that a month or so ago, but seem to have lost the code. For the Plouffe paper, it seems that individual digits are calculated by getting the value from the normal formula for all the digits (formula at the bottom of the second page). The hardest part about the Plouffe algorithm is that it requires relatively large bernoulli numbers, which are pretty computationally intense to calculate. There are other approaches to calculating the nth digit of pi, but they generally calculate the nth digit in some other base, like 16. It doesn't seem like there is currently an algorithm for the nth decimal digit of pi according to the internet: https://math.stackexchange.com/questions/4491628/is-there-a-formula-that-allows-you-calculate-the-n-th-decimal-digit-of-pi-withou