Jaki jest najszybszy sposób na uzyskanie wartości π?

Szukam najszybszego sposobu na uzyskanie wartości π jako osobistego wyzwania. Mówiąc dokładniej, korzystam ze sposobów, które nie wymagają użycia #definestałych, takich jak M_PIlub zakodowania na stałe liczby.

Poniższy program testuje różne znane mi sposoby. Wersja zestawu wbudowanego jest teoretycznie najszybszą opcją, choć oczywiście nie jest przenośna. Podałem go jako punkt odniesienia do porównania z innymi wersjami. W moich testach, z wbudowanymi 4 * atan(1)wersjami , wersja jest najszybsza na GCC 4.2, ponieważ automatycznie składa się atan(1)w stałą. Z -fno-builtinokreślony, atan2(0, -1)wersja jest najszybszy.

Oto główny program testujący ( 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;
}

A wbudowane elementy asemblera ( fldpi.c), które będą działać tylko w systemach x86 i x64:

double
fldpi()
{
    double pi;
    asm("fldpi" : "=t" (pi));
    return pi;
}

I skrypt kompilacji, który buduje wszystkie konfiguracje, które testuję ( 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

Oprócz testowania różnych flag kompilatora (porównałem również 32-bitowy z 64-bitowym, ponieważ różne są optymalizacje), próbowałem również zmienić kolejność testów. Mimo to atan2(0, -1)wersja wciąż wychodzi na wierzch za każdym razem.

Jak wspomniano, metoda Monte Carlo stosuje kilka świetnych koncepcji, ale najwyraźniej nie jest najszybsza, nie jest daleka, nie jest w żaden rozsądny sposób. Wszystko zależy również od tego, jakiej dokładności szukasz. Najszybszy, jaki znam, to ten z cyframi zakodowanymi na stałe. Patrząc na Pi i Pi [PDF] , istnieje wiele formuł.

Oto metoda, która szybko się zbiega - około 14 cyfr na iterację. PiFast , obecnie najszybsza aplikacja, używa tej formuły z FFT. Napiszę formułę, ponieważ kod jest prosty. Ta formuła została prawie znaleziona przez Ramanujana i odkryta przez Chudnovsky'ego . Tak właśnie wyliczył kilka miliardów cyfr tej liczby - więc nie można tego lekceważyć. Formuła przepełni się szybko, a ponieważ dzielimy czynniki, korzystne byłoby opóźnienie takich obliczeń w celu usunięcia warunków.

gdzie,

Poniżej znajduje się algorytm Brenta-Salamina . Wikipedia wskazuje, że gdy i b są „wystarczająco blisko”, wówczas (a + b) ² / 4T będzie przybliżeniem Õ. Nie jestem pewien, co oznacza „wystarczająco blisko”, ale z moich testów jedna iteracja otrzymała 2 cyfry, dwie dostały 7, a trzy miały 15, oczywiście jest to z podwójnymi, więc może mieć błąd na podstawie jego reprezentacji i prawdziwe obliczenie może być bardziej dokładne.

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)

Wreszcie, co powiesz na golfa pi (800 cyfr)? 160 znaków!

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

Naprawdę podoba mi się ten program, ponieważ przybliża on π patrząc na swój własny obszar.

IOCCC 1988: westley.c

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

Oto ogólny opis techniki obliczania liczby pi, której nauczyłem się w szkole średniej.

Dzielę się tym, ponieważ uważam, że jest to wystarczająco proste, aby każdy mógł je zapamiętać w nieskończoność, a ponadto uczy koncepcji metod „Monte-Carlo” - statystycznych metod uzyskiwania odpowiedzi, które nie wydają się od razu można wydedukować przez losowe procesy.

Narysuj kwadrat i wpisz kwadrant (ćwierć półkola) wewnątrz tego kwadratu (kwadrant o promieniu równym bokowi kwadratu, aby wypełniał jak najwięcej kwadratu)

Teraz rzuć strzałką w kwadrat i zapisz, gdzie wyląduje - to znaczy wybierz losowy punkt w dowolnym miejscu kwadratu. Oczywiście wylądował w kwadracie, ale czy jest w półkolu? Zapisz ten fakt.

Powtórz ten proces wiele razy - przekonasz się, że istnieje stosunek liczby punktów w półokręgu do całkowitej liczby rzuconych, nazwij ten stosunek x.

Ponieważ powierzchnia kwadratu wynosi r razy r, można wywnioskować, że powierzchnia półokręgu wynosi x razy r razy r (to znaczy x razy r do kwadratu). Stąd x razy 4 da ci pi.

