Pi and normality
David H. Bailey
Lawrence Berkeley National Laboratory (retired) andComputer Science Department, University of California, Davis

1 / 36

New York Times PiDay 2007 (March 14, 2007) crossword puzzle

2 / 36

Answers to crossword

3 / 36

I once received a strange fax
I In October 1992, I received this fax fromthe Simpsons TV show.
I They wanted the 40,000th digit of π.
I I faxed back the result: it is a “1.”I This was used in the Simpsons show,

dated 6 May 1993, “Marge in Chains.”

4 / 36

Smartphones and π
I Every smartphone or mobile phonecrucially relies on computations (e.g.,
the fast Fourier transform) thatinvolve π to resolve microwave signals.

I π appears in the fundamental equations of quantum mechanics, which are used todesign smartphone electronics. For example, Heisenberg’s uncertainty principle:
Z ∞ Z ∞
||f ||422
2
2
2
s |f (s)| ds
t |f (t)| dt ≥
16π 2
−∞
−∞
I π appears in the equations of general relativity, used in GPS:

Rµν −

Rgµ ν
8πG
= 4 Tµν2c
5 / 36

√“Exact” value of π = 17 − 8 3 = 3.1435935394 . . .
In November 2016, the IOSR Journal of Mathematics
published the following paper:
Abstract: It is believed that pi (π) is atranscendental number. In author’s opinion, it is
not the fact. The paper aims at showing that pi(π) is √an algebraic number with exact value
17 − 8 3. The derivation of this value is supported
by several geometrical constructions, arithmetic
calculations and use of some simple algebraic
formulae.

Another author√ has published 8 papers insisting that
π = (14 − 2)/4 = 3.1464466094 . . ., in supposedly“peer-reviewed” journals.
√

I “Exact value of pi π(17 − 8 3),” IOSR J. of Mathematics, vol. 12 (Nov.-Dec. 2016), http://www.iosrjournals.org/iosr-jm/papers/Vol12-issue6/Version-1/B1206010408.pdf.
6 / 36

Why compute π?
I Question: Do we need to know π to thousands or millions of digits in everyday

science and engineering?
Answer: No. 10–15 digits suffice for most scientific calculations.
I However, some research problems in mathematics and physics require hundreds or

thousands of digits.
I I have personally done computations that required π to 64,000-digit precision.
I Billions and even trillions of digits have been computed by mathematicians, inpart to explore the unanswered question “Are the digits of π ‘random’ ?”7 / 36

The first 1000 decimal digits of π
3.14159265358979323846264338327950288419716939937510582097494459230781
6406286208998628034825342117067982148086513282306647093844609550582231
7253594081284811174502841027019385211055596446229489549303819644288109
7566593344612847564823378678316527120190914564856692346034861045432664
8213393607260249141273724587006606315588174881520920962829254091715364
3678925903600113305305488204665213841469519415116094330572703657595919
5309218611738193261179310511854807446237996274956735188575272489122793
8183011949129833673362440656643086021394946395224737190702179860943702
7705392171762931767523846748184676694051320005681271452635608277857713
4275778960917363717872146844090122495343014654958537105079227968925892
3542019956112129021960864034418159813629774771309960518707211349999998
3729780499510597317328160963185950244594553469083026425223082533446850
3526193118817101000313783875288658753320838142061717766914730359825349
0428755468731159562863882353787593751957781857780532171226806613001927
876611195909216420198...
8 / 36

Pre-computer history of π calculations
Name
Archimedes
PtolemyLiu Hui
Aryabhata
Tsu Ch’ung Chi
Madhava
Al-Kashi
Romanus
Van Ceulen
Sharp and Halley
MachinStrassnitzky and Dase
Rutherford
W. Shanks
Ferguson (mechanical calculator)

Year
-250?
150?
265?
480?
480?
1400?
1429
1593
1615
1699
1706
1844
1853
1874
1947

Digits
3
3
5
5
7
13
14
15
35
71
100
200
440
∗
707
808

∗ Only the first 527 were correct.
9 / 36

