質問

個人的な課題として、π の値を最速で取得する方法を探しています。具体的には、を使わない方法を使っています。 #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) バージョンは今でも毎回トップになります。

役に立ちましたか?

解決

モンテカルロ法, 前述したように、これはいくつかの優れた概念を適用していますが、明らかに、最速ではなく、長期的な目標でも、合理的な尺度でもありません。また、それはすべて、どのような精度を求めるかによって異なります。私が知っている最速の π は、数字がハードコードされたものです。見つめている 円周率 そして 円周率[PDF], 公式がたくさんあります。

これは、反復ごとに約 14 桁で迅速に収束する方法です。 ピファスト, 現在最速のアプリケーションである は、FFT でこの式を使用します。コードは簡単なので、式を書きます。この公式はほぼ次のように発見されました。 ラマヌジャンとチュドノフスキーによって発見された. 。実際、これは彼が数十億桁の数値を計算した方法なので、無視できる方法ではありません。式はすぐにオーバーフローしてしまい、階乗を除算しているので、そのような計算を遅らせて項を削除する方が有利です。

enter image description here

enter image description here

どこ、

enter image description here

以下は、 ブレント・サラミンアルゴリズム. 。ウィキペディアにはそのときのことが記載されています ある そして b 「十分に近い」場合 (a + b)² / 4t はπの近似値になります。「十分に近い」が何を意味するのかはわかりませんが、テストによると、1 つの反復では 2 桁が得られ、2 つの反復では 7 が得られ、3 つの反復では 15 が得られました。もちろん、これは double を使用したものであるため、その表現に基づくエラーがある可能性があります。の 真実 計算はより正確になる可能性があります。

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);}

他のヒント

このプログラムは、独自の領域を見て π を近似するので、とても気に入っています。

IOCCC 1988 : ウェストリー.c

#define _ -F<00||--F-OO--;
int F=00,OO=00;main(){F_OO();printf("%1.3f\n",4.*-F/OO/OO);}F_OO()
{
            _-_-_-_
       _-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
        _-_-_-_-_-_-_-_
            _-_-_-_
}

ここでは私が高校で学んだ円周率の計算テクニックの概要を説明します。

私がこれを共有するのは、これが簡単で誰でも永久に覚えられると思うからです。さらに、「モンテカルロ」法の概念も教えてくれます。これは、すぐには答えが得られないと思われる答えに到達するための統計的手法です。ランダムなプロセスを通じて推定可能。

正方形を描き、その正方形の内側に四分円(半円の 4 分の 1)を内接します(四分円の半径が正方形の一辺に等しいため、正方形の大部分を埋めることができます)。

次に、ダーツを広場に投げて、どこに着地したかを記録します。つまり、広場内の任意の場所をランダムに選択します。もちろん四角の内側に着地しましたが、半円の内側でしょうか?この事実を記録してください。

このプロセスを何度も繰り返します。すると、半円内の点の数と投げられた総数の比があることがわかります。この比を x と呼びます。

正方形の面積は r × r であるため、半円の面積は x × r × r (つまり、x × r の 2 乗) であると推測できます。したがって、x に 4 を掛けると円周率が得られます。

これはすぐに使える方法ではありません。しかし、これはモンテカルロ法の良い例です。そして、周りを見回してみると、計算スキルの範囲外にある多くの問題が、そのような方法で解決できることに気づくかもしれません。

完全を期すため、C++ テンプレート バージョンは、最適化されたビルドの場合、コンパイル時に PI の近似値を計算し、単一の値にインライン化します。

#include <iostream>

template<int I>
struct sign
{
    enum {value = (I % 2) == 0 ? 1 : -1};
};

template<int I, int J>
struct pi_calc
{
    inline static double value ()
    {
        return (pi_calc<I-1, J>::value () + pi_calc<I-1, J+1>::value ()) / 2.0;
    }
};

template<int J>
struct pi_calc<0, J>
{
    inline static double value ()
    {
        return (sign<J>::value * 4.0) / (2.0 * J + 1.0) + pi_calc<0, J-1>::value ();
    }
};


template<>
struct pi_calc<0, 0>
{
    inline static double value ()
    {
        return 4.0;
    }
};

