CODEWITHSUNDEEP
c
Vesnog
Vesnog

Reputation: 785

Numerical Evaluation of Pi

I would like to evaluate Pi approximately by running the following code which fits a regular polygon of n sides inside a circle with unit diameter and calculates its perimeter using the function in the code. However the output after the 34th term is 0 when long double variable type is used or it increases without bounds when double variable type is used. How can I remedy this situation? Any suggestion or help is appreciated and welcome.

Thanks

P.S: Operating system: Ubuntu 12.04 LTS 32-bit, Compiler: GCC 4.6.3

#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#define increment 0.25
int main()
{
    int i = 0, k = 0, n[6] = {3, 6, 12, 24, 48, 96};
    double per[61] = {0}, per2[6] = {0};

    // Since the above algorithm is recursive we need to specify the perimeter for n = 3;
    per[3] = 0.5 * 3 * sqrtl(3);
    for(i = 3; i <= 60; i++)
    {
        per[i + 1] = powl(2, i) * sqrtl(2 * (1.0 - sqrtl(1.0 - (per[i] / powl(2, i)) * (per[i] / powl(2, i)))));
        printf("%d      %f \n", i, per[i]);
    }
    return 0;
    for(k = 0; k < 6; k++)
    {
        //p[k] = k
    }
}

Upvotes: 1

Views: 160

Answers (2)

chux
chux

Reputation: 153457

Some ideas:

Use y = (1.0 - x)*( 1.0 + x) instead of y = 1.0 - x*x. This helps with 1 stage of "subtraction of nearly equal values", but I am still stuck on the next 1.0 - sqrtl(y) as y approaches 1.0.

// per[i + 1] = powl(2, i) * sqrtl(2 * (1.0 - sqrtl(1.0 - (per[i] / powl(2, i)) * (per[i] / powl(2, i)))));
long double p = powl(2, i);
// per[i + 1] = p * sqrtl(2 * (1.0 - sqrtl(1.0 - (per[i] / p) * (per[i] / p))));
long double x = per[i] / p;
// per[i + 1] = p * sqrtl(2 * (1.0 - sqrtl(1.0 - x * x)));
// per[i + 1] = p * sqrtl(2 * (1.0 - sqrtl((1.0 - x)*(1.0 + x)) ));
long double y = (1.0 - x)*( 1.0 + x);
per[i + 1] = p * sqrtl(2 * (1.0 - sqrtl(y) ));

Change array size or for()

double per[61+1] = { 0 };  // Add 1 here
...
for (i = 3; i <= 60; i++) {
  ...
  per[i + 1] =

Following is a similar method for pi

unsigned n = 6;
double sine = 0.5;
double cosine = sqrt(0.75);
double pi = n*sine;
static const double mpi = 3.1415926535897932384626433832795;
do {
  sine = sqrt((1 - cosine)/2);
  cosine = sqrt((1 + cosine)/2);
  n *= 2;
  pi = n*sine;
  printf("%6u s:%.17e c:%.17e  pi:%.17e %%:%.6e\n", n, sine, cosine, pi, (pi-mpi)/mpi);
} while (n <500000);

Upvotes: 2

Sneftel
Sneftel

Reputation: 41474

Subtracting 1.0 from a nearly-1.0 number is leading to "catastrophic cancellation", where the relative error in a FP calculation skyrockets due to the loss of significant digits. Try evaluating pow(2, i) - (pow(2, i) - 1.0) for each i between 0 and 60 and you'll see what I mean.

The only real solution to this issue is reorganizing your equations to avoid subtracting nearly-equal nonzero quantities. For more details, see Acton, Real Computing Made Real, or Higham, Accuracy and Stability of Numerical Algorithms.

Upvotes: 1

Related Questions