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

Question

(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.

Answer 1

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)

Answer 2

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

Answer 3

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.