c - Pi 的数值评估

Question

我想通过运行以下代码来近似评估 Pi，该代码将 n 边的正多边形拟合在具有单位直径的圆内，并使用代码中的函数计算其周长。但是，当使用 long double 变量类型时，第 34 项之后的输出为 0，或者当使用 double 变量类型时，它会无限制地增加。我该如何补救这种情况？感谢并欢迎任何建议或帮助。

谢谢

P.S: 操作系统:Ubuntu 12.04 LTS 32-bit，编译器: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
    }
}

Answer 1

一些想法:

使用 y = (1.0 - x)*( 1.0 + x) 而不是 y = 1.0 - x*x。这有助于“减去几乎相等的值”的第一阶段，但我仍然停留在下一个 1.0 - sqrtl(y) 因为 y 接近 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) ));

更改数组大小或for()

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

---

下面是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);