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

Question

编辑：需求不够明确，他们只需要π的第n位而不是计算π的第n位，因此暴力方法符合要求。

我需要计算π的第n位，并尝试使用BBP公式，但遇到了困难。我打出的等式似乎没有正确地给出π。

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

我曾经用暴力方法成功地找到了π的值，但这种方法只有一定的精度，而且要找到第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;
}

- krizzo

---

你期望多少精度？这与你使用的类型支持的精度相比如何？算法的数字属性呢？减号总是意味着需要担心精度损失。 - dmckee --- ex-moderator kitten

我想计算圆周率，因为我们知道它要么正确，要么不正确（除了最后一位可能被四舍五入）。程序会提示用户需要多少个有效数字的圆周率，然后进行计算。据我所知，BBP公式将对0到无穷大的每个数字求和。每次都会增加一个圆周率数字。我将添加我的代码以帮助理解我的意图。 - krizzo

内置的浮点数表示仅支持6-7（32位）或15-16（64位）（可能还有17-18（80位））十进制数字的精度。要获得更多精度，您将需要使用某种任意精度包。互联网上有一份名为“计算机科学家应该了解的浮点运算知识”的文档。您需要阅读它。 - dmckee --- ex-moderator kitten

你是想要计算从第一位到第N位的所有数字吗？或者你只是想获得第N个二进制数字周围的一些数字？如果你想要所有数字，那么没有理由使用BBP公式，因为对于Pi来说有更快的公式。如果你想要从第一位到第N位的所有数字，你需要用到任意精度算术来完成它。 - Mysticial

BBP公式可以用来计算所有数字，但有更快的公式可以做到这一点。（参见Chudnovsky公式：http://en.wikipedia.org/wiki/Chudnovsky_algorithm）- BBP公式的威力在于它允许您直接计算Pi的第N个二进制数字，而无需计算它之前的数字。 - Mysticial

Answer 1

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

Chudnovsky公式是计算圆周率最快的已知公式，并且每个项收敛于14位数字。然而，实现它的效率非常困难。

由于这个公式的复杂性，除非你准备好全力以赴地使用任意精度算术，否则没有必要使用它来计算少于几千位数字的圆周率。

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

- Mysticial

---

这篇原始的BBP论文明确指出，他们的方法受益于不需要任意精度算术。那么，如何在这种情况下证明你所说的内容是正确的呢？ - JeremyKun

@JeremyKun BBP提取算法的问题在于，随着数字位数的增加，数值误差会以log(n)的速度累积。因此，如果您使用双精度浮点数，即53位精度，那么在求和超过2^53项之前，您就已经失去了所有结果，因为舍入误差太大。 - Mysticial

舍入误差+您希望提取的数字数量必须小于53位。所以，是的，如果您正在寻找一小部分不太远离小数点的数字，则可以避免所有任意精度算术。但是您将无法在任意偏移处提取100个数字。如果要使用BBP公式一次提取多个数字，则需要进行比字长更宽的算术运算。至少，您需要能够除以两个字长整数并保存长度超过一个字的结果。 - Mysticial

在所有情况下，您都可以避免执行大整数乘法 - 这被认为是任意精度算术中更困难的部分。我相信这就是论文所指的内容。 - Mysticial

通常，只需要一个数字N作为输入的问题（例如素性、分解等）是以输入字符串长度log(N)来评估的。这使得BBP成为一种指数算法。特别地，BBP并不排除使用π的数字作为加密安全的伪随机数生成器的可能性。 - JeremyKun

Answer 2

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

使用Newton公式要好得多：

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

它折叠为Horner方案：

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

因此，您可以将Pi写成位置系列，其中在每个小数位上使用不同的基数（n/(2n+1)），并且所有数字都等于2。显然它收敛，因为该基数小于1/2，因此要计算高达n个有效数字的Pi，您最多只需要log_2(10)*n个术语（N = 10*n/3+1是完美的内容）。

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

1.将所有元素乘以10。

2.重新计算每个元素[k]（从N到1）以使其“数字”小于分母（2*k+1），但同时还需要将商移动到左侧位置，因此：

q = element[k] / (2*k+1)； element[k] %= (2*k+1)； element[k-1] += q * k; //k是计数器，所以不要忘记乘。

3.取element[0]。它等于第一个数字的10倍，因此您需要输出element[0] / 10并存储

element[0] %= 10;

但是有一个诀窍：Newton公式最大可能位数（2*n）的最大和为2。因此，您可以从element[1]中获得高达19/10。当添加到element[0]（在步骤1中乘以10）时，您可以获得90 + 19 = 109。因此，输出的数字有时将是[10]。在这种情况下，您知道正确的数字是0，并且必须将1添加到先前输出的数字。

有两种解决此问题的方法：

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

2.将输出的数字放入结果数组中，因此如果出现[10]，您可以轻松地添加1。

在我的PC上，我可以在10秒内计算出Java的10,000位小数。复杂度为O(n^2)。

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

- Tomek

---

BBP不仅适用于十六进制数字，它本身就可以产生二进制数字。这可以轻松地映射到十六进制数字，因为16是2的幂，而10不是。您还可以使用BBP生成例如 pi 的第N个八进制数字。计算第n个十六进制数字只需要将[n，n + 3] 二进制位一起计算即可。 - Daniel Papasian

@DanielPapasian：当然，但试着切换到十进制，你会看到问题所在。那么，位流中任意位置的位数，如果不知道其前面的位数是什么以及它们对需要了解的十进制数字的影响，就无法帮助你太多。（你知道的，来自较低位置的进位等等。） - SasQ

Answer 3

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

David H. Bailey（Bailey-Borwein-Plouffe中的Bailey）编写了C和Fortran代码，使用BBP公式计算π的第n个十六进制数字。在具有IEEE 754双精度算术的机器上，从0开始计数，它的精度高达n≈1.18×107；即π＝（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+ 
/* David H. Bailey 2006-09-08 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
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];
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[])

- Daniel Trebbien