Computer-era π calculations
Name
Miyoshi and Kanada
Kanada-Yoshino-Tamura
GosperBaileyKanada et. al
Kanada and Tamura
Chudnovskys
Kanada and Takahashi
Kanada and Takahashi
Kanada-Ushiro-Kuroda
Takahashi
Bellard
Kondo and Yee
Trueb

Year
1981
1982
1985
1986
1987
1989
1994
1997
1999
2002
2009
2009
2010
2016

Digits2,000,036
16,777,206
17,526,200
29,360,111
134,217,700
1,073,741,799
4,044,000,000
51,539,600,000
206,158,430,000
1,241,100,000,000
2,576,980,377,524
2,699,999,990,000
5,000,000,000,000
22,459,157,718,361

If 22 trillion digits were printed in 12-point type, they would stretch nearly to Mars.
10 / 36

A random walk on the first 100 billion base-4 digits of πThis dataset can be explored online: http://gigapan.com/gigapans/106803

11 / 36

Some formulas for computing π

√
∞
2n
X
3 3
n
π=
− 24
4
(2n + 3)(2n − 1)42n+1
n=0(Newton, 1660)

∞
∞
X
X
(−1)n
(−1)n
π
=4
−
(Machin, 1730)
4(2n + 1)52n+1 n=0 (2n + 1)2392n+1
n=0√ ∞2 2 X (4n)!(1103 + 26390n)
1=(Ramanujan, 1930)π9801 n=0(n!)4 3964n√√
Set a0 = 6 − 4 2 and y0 = 2 − 1. Iterateyk+1 =

1 − (1 − yk4 )1/4
1 + (1 − yk4 )1/4

2
ak+1 = ak (1 + yk+1 )4 − 22k+3 (1 + yk+1 + yk+1).

Then ak converge quartically to 1/π: each iteration quadruples the number of correct digits.
20 iterations are sufficient to compute π to 2.9 trillion digits, provided all iterations are done tothis precision. (Jonathan and Peter Borwein, 1984)
12 / 36

How does one do arithmetic to extremely high precision?
Computing π or anything else to extremely high precision requires special software:
I High-precision numbers are stored as a sequence of computer words.
I Addition and subtraction are performed using relatively simple methods.
I Multiplication is performed using a fast Fourier transform, which is thousands oreven millions of times faster than conventional methods.
I Division and square roots are performed using Newton iterations, based on

multiplication and addition.
I Exponential and trigonometric functions are evaluated using special algorithms.

Software packages to perform these operations are readily available on the Internet, or
by using systems such as Mathematica, Maple or Sage.

13 / 36

The PSLQ integer relation algorithm
Given a vector (xn ) of real numbers, an integer relation algorithm finds integers (an )such thata1 x1 + a2 x2 + · · · + an xn = 0(to within the precision of the arithmetic being used), or else finds bounds within
which no relation can exist.
Helaman Ferguson’s PSLQ algorithm is the most widely used integer relation algorithm.
Integer relation detection (using PSLQ or any other algorithm) requires very highnumeric precision, both in the input data and in the operation of the algorithm.
1. H. R. P. Ferguson, D. H. Bailey and S. Arno, “Analysis of PSLQ, an integer relation finding
algorithm,” Mathematics of Computation, vol. 68, no. 225 (Jan 1999), 351–369.
2. D. H. Bailey and D. J. Broadhurst, “Parallel integer relation detection: Techniques andapplications,” Mathematics of Computation, vol. 70, no. 236 (Oct 2000), 1719–1736.
14 / 36

Helaman Ferguson’s “Umbilic Torus SC” sculpture at Stony Brook Univ.

15 / 36Computing binary digits of log(2) beginning at an arbitrary position
1996 result: Consider this well-known formula for log(2):log(2) =∞X1= 0.101100010111001000010111111101111101000111001111011 . . .2n2nn=1Note that the binary digits of log 2 beginning after position d can be written as{2d log 2}, where {·} denotes fractional part. Thus we can write:
( d
) ( ∞
)
X 2d−nX 2d−n{2d log(2)} =
+
n
n
n=1
n=d+1
( d
) ( ∞
)
X 2d−n mod n
X 2d−n
=
+
n
n
n=1

