我正在寻找最快的方法来获得π的值,作为一个个人挑战。更具体地说,我使用的方法不涉及使用#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)版本在每次测试中仍然名列前茅。


当前回答

实际上,有一本书专门介绍了π的快速计算方法:Jonathan和Peter Borwein写的《π和AGM》(在亚马逊上可以买到)。

我对年度股东大会和相关算法进行了相当多的研究:这非常有趣(尽管有时并非微不足道)。

请注意,要实现大多数现代算法来计算\pi,您将需要一个多精度算术库(GMP是一个很好的选择,尽管距离我上次使用它已经有一段时间了)。

最佳算法的时间复杂度为O(M(n)log(n)),其中M(n)是使用基于fft算法对两个n位整数(M(n)=O(n log(n) log(log(n)))相乘的时间复杂度,通常在计算\pi数字时需要fft算法,GMP中实现了该算法。

请注意,即使算法背后的数学可能并不简单,算法本身通常是几行伪代码,它们的实现通常非常简单(如果您选择不编写自己的多精度算术:-))。

其他回答

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

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

基本上是C版本的回形针优化器的答案,并且更加简化:

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

double calc_PI(int K) {
    static const int A = 545140134;
    static const int B = 13591409;
    static const int D = 640320;
    const double ID3 = 1.0 / ((double) D * (double) D * (double) D);
    double sum = 0.0;
    double b = sqrt(ID3);
    long long int p = 1;
    long long int a = B;
    sum += (double) p * (double) a * b;
    for (int k = 1; k < K; ++k) {
        a += A;
        b *= ID3;
        p *= (6 * k) * (6 * k - 1) * (6 * k - 2) * (6 * k - 3) * (6 * k - 4) * (6 * k - 5);
        p /= (3 * k) * (3 * k - 1) * (3 * k - 2) * k * k * k;
        p = -p;
        sum += (double) p * (double) a * b;
    }
    return 1.0 / (12 * sum);
}

int main() {
    for (int k = 1; k <= 5; ++k) {
        printf("k = %i, PI = %.16f\n", k, calc_PI(k));
    }
}

但为了更简化,这个算法采用Chudnovsky公式,如果你不太理解代码,我可以完全简化这个公式。

Summary: We will get a number from 1 to 5 and add it in to a function we will use to get PI. Then 3 numbers are given to you: 545140134 (A), 13591409 (B), 640320 (D). Then we will use D as a double multiplying itself 3 times into another double (ID3). We will then take the square root of ID3 into another double (b) and assign 2 numbers: 1 (p), the value of B (a). Take note that C is case-insensitive. Then a double (sum) will be created by multiplying the value's of p, a and b, all in doubles. Then a loop up until the number given for the function will start and add up A's value to a, b's value gets multiplied by ID3, p's value will be multiplied by multiple values that I hope you can understand and also gets divided by multiple values as well. The sum will add up by p, a and b once again and the loop will repeat until the value of the loop's number is greater or equal to 5. Later, the sum is multiplied by 12 and returned by the function giving us the result of PI.

好吧,这很长,但我想你会理解的……

在过去,由于字的大小很小,浮点运算很慢或者根本不存在,我们常常这样做:

/* Return approximation of n * PI; n is integer */
#define pi_times(n) (((n) * 22) / 7)

对于不需要很高精度的应用程序(例如电子游戏),这是非常快速和准确的。

Chudnovsky算法非常快如果你不介意做一个平方根和几个逆运算的话。它在2次迭代中收敛到两倍精度。

/*
    Chudnovsky algorithm for computing PI
*/

#include <iostream>
#include <cmath>
using namespace std;

double calc_PI(int K=2) {

    static const int A = 545140134;
    static const int B = 13591409;
    static const int D = 640320;

    const double ID3 = 1./ (double(D)*double(D)*double(D));

    double sum = 0.;
    double b   = sqrt(ID3);
    long long int p = 1;
    long long int a = B;

    sum += double(p) * double(a)* b;

    // 2 iterations enough for double convergence
    for (int k=1; k<K; ++k) {
        // A*k + B
        a += A;
        // update denominator
        b *= ID3;
        // p = (-1)^k 6k! / 3k! k!^3
        p *= (6*k)*(6*k-1)*(6*k-2)*(6*k-3)*(6*k-4)*(6*k-5);
        p /= (3*k)*(3*k-1)*(3*k-2) * k*k*k;
        p = -p;

        sum += double(p) * double(a)* b;
    }

    return 1./(12*sum);
}

int main() {

    cout.precision(16);
    cout.setf(ios::fixed);

    for (int k=1; k<=5; ++k) cout << "k = " << k << "   PI = " << calc_PI(k) << endl;

    return 0;
}

结果:

k = 1   PI = 3.1415926535897341
k = 2   PI = 3.1415926535897931
k = 3   PI = 3.1415926535897931
k = 4   PI = 3.1415926535897931
k = 5   PI = 3.1415926535897931

双打:

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

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

还有Seth,是3.141592653589793238463,不是64。