# Playing in Python with Pickover’s Premise per Pi

By Monty, 1st June 2015

Pi by J.Gabás Esteban

My favourite tweeter is Clifford A. Pickover. His tweets are delightful nuggets of math, physics & more. I have a copy of his fascinating The Physics Book. However, in these tweets, he claims that the string 44899 first occurs in pi at position 44899, counting from the first digit after the decimal point. I thought it might be interesting to check this, and to see if there are any more such ‘self-references’ in say, the first 100,000,000 digits of pi. There’s obviously 1 at position 1, which mathematicians would probably call ‘trivial’.

So my first task, is to find a way to compute pi to far more than the accuracy of math.pi in Python. A quick search of the web sent me to Nick Craig-Wood’s pages on computing pi in Python using the Chudnovsky algorithm. There I found pi_chudnovsky_bs_gmpy.py. It needs the gmpy2 library. I needed to change a linesqrtC = (10005*one_squared).sqrt()  (line 63) to sqrtC = gmpy2.isqrt(10005*one_squared).

XKCD: Duty Calls

I then wrote the following quick and dirty script, which probably isn’t as efficient or as beautiful as it could be – but it got the job done! The len(str(digits))-1 bit adds a few digits to allow for the length of the string being searched for. In the first 100,000,000 digits of pi, the numbers 1, 16470, 44899, 79873884 are the only ones which occur at the positions they index. I didn’t do any timing, but it only took a few minutes! When I first wrote this, I overlooked the word first in the claim; so then I added pistr.find(istr) to check where the found number first occurred. Only 1 satisfies the requirements. The program’s output is:

Computing Pi
Searching for Pickover numbers
1  first occurs at  1
16470  first occurs at  1602
44899  first occurs at  13714
79873884  first occurs at  46267046


As a ball-park check on this, go to 100,000 digits of Pi and use your browser’s find function to locate 44899. This number is slightly under one half of 100,000, but there’s a couple of occurrences of it well before halfway. Or better yet, go to The Pi-Search Page. I’ve replied to Pickover’s tweet and emailed him to let him know the claim isn’t correct. He is wrong on the Internet.

UPDATE: He kindly responded to my email, and agrees that the claim as worded is wrong.

#   Pi_ckover.py
#
#   http://www.craig-wood.com/nick/articles/pi-chudnovsky/

import pi_chudnovsky_bs as pi

digits = 100000000
print("Computing Pi")
pi = pi.pi_chudnovsky_bs(digits+len(str(digits))-1)     # compute pi
pistr = str(pi)                                         # convert to string

print("Searching for Pickover numbers")
for i in range(1, digits+1):        # search along from 1 to digits
istr = str(i)                       # convert i to string
l = len(istr)                       # get its length
if istr == pistr[i:i+l]:
print(istr," first occurs at ",pistr.find(istr))


### Further Information

Wikipedia:

The Chudnovsky algorithm is a fast method for calculating the digits of π. It was published by the Chudnovsky brothers in 1989, and was used in the world record calculations of 2.7 trillion digits of π in December 2009, 5 trillion digits in August 2010, 10 trillion digits in October 2011, 12.1 trillion digits in December 2013 and 22.4 trillion digits of π in November 2016.

The algorithm is based on the negated Heegner number ${\displaystyle d=-163}$, the j-function ${\displaystyle \scriptstyle j\left({\frac {1+{\sqrt {-163}}}{2}}\right)=-640320^{3}}$, and on the following rapidly convergent generalized hypergeometric series:

${\displaystyle {\frac {1}{\pi }}=12\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}\left(640320\right)^{3k+3/2}}}}$

For a high performance iterative implementation, this can be simplified to

${\displaystyle {\frac {(640320)^{3/2}}{12\pi }}={\frac {426880{\sqrt {10005}}}{\pi }}=\sum _{k=0}^{\infty }{\frac {(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}\left(-262537412640768000\right)^{k}}}}$

There are 3 big integer terms (the multinomial term Mk, the linear term Lk, and the exponential term Xk) that make up the series and π equals the constant C divided by the sum of the series, as below:

${\displaystyle \pi =C\left(\sum _{k=0}^{\infty }{\frac {M_{k}\cdot L_{k}}{X_{k}}}\right)^{-1}}$, where:
${\displaystyle C=426880{\sqrt {10005}},\quad \quad M_{k}={\frac {6k!}{(3k)!(k!)^{3}}},\quad \quad L_{k}=545140134k+13591409,\quad \quad X_{k}=(-262537412640768000)^{k}}$

The terms Mk, Lk, and Xk satisfy the following recurrences and can be computed as such:

{\displaystyle {\begin{alignedat}{4}L_{k+1}&=L_{k}+545140134\,\,&&{\textrm {where}}\,\,L_{0}&&=13591409\\[4pt]X_{k+1}&=X_{k}\cdot (-262537412640768000)&&{\textrm {where}}\,\,X_{0}&&=1\\[4pt]M_{k+1}&=M_{k}\cdot \left({\frac {(12k+2)(12k+6)(12k+10)}{(k+1)^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[4pt]\end{alignedat}}}

The computation of Mk can be further optimized by introducing an additional term Kk as follows:

{\displaystyle {\begin{alignedat}{4}K_{k+1}&=K_{k}+12\,\,&&{\textrm {where}}\,\,K_{0}&&=6\\[4pt]M_{k+1}&=M_{k}\cdot \left({\frac {K_{k}^{3}-16K_{k}}{(k+1)^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[12pt]\end{alignedat}}}

Note that

${\displaystyle e^{\pi {\sqrt {163}}}\approx 640320^{3}+743.99999999999925\dots }$ and
${\displaystyle 640320^{3}=262537412640768000}$
${\displaystyle 545140134=163\cdot 127\cdot 19\cdot 11\cdot 7\cdot 3^{2}\cdot 2}$
${\displaystyle 13591409=13\cdot 1045493}$

This identity is similar to some of Ramanujan's formulas involving π, and is an example of a Ramanujan–Sato series.

The time complexity of the algorithm is ${\displaystyle O(n\log(n)^{3})}$.