7

编辑：要求是模糊的，而不是计算 pi 的第 n 位数字，他们只是希望 pi 到第 n 位数字不超过浮点数限制，因此蛮力方式适用于要求。

我需要计算第 n 位的 PI，我想尝试使用BBP 公式，但遇到了困难。我输入的方程式似乎没有正确地给我 PI。

(1 / pow(16,n))((4 / (8 * n + 1)) - (2 / (8 * n + 4)) - (1 / (8 * n + 5)) - (1 / (8 * n + 6)))

我用蛮力找到 PI 的方法很成功，但这只是如此准确，而且很难找到第 n 个数字。

(4 - (4/3) + (4/5) - (4/7)...)

我想知道是否有人对如何做到这一点有更好的想法，或者可能帮助我的 BBP 方程解决我搞砸的问题？

谢谢你，
LF4

功能性但不是很准确，直到几次迭代，然后你不得不忽略最后几个。

#include <iostream>

using namespace std;

int main()
{
    int loop_num = 0;
    cout << "How many digits of pi do you want?: ";
    cin  >> loop_num;

    double my_pi = 4.0;
    bool add_check = false;
    int den = 3;
    for (int i = 0; i < loop_num; i++)
    {
        if (add_check)
        {
            my_pi += (4.0/den);
            add_check = false;
            den += 2;
        }
        else
        {
            my_pi -= (4.0/den);
            add_check = true;
            den += 2;
        }
    }
    cout << "Calculated PI is: " << my_pi << endl;
    system("pause");

    return 0;
}

我希望会是一个更好的程序。

#include <iostream>
#include <cmath>

using namespace std;

const double PI_BASE = 16.0;

int main()
{
    int loop_num = 0;
    cout << "How many digits of pi do you want?: ";
    cin  >> loop_num;

    double my_pi = 0.0;
    for (int i = 0; i <= loop_num; i++)
    {
        my_pi += ( 1.0 / pow(PI_BASE,i) )( (4.0 / (8.0 * i + 1.0)) -
                                           (2.0 / (8.0 * i + 4.0)) -
                                           (1.0 / (8.0 * i + 5.0)) -
                                           (1.0 / (8.0 * i + 6.0)) );
    }
    cout << "Calculated PI is: " << my_pi << endl;
    system("pause");

    return 0;
}
4

3 回答 3

4

无论您使用什么公式，您都需要任意精度的算术才能获得 16 位以上的数字。（因为“double”只有 16 位精度）。

Chudnovsky 公式是已知最快的计算 Pi 的公式，每项收敛到 14 位。但是，要有效地实施是极其困难的。

由于这个公式的复杂性，用它来计算 Pi 到几千位数是没有意义的。所以除非你准备好用任意精度算术全力以赴，否则不要使用它。

使用 GMP 库的 Chudnovsky 公式的一个很好的开源实现在这里： http: //gmplib.org/pi-with-gmp.html

于 2011-09-01T05:33:55.713 回答
4

当 BBP 公式主要用于计算 π 的任意十六进制数字时，您似乎正在尝试计算 π 的十进制数字。基本上，BBP 公式可以用来计算π 的第n个十六进制数字，而无需计算前面的数字，十六进制数字 0, 1, ..., n - 1。

David H. Bailey（Bailey-Borwein-Plouffe 的 Bailey）编写了C 和 Fortran 代码来使用 BBP 公式计算 π 的第n个十六进制数字。在具有 IEEE 754 双算术的机器上，从 0开始计算，精确到n  ≈ 1.18 × 10 7 ；即 π = (3.243F6A8...) 16所以当n = 3 时程序的输出以“F”开头：

位置 = 3
 分数 = 0.963509103793105
 十六进制数字 = F6A8885A30

我喜欢稍微修改 C 版本，以便n（在代码中命名id）可以被命令行参数覆盖：

--- piqpr8.c.orig 2011-10-08 14:54:46.840423000 -0400
+++ piqpr8.c 2011-10-08 15:04:41.524437000 -0400
@@ -14,14 +14,18 @@
 /* 大卫 H. 贝利 2006-09-08 */

 #include <stdio.h>
+#include <stdlib.h>
 #include <math.h>

-主要的（）
+int main(int argc, char *argv[])
 {
   双pid，s1，s2，s3，s4；
   双系列（int m，int n）；
   void ihex (double x, int m, char c[]);
   int id = 1000000;
+ if (argc == 2) {
+ id = atoi(argv[1]);
+ }
 #define NHX 16
   字符 chx[NHX];

@@ -36,6 +40,8 @@
   ihex (pid, NHX, chx);
   printf (" 位置 = %i\n 分数 = %.15f \n 十六进制数字 = %10.10s\n",
   ID，PID，CHX）；
+
+ 返回 EXIT_SUCCESS；
 }

 void ihex (double x, int nhx, char chx[])
于 2011-10-08T19:24:19.593 回答
4

BBP 公式不适合轻松查找第 n 个十进制数字，因为它很容易返回十六进制且仅返回十六进制数字。因此，要重新计算成小数，您需要收集所有十六进制数字。

使用牛顿公式要好得多：

Pi/2 = 1 + 1/3 + 1*2/3*5 + 1*2*3/3*5*7 + .... n!/(2n+1)！！+ ....

它崩溃到霍纳的计划：

Pi/2 = 1 + 1/3*(1 + 2/5*(1 + 3/7*(1 + ...... n/(2n+1)*(1) ......) ))

因此，您将 Pi 写为位置序列，其中在每个小数位置上使用不同的基数 (n/(2n+1))，并且所有数字都等于 2。它显然收敛，因为该基数小于 1/2 ，因此要计算最多 n 个有效十进制 dgit 的 Pi，您只需要 log_2(10)*n 项（N = 10*n/3+1 是完美的东西）。

您从 N 个整数元素的数组开始，所有元素都等于 2，然后重复 n 次，执行以下操作：

1.) 将所有元素乘以 10。

2.) 重新计算每个元素[k]（从 N 到 1）的“数字”小于分母 (2*k+1)，
但同时您需要将 qoutient 移动到左侧位置，所以:
q = 元素[k] / (2*k+1);
元素[k] %= (2*k+1);
元素[k-1] += q * k; //k 是计数器，所以不要忘记相乘。

3.) 取元素[0]。它等于 10 * 第一个数字，所以你需要输出 element[0] / 10 并存储
element[0] %= 10;

但是有一个线索：牛顿公式的最大可能数字 (2*n) 的最大和是 2。所以你可以从元素 [1] 中获得多达 19/10。添加到 element[0] 时（在步骤 1 中乘以 10。）您可以获得 90+19=109。所以有时会发生输出数字是[10]。在这种情况下，您知道正确的数字是 0，并且必须将 1 添加到先前输出的数字上。

有两种方法可以解决这个问题：

1.) 在计算下一位之前不要输出最后一位。此外，存储连续 9 的数量并将它们输出为 9 或 1 后跟零，具体取决于第一个非 9 位。

2.) 将输出的数字放入结果数组中，这样如果出现 [10] 就可以轻松加 1。

在我的 PC 上，我可以在 10 秒内（用 Java 计算）10,000 个十进制数字。复杂度为 O(n^2)。

element[k] 的值永远不会超过 12*k，因此在快速机器上使用 64 位长类型可以计算超过 10^15 位（非常强大的近似值）。

于 2012-10-03T11:20:22.157 回答