Bailey–Borwein–Plouffe 公式在 C++ 中的实现？

Question

编辑：要求很模糊，他们不计算 pi 的第 n 位数字，而是希望将 pi 计算到第 n 位数字而不超出浮点数限制，因此暴力方式可以满足要求。

我需要计算 PI 第 n 位数字，我想尝试使用 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;
}

Answer 1

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

Chudnovsky 公式是计算 Pi 最快的已知公式，每项收敛于 14 位数字。然而，有效实施却极其困难。

由于该公式的复杂性，将 Pi 计算为小于几千位是没有意义的。因此，除非您准备好全力以赴进行任意精度算术，否则不要使用它。

使用 GMP 库对 Chudnovsky Formula 进行良好的开源实现如下：http://gmplib .org/pi-with-gmp.html

Answer 2

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

David H. Bailey（Bailey–Borwein–Plouffe 的贝利）撰写了 C 和 Fortran 代码，使用 BBP 公式计算 π 的n^第位十六进制数字。在采用 IEEE 754 双精度算术的机器上，从 0 开始计数，精确度可达 n  ≈ 1.18 × 10⁷；即 π = (3.243F6A8...)₁₆ 所以当 n = 3 时程序的输出以“F”开头：

 position = 3
 fraction = 0.963509103793105
 hex digits =  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 @@
 /*  David H. Bailey     2006-09-08 */

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

-main()
+int main(int argc, char *argv[])
 {
   double pid, s1, s2, s3, s4;
   double series (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
   char chx[NHX];

@@ -36,6 +40,8 @@
   ihex (pid, NHX, chx);
   printf (" position = %i\n fraction = %.15f \n hex digits =  %10.10s\n",
   id, pid, chx);
+
+  return EXIT_SUCCESS;
 }

 void ihex (double x, int nhx, char chx[])

Answer 3

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，因此要计算 Pi 最多 n 个有效十进制数字，您只需要 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。因此，您可以从 element[1] 获得多达 19/10。当添加到 element[0]（在步骤 1 中乘以 10）时，您可以得到 90+19=109。所以有时输出的数字会是[10]。在这种情况下，您知道正确的数字是 0，并且必须将 1 添加到先前输出的数字上。

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

1.) 在计算下一位之前不输出最后一位数字。此外，存储连续 9 的数量，并根据第一个非 9 数字将它们输出为 9 或 1 后跟 0。

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

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

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