求π的最快方法是什么?

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

当前回答

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

2020-05-08 00:52:03

其他回答

从圆面积计算π:-)

<input id="range" type="range" min="10" max="960" value="10" step="50" oninput="calcPi()"> <br> <div id="cont"></div> <script> function generateCircle(width) { var c = width/2; var delta = 1.0; var str = ""; var xCount = 0; for (var x=0; x <= width; x++) { for (var y = 0; y <= width; y++) { var d = Math.sqrt((x-c)*(x-c) + (y-c)*(y-c)); if (d > (width-1)/2) { str += '.'; } else { xCount++; str += 'o'; } str += " " } str += "\n"; } var pi = (xCount * 4) / (width * width); return [str, pi]; } function calcPi() { var e = document.getElementById("cont"); var width = document.getElementById("range").value; e.innerHTML = "<h4>Generating circle...</h4>"; setTimeout(function() { var circ = generateCircle(width); e.innerHTML = "<pre>" + "π = " + circ[1].toFixed(2) + "\n" + circ[0] +"</pre>"; }, 200); } calcPi(); </script>

2017-06-03 17:13:12

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

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

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

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

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

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

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

2008-08-01 13:37:59

BBP公式允许你计算第n位数字-以2为基数(或16)-甚至不需要麻烦之前的n-1位:)

2008-08-29 09:22:16

如果你想计算π值的近似值(出于某种原因)，你应该尝试二进制提取算法。Bellard对BBP的改进给出了O(N²)中的PI。