n=d+1

We have inserted “modn” since were are only interested in the fractional part when
divided by n. Now note that the numerator 2d−n mod n can be calculated very rapidly
using the binary algorithm for exponentiation.
16 / 36

The binary algorithm for exponentiation
Problem: What is 317 mod 10? (i.e., what is the last decimal digit of 317 ?)Algorithm A:
317 = 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 × 3 = 129140163,
so answer = 3.
Algorithm B (faster): 317 = ((((32 )2 )2 )2 ) × 3 = 129140163, so answer = 3.
Algorithm C (fastest):
317 = ((((32 mod 10)2 mod 10)2 mod 10)2 mod 10) × 3 mod 10 = 3.
Note that in Algorithm C, we never have to deal with integers larger than 9 × 9 = 81,
so the entire operation can be performed very rapidly on a computer.

17 / 36

General BBP-type formulas
The same “trick” that was used for log(2) can be applied for any real constant α thatcan be written in the form
α=

∞
X
p(n)
b n q(n)n=0where p and q are integer polynomials, deg p < deg q, and q has no zeroes for n ≥ 0,or as a linear sum of such formulas.
What other well-known mathematical constants can be written by such a formula?
Can π be written in this form? None was known at the time (1996).

18 / 36

The BBP formula for π
In 1996, a PSLQ program discovered this new formula for π:


∞
X
4
1
2
1
1
π=
−
−
−
16k 8k + 1 8k + 4 8k + 5 8k + 6
k=0

Indeed, this formula permits one to compute base-16 (or binary) digits of π beginningat an arbitrary starting position. The proof is simple.
This is the first known instance of a computer program discovering a fundamentallynew formula for π.
BBP-type formulas (also discovered using PSLQ) are now known for numerous othermathematical constants.
Sadly, there is no such similar formula for base-10 digits of π.
I D. H. Bailey, P. B. Borwein and S. Plouffe, “On the rapid computation of various polylogarithmic

constants,” Mathematics of Computation, vol. 66 (Apr 1997), 903–913.

19 / 36

Some other BBP-type
formulas found using PSLQ

∞


216
72
54
9
144
−
−
−
+
(6k + 1)2
(6k + 2)2
(6k + 3)2
(6k + 4)2
(6k + 5)2
k=0

∞
405
81
27
2 X 1
243
−
−
−
π2 =
27
729k (12k + 1)2(12k + 2)2(12k + 4)2(27k + 5)2
k=0

9
9
5
1
72
−
−
−
+
−
(12k + 6)2(12k + 7)2(12k + 8)2(12k + 10)2
(12k + 11)2

∞
1 X 1
18432
69120
18432
25344
27648
π 2 log(2) =
−
+
+
+
32
4096k (24k + 2)3(24k + 3)3(24k + 4)3(24k + 6)3(24k + 8)3
π2 =

1X 1
8
64k

k=0

8640
1152
2880
288
1080
1728
+
+
+
+
+
+
(24k + 9)3(24k + 10)3
(24k + 12)3
(24k + 14)3
(24k + 15)3
(24k + 16)3

396
72
135
18
+
+
−
+
(24k + 18)3
(24k + 20)3
(24k + 21)3
(24k + 22)3

I David H. Bailey, “A compendium of BBP-type formulas for mathematical constants,” updated 15

Aug 2017, http://www.davidhbailey.com/dhbpapers/bbp-formulas.pdf.
20 / 36

The BBP formula for π in action
In July 2010, Tsz-Wo Sze used a variant of the BBP formula to compute the base-16 digits ofπ starting at position 500 trillion (corresponding to binary position
2 quadrillion). The run required 16 billion CPU-seconds of computing.
The result was checked by repeating the calculation to find digits starting at position 500
trillion + 1. The two results were (in base-16 digits):
0 E6C1294A ED40403F 56D2D764 026265BC A98511D0 FCFFAA10 F4D28B1B
E6C1294A ED40403F 56D2D764 026265BC A98511D0 FCFFAA10 F4D28B1B
Note that the results precisely overlap. The probability that two randomly generated 56-longstrings of base-16 digits perfectly agree is approximately 3.7 × 10−68 .
Philosophical question: What is more securely established?:
I The computational assertion that the 500 trillionth hex digit of pi is “0.”I A theorem whose proof is hundreds of pages long, relies on dozens of earlier results (anyone of which, if later found to be in error, would render the main theorem invalid), andwhich has been read in detail by only a handful of mathematicians worldwide?
In some cases, at least, computation is at least as compelling as formal proof.

