BLOKPI: Yet Another Program to Calculate Pi


Introduction
------------

This documentation explains how blokpi.c works and what you can do with 
it. Like the source code, this file is placed in the public domain by
its author, Jason Papadopoulos. Give it to anyone you want, but please
try and keep it bundled together with blokpi.c so that future generations
will have a chance to decipher what I did.

What's so special about this particular program? A couple of things.

First, it's extremely fast; on modern processors it averages three times
faster than all the other public domain calculators that add up series.
On a 200MHz Pentium it can add up a million digits of pi in under 2.5
hours!

Second, it's versatile: I've written it to compile portably on an ANSI
C compiler, and also included options to tailor performance to a specific
machine and FPU. For my favorite combo (gcc on a Pentium) there's some
inline assembly code for better performance. I also made it generic enough
to sum any series that has a special form, so you can use it to compute
other constants to monster decimal places, at breakneck speed.

Third, it's new. The formula used is an unusual one I haven't seen used
elsewhere, and the digit hunting tricks that make it so fast don't appear
anywhere else. I've seen lots of pi-crunching code over the years, and
it's always a treat to stumble onto something different.

I've verified that it produces correct results up to 20,000 digits of pi.


Configuring the Thing
---------------------

Basically, the only part of the code you have to worry about is at the
very beginning. You should tweak the #defines for your particular compiler
and system. Brief explanations follow:

BLOCKPI: This should be set to (at most) the number of doubles that
         will fit into half your cache. For a Pentium it should be
         large but less than 512; with MMX it can go up to 1024, but
         the exact number isn't critically important. Bigger=better.

BLOCKTERM: This should be set to (at most) the number of doubles that
         will fit into 1/8 of your cache. For a Pentium it should be
         large but less than 128; with MMX it can go up to 256.
         Bigger=better, but unlike BLOCKPI the exact number you choose
         is critical for the program to work correctly. See the code
         for details.

INTEL:  Intel FPU registers have 64 bits, and rounding a floating point
        value to the nearest whole number requires a little finesse. For
        any other machine, don't define this unless you know for sure
        that its FPU registers are 64 bits long.

