Why does this loop yield $\pi/8$?

Question

Was doing some benchmarks for the mpmath python library on my pc with randomly generated tests.

I found that one of those tests was returning a multiple of $\pi$ consistently. I report it here:

from mpmath import *
mp.dps=50

print "correct pi: ", pi

den=mpf("3")
a=0
res=0

for n in xrange(0,1000):
    res=res+1/den
    a=a+32
    den=den+a

print "result:     ", res*8

Result:

correct pi:  3.1415926535897932384626433832795028841971693993751
result:      3.1410926536210432286970258295579986968615976706527

The convergence rate is very low. 10000 terms only give accuracy to the 5th-4th decimal digit. Is it there somewhere a formula for pi that involves the number 32?

Answer 1

Per String's answer, continuing from his result we have \begin{align*} res &= \sum_{i=0}^{\infty}\frac{1}{16i(i+1)+3}\\ &= \sum_{i=0}^{\infty} \frac{1}{(4i+1)(4i+3)} \\ &= \frac{1}{2}\sum_{i=0}^{\infty} \left(\frac{1}{4i+1}-\frac{1}{4i+3}\right) \\ &= \frac{1}{2}\left(\frac{1}{1} - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots\right) \\ &= \frac{\pi}{8}. \end{align*} This derivation also explains why the series converges so slowly....this particular approximation converges as the alternating harmonic series converges to zero, which is very slow.

Answer 2

From your program one can infer that we must have $$ den_n=f(n)=32 T(n)+3 $$ Where $T(n)=1+2+...+n=\frac{n(n+1)}{2}$ is a so-called Triangular Number. Then you form the sum $$ \begin{align} res_n&=\sum_{i=0}^n\frac{1}{f(i)}\\ &=\sum_{i=0}^n\frac{1}{32 T(i)+3}\\ &=\sum_{i=0}^n\frac{1}{16i(i+1)+3}\\ \end{align} $$ So the question remains: is this series approximating $\pi/8$ and why ... Others have elaborated on that :)

Answer 3

Your algorithm is equivalent to the the sum: $$\pi'=8\left(\frac{1}{3}+\frac{1}{3+32}+\frac{1}{3+32+64}+\frac{1}{3+32+64+96}+\ldots\right)$$ Or: $$\pi'=\sum_{k=0}^\infty \frac{8}{3+16k(1+k)}= \frac{1}{2}\sum_{k=0}^\infty \frac{1}{(1/4+k)(3/4+k)}$$ These sums can be expressed as digamma functions, since: $$\sum_{k=0}^\infty \frac{1}{(a+k)(b+k)}=\frac{1}{b-a}\left(\psi(b)-\psi(a)\right)$$ Finally, we are left with: $$\psi(3/4)-\psi(1/4)$$ Which by the identity: $$\psi(1 - x) - \psi(x) = \pi\,\!\cot{ \left ( \pi x \right ) }$$ Is indeed equal to $\pi$.