21 / 36

BBP-type formulas for π 2
Whereas only base-2 (binary) BBP-type formulas exist for π, there are both binary(base-2) and ternary (base-3) formulas for π 2 , both discovered by PSLQ:
∞



16
24
8
6
1
−
−
−
+
(6k + 1)2 (6k + 2)2 (6k + 3)2 (6k + 4)2 (6k + 5)2
k=0

∞
243
405
81
27
2 X 1
2
−
−
−
π =
2
2
2
k
27
(12k + 1)(12k + 2)(12k + 4)(12k + 5)2
729
k=0

72
9
9
5
1
−
−
−
−
+
(12k + 6)2 (12k + 7)2 (12k + 8)2 (12k + 10)2 (12k + 11)2
9X 1
π =
8
64k
2



We decided to use these formulas to compute base-64 and base-729 digits of π 2 ,starting at position ten trillion.

22 / 36

BBP-type formula for Catalan’s constant
We also decided to calculate digits of Catalan’s constant:
∞
X
(−1)n
G=
= 0.91596559417722 . . .
(2n + 1)2n=0which is closely related to π 2 :
∞
1
π2 X
=
= 1.2337005501362 . . .
8(2n + 1)2n=0

We employed this formula, which we discovered using the PSLQ algorithm:

∞
36864
30720
30720
6144
1536
1 X 1
−
−
−
−
G=
2
2
2
2
k
4096
(24k + 2)(24k + 3)(24k + 4)(24k + 6)(24k + 7)2
4096
k=0

2304
768
480
384
1536
2304
+
+
+
+
+
+
2
2
2
2
2
(24k + 9)(24k + 10)(24k + 14)(24k + 15)(24k + 11)(24k + 12)2

24
120
36
48
6
+
−
−
+
−
.
(24k + 19)2 (24k + 20)2 (24k + 21)2 (24k + 22)2 (24k + 23)2
23 / 36

Andrew Mattingly, Glenn Wightwick, and the IBM BlueGene
For the actual computations, Jonathan Borwein and I turned to our colleagues Andrew
Mattingly and Glenn Wightwick at IBM Australia, who were willing to help with
programming and tuning. They received permission from IBM to use an IBM
BlueGene supercomputer for this purpose.

24 / 36

Our results — two double-checking runs each
1. Base-64 digits of π 2 beginning at position 10 trillion (a base-64 digit is a pair ofbase-8 digits):
75|60114505303236475724500005743262754530363052416350634|60114505303236475724500005743262754530363052416350634

2. Base-729 digits of π 2 beginning at position 10 trillion (a base-729 digits is a
triplet of base-9 digits):001|12264485064548583177111135210162856048323453468|12264485064548583177111135210162856048323453468

3. Base-4096 digits of Catalan’s constant beginning at position 10 trillion (abase-4096 digit is a quadruplet of base-8 digits):
0176|34705053774777051122613371620125257327217324522|34705053774777051122613371620125257327217324522

These runs required 22.1 billion CPU-seconds.
25 / 36

New calculation of base-16 digits of π
In December 2016, Daisuke Takahashi finished the computation of hexadecimal(base-16) digits of π beginning at position 100 quadrillion, or 1017 .
The run used Bellard’s formula (a variation of the BBP formula for π). Both the main
run and the verification run each required 320 hours on 512 nodes of a Fujitsu cluster
at the Joint Center for Advanced High Performance Computing (JCAHPC) in Japan.
The hexadecimal digits of π from position 1017 to 1017 + 15 are: A937EB59439E485E

26 / 36

