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

(Don't be alarmed by the title; this is a question about mathematics, not programming.)

In the documentation for the decimal module in the Python Standard Library, an example is given for computing the digits of $$\pi$$ to a given precision:

def pi():
"""Compute Pi to the current precision.

>>> print(pi())
3.141592653589793238462643383

"""
getcontext().prec += 2  # extra digits for intermediate steps
three = Decimal(3)      # substitute "three=3.0" for regular floats
lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
while s != lasts:
lasts = s
n, na = n+na, na+8
d, da = d+da, da+32
t = (t * n) / d
s += t
getcontext().prec -= 2
return +s               # unary plus applies the new precision


I was not able to find any reference for what formula or fact about $$\pi$$ this computation uses, hence this question.

Translating from code into more typical mathematical notation, and using some calculation and observation, this amounts to a formula for $$\pi$$ that begins like:

\begin{align}\pi &= 3+\frac{1}{8}+\frac{9}{640}+\frac{15}{7168}+\frac{35}{98304}+\frac{189}{2883584}+\frac{693}{54525952}+\frac{429}{167772160} + \dots\\ &= 3\left(1+\frac{1}{24}+\frac{1}{24}\frac{9}{80}+\frac{1}{24}\frac{9}{80}\frac{25}{168}+\frac{1}{24}\frac{9}{80}\frac{25}{168}\frac{49}{288}+\frac{1}{24}\frac{9}{80}\frac{25}{168}\frac{49}{288}\frac{81}{440}+\frac{1}{24}\frac{9}{80}\frac{25}{168}\frac{49}{288}\frac{81}{440}\frac{121}{624}+\frac{1}{24}\frac{9}{80}\frac{25}{168}\frac{49}{288}\frac{81}{440}\frac{121}{624}\frac{169}{840}+\dots\right) \end{align}

or, more compactly,

$$\pi = 3\left(1 + \sum_{n=1}^{\infty}\prod_{k=1}^{n}\frac{(2k-1)^2}{8k(2k+1)}\right)$$

Is this a well-known formula for $$\pi$$? How is it proved? How does it compare to other methods, in terms of how how quickly it converges, numerical stability issues, etc? At a glance I didn't see it on the Wikipedia page for List of formulae involving π or on the MathWorld page for Pi Formulas.

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

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)

• Thank you, great to hear from the original author of the code! Dec 8 '18 at 19:34
• Can you explain how you rearranged $\pi = 3+1/8+\cdots$ to $\pi= 3+1/24+\cdots$? Dec 12 '18 at 7:50
• @taritgoswami That was a typo, the latter should have been $pi = 3 + 3(1/24) + 3(1/24)(9/80) ...$ . It's fixed now. Thanks for noticing. Dec 13 '18 at 7:54

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

• Thanks, would you know anything about how it compares to other methods? E.g. I imagine it's better than the Leibniz formula for π which converges very slowly, and worse than the best methods. Dec 6 '18 at 19:56
• Leibniz formula for $\pi$ has very slow convergence. Formulas using power series based on inverse trigonometric functions (like the one above) converge much faster, but there are even faster algorithms such as Brent-Salamin (which doubles the number of correct digits in each iteration). There is a chronology here: en.wikipedia.org/wiki/Chronology_of_computation_of_%CF%80 Dec 6 '18 at 21:45
• Here is a very fast algorithm for computing $\pi$ and its implementation in python: en.wikipedia.org/wiki/Chudnovsky_algorithm Dec 6 '18 at 21:52
• @mlerma54: Note that the Leibniz series too can be viewed as being "based on inverse trigonometric functions" -- it corresponds to $4\arctan(1)$, with the series for the arctangent evaluated right at its radius of convergence. Dec 7 '18 at 2:15
• Yes, that is correct, Leibniz formula for $\pi$ can be seen as based on the Taylor series for $\arctan(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \cdots$, but evaluated at a particularly "bad" place ($x=1$), so it does not take advantage of the exponential convergence of the $n$th term that we see in the other power series such as $\arcsin(x)$ evaluated at $1/2$, and a few other such as 1706 Machin's $\frac{\pi}{4} = 4\arctan{\frac{1}{5}} - \arctan{\frac{1}{239}}$. Dec 7 '18 at 3:38

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.