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


3 Answers 3


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, 2018 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, 2018 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, 2018 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, 2018 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, 2018 at 3:38

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

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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .