Finding Pi to more decimal places

Hits: 178

Introduction

\pi, or the circle constant, is a number that is irrational. Irrational numbers consist of infinitely many decimal places which never repeat. This assures that you can never convert it into an exact fraction. But \pi is worse than that. It is also a transcendental number, meaning that it will never be the solution of a polynomial with integer coefficients and a finite number of terms. By the year 1400, the greatest advances in the accuracy of \pi was made by the Chinese, who were able to work it out to 6 decimal places.

Pi, known as the circle constant, or the Archimedian constant.
Circles.

You can come up with a lower and upper bound for \pi if you consider a circle of radius r, and take the area of an inscribed polygon of n sides whose corners touch the circles’ edge; and then the area of a circumscribed polygon, with the same number of sides, whose sides touch the circle’s edge. This was the tactic used by the Greek philosopher Archimedes about 2400 years ago, around 350 BC. Archimedes observed that the circle constant \pi can be had by dividing the circle’s area by the radius squared \left(\pi = \frac{A}{r^2}\right). But to get a circle area of any accuracy you had to have an accurate value for \pi. Since he could find the areas of polygons with much greater accuracy, Archimedes decided to circumscribe a polygon of n sides, whose sides touched the circle’s edge. But these areas were always too big. So, he also tried to find the area of the inscribed polygon of the same n sides. Of course, these areas were smaller than the circle’s area. He knew that \pi lay between these areas, or rather that: \frac{A_i}{r^2} < \pi < \frac{A_c}{r^2}, where A_i represents the area of the inscribed polygon and A_c represents the area of the circumscribed polygon. Another way of expressing this inequality may be more familiar to some: A_i < \pi r^2 < A_c. In fact, for a unit circle, this gets really simple: A_i < \pi < A_c, since r = 1 for a unit circle.

So, to increase the accuracy of \pi, all Archimedes had to do was increase the number of sides of the two polygons. You would achieve full accuracy for \pi if the number of sides goes to infinity, where you would then have: A_i = \pi r^2 = A_c. For over 2000 years, this was the holy grail for achieving \pi to perfect accuracy. Using this tedious method, Archimedes calculated polygon areas of up to 96 sides, so he was able to say \frac{223}{71} < \pi < \frac{22}{7}, or to at least estimate that \pi was about equal to 3.14 (decimals were also not known to the Greeks, so 3.14 is a modern estimate). Remeber that this was before computers or even before the discovery of irrational numbers. Square roots were known in his day, so his own calculations would involve fractions and nested square roots, which didn’t necesarily have integer or rational solutions, so an expression like 12\sqrt{6-\sqrt{3}} was left as it was. What was clear from his calculations, was that it was not possible to express \pi as a rational number, or an exact ratio of whole numbers. The good news is, the estimation of 3.14 is good enough for quick calculations today.

But of course, there was frustration in not being able to resolve \pi to a rational number, and so over the next 2000 years, across China, India, Persia, northern Africa and Europe, as we acquired more and more mathematical knowledge (including acquiring the decimal system), we could increase the number of sides of the polygon to a greater and greater number of sides. In France, François Viète used polygons of 393,216 sides in the year 1593, but could only achieve an accuracy of 10 decimals that way: 3.1415926536 (decimals were known by then). This would be a level of accuracy found on most inexpensive scientific calculators, and looked like a small prize given the number of sides in Viète’s polygon. But this was not state of the art in the early days of the Renaissance. Ludolph van Ceulen, a Dutchman in the late 1600s, who used a polygon of 2^{62} sides, or 4,611,686,018,427,387,904, or exceeding 4.6 quintillion sides. In these days before computers, or even electricity, this took van Ceulen 25 years as his lifetime achievement, but only generated \pi to 35 decimals: 3.141 592 653 589 793 238 462 643 383 279 502 88. There were those who surpassed this as well, but not by a lot.

The polynomial method of Newton, expanded to 12 terms, gives about the same accuracy as a scientific calculator, 3.141592654.
Sir Isaac Newton, or, what were you doing during the pandemic?

1666 was the year that Isaac Newton, a young Cambridge University undergrad student, was sitting at home, quarantining himself (as was everyone else) as the Bubonic Plague was raging through Europe. During his two years of social distancing, among his discoveries in optics, mechanics and Calculus was his method of computing the value of \pi. For this, the 23 year-old Newton used integration (which he called “fluxions”) on a polynomial with rational coefficients. He acquired accuracy to 9 decimals with only a 12-term polynomial. This was achieved by combining the circle formula with the binomial theorem. The circle formula for a unit circle is: x^2 + y^2 = 1

