Pi的数值计算
Pi的数值计算
提问于 2013-11-26 18:42:19
我想通过运行下面的代码来近似地计算Pi，它适合于具有单位直径的圆内n个正多边形，并使用代码中的函数计算其周长。然而，当使用长双变量类型时，第34项后的输出为0，或者使用双变量类型时，输出无界增加。我怎样才能补救这种情况？任何建议或帮助都会受到赞赏和欢迎。
谢谢
操作系统: 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
}
}回答 2
Stack Overflow用户
回答已采纳
发布于 2013-11-26 18:51:54
从接近- 1.0的数字中减去1.0将导致“灾难性的取消”，其中FP计算中的相对误差由于损失了重要的数字而急剧上升。试着对每个i在0到60之间评估pow(2, i) - (pow(2, i) - 1.0)，您就会明白我的意思。
唯一真正解决这个问题的办法是重组你的方程，以避免减去几乎相等的非零量。有关更多细节，请参见Acton，，精确性和数值算法的稳定性。
Stack Overflow用户
发布于 2013-11-26 20:30:28
一些想法：
使用y = (1.0 - x)*( 1.0 + x)而不是y = 1.0 - x*x。这有助于实现“几乎相等值的减法”的一个阶段，但随着1.0 - sqrtl(y)接近1.0，我仍然停留在下一个y上。
// 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);页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/20225326
复制相关文章
相似问题
腾讯云开发者