Are the digits of π random?
Given a positive integer b, a real number α is normal base b if every m-long string ofdigits appears in the base-b expansion of α with limiting frequency 1/b m . It can beshown that almost all real numbers are normal base b, for all bases b.
These constants are widely believed to be normal base b, for all bases b:
I π = 3.14159265358979323846 . . .
I e = 2.7182818284590452354 . . .

√I

2 = 1.4142135623730950488 . . .

I log(2) = 0.69314718055994530942 . . .
I Every irrational algebraic number (this conjecture is due to Borel).

But there are no proofs of normality for any of the above — not even for b = 2 andm = 1 (i.e., equal numbers of zeroes and ones in the binary expansion).
Until recently, normality proofs were known only for a few constants, such as
Champernowne’s constant = 0.12345678910111213141516 . . . (normal base 10).
27 / 36

One (very weak) result for algebraic numbers
If x is algebraic of degree d > 1, then its binary expansion through position n musthave at least Cn1/d 1-bits, for all sufficiently large n and for some C that depends on x.
√√
Simple case: The first n binary digits of 2 must have at least n one bits.
In this case, the result follows by noting that the
one-bit count of the product of two integers is lessthan or equal to the product of the one-bit countsof the two integers. The more general result above
requires a more sophisticated approach.

1

0

0

1

×
1

1
1
0

0
1
1

1
0
1

1
1
0

1
1
1

However, note that these results are still a far cry from even single-digit normality.
I D. H. Bailey, J. M. Borwein, R. E. Crandall and C. Pomerance, “On the Binary Expansions of Algebraic
Numbers,” Journal of Number Theory Bordeaux, vol. 16 (2004), pg. 487–518.

28 / 36

BBP formulas and normality
Consider a general BBP-type constant (i.e., a formula that permits the BBP “trick”):
α =

∞
X
p(n)
,
b n q(n)n=0where p and q are integer polynomials, deg p < deg q, and q has no zeroes for n ≥ 0.
Richard Crandall (deceased 2012) and I proved that α is normal base b if and only if
the sequence x0 = 0, and


p(n)
,
xn =
bxn−1 +q(n)is equidistributed in the unit interval. Brackets {·} denote fractional part, as before.
Here equidistributed means that the sequence visits each subinterval (c, d) withlimiting frequency d − c.
I D. H. Bailey and R. E. Crandall, “On the random character of fundamental constant

expansions,” Experimental Mathematics, vol. 10 (Jun 2001), 175–190.

29 / 36

Two specific examples: log(2) and πConsider the sequence x0 = 0 and1
xn = 2xn−1 +nThen log(2) is normal base 2 if and only if (xn ) is equidistributed in the unit interval.
Similarly, consider the sequence y0 = 0 and
120n2 − 89n + 16
yn = 16yn−1 +512n4 − 1024n3 + 712n2 − 206n + 21
Then π is normal base 16 (and hence normal base 2) if and only if (yn ) is
equidistributed in the unit interval.
Sadly, we have not yet been able to prove equidistribution for either sequence.
Curiously, the sequence (yn ), when mapped to the 16 divisions of the unit interval,
appears to generate, digit by digit, the entire base-16 expansion of π, error-free.
30 / 36

A class of provably normal constants
Crandall and I also proved that the following constant is normal base 2:∞
X1
α2,3 =
n3 23nn=1= 0.041883680831502985071252898624571682426096 . . .10= 0.000010101011100011100011100011110110100001 . . .2
This constant was proven normal by Stoneham in 1971, but we have extended this tothe case where (2, 3) are any pair (p, q) of relatively prime integers, and also to alarger, uncountably infinite class.
The original proof is difficult, but a subsequent proof using a “hot spot lemma” (viaergodic theory) is quite simple.
1. D. H. Bailey and R. E. Crandall, “Random generators and normal numbers,” Experimental
Mathematics, vol. 11 (2002), 527–546.
2. D. H. Bailey and M. Misiurewicz, “A strong hot spot theorem,” Proceedings of the American
Mathematical Society, vol. 134 (2006), 2495-2501.
31 / 36

