pete said:
Keith said:
Michael Mair wrote:
Michel Rouzic wrote:
[...]
What's wrong with the way I do my pi?
1) The C standard makes no guarantee whatsoever about the
precision of the math functions from the library. So, "your"
pi could be way off.
2) Your pi has to be recomputed at every function call.
oh, ok, well, I guess it's gonna be better
if I stop doing that and put
#define pi 3.1415926535897932 instead
That's 17 significant digits; long double often has more precision
than that.
Using 40 decimal digits will cover
anything with a mantissa up to 128
bits, which is enough for any hardware I've ever used. Use an 'L'
suffix to make sure you get the full precision.
It's very likely that you can get away with far less precision,
depending on the application, but using more digits than you'll ever
need isn't a burden, and it means one less thing to worry about if
your program misbehaves.
I got all the precision of a double right here:
/* BEGIN pi.c */
Subsequent testing indicates maybe I don't
got all the precision of a double.
I'm pretty close though.
My feeling is that since
(double)3.1415926535897932384626433832795028841971693993751
compares equal to
4 * atan(1)
on my system, then,
4 * atan(1)
is probably correct.
The return value of pi(), seems to be low by 2 * DBL_EPSILON.
I tried adding the positive half of an extra term
in each loop in pi3, but it didn't help.
/* BEGIN pi.c output */
PI = 3.1415926535897932384626433832795028841971693993751;
Pi = pi();
Pi2 = pi2();
Pi3 = pi3();
PI is 3.141593
Pi is 3.141593
Pi2 is 3.141593
Pi3 is 3.141593
PI - 4 * atan(1) is 0.000000e+000
Pi - 4 * atan(1) is -4.440892e-016
Pi2 - 4 * atan(1) is -4.440892e-016
Pi3 - 4 * atan(1) is -4.440892e-016
DBL_EPSILON * 2 is 4.440892e-016
Pi - Pi2 is 0.000000e+000
Pi - PI + DBL_EPSILON * 2 is 0.000000e+000
/* END pi.c output */
/* BEGIN pi.c */
#include <stdio.h>
#include <float.h>
#include <math.h>
double pi(void);
double pi2(void);
double pi3(void);
int main(void)
{
double Pi, Pi2, Pi3, PI;
PI = 3.1415926535897932384626433832795028841971693993751;
Pi = pi();
Pi2 = pi2();
Pi3 = pi3();
puts("/* BEGIN pi.c output */\n");
puts("PI = 3.1415926535897932384626"
"433832795028841971693993751;");
puts("Pi = pi();\nPi2 = pi2();\nPi3 = pi3();\n");
printf("PI is %f\n", PI);
printf("Pi is %f\n", Pi);
printf("Pi2 is %f\n", Pi2);
printf("Pi3 is %f\n", Pi3);
printf("PI - 4 * atan(1) is %e\n", PI - 4 * atan(1));
printf("Pi - 4 * atan(1) is %e\n", Pi - 4 * atan(1));
printf("Pi2 - 4 * atan(1) is %e\n", Pi2 - 4 * atan(1));
printf("Pi3 - 4 * atan(1) is %e\n", Pi3 - 4 * atan(1));
printf("DBL_EPSILON * 2 is %e\n", DBL_EPSILON * 2);
printf("Pi - Pi2 is %e\n", Pi - Pi2);
printf("Pi - PI + DBL_EPSILON * 2 is %e\n",
Pi - PI + DBL_EPSILON * 2);
puts("\n/* END pi.c output */");
return 0;
}
double pi(void)
{
double pi, b, numerator;
unsigned denominator;
pi = 0;
numerator = 1 / 5.0;
denominator = 1;
do {
b = numerator / denominator;
numerator /= 25;
denominator += 2;
b -= numerator / denominator;
numerator /= 25;
denominator += 2;
pi += b;
} while (b > DBL_EPSILON / 8);
pi *= 4;
numerator = 1 / 239.0;
denominator = 1;
do {
b = numerator / denominator;
numerator /= 57121;
denominator += 2;
b -= numerator / denominator;
numerator /= 57121;
denominator += 2;
pi -= b;
} while (b > DBL_EPSILON / 2);
return 4 * pi;
}
double pi2(void)
{
double pi, b, numerator;
unsigned denominator;
pi = 0;
numerator = 1 / 2.0;
denominator = 1;
do {
b = numerator / denominator;
numerator /= 4;
denominator += 2;
b -= numerator / denominator;
numerator /= 4;
denominator += 2;
pi += b;
} while (b > DBL_EPSILON / 4);
numerator = 1 / 3.0;
denominator = 1;
do {
b = numerator / denominator;
numerator /= 9;
denominator += 2;
b -= numerator / denominator;
numerator /= 9;
denominator += 2;
pi += b;
} while (b > DBL_EPSILON / 2);
return 4 * pi;
}
double pi3(void)
{
double pi, b, numerator;
unsigned denominator;
pi = 0;
numerator = 1 / 2.0;
denominator = 1;
do {
b = numerator / denominator;
numerator /= 4;
denominator += 2;
b -= numerator / denominator;
numerator /= 4;
denominator += 2;
pi += b;
} while (b > DBL_EPSILON / 4);
b = numerator / denominator;
pi += b;
numerator = 1 / 3.0;
denominator = 1;
do {
b = numerator / denominator;
numerator /= 9;
denominator += 2;
b -= numerator / denominator;
numerator /= 9;
denominator += 2;
pi += b;
} while (b > DBL_EPSILON / 2);
b = numerator / denominator;
pi += b;
return 4 * pi;
}
/* END pi.c */