Krdytkyu

(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\
&= 3left(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}+dotsright)
end{align}$$

or, more compactly,

$$pi = 3left(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:

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

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

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

$$6sin^{-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.