π, again

Yesterday was pi day and Stu put up a post with a formula for calculating the value of pi, so naturally I had to try it out. This one promised to be a little better than the traditional one. I had written a program to implement that algorithm a while back, and while it worked, it was very slow to close in on the value of pi, requiring hundreds of iterations. This formula closes in on pi in only 8 iterations. Well, it gets as close as it can. It is limited by the format of double precision floating point numbers.

/* Simple program to calculate the value of pi.
  Original formula:

           /   4        2        1        1    \
  16^-K * (  ------ - ------ - ------ - ------  )
           \ 8K + 1   8K + 4   8K + 5   8K + 6 /

  The caret (^) indicates raising to a power: 16^K means 16 to the power of K.
  Multiplying by 16^-K is the same as dividing by 16^K

  A double precision floating point number has 64 bits (8 bytes).
   1 bit for the sign
  11 bits for the exponent
  52 bits for the mantissa

  Three decimal digits require approximately 10 bits. 2^9 = 1 << 10 ~= 1024
  We were able to get 13 decimal digits, which works out to about 43 bits.
  The double angle brackets (<<) indicate a bitwise shift,
  in this case to the left.
*/

#include <stdio.h>
#include <math.h>

int main(void)
{
   int       k;
   int       n;
   double    pi;
   double    x;

   pi = 0;
   for (k=0; k<8; k++)
   {
       n = 8 * k;
       x = ( 4.0/(n + 1)
           - 2.0/(n + 4)
           - 1.0/(n + 5)
           - 1.0/(n + 6))
           / (1 << (4 * k));    // 16^K
       pi += x;
       printf("%16.12f    ", x);
       printf("%16.12f\r\n", pi);
   }
   return 0;
}


Results:
      x                   pi
3.133333333333      3.133333333333
0.008089133089      3.141422466422
0.000164923924      3.141587390347
0.000005067221      3.141592457567
0.000000187893      3.141592645460
0.000000007768      3.141592653228
0.000000000345      3.141592653573
0.000000000016      3.141592653589