BASE:  The base you want to work in. Good choices are 1.0e10 for Intel
       stuff, and 1.0e7 for everybody else. If you know for a fact that
       your FPU registers have more than 53 bits of precision you can
       set this higher, but unless the mantissa is 64 bits long you have
       to fiddle with the value of ROUND yourself (it's in the code).
       The higher you make BASE, the faster the program will run but the
       fewer terms you can add up. Details are in the code; you can
       always make it lower, and every power of ten will multiply by ten
       the number of terms you can sum up.

BASEDIGITS: the power of 10 in BASE

GCCDOS: compile assembly code optimized for DJGPP for DOS. The program
        will run about 20% faster if you do this.

Everything else allows you to change formulas at compile time; an expla-
nation appears a little later.


How it Works
------------

The idea is to take infinite series of a certain form and rearrange them
into another form, suitable for compuation. As an example, take the
usual series for "e":

           1    1      1        1 
   e = 1 + - + --- + ----- + ------- + ...
           2   2*3   2*3*4   2*3*4*5

The traditional method of adding up such a series is to have two arrays,
one to hold the sum (starts with 1) and the other to hold one term
(starts with 1/2). Divide term by 3, add to sum; divide term by 4, add to
sum; divide term by 5, add to sum; and so on.

There are advantages to this method (hopefully obvious): first, as the
computation progresses the array "term" becomes smaller and smaller
because you're dividing by integers over and over again. Thus the program
will speed up as time goes by. Second, you don't have to figure out in
advance how many terms to add up, just keep going until the term array is
all zeros (this is a real plus in the case of this series because its
convergenece is superlinear, i.e. faster than "x digits per term").

It also has disadvantages, though: every time you add on a term you have
to be careful that the elements of "sum" don't overflow, so either the
elements in the sum have to have lots of unnecessary extra digits, or
you "carry the one" with every term and pay a substantial loss in per-
formance.

There's another way to compute this series, though. It first appeared
officially in a paper by Rabinowitz and Wagon in the American Mathematical
Monthly of about 1995 but the idea is much older; the first mention of it
in the literature was from a little paper in an ACM journal in the 60s.

        1       1       1       1
e = 1 + - ( 1 + - ( 1 + - ( 1 + - ( 1 + ...   )))...)
        2       3       4       5

Adding this thing up is a completely different matter from the first
method. You have one array, initialized to 1. Divide by huge number,
add 1; divide by huge number minus 1, add one; divide by huge number
minus 2, add 1; and so on. After dividing by 2 and adding one the
result is e.

First the advantages compared to the standard way: the memory needed is
cut in half; you only need one array, for the answer. Second, instead of
a multiple-precision add (i.e. adding one array to another) you only have
a single-precision add between divisions (i.e. add 1). This cuts lots
of complexity from the code for a significant time savings.

The disadvantages are hopefully obvious too: the worst is that the array
being worked on *never gets smaller* so the computation will never speed
up. This amounts to a large increase in computation time.

Note that the nested method above involves less arithmetic per term than
before but requires full precision in every term. We'll see in a minute
for a simpler series whether time is wasted or saved from these two issues.

There are other, comparatively minor disadvantages to worry about too:
first, you have to start at the end and work backwards, meaning you have
to know in advance how many terms to add up. For the e series this is a
problem but for the series we'll look at the convergence rate is very
regular. Also, you get no correct digits of the answer at all, until the
very end when you get them all at once. This isn't a problem either since
you'd have to be nuts to print out partial results while the computation
is in progress.

For the curious, the series I use is

     8       1    2        3    4        5    6
pi = -- (6 - -- * -- (13 - -- * -- (20 - -- * -- (26 - ...
     15      21   11       39   17       51   23

first discovered by RW Gosper. See RW Gosper, "Acceleration of Series",
Memo #304, MIT AI Lab, 1973. It's way more complicated than the
series for e, but it converges at a nice regular 1.43 digits per term,
about as fast as Machin's arctangent formula. Incidentally, Gosper's
paper is wonderful, you'll want to spend hours with Mathematica after
reading it (that is, if you're very sick like me).


Changing Formulas
-----------------

Do you have another series you can rewrite in nested form? blokpi.c
will add it up too, if you change a few things at compile time.

Suppose you wanted to use another pi series I discovered:

     2       1  1       2   3        3    5
pi = - (5 - - * - (11 - - * -- (17 - -- * -- (23 - ...
     3       5  7       9   11       13   17

Here each term is 1/8 the size of the previous one, so it's slower
than the default series the prorgram comes with. You only have to
change a few things:

FRACTIONS_PER_TERM: How many fractions are in front of each parenthesis?

CORRECT_DIGITS_PER_TERM: This will be log10(8); for the default series
                         the program comes with, this is log10(27).

Finally, you have to describe the series in the stuff[] array. To do this,
you have to look at the first term in the series that has all its fractions.
Here it's the one with the 1/5 and the 1/7. Then the stuff[] array is

   { ignore,
     ignore,
     the number after the ( in the first term (11 here),
     how much that number increases with every term (6 here),
     numerator of the fraction in front (2),
     denominator of the fraction in front (3),

     numerator of first fraction (1),
     how much it increases with every term (1),
     denominator of first fraction (-5),
     how much it increases with every term (-4),

     numerator of second fraction (1),
     how much it increases with every term (2),
     denominator of second fraction (7),
     how much it increases with every term (4) };

If there are more fractions, just enter them in the same fashion. Note
that I made one of the fractions negative, because there are minus signs
in the series. Everything would be positive if there were all plus signs.
Also, if you had to add more fractions make sure stuff[] is declared big
enough to hold them.

If you think pi is passe' you can do a lot of good with this program. Here
are some series you can crunch away on:


           5      1   1   1      2    2   2       3    3   3
 Zeta(3) = - (1 - - * - * - (1 - -- * - * - ( 1 - -- * - * - (1 - ...     
           4      6   2   2      10   3   3       14   4   4

                   1   3       3   7       5    11
Lemniscate A = 1 + - * -- (4 + - * -- (7 + -- * -- (10 + ...
                   5   10      9   18      13   26

The first has been known for a while, the second I found. In both cases
the world record is only a few million digits, and a fast Pentium can
get that many in a few days with this program. See the CECM (the
Center for Experimental and Constructive Mathematics at Simon Fraser U)
for up-to-date reports on the present records.


Some Analysis
-------------

Here's a simpler series as another example, found by using Euler's
transformation used on Gregory's series ( pi = 4 - 4/3 + 4/5 - 4/7 + ...):

               1   1 2   1 2 3   1 2 3 4 
pi  =  2 ( 1 + - + -*- + -*-*- + -*-*-*- + ... )                   (1)
               3   3 5   3 5 7   3 5 7 9

           1       2       3       4
    =  2 + - ( 2 + - ( 2 + - ( 2 + - ( 2 + ...   )))...)           (2)
           3       5       7       9

Every term in (1) is about half the size of the previous term, so (1)
converges at at about log10(.5) = .301 correct digits per term. Like
the series for e, (2) saves one full-size addition.

You don't have to do full size addition, and pretend for the moment that
you don't have to propagate carries when multiplying by a fraction. Is
equation (2) going to be faster or slower than equation (1)? A back of
the envelope calculation shows the following:

time for (1)     1        extra work in (1)
------------ ~=  - ( 1 +  ----------------- )
time for (2)     2           work in (2)

For the Euler series the work in (2) involes multiplying by a single
fraction, and the extra work for (1) is a full size add and carry adjust.
Suppose the extra work comes out to 1/5 of the work in (2) (not unrea-
sonable, since division takes a long time); then equation (1) will finish
in only 3/5 the time that (2) needs! If all the work you do for one term
can carry over to the next term the standard method of summing series
runs almost twice as fast as the nested method!

Do you have to carry any 1's in formulas like (2)? Suppose you had to
multiply the array by the fraction num/den; in C the code would look
like this:

   remainder = 0;
   for(i=0; i<size; i++) {
      dividend = remainder * BASE + num * array[i];
      array[i] = dividend / den;
      remainder = dividend % den;
   }

where BASE is some big number like 10000 or 65536. In effect this
code multiplies each array element by num and divides it by den; the
question is whether array[i] and remainder will always be < BASE. If
that's the case no carries need ever be propagated.

First, you can assume that the old remainder is always < BASE, because
the denominator is < BASE. Suppose remainder has r bits, BASE has
B bits, and denom has d bits. Then in the worst case dividend will have
r + B + 1 bits (assuming num < den, which almost always is the case).
The new remainder is safe (must be < denom, denom < BASE) but the new
quotient has r + B + 1 - d bits. The worst case occurs when remainder
and denom are about the same size; then the quotient will have B + 1
bits at most. The next time around, this B+1 bit quantity is going to
be multiplied by the next num, and if num is too big the quotient will
have B+2 bits, and so on until overflow occurs.

All this messy math amounts to the following: if the numerators you use
are less than half the size of the base you're working in (i.e. have
B-1 bits or less), you will never get overflow and thus never have to
propagate carries. You *do* have to deal with numbers one bit larger than
you expected, however, and so before you print anything out you should
perform all the carries and borrows.


So...back to pi formulas. If the nested method runs so much slower, what
good is it? Firstly, there are other pi formulas that can be nested in
the same way, and they involve work that doesn't carry over from one
term to the next. In that case the ratio becomes more favorable for nested
methods. Second, the analysis above only counts arithmetic operations; on
modern processors memory performance is just as important, if not more so,
and here nested summation has a distinct advantage.

The biggest advantage nested methods have, however, is that they're pre-
dictable. You know exactly how many fractions to apply, and how big the
array is. This lets you rearrange the computations in a way that modern
processors can execute quickly.


Rearranging Things
------------------

Believe it or not, there's actually quite a lot you can do to speed up
the inner loop as given above. Here it is again:

   remainder = 0;
   for(i=0; i<size; i++) {
      dividend = remainder * BASE + num * array[i];
      array[i] = dividend / den;
      remainder = dividend % den;
   }

This will multiply an array by a fraction num/den, and needs about 60
cycles per iteration on a Pentium. Rewriting the same thing in assembly
language won't change much; the loop will have to be changed completely.

The most expensive element by far is the division; on a Pentium an
integer division takes 4.5 times as long as an integer multiply, and a
multiply takes 10-20 times as long as integer addition. Other machines
aren't all that much better, because division is fundamentally slow.
Actually, the above has TWO divides, one for the / and the second for
the %, but a good optimizing compiler will recognize that a single division
can provide both answers.

There's no way to speed up a division, other than to get rid of it.
One can instead multiply by a reciprocal, but that requires floating
point math, and there's no remainder operation built in. The code is now:

   remainder = 0;
   recip = 1.0 / den;

   for(i=0; i<size; i++) {
      dividend = (double)remainder * BASE + (double)num * array[i];
      quotient = (long)(dividend * recip + .5);
      array[i] = quotient;
      remainder = (long)(dividend - (double)den * quotient);
   }

This first version trusts the C compiler to perform all the conversions
from integer to floating point to integer. On a Pentium, floating point
multiplies are three times faster than integer multiplies, but the speed
at which this code runs depends very heavily on how the instructions
are ordered. Also, loading an integer into the FPU takes a while, and
storing floating point numbers into an integer takes even longer.
The one major advantage of this method is that dividend can have 15
digits in it, so that BASE can be much larger than the integer version
and you can do more work at a time. With BASE = 10^9 the floating point
version above is over 2.1 times faster than the integer version.

You can get rid of most of the conversion problems by *always* using
floating point, and never converting to integers. THIS IS NOT DANGEROUS!
FPU registers act just like integer registers, as long as you never
overflow them. This is extra nice on Intel machines because the FPU
registers are twice as big as the integer ones, allowing BASE to be
even larger than before. The only problem is that you still have to round
the quotient to an integer, i.e. get rid of the fractional part somehow.
There's an old trick to doing this, much faster than some floor() function:
add a magic number, and then subtract it off. It's that simple, but
explaining why it works requires technical details about floating point
math that I don't want to get into. Anyway, here's the new code:

   recip = 1.0 / den;
   f0 = 0.0;

   for(i=0; i<size; i++) {
      f0 *= BASE;
      f1 = array[i];
      f1 += num;
      f0 += f1;                /* remainder * BASE + num * array[i]; */
      f1 = f0;
      f1 *= recip;
      f1 += ROUND;
      f1 -= ROUND;             /* compute quotient */
      array[i] = f1;
      f1 *= den;
      f0 -= f1;                /* compute remainder */
   }

f0 and f1 are register doubles, to make sure you get the extended
precision of FPU registers. I wrote everything out explicitly, to
minimize unnecessary loads and stores. Each iteration of the loop
requires about 28 cycles on a Pentium; this is 1.7 times
faster than the first floating point version, for a total factor of
3.6 better than the original integer code!


Even this can be improved on! Floating point on just about any
modern processor is "pipelined", i.e. one operation may take a while,
but the FPU can do *something else* while waiting for an operation
to finish.

The best analogy is what you do when you have 10 loads of laundry to
do but only one washer and dryer. You *could* start the first load
of laundry, wait till it finishes in the washer and dryer, then start
the second load of laundry, etc. However, it's a much better idea to
start the washer on the first load, and when you stick it in the dryer
immediately stick the second load in the washer. Once you get the
pipeline going you finish one load of laundry every cycle, in effect
getting twice the performance out of the equipment you have.

Floating point pipelining works the same way. If there are several
things you can do (i.e. loads, stores, adds, multiplies) that have
nothing to do with each other, then a modern FPU can do two or even
three things at a time. The FPU on a Pentium is much more finicky than
this, and there are a host of rules and exceptions, but the bottom line
is that multiplies and adds may take three cycles to finish, but two
of those three cycles can be overlapped with other stuff, making arithmetic
seem to work three times faster.

Unfortunately, for our fraction code, there are dependencies everywhere.
Each operation depends on the result of the last operation, so almost
nothing can overlap to save time. The only way around that is to multiply
by MORE THAN ONE FRACTION AT A TIME.

You have to be careful, though. It won't change anything to cut and
paste the code above three or four times, the whole thing has to be
reformulated yet again. Pretend this is your array:

        *   *   *   *   *   *

Multiplying it by one fraction will leave a trail of modified array
elements, plus the one that's being worked on (the &).

        X   X   X   &   *   *

Any of the three X's can be worked on at the same time, but you can't
touch the & because the computation is still in progress. Thus you
can multiply by another fraction at the same time, as long as that fraction
(the &2) stays at the heels of the first one.

       XX   XX   &2   &   *   *

and you can repeat the process with more fractions on the XX's. The code
as written uses three of them; fraction 1 multiplies x[i] while fraction
2 multiplies x[i-1] while fraction 3 multiplies x[i-2], and since each
multiplication has nothing to do with the other two, they can pipeline
perfectly. The code also includes an assembly-coded version that multi-
plies by four fractions at a time.

The savings are dramatic. The un-pipelined loop you've seen requires 28
cycles per iteration, but the assembly-coded version needs 50 cycles for
four of them, for a net time of 12.5 cycles per multiply per fraction!
This is almost 2.3 times faster than the basic FPU code. With all the
speedups and changes the code to multiply by fractions runs more than
EIGHT TIMES FASTER than what it started out as. Assembly code isn't vital
for that much speed, either: the assembly version is only 20% faster than
a version carefully written in C.

There's one final hurdle: you can only use pipelined multiplication
if the array is pipelined to begin with. To get it that way you have to
use basic multiplication; multiply two elements by the first fraction,
then one element by the second fraction. From that point on the array
is pipelined, and you can multiply by three fractions at a time. If
the array involved is large, the overhead of using slow multiplies the
first few times is negligible. In fact, the only real problem is that
the code to do this is painful to debug.


Block Decomposition
--------------------

Okay, so the basic unit of work has been sped up tremendously from before.
Can more time be saved by rethinking the computation somehow?

Here's one way: think what happens when you have a really long array, and
you want to multiply it by a bunch of fractions (pipelined or not). The
first thing I'd think of is to just burn through the whole array from left
to right, once for each fraction or group of fractions.

Once upon a time it didn't matter, but now it does! Modern processors
may be capable of doing arithmetic very fast, but memory performance
hasn't made the speed gains that arithmetic has made. If you only work
with small chunks of data, and use them over and over again before
switching to the next chunk of data, memory performance can improve
dramatically and programs that make heavy use of memory can be greatly
sped up.

Unfortunately, really big arrays require really big chunks of memory,
so racing through a really big array over and over again will not
re-use any data (by the time you get to the end of the array your
computer will have thrown away the chunks associated with the beginning
of the array, all in an effort to keep up with you). The solution?
Break up your big array into convenient-size chunks and do as much work
as possible on each chunk before moving on to the next one. Algorithms
that do this are called "block algorithms", and are very common in
high-performance environments.

Here's an example to put things in focus. Suppose your array is 12
elements long and you want to multiply it 1/3, 2/5, and 3/7. Suppose
also that your computer can only efficiently handle arrays 4 elements
long; the first time an array element is accessed there's a long delay
(a "cache miss") but later on if you want that same element again,
and your computer has accessed it recently, you can get it very quickly
(a "cache hit").

The basic way to multiply by three fractions is to take your array,
multiply every element by 1/3, then every element by 2/5, then every
element by 3/7. For the multiply by 1/3,

array:    *    *    *    *    *    *    *    *    *    *    *    *

x 1/3    miss miss miss miss miss miss miss miss miss miss miss miss

Your computer hasn't seen any of the array elements, so you get a cache
miss every time, and this first multiply will seem to take forever. Now
when you multiply by 2/5 your computer knows about the entire array,
but the only part that it can remember and access quickly is the last
four elements, since those are the last ones you used. Hence when you loop
around to the beginning of the array your computer has to rediscover it
again. In short, when your array is too big and you use this basic method,
you'll get a cache miss every last time.

Making a block algorithm out of this process is easy. You process the
array in chunks that your computer can remember easily.

array:    *    *    *    *    *    *    *    *    *    *    *    *

x 1/3    miss miss miss miss
x 2/5    hit  hit  hit  hit
x 3/7    hit  hit  hit  hit
               
                   x 1/3    miss miss miss miss
                   x 2/5    hit  hit  hit  hit
                   x 3/7    hit  hit  hit  hit

                                        x 1/3    miss miss miss miss
                                        x 2/5    hit  hit  hit  hit
                                        x 3/7    hit  hit  hit  hit

Once again, if an array element hasn't been seen before you suffer a
cache miss. However, once you've processed a few elements in a row,
instead of forgetting about them and moving on, you go back and multiply
those same elements by a different fraction. This scheme cuts the number
of cache misses dramatically, especially if you can afford to multiply
by many fractions before moving on. I'll call this "block summation",
for lack of anything better.

However, there are several gotchas to doing this. In order to multiply
x[i] by 1/3, you need the value of x[i] and also the remainder from multi-
plying x[i-1] by 1/3. Hence, for each of 1/3, 2/5, and 3/7 you have to
store some information, so that when you switch blocks you can pick up
the computation where you left off. This was never a problem before, but
if you want efficient memory use you have to store things. The best way
to do this is to have an array of terms, each element of which can store
all the info for a specific fraction.

But this is another array! It follows the same rules about accessing
memory that we had to deal with before. You're now dealing with multiplying
the array x[] by a block of fractions, and doing the multiplication on a
block of x[] at a time. The code for this is fiendishly complex, because
you have to deal with pieces of blocks. Not only that, but we're still
using pipelined methods to multiply by several fractions at once! This
is the reason I'm rewriting infinite series in nested form: the size of
the array worked on is always the same, so you only need to chop the
problem into blocks once.

How much time does it save? Less than I thought. Breaking up the problem
into blocks and using block summation speeds the whole program up by
about 10%, but without an external cache the speedup becomes 20%. Block
summation is the reason the code as is looks so messy, and if it only
speeds things up a little then it may be a good idea to abandon it (and
the slowdown inherent in nested series) entirely.

The End (for now)
-----------------