template<int I>
struct pi
{
    inline static double value ()
    {
        return pi_calc<I, I>::value ();
    }
};

int main ()
{
    std::cout.precision (12);

    const double pi_value = pi<10>::value ();

    std::cout << "pi ~ " << pi_value << std::endl;

    return 0;
}

I > 10 の場合、最適化されたビルドは、最適化されていない実行と同様に遅くなる可能性があることに注意してください。12 回の反復では、value() への呼び出しは約 80,000 回あると思います (メモ化がない場合)。

実際、(とりわけ)次のことに特化した本が 1 冊あります。 速い \pi の計算方法:「Pi と AGM」、ジョナサン ボーワインとピーター ボーワイン (アマゾンで入手可能).

私は AGM と関連アルゴリズムをかなり勉強しました。それは非常に興味深いものです(ただし、時には自明ではない場合もあります)。

\pi を計算する最新のアルゴリズムを実装するには、多精度算術ライブラリ (GMP 最後に使用してからかなり時間が経ちましたが、これは非常に良い選択です)。

最良のアルゴリズムの時間計算量は O(M(n)log(n)) です。ここで、M(n) は 2 つの n ビット整数の乗算の時間計算量です (M(n)=O(n) log(n) log(log(n))) FFT ベースのアルゴリズムを使用します。通常、\pi の桁を計算するときに必要であり、そのようなアルゴリズムは GMP で実装されています)。

アルゴリズムの背後にある数学は簡単ではないかもしれませんが、アルゴリズム自体は通常数行の疑似コードであり、その実装は通常非常に簡単であることに注意してください (独自の多精度算術演算を作成しないことを選択した場合:-) 。

以下の回答 正確には、最小限のコンピューティング労力で、可能な限り最速の方法でこれを行う方法. 。たとえその答えが気に入らないとしても、それが確かに PI の値を取得する最も早い方法であることを認めなければなりません。

最速 Pi の値を取得する方法は次のとおりです。

1)お気に入りのプログラミング言語を選択しました2)数学ライブラリ3)そして、Piがすでに定義されていることを見つけます - 使用の準備ができています!

手元に数学ライブラリがない場合に備えて..

2 番目に速い 方法(より普遍的な解決策)は次のとおりです。

インターネットで Pi を調べます。例:ここ:

http://www.eveandersson.com/pi/digits/1000000 (100万桁..浮動小数点の精度はどれくらいですか?)

またはここ:

http://3.141592653589793238462643383279502884197169399375105820974944592.com/

またはここ:

http://en.wikipedia.org/wiki/Pi

使用したい精度の演算に必要な桁を見つけるのは非常に速く、定数を定義することで貴重な CPU 時間を無駄にしないようにできます。

これは部分的にユーモアのある答えであるだけでなく、実際に、もし誰かが実際のアプリケーションで Pi の値を計算してみてくれたら...それはかなりの CPU 時間の無駄になりますね。少なくとも、これを再計算するための実際のアプリケーションは見当たりません。

モデレーター様:OPが次のように尋ねたことに注意してください。「PI値を取得する最速の方法」

BBP式 最初に前の n-1 桁を気にすることなく、n 番目の桁を基数 2 (または 16) で計算できます:)

pi を定数として定義する代わりに、私は常に次のように使用します。 acos(-1).

完全を期すためにここにあるはずのこれを見つけました。

Piet で PI を計算する

プログラムを大きくするほど精度が向上するという、かなり優れた特性を持っています。

ここ言語自体についての洞察

もし この記事 が true の場合、 ベラードのアルゴリズム が作成したものは、利用可能なものの中で最も高速なものの 1 つになる可能性があります。彼はデスクトップ PC を使用して 2.7 兆桁の円周率を作成しました。

...そして彼は彼の作品を出版しました ここで働きます

頑張れベラード、あなたはパイオニアです!

http://www.theregister.co.uk/2010/01/06/very_long_pi/

これは「古典的な」方法であり、実装が非常に簡単です。この実装は Python (それほど高速ではない言語) で行われます。

from math import pi
from time import time


precision = 10**6 # higher value -> higher precision
                  # lower  value -> higher speed

