Calculating PI, with computer or pencil
Over the centuries a number of people
have wasted their lives calculating ever more useless
digits of PI by hand.
Computers have made the job of calculating ever more useless
digits much easier  the current record is approximately
6.4 billion (you can download
a very fast Windows based program for calculating millions
of digits). Update: as of December 28, 2013, the record is 12.1 trillion digits of pi.
In the course of completing my high precision number class for
Fractal eXtreme
I had to create high precision versions of
functions like exp, log, sin, cos, atan, sqrt, etc. I
wanted to test these, and calculating PI to a few thousand
digits seemed like a good way, so I scanned around for an
appropriate formula for PI.
A little trigonometric thought will show that, since the
tangent of a forty five degree angle is 1 (since sin
and cos are equal at that point) then atan(1)
is equal to forty five degrees, or PI/4 radians.
No problem, a naive implementor might think. I'll just
calculate atan(1), perhaps using the handy Taylor series.
The Taylor series for atan(x) is:
atan(x) = x  x^3/3 + x^5/5  x^7/7 + x^9/9...*
therefore, simplifying for the case x equals 1:
PI/4 = atan(1) = 1  1/3 + 1/5  1/7 + 1/9...
Trouble is, the Taylor series for atan converges very
slowly for values of x near one, and it converges
extremely slowly when x is equal to one.
My guess would be that if you wanted PI accurate to
ten digits using this method, you would probably have
to sum up at least a few billion terms.
I would guess
that this sequence is not terribly stable, so you
would have to do all your math to quite high
precision, just to produce a very small number of
digits.
Not a good method.
A much better method is available however. It turns out
that:
PI/4 = atan(1) = 4 * atan(1/5)  atan(1/239)**
This is much more promising, because the Taylor series
for atan(x) actually converges at a reasonable rate
when x is equal to 1/5 or 1/239. Since each successive
term is reduced by a factor of x squared, that means that
with each term of atan(1/5) we gain about 1.4 digits.
atan(1/239) converges much faster,
with each term gaining about 4.75
digits. These swiftly diminishing terms have two
advantages. Number one is that we need to do a
relatively small number of terms to calculate PI to
thousands of digits. Number two is that the calculations
are quite stable. If a particular term is, say, 1e50,
then we know that all subsequent terms, and the sum of
all subsequent terms, will be less than 1e50. Therefore
we know that at that point we have calculated atan
to fifty digits of accuracy.
If we have a C++ high precision number class called InfPrec
then a generic method for calculating atan(x) might
look like this:
InfPrec atan(const InfPrec& x)
{
InfPrec Result = x;
InfPrec XSquared = x * x;
InfPrec Term = x;
int Divisor = 1;
while (Term != 0)
{
Divisor += 2;
Term *= XSquared;
Result = Term / Divisor;
Divisor += 2;
Term *= XSquared;
Result += Term / Divisor;
}
return Result;
}
Therefore, our high precision number class needs to support
multiplying two high precision numbers together, dividing
a high precision number by an integer, comparing a high
precision number to zero, and adding and
subtracting high precision numbers.
Unseen, but equally important, we need a way to initialize
x with the values 1/5 and 1/239. This is probably best done
by initializing the high precision numbers to one, and then
dividing by an integer. We also probably want a way to print
the results, preferably in decimal. This requires having
a way to multiply a high precision number by an integer, so
that we can repetitively multiply by ten and strip off the
integer portion for printing.
So, calculating atan, and by extension PI, isn't too difficult
at all.
However this method is both harder and slower than it needs to
be. The hardest and slowest part of our algorithm is muliplying
two high precision numbers together. This is mildly messy to
write, especially in a high level language, and (unless
you employ very sophisticated algorithms***)
it is an O(n^2) algorithm. O(n^2) is just a fancy way of saying
that doubling the number of digits will quadruple your calculation
time.
We are using a generic atan routine. This is a handy thing
to have, but overkill for this situation. We know that x will
always be either 1/5 or 1/239. Therefore, instead of multiplying
by a high precision number, we could simply divide by
an integer  in this case, dividing by either 5^2 or 239^2.
Normally I would not advocate replacing multiplication with division 
it goes against everything we are taught in optimization school.
But when the multiplication is going to involve millions of
operations, and the division is going to involve just a few
thousand, it can be a huge win.
Here's a special version of atan that handles the case that
we need much more efficiently.
InfPrec ataninvint(int x)
{
InfPrec Result = InfPrec(1) / x;
int XSquared = x * x;
int Divisor = 1;
InfPrec Term = Result;
InfPrec TempTerm;
while (Term != 0)
{
Divisor += 2;
Term /= XSquared;
Result = Term / Divisor;
Divisor += 2;
Term /= XSquared;
Result += Term / Divisor;
}
return Result;
}
So now our high precision number class just needs to support
multiplying and dividing a high precision number by an integer,
comparing a high precision number to zero, and adding and
subtracting high precision numbers. These are all relatively
easy and fast operations, even in a high level language.
With careful coding of these routines, a desktop computer can
now calculate PI accurate to millions of digits! The code
to calculate it is simply:
InfPrec PI = (ataninvint(5) * 4  ataninvint(239)) * 4;
Pretty darned simple. On a 233Mhz Pentium II I was able to
calculate one million digits of PI in just under six hours.
What I think is really cool about this is how easily this
algorithm works for manual calculation of PI. With a few
sheets of paper, some patience, and a vague memory of
long division, you can calculate PI to quite a high
degree of accuracy. I estimate that manually calculating
PI to thirty six digits would probably take a couple of
days. Less than
three hundred years ago****, that would have been enough for
a world record, beating the efforts of a number of people
who devoted much more time for worse results. They weren't
using such efficient methods, because they hadn't yet been
discovered.
When calculating PI by hand it makes sense to premultiply
the two atan calls. This allows us to read off digits
more easily, and makes the convergence more obvious.
So, our formula is:
PI = 4 * (4 * atan(1/5)  atan(1/239))
or:
PI = 16 * atan(1/5)  4 * atan(1/239)
Plugging this into our atan formula we get:
16 * atan(1/5) = 16/5  16/3*5^3 + 16/5*5^5  16/7*5^7 + 16/9*5^9...
and
4 * atan(1/239) = 4/239  4/3*239^3 + 16/5*239^5  16/7*239^7 + 16/9*239^9...
Calculating 16/5^n is particularly easy for people accustomed to
dealing with powers of two. Since 1/5 equals 2/10, (1/5)^n can be
rewritten as 2^n/10^n. In other words, a power of two, shifted right
n digits. Multiplying by the additional sixteen makes our final
calculation equivalent to 2^(n+4)/10^n.
Calculating the terms of atan(1/239) is a bit harder, because
dividing by 239^2 is a bit messy. But, you need comparatively much
fewer terms of atan(1/239) than of atan(1/5).
In the following tables, underlined digits represent digits
which repeat. Noticing these repeating digits simplify the calculations
tremendously because they frequently allow you to get an arbitrary
number of digits of accuracy for a fixed amount of work.
16*atan(1/5)
Sign  n  16/5^n (Term)  16/n*5^n (Term / n) 
+  1  3.2  3.2 
  3  0.128  0.0426 