Nie jest to szybka metoda. Ale to dobry przykład metody Monte Carlo. A jeśli się rozejrzysz, możesz odkryć, że wiele problemów poza umiejętnościami obliczeniowymi można rozwiązać za pomocą takich metod.

W trosce o kompletność, wersja szablonu C ++, która dla zoptymalizowanej kompilacji oblicza aproksymację PI w czasie kompilacji i wstawia do pojedynczej wartości.

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

Uwaga dla I> 10, zoptymalizowane kompilacje mogą być powolne, podobnie w przypadku niezoptymalizowanych przebiegów. Wydaje mi się, że dla 12 iteracji istnieje około 80 000 wywołań wartości () (przy braku zapamiętywania).

W rzeczywistości istnieje cała książka poświęcona (między innymi) szybkim metodom obliczania \ pi: „Pi and the AGM” autorstwa Jonathana i Petera Borweina ( dostępna na Amazon ).

Studiowałem AGM i pokrewne algorytmy: jest to dość interesujące (choć czasem nietrywialne).

Pamiętaj, że aby zaimplementować większość nowoczesnych algorytmów do obliczania \ pi, potrzebujesz biblioteki arytmetycznej z wieloma dokładnościami ( GMP jest całkiem dobrym wyborem, chociaż minęło trochę czasu, odkąd go ostatnio używałem).

Złożoność czasowa najlepszych algorytmów wyrażona jest w O (M (n) log (n)), gdzie M (n) jest złożonością czasową mnożenia dwóch liczb całkowitych n-bitowych (M (n) = O (n log (n) log (log (n))) przy użyciu algorytmów opartych na FFT, które są zwykle potrzebne przy obliczaniu cyfr \ pi, a taki algorytm jest implementowany w GMP).

Zauważ, że chociaż matematyka stojąca za algorytmami może nie być trywialna, same algorytmy są zwykle kilkoma liniami pseudokodu, a ich implementacja jest zwykle bardzo prosta (jeśli nie zdecydujesz się pisać własnej arytmetyki wielrecyzyjnej :-)).

Poniższe odpowiedzi dokładnie wyjaśniają, jak to zrobić w najszybszy możliwy sposób - przy minimalnym wysiłku obliczeniowym . Nawet jeśli nie podoba ci się odpowiedź, musisz przyznać, że jest to rzeczywiście najszybszy sposób na uzyskanie wartości PI.

NAJSZYBSZY sposób, aby uzyskać wartość PI:

1) wybierz swój ulubiony język programowania 2) załaduj bibliotekę Math 3) i przekonaj się, że Pi jest już tam zdefiniowany - gotowy do użycia!

Jeśli nie masz pod ręką biblioteki Math ...

The Drugim najszybciej sposób (bardziej uniwersalne rozwiązanie) to:

wyszukaj Pi w Internecie, np. tutaj:

http://www.eveandersson.com/pi/digits/1000000 (1 milion cyfr .. jaka jest twoja precyzja zmiennoprzecinkowa?)

lub tu:

http://3.141592653589793238462643383279502884197169399375105820974944592.com/

lub tu:

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

Znalezienie potrzebnych cyfr dla dowolnej precyzyjnej arytmetyki jest bardzo szybkie, a definiując stałą, możesz mieć pewność, że nie zmarnujesz cennego czasu procesora.

To nie tylko częściowo humorystyczna odpowiedź, ale w rzeczywistości, jeśli ktoś poszedłby naprzód i obliczyłby wartość Pi w prawdziwej aplikacji ... byłoby to dość marnowaniem czasu procesora, prawda? Przynajmniej nie widzę prawdziwej aplikacji do próby ponownego obliczenia tego.

Drogi Moderatorze: pamiętaj, że PO zapytał: „Najszybszy sposób na uzyskanie wartości PI”

Wzór BBP pozwala obliczyć n-cyfry - w podstawie (2 lub 16), - bez konieczności nawet przejmować z wcześniejszych n-1 cyframi pierwszy :)

Zamiast definiować pi jako stałą, zawsze używam acos(-1).

Jest to „klasyczna” metoda, bardzo łatwa do wdrożenia. Ta implementacja w Pythonie (nie w najszybszym języku) robi to:

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 "Constant pi: %.40f" % pi
print "Difference: %.40f" % abs(calc-pi)
print "Time elapsed: %s" % repr(t)

Więcej informacji znajdziesz tutaj .

W każdym razie najszybszym sposobem uzyskania dokładnej wartości liczby pi w pythonie jest:

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

Oto źródło dla metody gmpy pi, nie sądzę, aby kod był tak przydatny jak komentarz w tym przypadku:

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

EDYCJA: Miałem pewne problemy z wycinaniem i wklejaniem oraz wcięciem, możesz znaleźć źródło tutaj .