Solving for y in terms of x, we get y = \pm\sqrt{1-x^2}. It would be a good idea to find the area of the circle going from 0 to 1 (a quarter circle), so we would only need the “positive” version of this formula: y = \sqrt{1-x^2}. The integral would look like:

[
\int^1_0 \left(1 – x^2\right)^{\frac{1}{2}} dx
]

But because this only gets us the area of a quarter of the unit circle, we need to multiply the integral by 4 to get the exact answer \pi:

[
4\int^1_0 \left(1 – x^2\right)^{\frac{1}{2}} dx = \pi
]

While we have a good idea that this will get us \pi, it doesn’t yet yield any decimals, since \left(1 - x^2\right)^{\frac{1}{2}} must be expanded into a polynomial. For that, you have to break some rules about the Binomial Theorem.

The binomial theorem applies to the expansions of binomials such as: \left(1 + x\right)^n, where it is normally understood that n is a whole number greater than or equal to 1 (or n \ge 1 | n \in \mathbb{Z}). The coefficients of a general binomial expasion go as follows:

[
\left(1 + x\right)^n = {n \choose 0}+{n \choose 1}x+{n \choose 2}x^2 + … + {n \choose k}x^k + … + {n \choose n}x^n
]

Squaring x makes the expression (1 + x^2)^n. Its expansion looks like this:

[
\left(1 + x^2\right)^n = {n \choose 0}+{n \choose 1}x^2+{n \choose 2}x^4 + … + {n \choose k}x^{2k} + … + {n \choose n}x^{2n}
]

Now turning this into (1-x^2)^n, where we subtract instead of add, gives us a similar polynomial, but this time with alternating plus and minus signs:

[
\left(1 – x^2\right)^n = {n \choose 0}-{n \choose 1}x^2+{n \choose 2}x^4 – … \pm {n \choose k}x^{2k} \pm … \pm {n \choose n}x^{2n}
]

So, the strange thing Newton did with this is that he broke the rule that n be a natural number. Instead, for the purpose of finding the area of a unit circle (that is, pi, to as many decimals as possible), he had to go back to the formula for combinations, that is, {n \choose r} = \frac{n!}{(n-r)!r!} to see what would happen if n = 1/2. The reason for this is because it fit in with the circle formula of (1-x^2)^{1/2}.

n!, called “n-factorial”, is normally the whole number n, multiplied by all of its integer predecessors down to 1, as in: 3! = 3\times 2\times 1 = 6. You stop multiplying when you reach 1. But what if n = 1/2? For normal counting numbers, subtracting 1 gets you to the next lower whole number. At some point, you will reach 1 if you keep subtracting 1. But if n = 1/2, then subtracting 1 will always yield a fraction, and you never reach 1. In fact, you proceed forever into the negative numbers. Thus, a number like 1/2, when the formula is applied, leads to: \frac{1}{2}\times \frac{-1}{2}\times \frac{-3}{2}\times \frac{-5}{2}\times ... going forever. Is this useful?

It turns out, it does have  a use when placed in the formula for n \choose r, because, for example, the first term is 1, due to n \choose 0. The next term has the coefficient:

[
{n \choose 1} = \frac{n(n-1)(n-2)…(3)(2)(1)}{(n-1)(n-2)…(3)(2)(1)\times 1} = n
]

So, after all that, if n = 1/2, then this coefficient becomes 1/2 so that the second term is \frac{-x^2}{2}. The third term has the coefficient:

[
{n \choose 2} = \frac{n(n-1)(n-2)…(3)(2)(1)}{(n-2)(n-3)…(3)(2)(1)\times 2} = \frac{n(n-1)}{2} = \frac{n^2 – n}{2}
]

So, for n = 1/2:

[
\frac{\left(\frac{1}{2}\right)^2-\frac{1}{2}}{2} = \frac{\frac{1}{4}-\frac{1}{2}}{2} = \frac{-1}{4}\times\frac{1}{2}=\frac{-1}{8}
]

The first 12 terms in the infinite expansion are given in the above illustration at the start of this section.

How much of \pi could we get “on the cheap” by performing integration on only the first 7 terms?

\begin{align}
\pi &= 4\int_0^1 \left( 1-\frac{x^2}{2}+\frac{x^4}{8}-\frac{x^6}{16}+\frac{15x^8}{384}-… \right) dx \
&= 4\left( x-\frac{x^3}{6}-\frac{x^5}{40}-\frac{x^7}{112}-\frac{5x^9}{1152}-\frac{7x^{11}}{2816}-\frac{7x^{13}}{1024}-… \right) \Big|_0^1 \
&= 4\left[1-\frac{1}{6}-\frac{1}{40}-\frac{1}{112}-\frac{5}{1152}-\frac{7}{2816}-\frac{7}{1024}-…\right] \
&= 3.14297
\end{align}

