计算圆周率的数字

4
我使用GMP库和C++编写了高斯-勒让德算法的实现来计算pi的位数。它有正确的输出,但问题在于我不知道输出何时“变坏”,因为我必须在代码中指定精度。以下是使用64位精度的输出:3.141592653589793238*35*,最后两个数字是不正确的。我的问题是,如果我想要n位pi,需要多少位精度b和算法迭代次数i?谢谢。

是的,可以去http://chat.stackoverflow.com/rooms/10/loungec与Mysticial联系。如果有一个类似的方法可以联系Dietmar Kühl关于iostream问题就好了... - Mooing Duck
1个回答

10

高斯-勒让德算法(又称AGM算法)需要始终保持完全精度。

与牛顿迭代法不同,AGM迭代没有自我修正功能。因此,您需要从一开始就保持完全精度。此外,您需要额外的保护位数。

我的问题是,如果我想要π的n位小数,需要多少位精度b

首先,您需要将n(十进制)位转换为b位二进制位。所以是:

        log(10)
b = n ---------- + epsilon
        log(2)

其中epsilon是保护位的数量。需要多少取决于实现和舍入行为,但通常对于任何迭代次数,100位足够了。

至于需要多少迭代,这可以通过经验数据找到。

以下是我编写的一个小应用程序使用 Gauss-Legendre 算法计算 100 百万位数的圆周率的输出:

================================================================
Computing pi to 100000000 digits:
Threads: 8

Starting AGM...         1.394965 seconds
Starting Iteration 0...    -3 (error in decimal digits)
Starting Iteration 1...    -9
Starting Iteration 2...    -20
Starting Iteration 3...    -42
Starting Iteration 4...    -85
Starting Iteration 5...    -173
Starting Iteration 6...    -347
Starting Iteration 7...    -696
Starting Iteration 8...    -1395
Starting Iteration 9...    -2792
Starting Iteration 10...    -5586
Starting Iteration 11...    -11175
Starting Iteration 12...    -22352
Starting Iteration 13...    -44706
Starting Iteration 14...    -89414
Starting Iteration 15...    -178829
Starting Iteration 16...    -357661
Starting Iteration 17...    -715324
Starting Iteration 18...    -1430650
Starting Iteration 19...    -2861302
Starting Iteration 20...    -5722607
Starting Iteration 21...    -11445216
Starting Iteration 22...    -22890435
Starting Iteration 23...    -45780871
Starting Iteration 24...    -91561745
Starting Iteration 25...    -183123492
AGM:                    118.796792 seconds
Finishing...            3.033239   seconds
Radix Conversion...     2.924844   seconds

Total Wall Time:        126.151086 seconds

CPU Utilization:        495.871%
CPU Efficiency:         61.984%

Writing to "pi.txt"...  Done

因此,25次迭代足以处理1.83亿位数字。同样,每次迭代都会翻倍,因此您可以运用一些基本的对数数学知识来计算需要多少次迭代。


谢谢,我已经有一个基本的实现,但是它非常慢。我做了明显的改进,避免在每次迭代中重新计算阶乘,但整体代码仍然相当慢。例如,10000位数字需要8秒钟。如果有人想看一下我的代码,这是它的链接:http://pastebin.com/9Bs2WHRV 编辑:抱歉,那是我的Chudnovsky实现,Gauss-Legendre实现还可以。 - Sam Kennedy
是的,Chudnovsky算法是一种不同的动物。它更难实现,因为您需要将其与二分法结合使用。但是,如果正确执行,它非常快速。 - Mysticial
错误是如何计算的?这8个单独的线程都做什么?平方根操作非常昂贵,所以我可以理解为什么要将其放在单独的线程中,但其他操作所需时间相对较短。 - Sam Kennedy
错误计算为前一次迭代的差值的两倍。线程不相关。该应用程序使用与y-cruncher相同的后端,该后端在大量乘法中进行并行化处理。它根本不使用GMP库。实际上,在顶层没有进行任何并行处理。 - Mysticial
起始错误是什么?看起来像是2点什么,但我找不到任何地方提到它。 - Sam Kennedy
我不知道确切的数字——因为它既不是2也不是3(介于两者之间)。建议使用后面的迭代作为起点来确定需要多少次迭代。在我的具体实现中,我只需循环直到与上一次迭代的差异超过最终精度的一半。因此,在这种意义上,您甚至不需要知道需要多少次迭代。 - Mysticial

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接