Jeśli przez najszybszy masz na myśli najszybszy wpisanie kodu, oto rozwiązanie w golfscript :

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

Jeśli chcesz użyć przybliżenia, 355 / 113jest dobre dla 6 cyfr dziesiętnych i ma tę dodatkową zaletę, że można go używać z wyrażeniami całkowitymi. Nie jest to obecnie tak ważne, ponieważ „zmiennoprzecinkowy koprocesor matematyczny” przestał mieć jakiekolwiek znaczenie, ale kiedyś był dość ważny.

Użyj formuły podobnej do Machina

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.

Zaimplementowane w programie, na przykład:

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

Z podwójnymi:

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

Będzie to dokładne z dokładnością do 14 miejsc po przecinku, wystarczające do wypełnienia podwójnej (niedokładność wynika prawdopodobnie z tego, że reszta miejsc po przecinku w stycznych łuku jest obcięta).

Także Set, to 3.14159265358979323846 3 , a nie 64.

Pi ma dokładnie 3! [Prof. Frink (Simpsons)]

Żart, ale tutaj jest w C # (wymagany .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;
        }
    }
}

Oblicz PI w czasie kompilacji za pomocą D.

(Skopiowane z 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))) );

Ta wersja (w Delphi) nie jest niczym specjalnym, ale jest co najmniej szybsza niż wersja, którą Nick Hodge opublikował na swoim blogu :). Na mojej maszynie wykonanie miliarda iteracji zajmuje około 16 sekund, co daje wartość 3,14159265 25879 (dokładna część jest pogrubiona).

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.

W dawnych czasach, z małymi rozmiarami słów i powolnymi lub nieistniejącymi operacjami zmiennoprzecinkowymi, robiliśmy takie rzeczy:

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

W przypadku aplikacji, które nie wymagają dużej precyzji (na przykład gry wideo), jest to bardzo szybkie i wystarczająco dokładne.

Jeśli chcesz obliczyć przybliżoną wartość π (z jakiegoś powodu), powinieneś wypróbować binarny algorytm ekstrakcji. Bellard jest poprawa BBP daje nie PI w O (N ^ 2).

Jeśli chcesz uzyskać przybliżenie wartości π do wykonania obliczeń, to:

PI = 3.141592654

To prawda, to tylko przybliżenie i nie do końca dokładne. Jest wyłączony o nieco więcej niż 0,00000000004102. (cztery dziesięć-trillionths około ⁴ / _10000000000 ).

Jeśli chcesz robić matematykę z π, zdobądź ołówek i papier lub pakiet algebry komputerowej i użyj dokładnej wartości π, π.

Jeśli naprawdę potrzebujesz formuły, ta jest fajna:

π = - i ln (-1)

Metoda Brenta opublikowana powyżej przez Chrisa jest bardzo dobra; Brent ogólnie jest gigantem w dziedzinie arytmetyki o dowolnej precyzji.

Jeśli wszystko, czego chcesz, to N-ta cyfra, znana formuła BBP jest przydatna w trybie szesnastkowym

Obliczanie π na podstawie okręgu :-)

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

Algorytm Chudnovsky'ego jest dość szybki, jeśli nie masz nic przeciwko wykonaniu pierwiastka kwadratowego i kilku odwrotności. Zbiega się do podwójnej precyzji w zaledwie 2 iteracjach.

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

Wyniki:

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

Lepsze podejście

Aby uzyskać wynik standardowych stałych, takich jak pi lub standardowe pojęcia, powinniśmy najpierw przejść do metod wbudowanych dostępnych w języku, którego używasz. Zwróci wartość w najszybszy i najlepszy sposób. Używam Pythona, aby uruchomić najszybszy sposób na uzyskanie wartości pi.

• zmienna pi z biblioteki matematycznej . Biblioteka matematyczna przechowuje zmienną pi jako stałą.

math_pi.py

import math
print math.pi

Uruchom skrypt za pomocą narzędzia Linux-a /usr/bin/time -v python math_pi.py

Wynik:

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

• Użyj metody matematycznej arc cos

acos_pi.py

import math
print math.acos(-1)

Uruchom skrypt za pomocą narzędzia Linux-a /usr/bin/time -v python acos_pi.py

Wynik:

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

• użyj formuły BBP

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

Uruchom skrypt za pomocą narzędzia Linux-a /usr/bin/time -v python bbp_pi.py

Wynik:

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

Zatem najlepszym sposobem jest użycie wbudowanych metod dostarczanych przez język, ponieważ są one najszybsze i najlepsze, aby uzyskać wynik. W pythonie użyj math.pi