我在C++精度问题中使用Chudnovsky计算圆周率的程序

我的代码：

#include <iostream>
#include <iomanip>
#include <cmath>
long double fac(long double num) {
long double result = 1.0;
for (long double i=2.0; i<num; i++)
result *= i;
return result;
}
int main() {
using namespace std;
long double pi=0.0;
for (long double k = 0.0; k < 10.0; k++) {
pi += (pow(-1.0,k) * fac(6.0 * k) * (13591409.0 + (545140134.0 * k))) 
/ (fac(3.0 * k) * pow(fac(k), 3.0) * pow(640320.0, 3.0 * k + 3.0/2.0));
}
pi *= 12.0;
cout << setprecision(100) << 1.0 / pi << endl;
return 0;
}

我的输出：

3.1415926535897637228433865175247774459421634674072265625

此输出的问题是它输出的是56位数字，而不是100位；我该怎么解决？

首先您的阶乘是错误的，循环应该是for (long double i=2.0; i<=num; i++)而不是i<num！！！

如注释中所述，double最多只能容纳~16位数字，因此您的100位数字不适用于此方法。要解决这个问题，有两种方法：

使用高精度数据类型

这方面有libs，或者你可以自己实现它——你只需要一些基本的操作。请注意，要表示100位数字，您至少需要
```
ceil(100 digits/log10(2)) = 333 bits
```
尾数或不动点整数，而double只有53
```
53*log10(2) = 15.954589770191003346328161420398 digits
```

使用不同的PI计算方法

对于任意精度，我建议使用BPP。但是，如果你只想要100个数字，你可以在字符串上使用这样的简单泰勒级数(不需要任何高精度数据类型或FPU(：

//The following 160 character C program, written by Dik T. Winter at CWI, computes pi  to 800 decimal digits. 
int a=10000,b=0,c=2800,d=0,e=0,f[2801],g=0;main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}

除了明显的精度限制从性能和精度两个方面来看，您的实现都非常糟糕，这就是为什么在非常低的k迭代中达到双倍精度限制时，您会更快地失去精度。如果你重写迭代，使子结果尽可能小(就尾数而言(，并且不使用太多不必要的计算，这里有几个提示：

为什么要一次又一次地计算相同的阶乘

循环中有k!，其中k是递增的，为什么不将k乘以某个包含实际阶乘的变量呢？例如：
```
//for (    k=0;k<10;k++){              ... fac(k) ... }
for (f=1,k=0;k<10;k++){ if (k) f*=k; ... f      ... }
```
为什么要一次又一次地用阶乘除法

如果你想一想，那么如果(a>b)，你可以计算这个：
```
a! / b! = (1*2*3*4*...*b*...*a) / (1*2*3*4*...*b)
a! / b! = (b+1)*(b+2)*...*(a)
```
我根本不会为此使用pow

CCD_ 11是"；非常复杂"；造成进一步精度和性能损失的函数，例如pow(-1.0,k)，可以这样做：
```
//for (     k=0;k<10;k++){       ... pow(-1.0,k) ... }
for (s=+1,k=0;k<10;k++){ s=-s; ... s           ... }
```
pow(640320.0, 3.0 * k + 3.0/2.0))也可以用阶乘的方法计算，pow(fac(k), 3.0)可以用3倍乘以变量fac(k)来代替。。。
热pow(640320.0, 3.0 * k + 3.0/2.0)甚至生长出(6k)!

所以你可以把它除以它，使子结果更小。。。

这些简单的调整将大大提高精度，因为你会溢出双倍精度，因为子结果将比原始结果小得多，因为阶乘往往增长非常快

将所有这些放在一起会导致以下情况：

double pi_Chudnovsky() // no pow,fac lower subresult
{                   // https://en.wikipedia.org/wiki/Chudnovsky_algorithm
double pi,s,f,f3,k,k3,k6,p,dp,q,r;
for (pi=0.0,s=1.0,f=f3=1,k=k3=k6=0.0,p=640320.0,dp=p*p*p,p*=sqrt(p),r=13591409.0;k<27.0;k++,s=-s)
{
if (k)  // f=k!,  f3=(3k)!, p=pow(640320.0,3k+1.5)*(3k)!/(6k)!, r=13591409.0+(545140134.0*k)
{
p*=dp; r+=545140134.0;
f*=k; k3++; f3*=k3; k6++; p/=k6; p*=k3;
k3++; f3*=k3; k6++; p/=k6; p*=k3;
k3++; f3*=k3; k6++; p/=k6; p*=k3;
k6++; p/=k6;
k6++; p/=k6;
k6++; p/=k6;
}
q=s*r; q/=f; q/=f; q/=f; q/=p; pi+=q;
}
return 1.0/(pi*12.0);
}

正如您所看到的，k会上升到27，而您的天真方法在溢出之前的64位双打中只能上升到18。然而，结果与双尾数在2次迭代后饱和相同。。。

由于以下代码，我感到很高兴：(

/*
I have compiled using cygwin
change "iostream...using namespace std" OR iostream.h based on your compiler at related OS.
*/
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
long double fac(long double num)
{
long double result = 1.0;
for (long double i=2.0; num > i; ++i)
{
result *= i;
}
return result;
}
int main()
{
long double pi=0.0;
for (long double k = 0.0; 10.0 > k; ++k)
{
pi += (pow(-1.0,k) * fac(6.0 * k) * (13591409.0 + (545140134.0 * k)))
/ (fac(3.0 * k) * pow(fac(k), 3.0) * pow(640320.0, 3.0 * k + 3.0/2.0));
}
pi *= 12.0;
cout << "BEFORE USING setprecision VALUE OF DEFAULT PRECISION " << cout.precision() << "n";
cout << setprecision(100) << 1.0 / pi << endl;
cout << "AFTER  USING setprecision VALUE OF CURRENT PRECISION   WITHOUT USING fixed " << cout.precision() << "n";
cout << fixed;
cout << "AFTER  USING setprecision VALUE OF CURRENT PRECISION   USING fixed " << cout.precision() << "n";
cout << "USING fixed PREVENT THE EARTH'S ROUNDING OFF INSIDE OUR UNIVERSE :)n";
cout << setprecision(100) << 1.0 / pi << endl;
return 0;
}
/*
$ # Sample output:
$ g++  73256565.cpp -o ./a.out;./a.out
$ ./a.out
BEFORE USING setprecision VALUE OF DEFAULT PRECISION 6
3.14159265358976372457810999350158454035408794879913330078125
AFTER   USING setprecision VALUE OF CURRENT PRECISION   WITHOUT USING fixed 100
AFTER   USING setprecision VALUE OF CURRENT PRECISION   USING fixed 100
USING fixed PREVENT THE EARTH'S ROUNDING OFF INSIDE OUR UNIVERSE :)
3.1415926535897637245781099935015845403540879487991333007812500000000000000000000000000000000000000000
*/

相关内容

最新更新

热门标签：