什么是获取π值最快的方法？

Question

什么是获取π值最快的方法？

performancealgorithmlanguage-agnosticunixpi

353

我正在寻找最快的方法来获取π的值，这是一个个人挑战。更具体地说，我正在使用不涉及使用#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)版本最优。

- C. K. Young

2

你为什么认为使用atan(1)与使用M_PI不同？如果你只使用算术运算，我可以理解你为什么这样做，但是对于atan，我看不出有什么意义。 - erikkallen

11

为什么你不想使用常量？比如由库定义或自己定义的常量？计算圆周率是浪费CPU周期，因为该问题已经被反复解决，并得出了比日常计算所需更多的有效数字。 - Tilo

1

除了预先计算常数π之外，只有一种解决方案更快：预先计算公式中出现的所有值，例如当需要周长时，您可以预先计算2*PI而不是在运行时每次将PI乘以2。 - ern0

5

在我所说的英语方言中，“optimisation”的拼写是带有“s”而不是“z”的（顺便提一下，“z”读作“zed”，而不是“zee” ;-))。如果您查看审查历史记录，您会发现这并不是我第一次撤销这种编辑。 - C. K. Young

3

@Pacerier 请见http://en.wiktionary.org/wiki/boggle 和 http://en.wiktionary.org/wiki/mindboggling。 - C. K. Young

显示剩余12条评论

24个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- NihilistDandy · Answer 1

使用 Machin 式公式。

176 * arctan (1/57) + 28 * arctan (1/239) - 48 * arctan (1/682) + 96 * arctan(1/12943) 

[; \left( 176 \arctan \frac{1}{57} + 28 \arctan \frac{1}{239} - 48 \arctan \frac{1}{682} + 96 \arctan \frac{1}{12943}\right) ;], for you TeX the World people.

例如，在Scheme中实现：

(+ (- (+ (* 176 (atan (/ 1 57))) (* 28 (atan (/ 1 239)))) (* 48 (atan (/ 1 682)))) (* 96 (atan (/ 1 12943))))

- Daniel C. Sobral · Answer 2

如果您愿意使用近似值，355 / 113 可以保留6位小数，并且具有与整数表达式一起使用的优点。这在当今并不像以前那样重要，因为"浮点数协处理器"已经没有了任何意义，但它曾经非常重要。

- Jason Delife · Answer 3

派的确是3！[弗林克教授（辛普森一家）]

开个玩笑，这里有一个用C#编写的笑话（需要.NET框架）。

using System;
using System.Text;

class Program {
    static void Main(string[] args) {
        int Digits = 100;

        BigNumber x = new BigNumber(Digits);
        BigNumber y = new BigNumber(Digits);
        x.ArcTan(16, 5);
        y.ArcTan(4, 239);
        x.Subtract(y);
        string pi = x.ToString();
        Console.WriteLine(pi);
    }
}

public class BigNumber {
    private UInt32[] number;
    private int size;
    private int maxDigits;

    public BigNumber(int maxDigits) {
        this.maxDigits = maxDigits;
        this.size = (int)Math.Ceiling((float)maxDigits * 0.104) + 2;
        number = new UInt32[size];
    }
    public BigNumber(int maxDigits, UInt32 intPart)
        : this(maxDigits) {
        number[0] = intPart;
        for (int i = 1; i < size; i++) {
            number[i] = 0;
        }
    }
    private void VerifySameSize(BigNumber value) {
        if (Object.ReferenceEquals(this, value))
            throw new Exception("BigNumbers cannot operate on themselves");
        if (value.size != this.size)
            throw new Exception("BigNumbers must have the same size");
    }