7 polynomial terms is all it takes to get the same accuracy Archimedes once had when he attempted to average out the areas of a 96-sided polygon both inscribed and circumscribed around a circle.

Integration will measure the area of the part of the circle from x=0 to x=0.5, down to the x axis. This shape can be thought of as a 30 degree sector combined with a right triangle.

There is an even better calculus-based solution, which has greater accuracy in fewer terms. What you want to do is to calculate the area under the same curve over 0 to 1/2 instead of 0 to 1. Substituting x = 1/2 allows each term to shrink faster, and thus to obtain accurate results with less work. Newton knew that on a unit circle centered at the origin, x = 1/2 is actually \cos(60^\circ). The y-value, by extension, is \sin(60^\circ) = \frac{\sqrt{3}}{2}. The area he is calculating describes a sector of the circle 30^\circ  from the perpendicular (this is 1/12th of the total area of the circle, and so its area is \pi/12). Added to this sector is a right triangle whose base is 1/2 and whose height is \frac{\sqrt{3}}{2}. The triangle would have an area of \frac{\sqrt{3}}{8} by the area formula A=\frac{1}{2}bh. It would be useful to get only the sector and subtract out the right triangle from the integral. So we will illustrate the first 5 terms of what Newton did:

\begin{align}
\pi &= 12 \left(\int_0^{\frac{1}{2}} (1-x^2)^{1/2} dx – \frac{\sqrt{3}}{8}\right) \
&= 12\int_0^{\frac{1}{2}} (1-x^2)^{1/2} dx – \frac{12\sqrt{3}}{8} \
&=12\left( x-\frac{x^3}{6}-\frac{x^5}{40}-\frac{x^7}{112}-\frac{5x^9}{1152}\right)\Big|_0^{\frac{1}{2}} – \frac{12\sqrt{3}}{8} \
&= 12\left(\frac{1}{2}-\frac{1}{48}-\frac{1}{1280}-\frac{1}{14336}-\frac{5}{589284}\right)-\frac{12\sqrt{3}}{8} \
&= 3.14161
\end{align}

That is an accuracy within \frac{1}{50,000}, which is impressive for a mere 5 terms. If you wanted more accuracy, just add more terms. Applying 7 terms as we did above would give us \pi\approx 3.141585, or an accuacy of within \frac{7}{1,000,000}.

Modern \pi calculations reflect the power of computers rather than the power of human calculation

Not much progress was observed after that until the invention of computers. So by 1949, the invention of the ENIAC computer gave us 1,120 decimal places, taking 70 hours. By 1973, the CDC-7600, or Cray computer was the first device to give more than 1 million decimal places in 23.3 hours. August 1989 was when the 1 billion digit barrier was broken, using an IBM mainframe. \pi became known to 1 trillion digits in 2002, by a Japanese team using a Hitachi supercomputer, after taking 25 days. Since then the orders of magnitude and computational time appeared to have plateaued somewhat. The latest world record, according to The Guiness Book of World Records, was in November 2020, when 50 trillion digits were found over a period of 8 months of computation. If printed in a normal sized font on letter-sized paper, it would require nearly 28 billion pages to contain the constant at 1800 characters per page. The supercomputer used 4 Intel Xeon processors (15 cores per CPU) running at 2.5GHz each, had 320 GB DDR3 RAM, and 336 terabytes across 60 hard drives, most of them used for computation.

As of 2021, a claim has been made of 62 trillion digits in a little over 3 months. The Swiss-based computer was able to find more digits, faster because of slightly more advanced hardware: Two AMD Epyc 7542 32-core processors at 2.9 GHz each, 1 terabyte of RAM, and 510 terabytes of storage across 38 hard drives, most used for swapping data with the RAM. Their “book” would be over 34 billion pages long if the digits were printed out. It was never made clear what kind of RAM was used. It is assumed they used DDR4.

The current computation of \pi to ever more digits has become a quasi-annual event; an opportunity for computer companies to flex their technological muscle and promote their brand. These impressive calculations are not so much human achievements in the sense of Newton or van Ceulen, but are really achievements of computer engineers and software programmers.

But how many digits do we really need? Even to express the diameter of the universe to within a Planck length (1.6162\times 10^{-35} meters), we would only require 62 digits of \pi.

This article is based on a video by YouTube vlogger Veritasum, plus some additional sources, such as the University of St. Andrews department of Math History.

What are your thoughts?