Calculate Pi with C++
Calculate Pi with C++
March 22, 2021
Keith Greiner
This essay is about the calculation of pi, also known as "π" with code examples in C++ XCode for Mac. The programs may be adapted to other versions of C++, but some, including C++ in Microsoft Visual Studio 2010 have limits on the number of digits that are displayed. A work-around for the limits of precision may be found in the spigot functions that are described below.
This essay has the following sections:
1. Introduction
2. Observations Regarding Number Precision
3. Pi Calculations Based on Geometry
4. Pi Calculations Based on Infinite Series
5. Pi Calculations Based on Spigot Functions
6. Pi Calculations Based on Ratios
7. Two Distributions of Pi Digits
1. Introduction
In 1873, William Shanks completed his calculation of 707 places of pi. Given the complexity of the calculations, that was an amazing accomplishment. Shanks completed the calculations by hand with no calculator and no computer. One source suggested the numbers were published in three installments over a 20 year period. Always attentive to accuracy, Shanks worked for months to calculate additional digits each morning and then would check his work each afternoon. Some of his calculations were checked by mathematician William Rutherford. Shanks work was considered the standard until, 1944 when D. F. Ferguson used a mechanical calculator to check the numbers and discovered Shanks made a mistake when calculating the 527th digit. So basically, Shanks was nearly 75% correct. According to Brian Hayes, the first of several errors appears to be the clerical omission of a zero, perhaps after a period of time when the project had been set aside. (Berggren, Lennart, & Jonathon and Peter Borwein. Pi: A Source Book. New York, Springer-Verlag, 1997, page 627, and http://sites.math.rutgers.edu/~cherlin/History/Papers2000/wilson.html. and Hayes, Brian, Pencil, Paper, and Pi, American Scientist, https://www.americanscientist.org/article/pencil-paper-and-pi)
Shanks calculations were based on a pi formula discovered by John Manchin in 1706. Manchin’s formula uses the arc-tangent of two ratios.
Today, the calculation of pi, using the Manchin formula, is a simple task. In C++, some workable code is….
If you run this program, however, the best result you will obtain will have only 16 correct digits including the 3 at the beginning, but not the decimal point “.” . So the mathematicians from the 18th and 19th centuries, working for years by hand, calculated pi to greater accuracy than the computer you are using right now. That’s because the designers of computer math systems decided that today’s calculations don’t need the number-of-digits level of accuracy calculated by the likes of Manchin and Shanks. The modern mathematicians and computer engineers determined that numbers like 3.402823466e+38 are best represented by just a few digits and a large exponent, which means you don't get all the digits to the right of the decimal but you do get the essential answer. For many math problems that will do. When calculating pi and other constants, however, we find that floating point system is far far less than what is needed. It is also why computer systems store a value of pi as a constant that has its roots in those 19th century calculations by Shanks. In C++ that constant can be obtained by invoking M_PI. Now it is possible to calculate pi to greater digits than ever imagined to be possible in the last century, but to do so, you need a specialized big-number add-in or one of the spigot formulae that are described below. Without those, you’re better off working to calculate the value by hand. For those of you who are designing future computer systems and compilers, how about you include functions and class methods for big math problems like the ones shown below.
This following sections explore some of the programs that can be used, today, to calculate pi. The discussion focuses on programs written in C++ rather than commercially available programs, because C++ allows you to observe details of how the formulae can translate to the mathematical steps used to accomplish those calculations. C++ also allows you to customize the output as you use computational experimentation to better understand the behavior of numbers.
2. Observations Regarding Number Precision
When evaluating a calculation of pi, I used three sources:
1) A list of one million digits of pi that were posted on an MIT.edu website. The post has since disappeared from Google and Bing search results.
2) A list of 10,000 digits from the University of Utah at http://www.math.utah.edu/~alfeld/math/pi.html available on June 2, 2019. My reformatted text file of the 10,000 numbers may be found at this link.
3) The constant M_PI encoded into the C++ compiler and accessible via one of two options. With C++ in XCode, the following include statement is required.
#include <cmath>.
In Visual Studio 2010, M_PI is available with the following two statements;
#define_USE_MATH_DEFINES
#include<math.h> .
The first 50 digits of each of the three sources are shown below. In C++ for Visual Studio 2010, the 16th digit is the last digit showing a number. All digits to the right of the 16th are set to zero.
MIT = 3.14159265358979311599796346854418516159057617187500
Utah = 3.14159265358979323846264338327950288419716939937510
M_PI = 3.14159265358979311599796346854418516159057617187500
The values of pi that are returned by the programs that are presented with this essay, can depend on which C++ compiler is used. Below is a comparison of output from the same program run first on a Mac, with an XCode compiler, and second, with a Visual Studio 2010 C++ compiler. The program is CPP_Test_Long_Double. The program presents numbers as "long double" data types and "string" data types. The first long double, testValue, is defined to be 0.000000000000000000000000000001which is a one in the 30thposition to the right of the decimal. The number is then printed with the following statement:
printf("testValue = %2.60Lf \n",testValue);
The statement presents a total of 60 digits to the right of the decimal. When we look at the 60 digits, we see that on the Mac, there are extraneous values beginning at position 47, while on the Visual Studio 2010 C++ output there is an extraneous 1 at position 46.
When the testValue is multiplied by 2, so are the extraneous values. In the second section, the values of pi are defined as strings. In the third section, they are converted to long double values and then printed with statements that, again, have a large number of digits to the right of the decimal. On the Mac, the string’s number values completely fill the space specified in the printf statement. In Visual Studio, however, the numbers are all truncated at the 16thdigit, and the remaining characters are set to zero. Below are lists of output from the same program run first in Mac XCode, and second in Visual Studio C++ 2010.
CPP_Test_Long_Double.xcodeproj with Mac XCode
Long Double
testValue = 0.000000000000000000000000000001000000000000000083336420607586
testValue multiplied by 2
testValue = 0.000000000000000000000000000002000000000000000166672841215172
String
FIT no mod= 3.14159265358979323846264338327950288419716939937510582097494459230781
MIT = 3.14159265358979323846264338327950288419716939937510582097494459230781
Utah = 3.14159265358979323846264338327950288419716939937510582097494459230781
Long Double
CPP_PI = 3.141592653589793115997963468544185161590576171875000000000000
M_PI = 3.141592653589793115997963468544185161590576171875000000000000
Long Double
LD_FIT = 3.14159265358979311599796346854418516159057617187500
LD_MIT = 3.14159265358979311599796346854418516159057617187500
LD_Utah = 3.14159265358979311599796346854418516159057617187500
Done
Program ended with exit code: 0
CPP_VS_Test_Long_Double with Visual Studio 2010
Long Double
testValue = 0.000000000000000000000000000001000000000000000100000000000000
testValue multiplied by 2
testValue = 0.000000000000000000000000000002000000000000000200000000000000
String
FIT no mod= 3.14159265358979323846264338327950288419716939937510582097494459230781
MIT = 3.14159265358979323846264338327950288419716939937510582097494459230781
Utah = 3.14159265358979323846264338327950288419716939937510582097494459230781
Long Double
CPP_PI = 3.141592653589793100000000000000000000000000000000000000000000
M_PI = 3.141592653589793100000000000000000000000000000000000000000000
Long Double
LD_FIT = 3.14159265358979310000000000000000000000000000000000
LD_MIT = 3.14159265358979310000000000000000000000000000000000
LD_Utah = 3.14159265358979310000000000000000000000000000000000
Done
Program ended with exit code: 0
3. Pi Calculations Based on Geometry
The great advancement that Archimedes contributed to our modern calculation and use of pi, was his application of a generalized formula. Knowing Pythagoras’ formulas for the three sides of a right triangle, and knowing the size of an angle that originates at the center of a circle, Archimedes reasoned that he could slice up a circle into smaller and smaller triangular parts having a side that lies along the circle, and a center angle that is also the center the of the circle. He could keep making the circle-parts smaller and smaller until the sums of the side that lies along the corresponding circle, closely estimated the circumference of a circle. It sounds complicated, but actually it is straight-forward when the calculations are done.
Today, students of high school geometry should be able to do this calculation with the aid of a calculator or computer. However, in circa 200 BC, it was a challenging calculation and a great leap forward for mathematics. Archimedes appears be the first to do it in that manner. Other mathematicians already estimated pi, but Archimedes’ method was a major step forward as an application of Pythagoras’ formulas.
Although Archimedes calculations were highly advanced for their time, they lacked the precision of today’s decimal system. As a result Archimedes was able to find only a range of estimates. He concluded that the actual value of pi was somewhere between a high of 3 1/7 (same as 22/7) and a low of 3 10/71 (same as 223/71). Expressed in our modern decimal system, the high would be 3.14286 and the low would be 3.14084. The two are correct to two digits to the right of the decimal. The average of these high and low values is 3.14185, which is correct to three digits to the right of the decimal.
The project described in this essay started with an attempt to replicate Archimedes calculations, using a computer and the C++ programming language, and continued through several iterations of calculations. And so we begin with the program that somewhat replicates Archimedes calculations.
CPP_Calculate_Pi_from_Pythagorean was written for C++ in XCode 10.1. The program calculates pi using Pythagorean formulas. It uses a similar method to the one used by Archimedes with the advantage of the modern-day algebraic terms, the modern day decimal system, a computer, and C++. The program stores each digit of pi in a vector of class objects. The program uses two functions that replicate Archimedes’ method. Each function calculates the circumference of a circle that has a radius of 1 and a diameter of 2. As the circle is divided into smaller and smaller parts, the accuracy of pi becomes more and more precise. The two functions are;
1) circle_inscribed_in_polygon(std::vector<pi_Data0> & pi_List ), and,
2) polygon_inscribed_in_circle(std::vector<pi_Data0> & pi_List)
It is easier to calculate the circumference of a polygon because the side that lines up with the circumference is made of many short, straight lines. As the number of sides increases, the length of each side decreases, and the polygon looks more and more like a circle. As the number of sides increases the side of each polygon that lies closest to the circle becomes indistinguishable from the circle. When we sum the polygon pieces, we get a better and better estimate of the actual circumference and the best possible estimate of pi. The idea of basing a sum of increasingly smaller segments of a circle tells us the value of pi and foreshadows Newton’s and Leibniz’s calculus,
In this program, the number of polygon sides, doubles with each iteration. After 21 iterations, there are 2,097,152sides to the polygons and the two polygons (inside and outside the circle) yield the same pi value to the third digit to the rightof the decimal. With greater numbers of iterations, the values returned by the two functions soon become indeterminate and non-numbers. There is a good explanation of the method at http://www.physicsinsights.org/pi_from_pythagoras-1.htmland another at https://itech.fgcu.edu/faculty/clindsey/mhf4404/archimedes/archimedes.html.
Following are estimates of pi that came from my implementation of the C++ program. The high value comes from the circle outside a polygon, and the low value comes from the polygon inside a circle. The polygon inside a circle agrees with the C++ constant, M_PI to the seventh digit to the right of the decimal.
Archimedes Estimates from Geometrical Calculations
Iteration = 20
Circle Sides = 2,097,152
High: Circle in Polygon = 3.14159268140792846680
Low: Polygon in Circle = 3.14159266429396318867
Average = 3.14159267285094582784
M_PI = 3.141592653589793115997963468544
While those calculations use Archimedes’ geometric method, he, of course, did not have access to the same decimal mathematics that we have, today. If we limit ourselves to the ratios that Archimedes reported, (as input) we have the fractional values shown below. The high and low values agree with the C++ constant, M_PI, to two digits to the right of the decimal. The average agrees with M_PI to three digits to the right of the decimal.
Archimedes Estimates from Ratios Converted to Decimal Form
High 3+(10/71) = 3.14285714285714285724
Low 3+(1/7) = 3.14084507042253521119
Average = 3.14185110663983903421
M_PI = 3.141592653589793115997963468544
The program CPP_Calculate_Pi_from_Pythagorean helps us understand that the decimal notation system as applied in the C++ language has limits. As the number of sides increases beyond 2,097,152, the calculations start to overflow, while the number of correct digits remains at only three to five digits to the right of the decimal. That is hardly sufficient. The problem is an outcome of the long double numbering system used in Microsoft Windows and Apple Mac operating systems. Those who wish to calculate to greater precision need to use a programming tool that allows fractional numbers to be calculated to greater precision. I chose not to use a C++ library that provides increased floating-point precision because the download page included warnings about how the library might have been compromised by nefarious operators. It would be appropriate of the providers of C++ compilers would include increased floating point precision libraries in their standard compilers. Interestingly, the floating-point limit for modern computers was not a limitation for people who calculated pi by hand. Apparently in 1874 William Shanks calculated pi to 707 digits (by hand). Luckily for him, the mistake at digit 527 wasn’t discovered until 1944.
4. Pi Calculations Based on Infinite Series
Between 1340 AD and 1425 AD Madhava of Sangamagrama (an Indian mathematician) is said to have discovered that pi can be calculated from an infinite series. While it is said that Madhava taught the method in the 15thcentury, it is also said that Leibniz and Gregory also discovered the formula in the 17thcentury. Leibniz’s and Gregory’s use of an infinite series takes a step away from the geometrical calculation of Archimedes, but allows for additional floating point accuracy: up to a point. Most of the modern pi calculations are applications of the infinite series and not the geometric.
Following is the Madhava, Leibniz, Gregory formula. Note that the exponent of the (-1), is (n-1). The purpose of the exponent is to change the sign of the ratio so that ratios are alternatively added and subtracted.
Operationally, the formula becomes the following, and the C++ code that is also shown below. Notice in C++ that the range of n is from 0 to infinity. That is the same as 1 to infinity when the exponent of (-1) is rewritten as n-1 because n-1 is zero for the value n=1.
The C++ code is…
E_Val0 = 0;
for ( long long i = 0; i < n; i++)
{
numerator = pow( -1.0, (double) i );
denominator = (( 2.* (double) i ) + 1.0);
E_Val0 += numerator/denominator;
pi = E_Val0 * 4;
} // End for ( long long i ….
My interpretation of the formula is presented in program CPP_Calculate_Pi_Madhava_Leibniz_Gregory and yields a pi value shown below, after a million iterations of n. Below the estimate, is the C++ constant for pi that is in the #include <cmath> library as M_PI and below that is a value published on an MIT site. (The file pi_dec_1m.txt was found on an MIT web site when this project started, but has since been absent from Google and Bing searches of the MIT site). The three numbers have agreement to the 16thplace to the right of the decimal.
pi est. = 3.14159165358979318525417534502253147365991026163101
M_PI = 3.14159265358979311599796346854418516159057617187500
MIT = 3.14159265358979311599796346854418516159057617187500
If the program is run in Visual Studio 2010, Visual C++ The output will be truncated at the 16thdigit and all numbers to the right of the 16thdigit will be presented as zeros. Following is an example of the Windows output.
pi est. = 3.14159165358979310000000000000000000000000000000000
XCode C++ 10.1 displays a higher level of precision and you can see as many digits past the decimal as you can see in your output using the statement: printf("pi est. = %2.50Lf \n",pi);. Still, the estimate of pi matches other sources to only 16 digits. Which is correct? Perhaps the better way to think of it is, that neither is either right or wrong: they’re just different, and reflective of the capabilities of the software being used.
5. Pi Calculations Based on Spigot Functions
There is a substantial limitation of the pi calculations, shown above, that are based on formulas for geometry, ratios, and infinite series. All are limited by the number of digits to the right of the decimal that can be manipulated by the numbering systems in today’s computers. C++ compilers, for example, perform calculations based on the IEEE-754 standard for floating point arithmetic, and that limits the number of digits to the right of the decimal. When data are stored as the largest available floating point number, a “long double,” the largest possible decimal will have 16 digits and when there are “9”s in each digit the largest decimal number is 0.9999999999999999. When pi is represented this way, the most precise number that is possible becomes “3.141592653589793100000000000000000000000000000000000000000000” as represented in the Windows Visual Studio 2010 C++ output. The last two most precise numbers are 31 and the remaining digits are zeros. On a Mac C++ program, the output that I obtained had correct values to 48 digits as shown below. 3.141592653589793115997963468544185161590576171875000000000000. However, this was deceiving. While the number can be represented to 48 digits and even to several hundred digits, the limit of precision that one can count on still has a 16 digit limit of precision.
The work-around for this limitation is to use a function that calculates each digit individually. Such functions are known as “spigot” functions, after a name applied by Rabinowitz and Wagon in the presentation of their method. Their paper, “A Spigot Algorithm for the Digits of Pi,” may be found at http://stanleyrabinowitz.com/bibliography/spigot.pdf and https://www.maa.org/sites/default/files/pdf/pubs/amm_supplements/Monthly_Reference_12.pdf.A C program that applies the Rabinowitz and Wagon algorithm may be found on the Florida Institute of Technology website at https://cs.fit.edu/~mmahoney/cse3101/pi.c.txt. The C++ program, CPP_Calculate_Pi_Rabinowitz_and_Wagon_FIT_Modified, is an adaptation of the Florida Institute of Technology (FIT) program. For the adaptation, I added a loop that calculates values of pi with various values of “n”, and a timer that reports the amount of time required for calculations. The output was placed into a C++ vector and was compared, digit by digit, to digits from the MIT and Utah standards. If the “n” value was set to 20, the program generated 70 numbers, and the first 10 were exact matches with the standards. If “n” was set to 40,000 (entered as 40000) the program generated 159,990 numbers and the first 12,049 were exact matches. At 40,000 the proportion of exact matches to the numbers generated was 0.075 (7.5%). For 12,049 matches, the program required 65.558 seconds of processing time. The actual calculation time varies.
Following is output from an “n” setting of 200. The list has a total of 790 characters with the first 66 digits highlighted.
3.14159265358979323846264338327950288419716939937510582097494459269133196261242295273115610496445264996640759087874802953643203410319967987203465827380435527187072554703130495611202643885282529529528900882028602563625060782238211501502844049253713154412304570667918519410130313280277090426420982178837069996183307061842089531760771665136645899946634838758952801283372308241032329612043444642807275738092153837627880468657321438910277203550447650364045798198305260657060135096927439674297560624338590559950654689639911261367184150180139420795373982097722126453300173810946917962611182912917196863448096859515381498031949770742924107637490859490376004471760585006291487955159685693488616304421140821896361789742061065540859013402297800248466707077399081763900584977497541643749998434263138140
After 30 runs of the program with n = 40,000
The following summary table was produced to show the number of seconds required to complete the generation of 159,990 digits of Pi with 12,049 matching one of the standard lists.
Another spigot program by David Bailey, Peter Borwein, and Simon Plouffe returns hex values of pi, with each selection independent of previous calculations. I attempted to convert the hex output to decimal values, without success beyond only 15 digits. I am confident my conversion function is correct, and the solution remains an unknown. For further information, try this link: davidhbailey.com
A third spigot program of note is attributed to mathematician Dik T. Winter. This program is published in two forms: first is a very brief obfuscated code that produces 800 digits of pi, and second is an expanded C/C++ form that I used to produce 16,277 digits. My comments and the Winter code may be found at this link.
6. Pi Calculations Based on Ratios
There are so many individual mathematics calculations involved in pi that it becomes a cumbersome and time-consuming task, unless one is using a computer. So it is no wonder that, over the centuries, people have attempted to find shortcuts for the calculation. Some of the shortcuts are simple division problems that have a never-ending result. Below are a few of them.
In about 15 B.C., Vitruvius calculated pi as the ratio of 25/8 = 3.25. Ptolemy suggested 377/120 = 3.14166667.
Some other ratios include the following:
256 / 81 = 3.1604938272
339 / 108 = 3.1388888889
377 / 120 = 3.1416666667
3927 / 1250 = 3.1416000000
355 / 113 = 3.1415929204
(2143 / 22)^0.25 = 3.1415926526
16 / 5 = 3.200 (Indiana General Assembly, bill number 246 of the 1897 Assembly)
And of course, there are Archimedes' ratios that define the low and high points of a range for pi.
3 1/7 = 22/7 = 3.14285714285714
3 10/71 = 223/71 = 3.14084507042254
Some of these values are calculated in program CPP_Calculate_Pi_with_Historical_Ratios.
7. Two Distributions of Pi Digits
Below is a table showing the distribution of digits of Pi from the Utah list of 10,000 digits and the MIT list of one million (actually, I found one less than a million) digits of pi.