I started off calling this project "Calculating Pi" but soon realised that I needed to rename it "Estimating Pi". Pi is an irrational number starting off 3.14159 and then carrying on for an infinite number of digits with no pattern which anybody has ever discovered. Therefore it cannot be calculated as such, just estimated to (in principle) any number of digits.
So far π has been calculated to 22,459,157,718,361 digits (nearly twenty two and a half trillion) which, by coincidence, is the number of different ways there are of calculating it.
OK, that's a slight exaggeration but there are many ways of getting the job done, ranging from the ancient 22/7 to recently discovered fast but complex methods such as the Chudnovsky Algorithm. These generally calculate all the digits up to a certain point but there are also a few so-called spigot algorithms which calculate a given digit without calculating all the preceding ones first.
There is plenty of information around regarding the number of digits of π which are actually necessary for practical purposes. This caught my eye though from NASA. You can calculate the circumference of the Universe to the accuracy of a hydrogen atom using just 40dp which is the number reached over three centuries ago! And you only need 62 digits to calculate the circumference to the Planck length, the smallest meaningful unit of length which you can think of as a quantum of space.
Trillions of digits aren't therefore useful in their own right but the effort put in to their calculation no doubt benefits mathematics and computer science.
In this project I will code a few of the simpler methods to give a decidedly non-rigorous introduction to what is actually a vast topic. If you want to do something rather more serious than playing with my code you can download the application called y-cruncher used to break the record here.
The long double Type
All the calculations in this project will be carried out using the long double type. Exactly how much precision that gives is dependent on the compiler, but GCC on x86 uses 80 bits, although if you run sizeof on the type it will give 16 (128 bits) to use a multiple of the data size. (On a 32 bit machine 32 x 3 = 96 bits might be used.) I will start off the code with a simple function to output a few interesting facts about long double before getting stuck into the π-calculating stuff. When we get to that you'll see the maximum guaranteed precision is 18dp so that is what I'll use as the gold standard, and we will see how close to that we can get using various algorithms within the limits of the data type used.
M_PI
There is a common but non-standard #define called M_PI which you might have in your math.h file. If you don't take a look here to fix the problem. Needless to say, if you actually need to use π don't mess around calculating it, just make sure it is #defined somewhere. That's what Charles Babbage did - this punch card is his equivalent of #define M_PI.
For the full story go here
Starting to Code
Create a new folder and within it create a single file called estimatingpi.c, which will contain all the code for this project. You can download the code as a zip or clone/download from Github if you prefer.
Source Code Links
The main Function
Open estimatingpi.c and type or paste this code.
estimatingpi.c part 1
#include<stdio.h> #include<stdlib.h> #include<math.h> #include<string.h> #include<float.h> #define PI_STRING "3.141592653589793238" #define RED "\x1B[31m" #define GRN "\x1B[32m" #define RESET "\x1B[0m" //-------------------------------------------------------- // FUNCTION PROTOTYPES //-------------------------------------------------------- void long_double_info(); void print_as_text(long double pi); void fractions(); void francois_viete(); void john_wallis(); void john_machin(); void gregory_leibniz(); void nilakantha(); //-------------------------------------------------------- // FUNCTION main //-------------------------------------------------------- int main(int argc, char* argv[]) { puts("-----------------"); puts("| codedrome.com |"); puts("| Estimating Pi |"); puts("-----------------\n"); //long_double_info(); //fractions(); //francois_viete(); //john_wallis(); //john_machin(); //gregory_leibniz(); //nilakantha(); return EXIT_SUCCESS; }
After the #includes I have #defined π as a string to 18dp, which will be used later to check results digit by digit. (The 19th digit is 4 therefore the 18th doesn't need to be rounded up.) There are also a few #defines for printing in colour which I'll use later.
The rest is very straightforward - just some function prototypes and calls to those functions in main, all of which are commented out at the moment so we can implement and run them one at a time.
Info on The long double Type
The following function (which you can skip if you find it too boring!) just prints out a few pieces of information from float.h on the long double type.
estimatingpi.c part 2
//-------------------------------------------------------- // FUNCTION long_double_info //-------------------------------------------------------- void long_double_info() { printf("sizeof(long double) %ld\n", sizeof(long double)); printf("LDBL_MIN %Lg\n", LDBL_MIN); printf("LDBL_MAX %Lg\n", LDBL_MAX); printf("LDBL_DIG %d\n", LDBL_DIG); }
Uncomment long_double_info in main, and compile and run with this...
Compile and Run
gcc estimatingpi.c -std=c11 -lm -o estimatingpi
./estimatingpi
...which should give you something like this. Old or esoteric hardware or compilers might differ.
Program Output
-----------------
| codedrome.com |
| Estimating Pi |
-----------------
sizeof(long double) 16
LDBL_MIN 3.3621e-4932
LDBL_MAX 1.18973e+4932
LDBL_DIG 18
Printing Pi as Text
A quick utility function before we get into actually estimating pi, which takes a long double and prints it out below the definitive value for comparison, with matching digits in green and non-matching ones in red.
estimatingpi.c part 3
//-------------------------------------------------------- // FUNCTION print_as_text //-------------------------------------------------------- void print_as_text(long double pi) { char pi_string[21]; sprintf(pi_string, "%1.18Lf", pi); printf("Definitive %s\nCalculated ", PI_STRING); for(int i = 0, l = strlen(pi_string); i < l; i++) { if(pi_string[i] == PI_STRING[i]) { printf(GRN "%c" RESET, pi_string[i]); } else { printf(RED "%c" RESET, pi_string[i]); } } puts("\n"); }
Firstly we sprintf the long double to a char array, print the string value #defined earlier, and then iterate the digits of the calculated value. These are printed in green or red depending on whether they are correct.
Estimating Pi Using Fractions
The first of the six methods we will use is the simplest and just consists of dividing one number by another, or to look at it another way converting fractions to decimals
estimatingpi.c part 4
//-------------------------------------------------------- // FUNCTION fractions //-------------------------------------------------------- void fractions() { long double pi; pi = 22.0L / 7.0L; puts("22/7\n===="); print_as_text(pi); pi = 333.0L / 106.0L; puts("333/106\n======="); print_as_text(pi); pi = 355.0L / 113.0L; puts("355/113\n======="); print_as_text(pi); pi = 52163.0L / 16604.0L; puts("52163/16604\n==========="); print_as_text(pi); pi = 103993.0L / 33102.0L; puts("103993/33102\n============"); print_as_text(pi); pi = 245850922.0L / 78256779.0L; puts("245850922/78256779\n=================="); print_as_text(pi); }
The variable pi is set to each of six different fractions and then printed using the print_as_text function. Uncomment fractions in main and compile/run again to get this.
Program Output
22/7 ==== Definitive 3.141592653589793238 Calculated 3.142857142857142857 333/106 ======= Definitive 3.141592653589793238 Calculated 3.141509433962264151 355/113 ======= Definitive 3.141592653589793238 Calculated 3.141592920353982301 52163/16604 =========== Definitive 3.141592653589793238 Calculated 3.141592387376535774 103993/33102 ============ Definitive 3.141592653589793238 Calculated 3.141592653011902604 245850922/78256779 ================== Definitive 3.141592653589793238 Calculated 3.141592653589793160
As you can see the longer fractions give us ever increasing accuracy. The second fraction looks odd as it has a seemingly correct digit after a number of wrong ones but that it just a coincidence. It is quite impressive that pi can be calculated to such accuracy with a simple long division, although I suspect the later fractions were discovered after pi had already been calculated to that accuracy using more complex methods.
Francois Viete
The following formula was discovered by French mathematician Francois Viete in 1593.
These are just the first three terms of an infinite series and already grow to somewhat monstrous proportions. However, if you look closely you can see that each numerator (the expression on the top of the fraction) is just the square root of 2 plus the previous numerator. This even applies to the first term if you assume the previous numerator is 0, so let's use that insight to implement Monsieur Viete's infinite product.
estimatingpi.c part 5
//-------------------------------------------------------- // FUNCTION francois_viete //-------------------------------------------------------- void francois_viete() { puts("Francois Viete\n=============="); const int iterations = 28; long double numerator = 0.0L; long double pi = 1.0L; for(int i = 1; i <= iterations; i++) { numerator = sqrt(2.0L + numerator); pi*= (numerator / 2.0L); } pi = (1.0L / pi) * 2.0L; print_as_text(pi); }
The iterations variable specifies how many terms we will calculate, and the numerator is initialized to 0 as mentioned above. The pi variable is initialised to 1 so we can multiply the value of the first term by it without needing any special case for the first term.
Then we enter a for loop, setting the numerator to the square root of 2 + itself as described above, after which pi is multiplied by the evaluated term. Again we do not need a special case for the first term as pi is initialised to 1.
You probably noticed that the formula does not give us π but 2/π so after the loop we need to calculate pi itself and than call print_as_text. Uncomment francois_viete in main and then build and run.
Program Output
Francois Viete ============== Definitive 3.141592653589793238 Calculated 3.141592653589794511
Just 28 terms (I guess calculable by hand four hundred years ago with care and patience) give us 14dp of accuracy. Not bad!
John Wallis
English mathematician John Wallis came up with this in 1655, another infinite product (ie. each term is multiplied by the previous), this one giving us π/2 rather than 2/π.
It's easy to see what is going on here, the denominator and numerator are alternately incremented by 2. In the code, if the term is odd we'll increment the denominator, and if it's even we'll increment the numerator.
estimatingpi.c part 6
//-------------------------------------------------------- // FUNCTION john_wallis //-------------------------------------------------------- void john_wallis() { puts("John Wallis\n==========="); int iterations = 1000000; long double numerator = 2.0L; long double denominator = 1.0L; long double pi = 1.0L; for(int i = 1; i <= iterations; i++) { pi*= (numerator / denominator); ((i%2) == 1) ? (denominator+= 2.0L) : (numerator+= 2.0L); } pi*= 2.0L; print_as_text(pi); }
The number of iterations is set to 1,000,000 (yes, really) and the numerator and denominator are set to 2 and 1 respectively, as per the first term. The pi variable is again set to 1 so we can multiply the first term by it without a special case.
The for loop is simple, it just multiplies pi by the next term, incrementing the numerator or denominator depending on whether the current term is odd or even.
So far we only have half of pi so after the loop it's multiplied by 2 and printed. Uncomment john_wallis in main, compile and run.
Program Output
John Wallis =========== Definitive 3.141592653589793238 Calculated 3.141591082795430010
A million terms, waaaaaay beyond what Mr Wallis could have calculated, and we only get 5dp of accuracy. I described Francois Viete's effort as "not bad" but I have to describe John Wallis's as "not good".
I suspected this might be because the later terms choked off at the limit of the long double data type's accuracy. However, I printed every 100,000th term but it seems this is not the case.
John Machin
This is English astronomer John Machin's contribution from 1706, which he used to calculate 100 digits. Others used his formula to calculate many more digits by hand and it was the most used until the 20th century.
It looks simple - no infinite series, just one calculation. However, the complexity is hidden in the use of arctan although we have math.h whereas Professor Machin did not.
estimatingpi.c part 7
//-------------------------------------------------------- // FUNCTION john_machin //-------------------------------------------------------- void john_machin() { puts("John Machin\n==========="); long double pi = (4.0L * atanl(1.0L/5.0L) - atanl(1.0L/239.0L)) * 4.0L; print_as_text(pi); }
This function is very straightforward, just a translation of the formula into C. Note the use of atanl, the long double version of atan. Also, the formula only gives us a measly quarter of π so we need to multiply it by 4 at the end.
Uncomment john_machin in main and run the program.
Program Output
John Machin =========== Definitive 3.141592653589793238 Calculated 3.141592653589793239
This looks impressive but of course the accuracy is down to the arctangent we have available. I am not sure how this would have been calculated in Machin's time or how much accuracy he would have achieved.
Gregory-Leibniz
Next we have the Gregory-Leibniz series. This is very simple, it gives us π rather than a fraction or multiple of it and the numerator is constant. All we need to do is add 2 to the denominator each time and flip between subtraction and addition.
estimatingpi.c part 8
//-------------------------------------------------------- // FUNCTION gregory_leibniz //-------------------------------------------------------- void gregory_leibniz() { puts("Gregory-Leibniz\n==============="); int iterations = 400000; long double denominator = 1.0L; long double multiplier = 1.0L; long double pi = (4.0L / denominator); for(int i = 2; i <= iterations; i++) { pi += ( (4.0L / (denominator += 2.0L)) * (multiplier *= -1.0L) ); } print_as_text(pi); }
The denominator variable is first set to 1, and we also have a multiplier which will alternate between 1 and -1 to allow for the alternating addition and subtraction. Pi is initialised to the first term.
The for loop starts at the second term, taking the present value of pi and adding the next term. Within the evaluation of the term the denominator is incremented, and the value multiplied by 1 or -1 to allow for the alternating addition and subtraction. The multiplier is also "flipped" within the expression.
Uncomment the function in main and run again.
Program Output
Gregory-Leibniz =============== Definitive 3.141592653589793238 Calculated 3.141590153589793259
I have to admit the result is a bit baffling. The second batch of accurate digits cannot be chance as with 333/106, but why the wildly inaccurate two digits before them?
Nilakantha
Finally another infinite sum, this one from the 15th century Indian mathematician Nilakantha Somayaji, which has some obvious similarities with the Gregory-Leibniz series.
So we start off with 3, have a constant numerator of 4, and calculate the denominator by multiplying three numbers. We also alternate between addition and subtraction.
estimatingpi.c part 9
//-------------------------------------------------------- // FUNCTION nilakantha //-------------------------------------------------------- void nilakantha() { puts("Nilakantha\n========="); int iterations = 1000000; long double multiplier = -1.0L; long double start_denominator = 2.0L; long double pi = 3.0L; for(int i = 1; i <= iterations; i++) { pi += ( (4.0L / (start_denominator * (start_denominator + 1.0L) * (start_denominator + 2.0L)) ) * (multiplier *= -1.0L)); start_denominator += 2.0L; } print_as_text(pi); }
Again we have a multiplier to alternate between addition and subtraction. The start_denominator variable represents the first of the three numbers which are multiplied to get the denominator, and will be incremented by 2 on each iteration. We then start of pi at 3.
Within the loop the numerator 4 is divided by the denominator which is calculated by multiplying start_denominator by its two consective numbers, this then being multiplied by 1 or -1 alternately. Lastly within the loop 2 is added to start_denominator.
Uncomment nilakantha in main and compile/run.
Program Output
Nilakantha ========= Definitive 3.141592653589793238 Calculated 3.141592653589793217
This is pretty good, but even 50 iterations gets us 5dp of accuracy. Remember that each additional digit only gives us at most a 10% increase in accuracy.