Calculating PI, Gauss - Legendre method

[bobsmith1492]

Ok, so I thought Super PI was cool and wanted to make my own. So, I broke out the ol' C programming and fired up a simple program. However, it doesn't work.

Here's the algorithm:
http://cage.rug.ac.be/~hvernaev/Gauss-L.html

I was wondering if the limitations on the double number type could cause unintelligible results, or if maybe someone could spot a blatant error in my code?

- Thanks in advance

Here's my code by the way:

/* Program to calculate PI */

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

int main(void)
{
double a=1, b, t=0.25, x=1, y;
int its, i;

b = 1/sqrt(2);

printf("\n Please enter the number of iterations: >");
scanf("%u", &its);

for(i=0; i<its; i++)
{
y = a;
a = (a + b) / 2;
b = sqrt(b * y);
t = t - x*(y - a)*(y-a);
x = 2*x;
}

printf("\nPI = %u\n\n", ((a+b)*(a+b)/(4*t)));


return(0);
}

 

[TuxDave]

So.. uh... what part doesn't work? If you do 4 iterations, are you getting roughly 16 digits of accuracy?

 

[ghackmann]

printf("\nPI = %u\n\n", ((a+b)*(a+b)/(4*t)));

Last time I checked, pi wasn't an integer.

 

[bobsmith1492]

Ha; thanks ghackmann

Yeah, it works now.

Before it was giving me some random integer.

The next question is: this method gives only up to 6 decimals. Is that the limitation of the double number format? If so, what can be done to create a more accurate value? Chain more doubles together? But how can you split up a decimal number....

 

[ghackmann]

printf() only prints floating-point numbers to 6 decimal places by default. You can override this by using the formatting string %.[x]f instead of %f, where [x] is the number of decimal places you want printed.

Be warned that floating-point numbers are pretty inaccurate, so the more decimal places you go out, the less likely it'll be right. As far as I know, C and C++ don't have a built-in mechanism for representing decimals of arbitrary length and accuracy. (Somebody's inevitably built a third-party library to do it, but I don't know of any off the top of my head.) If you're willing to use a different language, Java has the BigDecimal class that should fit the bill.