What is the formula for pi used in the Python decimal library?

That is the Taylor series of $\arcsin(x)$ at $x=1/2$ (times 6).

This approximation for $\pi$ is attributed to Issac Newton:

• https://loresayer.com/2016/03/14/pi-infinite-sum-approximation/
• http://www.geom.uiuc.edu/~huberty/math5337/groupe/expresspi.html
• http://www.pi314.net/eng/newton.php

When I wrote that code shown in the Python docs, I got the formula came from p.53 in "The Joy of π". Of the many formulas listed, it was the first that:

1. converged quickly,
2. was short,
3. was something I understood well-enough to derive by hand, and
4. could be implemented using cheap operations: several additions with only a single multiply and single divide for each term. This allowed the estimate of $\pi$ to be easily be written as an efficient function using Python's floats, or with the decimal module, or with Python's multi-precision integers.

The formula solves for π in the equation $sin(\pi/6)=\frac{1}{2}$.

WolframAlpha gives the Maclaurin series for $6 \arcsin{(x)}$ as:

$$6 \arcsin{(x)} = 6 x + x^{3} + \frac{9 x^{5}}{20} + \frac{15 x^{7}}{56} + \frac{35 x^{9}}{192} + \dots $$

Evaluating the series at $x = \frac{1}{2}$ gives:

$$ \pi \approx 3+3 \frac{1}{24}+3 \frac{1}{24}\frac{9}{80}+3 \frac{1}{24}\frac{9}{80}\frac{25}{168}+\dots + \frac{(2k+1)^2}{16k^2+40k+24} + \dots\\ $$

From there, I used finite differences, to incrementally compute the numerators and denominators. The numerator differences were 8, 16, 24, ..., hence the numerator adjustment na+8 in the code. The denominator differences were 56, 88, 120, ..., hence the denominator adjustment da+32 in the code:

 1     9    25    49    numerators
    8    16    24       1st differences
       8     8          2nd differences


24    80   168   288    denominator 
   56    88   120       1st differences
      32    32          2nd differences

Here is the original code I wrote back in 1999 using Python's multi-precision integers (this predates the decimal module):

def pi(places=10):
    "Computes pi to given number of decimal places"
    # From p.53 in "The Joy of Pi".  sin(pi/6) = 1/2
    # 3 + 3*(1/24) + 3*(1/24)*(9/80) + 3*(1/24)*(9/80)*(25/168)
    # The numerators 1, 9, 25, ... are given by  (2x + 1) ^ 2
    # The denominators 24, 80, 168 are given by 16x^2 +40x + 24
    extra = 8
    one = 10 ** (places+extra)
    t, c, n, na, d, da = 3*one, 3*one, 1, 0, 0, 24
    while t > 1:
        n, na, d, da  = n+na, na+8, d+da, da+32
        t = t * n // d
        c += t
    return c // (10 ** extra)

It just computing $\pi = 6\sin^{-1}\left(\frac12\right)$ using the Taylor series expansion of arcsine. For reference,

$$6\sin^{-1}\frac{t}{2} = 3t+\frac{t^3}{8}+\frac{9t^5}{640}+\frac{15t^7}{7168}+\frac{35 t^9}{98304} + \cdots$$ and compare the coefficients with what you get.