    public void Add(BigNumber value) {
        VerifySameSize(value);

        int index = size - 1;
        while (index >= 0 && value.number[index] == 0)
            index--;

        UInt32 carry = 0;
        while (index >= 0) {
            UInt64 result = (UInt64)number[index] +
                            value.number[index] + carry;
            number[index] = (UInt32)result;
            if (result >= 0x100000000U)
                carry = 1;
            else
                carry = 0;
            index--;
        }
    }
    public void Subtract(BigNumber value) {
        VerifySameSize(value);

        int index = size - 1;
        while (index >= 0 && value.number[index] == 0)
            index--;

        UInt32 borrow = 0;
        while (index >= 0) {
            UInt64 result = 0x100000000U + (UInt64)number[index] -
                            value.number[index] - borrow;
            number[index] = (UInt32)result;
            if (result >= 0x100000000U)
                borrow = 0;
            else
                borrow = 1;
            index--;
        }
    }
    public void Multiply(UInt32 value) {
        int index = size - 1;
        while (index >= 0 && number[index] == 0)
            index--;

        UInt32 carry = 0;
        while (index >= 0) {
            UInt64 result = (UInt64)number[index] * value + carry;
            number[index] = (UInt32)result;
            carry = (UInt32)(result >> 32);
            index--;
        }
    }
    public void Divide(UInt32 value) {
        int index = 0;
        while (index < size && number[index] == 0)
            index++;

        UInt32 carry = 0;
        while (index < size) {
            UInt64 result = number[index] + ((UInt64)carry << 32);
            number[index] = (UInt32)(result / (UInt64)value);
            carry = (UInt32)(result % (UInt64)value);
            index++;
        }
    }
    public void Assign(BigNumber value) {
        VerifySameSize(value);
        for (int i = 0; i < size; i++) {
            number[i] = value.number[i];
        }
    }

    public override string ToString() {
        BigNumber temp = new BigNumber(maxDigits);
        temp.Assign(this);

        StringBuilder sb = new StringBuilder();
        sb.Append(temp.number[0]);
        sb.Append(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.CurrencyDecimalSeparator);

        int digitCount = 0;
        while (digitCount < maxDigits) {
            temp.number[0] = 0;
            temp.Multiply(100000);
            sb.AppendFormat("{0:D5}", temp.number[0]);
            digitCount += 5;
        }

        return sb.ToString();
    }
    public bool IsZero() {
        foreach (UInt32 item in number) {
            if (item != 0)
                return false;
        }
        return true;
    }

    public void ArcTan(UInt32 multiplicand, UInt32 reciprocal) {
        BigNumber X = new BigNumber(maxDigits, multiplicand);
        X.Divide(reciprocal);
        reciprocal *= reciprocal;

        this.Assign(X);

        BigNumber term = new BigNumber(maxDigits);
        UInt32 divisor = 1;
        bool subtractTerm = true;
        while (true) {
            X.Divide(reciprocal);
            term.Assign(X);
            divisor += 2;
            term.Divide(divisor);
            if (term.IsZero())
                break;

            if (subtractTerm)
                this.Subtract(term);
            else
                this.Add(term);
            subtractTerm = !subtractTerm;
        }
    }
}

- qwerty01 · Answer 4

使用双精度浮点数：

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

这将精确到小数点后14位，足以填满一个double（不准确可能是因为弧正切函数中的其余小数被截断了）。

此外，Seth，它是3.141592653589793238463，而不是64。

- Brad Gilbert · Answer 5

使用D语言在编译时计算圆周率。

（摘自DSource.org）

/** Calculate pi at compile time
 *
 * Compile with dmd -c pi.d
 */
module calcpi;

import meta.math;
import meta.conv;

/** real evaluateSeries!(real x, real metafunction!(real y, int n) term)
 *
 * Evaluate a power series at compile time.
 *
 * Given a metafunction of the form
 *  real term!(real y, int n),
 * which gives the nth term of a convergent series at the point y
 * (where the first term is n==1), and a real number x,
 * this metafunction calculates the infinite sum at the point x
 * by adding terms until the sum doesn't change any more.
 */
template evaluateSeries(real x, alias term, int n=1, real sumsofar=0.0)
{
  static if (n>1 && sumsofar == sumsofar + term!(x, n+1)) {
     const real evaluateSeries = sumsofar;
  } else {
     const real evaluateSeries = evaluateSeries!(x, term, n+1, sumsofar + term!(x, n));
  }
}