A (weak) hot spot theorem
The previous result on the normality of Stoneham numbers can now be proven much
more easily using this result from ergodic theory:
The (weak) hot spot theorem: Given the real constant α, if there exists some B suchthat for every subinterval [c, d) of [0, 1),#0≤j≤m {b j α} ∈ [c, d)≤Blim supm(d − c)m≥1then α is b-normal.
In other words, if α is not b-normal, then:
I There is some interval [c1 , d1 ) that is visited 10 times too often by shifts of thebase-b expansion of alpha;
I The is some other interval [c2 , d2 ) that is visited 100 times too often;
I There is some other interval [c3 , d3 ) that is visited 1000 times too often; etc.
However, one cannot conclude that these intervals are necessarily nested.
I L. Kuipers and H. Niederreiter, Uniform Distribution of Sequences, Dover, 1974, 77.

32 / 36

A strong hot spot theorem
In 2006 Michal Misiurewicz and I proved a stronger version of this result, using
methods of ergodic theory:
The strong hot spot theorem: Let 0.x1 x2 . . . xn be the base-b expansion of x out toposition n, so that [cn (x), dn (x)) = [0.x1 x2 . . . xn , 0.x1 x2 . . . (xn + 1)) is the n-long digit
interval containing x. If for every x in [0, 1),#0≤j<m {b j α} ∈ [cn (x), dn (x))<∞
lim inf lim sup
n≥1
mb −n
m≥1
then α is b-normal.
In other words: If α is not b-normal, then there is at least one “hot spot,” namely
some x ∈ [0, 1) such that shifts of the base-b expansion of α visit all sufficiently small
digit neighborhoods of x too often, by an arbitrarily large factor.
Conversely, if one can establish that there is no such “hot spot,” then α is b-normal.
I D. H. Bailey and M. Misiurewicz, “A strong hot spot theorem,” Proceedings of the American

Mathematical Society, vol. 134 (2006), 2495-2501.

33 / 36

Binary digits of α2,3 versus base-6 digits of α2,30.
0000101010111000111000111000111101101000
0100101111011010000100101111011010000100
1011111100110101101110100111100000011001
0100100010110000111111001101011011101001
1110000001100101001000101100001111110011
0101101110100111100000011001010010001011
0000111111101111001001001101111101010111
0111000010111001011010100110011100111110
0010100000001000011011011001000001010100
0100011110100011010010101100110001100000
1110101111111011110010010011011111010101
1101110000101110010110101001100111001111
1000101000000010000110110110010000010101
0001000111101000110100101011001100011000
0011101011111110111100100100110111110101
0111011100001011100101101010011001110011
1110001010000000100001101101100100000101
0100010001111010001101001010110011000110
0000111010111111111010011000011001111111
0001111101000000111101110011100010001001

0.
0130140430003334251130502130000001243555
0454322330115002435253205513523435410104
3000000000000000051411300540405554553031
4425043343510124134523511251421251345055
0354501505352205204434045215150510241155
2500425130051124454001044131150032420303
2130000000000000000000000000000000000000
0000014212034311121452013525445342113412
2402205253010542044235524110554150155204
3504145554003101453030335320025343404013
0124010445325434350214202043241502555510
1004043300045544114501031331451151014451
4123443342341240055131333504542353055315
1153501533452435450250055521453054234342
1530350125024205404135451231323245353031
5345523041150201542421211452015422225343
4034045053012332553444044310333244533214
1415014233454542412432031253400501341502
4551440430000000000000000000000000000000
0000000000000000000000000000000000000000 34 / 36

Binary digits of α2,3 versus base-6 digits of α2,3 : Random walks

I F. J. Aragon Artacho, D. H. Bailey, J. M. Borwein and P. B. Borwein, “Walking on real

numbers,” Mathematical Intelligencer, vol. 35 (2013), 42-60.
35 / 36

Summary
The study of π is a paradigm of experimental mathematics research:
I π is highly suitable for computational exploration.
I π can be explored visually as well as numerically.
I π has great public appeal, from grade school to serious scientific research.
I The question of the normality of π continues to fascinate (and frustrate!)

research mathematicians.

Thanks! This talk is available at:
http://www.davidhbailey.com/dhbpapers/dhb-pi-2017.pdf

36 / 36