Pi with arbitrary precision using Bernoulli numbers (the slow way)

pi_arb.py

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


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


if __name__ == "__main__":
    precision = 200
    print(p := pi(precision))
    print(b := bbp(precision))
    print(timeit(lambda: pi(precision), number=10))
    print(timeit(lambda: bbp(precision), number=10))
    print(abs(p-b))

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303817
0.20411962500656955
0.010316250001778826
2E-199