Pi的数值评估

Question

我想通过运行下面的代码来估算Pi，该代码适合具有单位直径的圆的n边正多边形，并使用代码中的函数计算其周长。但是，当使用长双变量类型时，第34项后的输出为0，或者当使用双变量类型时，其增加无界限。我该如何纠正这种情况？任何建议或帮助表示赞赏和欢迎。

由于

PS：操作系统：Ubuntu的12.04 LTS 32位，编译器：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

减去从一个几乎1.0 1.0数是导致 “灾难性消除”， FP计算中的相对误差由于丢失有效数字而突然增大。尝试评估pow(2, i) - (pow(2, i) - 1.0)为0到60之间的每个我，你会明白我的意思。

这个问题的唯一真正的解决方案是重新组织你的方程，以避免减去几乎相等的非零量。有关更多详细信息，请参阅Acton，Real Computing Made Real或Higham，数字算法的准确性和稳定性。

Answer 2

一些想法：

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