t = time()

calc = 0
for k in xrange(0, precision):
    calc += ((-1)**k) / (2*k+1.)
calc *= 4. # this is just a little optimization

t = time()-t

print "Calculated: %.40f" % calc
print "Costant pi: %.40f" % pi
print "Difference: %.40f" % abs(calc-pi)
print "Time elapsed: %s" % repr(t)

さらに詳しい情報を見つけることができます ここ.

とにかく、Python で正確な pi 値を必要なだけ取得する最速の方法は次のとおりです。

from gmpy import pi
print pi(3000) # the rule is the same as 
               # the precision on the previous code

これは gmpy pi メソッドのソースの一部ですが、この場合のコードはコメントほど役に立たないと思います。

static char doc_pi[]="\
pi(n): returns pi with n bits of precision in an mpf object\n\
";

/* This function was originally from netlib, package bmp, by
 * Richard P. Brent. Paulo Cesar Pereira de Andrade converted
 * it to C and used it in his LISP interpreter.
 *
 * Original comments:
 * 
 *   sets mp pi = 3.14159... to the available precision.
 *   uses the gauss-legendre algorithm.
 *   this method requires time o(ln(t)m(t)), so it is slower
 *   than mppi if m(t) = o(t**2), but would be faster for
 *   large t if a faster multiplication algorithm were used
 *   (see comments in mpmul).
 *   for a description of the method, see - multiple-precision
 *   zero-finding and the complexity of elementary function
 *   evaluation (by r. p. brent), in analytic computational
 *   complexity (edited by j. f. traub), academic press, 1976, 151-176.
 *   rounding options not implemented, no guard digits used.
*/
static PyObject *
Pygmpy_pi(PyObject *self, PyObject *args)
{
    PympfObject *pi;
    int precision;
    mpf_t r_i2, r_i3, r_i4;
    mpf_t ix;

    ONE_ARG("pi", "i", &precision);
    if(!(pi = Pympf_new(precision))) {
        return NULL;
    }

    mpf_set_si(pi->f, 1);

    mpf_init(ix);
    mpf_set_ui(ix, 1);

    mpf_init2(r_i2, precision);

    mpf_init2(r_i3, precision);
    mpf_set_d(r_i3, 0.25);

    mpf_init2(r_i4, precision);
    mpf_set_d(r_i4, 0.5);
    mpf_sqrt(r_i4, r_i4);

    for (;;) {
        mpf_set(r_i2, pi->f);
        mpf_add(pi->f, pi->f, r_i4);
        mpf_div_ui(pi->f, pi->f, 2);
        mpf_mul(r_i4, r_i2, r_i4);
        mpf_sub(r_i2, pi->f, r_i2);
        mpf_mul(r_i2, r_i2, r_i2);
        mpf_mul(r_i2, r_i2, ix);
        mpf_sub(r_i3, r_i3, r_i2);
        mpf_sqrt(r_i4, r_i4);
        mpf_mul_ui(ix, ix, 2);
        /* Check for convergence */
        if (!(mpf_cmp_si(r_i2, 0) && 
              mpf_get_prec(r_i2) >= (unsigned)precision)) {
            mpf_mul(pi->f, pi->f, r_i4);
            mpf_div(pi->f, pi->f, r_i3);
            break;
        }
    }

    mpf_clear(ix);
    mpf_clear(r_i2);
    mpf_clear(r_i3);
    mpf_clear(r_i4);

    return (PyObject*)pi;
}

編集: カットアンドペーストと識別に問題がありましたが、とにかくソースを見つけることができます ここ.

最速とは、コードを入力するのが最も速いことを意味する場合は、次のとおりです。 ゴルフスクリプト 解決:

