我正在寻找最快的方法来获得π的值,作为一个个人挑战。更具体地说,我使用的方法不涉及使用#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;
}

内联汇编的东西(fldpi.c)只适用于x86和x64系统:

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)版本在每次测试中仍然名列前茅。


当前回答

蒙特卡罗方法,如前所述,应用了一些伟大的概念,但很明显,它不是最快的,不是从任何合理的标准来看。此外,这完全取决于你想要什么样的准确性。我所知道的最快的π是数字硬编码的π。看看圆周率和圆周率,有很多公式。

Here is a method that converges quickly — about 14 digits per iteration. PiFast, the current fastest application, uses this formula with the FFT. I'll just write the formula, since the code is straightforward. This formula was almost found by Ramanujan and discovered by Chudnovsky. It is actually how he calculated several billion digits of the number — so it isn't a method to disregard. The formula will overflow quickly and, since we are dividing factorials, it would be advantageous then to delay such calculations to remove terms.

在那里,

下面是Brent-Salamin算法。维基百科提到,当a和b“足够接近”时,(a + b)²/ 4t将是π的近似值。我不确定“足够接近”是什么意思,但从我的测试来看,一次迭代得到2位数字,两次得到7位,3次得到15位,当然这是双精度,所以它可能会有一个基于它的表示的错误,真实的计算可能会更准确。

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)

最后,来点圆周率高尔夫(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);}

其他回答

我总是使用acos(-1),而不是将π定义为常数。

蒙特卡罗方法,如前所述,应用了一些伟大的概念,但很明显,它不是最快的,不是从任何合理的标准来看。此外,这完全取决于你想要什么样的准确性。我所知道的最快的π是数字硬编码的π。看看圆周率和圆周率,有很多公式。

Here is a method that converges quickly — about 14 digits per iteration. PiFast, the current fastest application, uses this formula with the FFT. I'll just write the formula, since the code is straightforward. This formula was almost found by Ramanujan and discovered by Chudnovsky. It is actually how he calculated several billion digits of the number — so it isn't a method to disregard. The formula will overflow quickly and, since we are dividing factorials, it would be advantageous then to delay such calculations to remove terms.

在那里,

下面是Brent-Salamin算法。维基百科提到,当a和b“足够接近”时,(a + b)²/ 4t将是π的近似值。我不确定“足够接近”是什么意思,但从我的测试来看,一次迭代得到2位数字,两次得到7位,3次得到15位,当然这是双精度,所以它可能会有一个基于它的表示的错误,真实的计算可能会更准确。

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)

最后,来点圆周率高尔夫(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);}

我认为圆周率的值是圆的周长和半径之比。

它可以通过常规的数学计算简单地实现

下面是我在高中时学过的计算圆周率的技巧。

我之所以分享它,是因为我认为它足够简单,任何人都可以无限期地记住它,而且它教会了你“蒙特卡罗”方法的概念——这是一种统计方法,可以得到答案,这些答案不会立即通过随机过程演绎出来。

画一个正方形,在这个正方形内画一个象限(半圆的四分之一)(一个半径等于正方形边的象限,这样它就能尽可能多地填充正方形)

现在向正方形投掷飞镖,并记录飞镖落在何处——也就是说,在正方形内任意选择一个点。当然,它落在了正方形内部,但它落在半圆内部吗?记录这个事实。

重复此过程多次,你会发现半圆内的点数量与抛出的总数量之比为x。

由于正方形的面积是r乘以r,可以推导出半圆的面积是x乘以r乘以r(即x乘以r的平方)。因此x乘以4会得到。

这不是一个快速使用的方法。但这是蒙特卡罗方法的一个很好的例子。如果你环顾四周,你可能会发现许多超出你计算能力的问题都可以用这种方法来解决。

双打:

4.0 * (4.0 * Math.Atan(0.2) - Math.Atan(1.0 / 239.0))

这将精确到小数点后14位,足以填充双精度(不准确可能是因为弧切线中的其余小数被截断了)。

还有Seth,是3.141592653589793238463,不是64。