pete wrote:

>

> pete wrote:

> >

> > Keith Thompson wrote:

> > >

> > > "Michel Rouzic" <(E-Mail Removed)> writes:

> > > > 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.1415926535897932384626433832795028841971 693993751

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.141592653589793238462643383279502884197169399375 1;

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.141592653589793238462643383279502884197169399375 1;

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 /

;

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 */

--

pete