¿Cuál es la fórmula de pi utilizada en la biblioteca decimal de Python?

Question

(No te alarmes por el título; esta es una pregunta sobre matemáticas, no sobre programación).

En la documentación del decimal de la Biblioteca Estándar de Python, un ejemplo se da para calcular los dígitos de $\pi$ con una precisión determinada:

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

No he podido encontrar ninguna referencia sobre qué fórmula o hecho sobre $\pi$ este cómputo utiliza, de ahí esta pregunta.

Traduciendo el código a una notación matemática más típica, y utilizando algunos cálculos y observaciones, esto equivale a una fórmula para $\pi$ que comienza como:

$$\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}$$

o, de forma más compacta,

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

¿Es ésta una fórmula conocida para $\pi$ ? ¿Cómo se demuestra? ¿Cómo se compara con otros métodos, en términos de rapidez de convergencia, problemas de estabilidad numérica, etc.? De un vistazo no lo he visto en la página de Wikipedia para Lista de fórmulas que implican o en la página de MathWorld para Fórmulas Pi .

Answer 1

Esta es la serie Taylor de $\arcsin(x)$ en $x=1/2$ (tiempos 6).

Answer 2

Esta aproximación para $\pi$ se atribuye a 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

Cuando escribí el código que se muestra en los documentos de Python, obtuve la fórmula de la página 53 de " La alegría de ". De las muchas fórmulas enumeradas, fue la primera que:

1. convergieron rápidamente,
2. era corto,
3. era algo que entendía lo suficientemente bien como para derivarlo a mano, y
4. podría implementarse utilizando operaciones baratas: varias sumas con una sola multiplicación y una sola división para cada término. Esto permitió la estimación de $\pi$ se puede escribir fácilmente como una función eficiente utilizando los flotadores de Python, o con el módulo decimal, o con los enteros de precisión múltiple de Python.

La fórmula se resuelve en la ecuación $sin(\pi/6)=\frac{1}{2}$ .

WolframAlpha da la serie Maclaurin para $6 \arcsin{(x)}$ como:

$$6 \arcsin{(x)} = 6 x + x^{3} + \frac{9 x^{5}}{20} + \frac{15 x^{7}}{56} + \frac{35 x^{9}}{192} + \dots $$

La evaluación de la serie en $x = \frac{1}{2}$ da:

$$ \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\\ $$

A partir de ahí, utilicé diferencias finitas para calcular los numeradores y denominadores de forma incremental. Las diferencias del numerador fueron 8, 16, 24, ... por lo que el ajuste del numerador na+8 en el código. Las diferencias del denominador fueron 56, 88, 120, ... de ahí el ajuste del denominador da+32 en el código:

 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

Aquí está el código original que escribí en 1999 utilizando los enteros de precisión múltiple de Python (esto es anterior al módulo decimal):

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 3

Es sólo la computación $\pi = 6\sin^{-1}\left(\frac12\right)$ utilizando la expansión en serie de Taylor del arcoseno. Como referencia,

$$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$$ y compara los coeficientes con lo que obtienes.