+  5  0.00512  0.001024 
  7  0.0002048  0.0000292571428 
+  9  0.000008192  0.0000009102 
  11  0.00000032768  0.0000000297890 
+  13  0.0000000131072  0.000000001008246 
  15  0.000000000524288  0.000000000034952 
+  17  0.00000000002097152  0.000000000001233 
  19  0.0000000000008388608  0.000000000000044 
4*atan(1/239)
Sign  n  4/239^n (Term)  4/n*239^n (Term / n) 
+  1  0.0167364016736401673  0.0167364 
  3  0.0000002929991014450  0.000000097666367 
+  5  0.0000000000051294462  0.000000000001025 
Since addition is easier than subtraction, we're going
to try to do as little subtraction as possible. Rather
than subtracting the terms of atan(1/239), we simply
negate them all so that we can then add them. Our
answer is then the sum of all of the terms.
Half of those terms are still negative, so we won't be
able to avoid all subtractions, but if we sum up all of
the positive terms together, and all of the negative
terms together, then we just have to subtract the
sum of the negatives from the sum of the positives.
+3.200000000000000
+0.001024000000000
+0.000000910222222
+0.000000001008246
+0.000000000001233
+0.000000097666367

+3.201025008898068
0.042666666666666
0.000029257142857
0.000000029789090
0.000000000034952
0.000000000000044
0.016736401673640
0.000000000001025

0.059432355308274
And the answer is...
PI = 3.201025008898068  0.059432355308274 = 3.141592653589794
The last digit should actually be a three rather than a four,
but we have still managed to calculate PI to 14 decimal places!
You can
download (27K)
a sample Win32 program, with full, portable source code
and generate a few thousand digits of PI yourself.
Thanks for reading. I hope you found this interesting. While you're
here, take a look at Fractal eXtreme, which makes
use of a much more heavily optimized and extended version of this
InfPrec class. This allows fast and easy zooming in on many different
types of fractals to
virtually infinite magnification.
You can look
here
for more information on historical methods for calculating PI.
Drop me a line if you have any comments.
.Bruce Dawson.
Last modified: August 8, 1997
* Black magic #1.
You'll just have to take it on faith that this series does actually
sum up to atan(x). I do.
** Black magic #2.
This formula (atan(1) equals 4 * atan(1/5)  atan(1/239)) is called the Machin Formula.
You can read more about it, including how to easily derive it, on this page.
*** Although it is intuitively obvious that multiplication of high
precision numbers must take n squared operations, it turns out that
there are algorithms that improve the performance to take a number
of operations proportional to n. There is considerable complication
and overhead, but on numbers with millions of digits this can make
multiplication run much much faster than the obvious method.
Luckily we don't need this technique. For more information, read
Knuth's excellent Art of Computer Programming.
**** Van Ceulen calculated PI to 35 digits c. 1600, apparently
using Archimedes' method of inscribed and circumsribed polygons.
In 1699 Sharp used PI/6 equals atan(30) to calculate 71 digits.
Two years later Machin discovered and used the method described
here and calculated 100 digits.
Why Buy FX? 
Download Area 
How to Buy FX 
The Gallery
Fractal Theory 
Comments Area 
Company Profile 
Tips & Tricks
Main Page 
Links 
Send Mail
