performance - 获得 π 值的最快方法是什么？

Question

作为个人挑战，我正在寻找获得 π 值的最快方法。更具体地说，我使用的方法不涉及使用 #define 常量(例如 M_PI)，或者将数字硬编码到其中。

下面的程序测试了我所知道的各种方法。理论上，内联汇编版本是最快的选择，尽管显然不可移植。我将其作为与其他版本进行比较的基准。在我的测试中，使用内置函数，4 * atan(1) 版本在 GCC 4.2 上最快，因为它自动将 atan(1) 折叠成常量.指定 -fno-builtin 后，atan2(0, -1) 版本最快。

这里是主要的测试程序(pitimes.c):

#include <math.h>
#include <stdio.h>
#include <time.h>

#define ITERS 10000000
#define TESTWITH(x) {                                                       \
    diff = 0.0;                                                             \
    time1 = clock();                                                        \
    for (i = 0; i < ITERS; ++i)                                             \
        diff += (x) - M_PI;                                                 \
    time2 = clock();                                                        \
    printf("%s\t=> %e, time => %f\n", #x, diff, diffclock(time2, time1));   \
}

static inline double
diffclock(clock_t time1, clock_t time0)
{
    return (double) (time1 - time0) / CLOCKS_PER_SEC;
}

int
main()
{
    int i;
    clock_t time1, time2;
    double diff;

    /* Warmup. The atan2 case catches GCC's atan folding (which would
     * optimise the ``4 * atan(1) - M_PI'' to a no-op), if -fno-builtin
     * is not used. */
    TESTWITH(4 * atan(1))
    TESTWITH(4 * atan2(1, 1))

#if defined(__GNUC__) && (defined(__i386__) || defined(__amd64__))
    extern double fldpi();
    TESTWITH(fldpi())
#endif

    /* Actual tests start here. */
    TESTWITH(atan2(0, -1))
    TESTWITH(acos(-1))
    TESTWITH(2 * asin(1))
    TESTWITH(4 * atan2(1, 1))
    TESTWITH(4 * atan(1))

    return 0;
}

以及仅适用于 x86 和 x64 系统的内联汇编内容 (fldpi.c):

double
fldpi()
{
    double pi;
    asm("fldpi" : "=t" (pi));
    return pi;
}

以及构建我正在测试的所有配置的构建脚本 (build.sh):

#!/bin/sh
gcc -O3 -Wall -c           -m32 -o fldpi-32.o fldpi.c
gcc -O3 -Wall -c           -m64 -o fldpi-64.o fldpi.c

gcc -O3 -Wall -ffast-math  -m32 -o pitimes1-32 pitimes.c fldpi-32.o
gcc -O3 -Wall              -m32 -o pitimes2-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -fno-builtin -m32 -o pitimes3-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -ffast-math  -m64 -o pitimes1-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall              -m64 -o pitimes2-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall -fno-builtin -m64 -o pitimes3-64 pitimes.c fldpi-64.o -lm

除了在各种编译器标志之间进行测试(我也比较了 32 位和 64 位，因为优化不同)，我还尝试切换测试的顺序。但是，atan2(0, -1) 版本仍然每次都名列前茅。

Answer 1

Monte Carlo method ，如前所述，应用了一些伟大的概念，但显然，它不是最快的，不是远射，也不是任何合理的衡量标准。此外，这完全取决于您正在寻找什么样的准确性。我所知道的最快的 π 是数字硬编码的那个。看着Pi和 Pi[PDF] , 有很多公式。

这是一种快速收敛的方法——每次迭代大约 14 个数字。 PiFast ，当前最快的应用程序，将此公式与 FFT 结合使用。我只写公式，因为代码很简单。 Ramanujan and discovered by Chudnovsky差点找到这个公式.这实际上是他如何计算出数十亿位数的数字——所以这不是一种可以忽视的方法。公式会很快溢出，并且由于我们要划分阶乘，因此延迟此类计算以删除项将是有利的。

在哪里，

下面是Brent–Salamin algorithm .维基百科提到，当 a 和 b “足够接近”时，(a + b)²/4t 将是 π 的近似值。我不确定“足够接近”是什么意思，但根据我的测试，一次迭代得到 2 位数字，两次迭代得到 7 位数字，三次迭代得到 15 位数字，当然这是 double ，所以它可能有基于其表示的错误和true 计算可能更准确。

let pi_2 iters =
    let rec loop_ a b t p i =
        if i = 0 then a,b,t,p
        else
            let a_n = (a +. b) /. 2.0 
            and b_n = sqrt (a*.b)
            and p_n = 2.0 *. p in
            let t_n = t -. (p *. (a -. a_n) *. (a -. a_n)) in
            loop_ a_n b_n t_n p_n (i - 1)
    in 
    let a,b,t,p = loop_ (1.0) (1.0 /. (sqrt 2.0)) (1.0/.4.0) (1.0) iters in
    (a +. b) *. (a +. b) /. (4.0 *. t)

最后，来点 pi golf(800 位)怎么样？ 160 个字符!

int a=10000,b,c=2800,d,e,f[2801],g;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);}