/*** Calculate atan(x) at compile time.
 *
 * Uses the Maclaurin formula
 *  atan(z) = z - z^3/3 + Z^5/5 - Z^7/7 + ...
 */
template atan(real z)
{
    const real atan = evaluateSeries!(z, atanTerm);
}

template atanTerm(real x, int n)
{
    const real atanTerm =  (n & 1 ? 1 : -1) * pow!(x, 2*n-1)/(2*n-1);
}

/// Machin's formula for pi
/// pi/4 = 4 atan(1/5) - atan(1/239).
pragma(msg, "PI = " ~ fcvt!(4.0 * (4*atan!(1/5.0) - atan!(1/239.0))) );

- kogus · Answer 6

这个（使用Delphi编写的）版本并不是特别出色，但至少比Nick Hodge在博客上发布的版本更快：）。在我的计算机上，进行十亿次迭代需要大约16秒钟，得到一个值为3.1415926525879（精确的部分用粗体表示）。

program calcpi;

{$APPTYPE CONSOLE}

uses
  SysUtils;

var
  start, finish: TDateTime;

function CalculatePi(iterations: integer): double;
var
  numerator, denominator, i: integer;
  sum: double;
begin
  {
  PI may be approximated with this formula:
  4 * (1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 .......)
  //}
  numerator := 1;
  denominator := 1;
  sum := 0;
  for i := 1 to iterations do begin
    sum := sum + (numerator/denominator);
    denominator := denominator + 2;
    numerator := -numerator;
  end;
  Result := 4 * sum;
end;

begin
  try
    start := Now;
    WriteLn(FloatToStr(CalculatePi(StrToInt(ParamStr(1)))));
    finish := Now;
    WriteLn('Seconds:' + FormatDateTime('hh:mm:ss.zz',finish-start));
  except
    on E:Exception do
      Writeln(E.Classname, ': ', E.Message);
  end;
end.

- Seth · Answer 7

如果您想要计算π的近似值（出于某种原因），您应该尝试使用二进制提取算法。 Bellard 对 BBP 的改进使得PI的复杂度为O(N^2)。

如果您想要获取π的近似值以进行计算，那么：

PI = 3.141592654

当然，这只是一个近似值，不完全准确。它的偏差超过0.00000000004102（四万亿分之四），大约是⁴/_{10,000,000,000}。

如果您想使用π进行数学计算，则需要准备铅笔和纸或计算机代数软件，并使用π的精确值π。

如果您真的需要一个公式，这个很有趣：

π = -i ln(-1)

- Kristopher Johnson · Answer 8

回到早期，由于字长较小且浮点运算速度缓慢或不存在，我们过去常常这样操作：

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

对于不需要太高精度的应用程序（例如视频游戏），这非常快速并且足够准确。

- user14273431 · Answer 9

基本上是纸夹优化器答案的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公式，如果您不太理解代码，我可以完全简化它。

总结：我们将从1到5得到一个数字，并将其添加到我们将用于获取PI的函数中。然后向您提供3个数字：545140134（A），13591409（B），640320（D）。然后，我们将使用D作为double将其自身乘以3次成为另一个double（ID3）。然后我们将ID3的平方根带入另一个double（b）中，并分配2个数字：1（p）和B的值（a）。请注意，C不区分大小写。然后通过将p，a和b的值全部乘以double来创建一个double（sum）。然后开始循环直到给定函数的数字，将A的值加起来，将b的值乘以ID3，将p的值乘以多个值并除以多个值。总和将再次由p，a和b相加，循环将重复，直到循环的数字的值大于或等于5。稍后，将12乘以总和，并通过函数返回，从而给出PI的结果。

好的，那太长了，但我想你会理解的...

- Rajanand · Answer 10

我认为圆周长和半径之间的比值是圆周率的值。它可以通过简单的数学计算得出。