/* 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
Pages, some stolen, some original
▼
Sunday, March 15, 2009
π, 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.
No comments:
Post a Comment