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

    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)
  • 8
    $\begingroup$ Thank you, great to hear from the original author of the code! $\endgroup$ Dec 8 '18 at 19:34
  • $\begingroup$ Can you explain how you rearranged $\pi = 3+1/8+\cdots $ to $\pi= 3+1/24+\cdots$? $\endgroup$ Dec 12 '18 at 7:50
  • 1
    $\begingroup$ @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. $\endgroup$ Dec 13 '18 at 7:54

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

  • $\begingroup$ 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. $\endgroup$ Dec 6 '18 at 19:56
  • 4
    $\begingroup$ 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 $\endgroup$
    – mlerma54
    Dec 6 '18 at 21:45
  • 2
    $\begingroup$ Here is a very fast algorithm for computing $\pi$ and its implementation in python: en.wikipedia.org/wiki/Chudnovsky_algorithm $\endgroup$
    – mlerma54
    Dec 6 '18 at 21:52
  • $\begingroup$ @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. $\endgroup$ Dec 7 '18 at 2:15
  • $\begingroup$ 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}}$. $\endgroup$
    – mlerma54
    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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.