;''6666,-2%{2+.2/@*\/10.3??2*+}*`1000<~\;

Machine のような公式を使用する

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))))

近似値を使用してもよい場合は、 355 / 113 は 6 桁の 10 進数に適しており、整数式にも使用できるという利点もあります。「浮動小数点演算コプロセッサ」が意味を持たなくなったため、最近ではそれほど重要ではなくなりましたが、かつては非常に重要でした。

ダブルスの場合:

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

これは、小数点以下 14 桁まで正確で、double を埋めるのに十分です (不正確なのは、おそらく逆正接の小数点以下が切り捨てられるためです)。

あとセス、3.14159265358979323846です3, 、64ではありません。

コンパイル時に D を使用して PI を計算します。

( からコピー 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))) );

円周率はちょうど 3 です![教授]フリンク(シンプソンズ)】

冗談ですが、これは C# でのものです (.NET-Framework が必要です)。

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;
        }
    }
}

このバージョン (Delphi 版) は特別なものではありませんが、少なくとも以下よりも高速です。 Nick Hodge がブログに投稿したバージョン :)。私のマシンでは、10 億回の反復を実行するのに約 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.

昔は、ワード サイズが小さく、浮動小数点演算が遅かったり存在しなかったりして、次のようなことを行っていました。

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

あまり精度を必要としないアプリケーション (ビデオ ゲームなど) の場合、これは非常に高速であり、十分な精度です。

あなたがしたい場合は 計算する (何らかの理由で) π の値の近似値を取得するには、バイナリ抽出アルゴリズムを試してください。 ベラードさん の改善 BBP O(N^2) で PI を与えます。


あなたがしたい場合は 得る 計算を行うための π の値の近似値を使用すると、次のようになります。

PI = 3.141592654

確かに、これは単なる近似値であり、完全に正確ではありません。0.00000000004102 より少しだけずれています。(10 兆分の 4、約 4/10,000,000,000).


やりたいなら 数学 π を使用する場合は、紙と鉛筆、またはコンピューター代数パッケージを用意して、π の正確な値 π を使用します。

本当に数式が必要な場合は、これが楽しいです。

π = - ln(-1)

Chris によって投稿された Brent の方法は非常に優れています。ブレントは一般に、任意精度演算の分野では巨人です。

N 番目の数字だけが必要な場合は、有名なBBP式16進数で役に立ちます

円の面積から π を計算する :-)

<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 += "&nbsp;" 
        }
        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>

より良いアプローチ

次のような標準定数の出力を取得するには 円周率 または標準的な概念を使用する場合は、まず、使用している言語で利用可能な組み込みメソッドを使用する必要があります。最速かつ最善の方法で値を返します。Pythonを使用して値piを取得する最速の方法を取得しています

  • 数学ライブラリの pi 変数. 。数学ライブラリは変数 pi を定数として保存します。

math_pi.py

import math
print math.pi

Linuxのtimeユーティリティでスクリプトを実行します。 /usr/bin/time -v python math_pi.py

出力:

Command being timed: "python math_pi.py"
User time (seconds): 0.01
System time (seconds): 0.01
Percent of CPU this job got: 91%
Elapsed (wall clock) time (h:mm:ss or m:ss): 0:00.03
  • arc cos 法による計算を使用する

acos_pi.py

import math
print math.acos(-1)

Linuxのtimeユーティリティでスクリプトを実行します。 /usr/bin/time -v python acos_pi.py

出力:

Command being timed: "python acos_pi.py"
User time (seconds): 0.02
System time (seconds): 0.01
Percent of CPU this job got: 94%
Elapsed (wall clock) time (h:mm:ss or m:ss): 0:00.03

bbp_pi.py

from decimal import Decimal, getcontext
getcontext().prec=100
print sum(1/Decimal(16)**k * 
          (Decimal(4)/(8*k+1) - 
           Decimal(2)/(8*k+4) - 
           Decimal(1)/(8*k+5) -
           Decimal(1)/(8*k+6)) for k in range(100))

Linuxのtimeユーティリティでスクリプトを実行します。 /usr/bin/time -v python bbp_pi.py

出力:

Command being timed: "python c.py"
User time (seconds): 0.05
System time (seconds): 0.01
Percent of CPU this job got: 98%
Elapsed (wall clock) time (h:mm:ss or m:ss): 0:00.06

したがって、最良の方法は、言語によって提供される組み込みメソッドを使用することです。これは、出力を取得するのに最も速く最適であるためです。Pythonではmath.piを使用します

ライセンス: CC-BY-SA帰属
所属していません StackOverflow
scroll top