Wie kann man den Wert von π am schnellsten ermitteln?

Question

Ich suche nach dem schnellsten Weg, den Wert von zu erhalten, als persönliche Herausforderung. Genauer gesagt, suche ich nach Wegen, die nicht die Verwendung von #define Konstanten wie M_PI oder die Nummer fest eintippen.

Das folgende Programm testet die verschiedenen mir bekannten Möglichkeiten. Die Inline-Assembly-Version ist theoretisch die schnellste Option, obwohl sie natürlich nicht portabel ist. Ich habe sie als Basis für den Vergleich mit den anderen Versionen aufgenommen. In meinen Tests, mit eingebauten Komponenten, war die 4 * atan(1) Version ist auf GCC 4.2 am schnellsten, weil sie die automatische Faltung der atan(1) in eine Konstante. Mit -fno-builtin angegeben, die atan2(0, -1) Version ist am schnellsten.

Hier ist das Haupttestprogramm ( 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;
}

Und die Inline-Assembly-Sachen ( fldpi.c ), die nur für x86- und x64-Systeme funktioniert:

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

Und ein Build-Skript, das alle Konfigurationen erstellt, die ich teste ( 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

Abgesehen vom Testen zwischen verschiedenen Compiler-Flags (ich habe auch 32-Bit mit 64-Bit verglichen, weil die Optimierungen unterschiedlich sind), habe ich auch versucht, die Reihenfolge der Tests umzukehren. Aber immer noch ist die atan2(0, -1) Version immer noch jedes Mal die Nase vorn hat.

Answer 1

Les Monte-Carlo-Verfahren wendet, wie erwähnt, einige großartige Konzepte an, aber es ist eindeutig nicht das schnellste, bei weitem nicht, bei keinem vernünftigen Maß. Außerdem hängt alles davon ab, welche Art von Genauigkeit Sie suchen. Die schnellste mir bekannte Lösung ist die mit den fest kodierten Ziffern. Blick auf Pi y Pi[PDF] gibt es eine Vielzahl von Formeln.

Hier ist eine Methode, die schnell konvergiert - etwa 14 Ziffern pro Iteration. PiFast , die derzeit schnellste Anwendung, verwendet diese Formel mit der FFT. Ich werde nur die Formel schreiben, da der Code einfach ist. Diese Formel wurde fast gefunden von Ramanujan und entdeckt von Chudnovsky . Auf diese Weise hat er mehrere Milliarden Ziffern der Zahl berechnet - es ist also keine Methode, die man außer Acht lassen sollte. Die Formel wird schnell überlaufen, und da wir durch Faktoren dividieren, wäre es von Vorteil, solche Berechnungen zu verzögern, um Terme zu entfernen.

wo,

Nachstehend finden Sie die Brent-Salamin-Algorithmus . Wikipedia erwähnt, dass, wenn a y b "nahe genug" sind, dann (a + b)² / 4t Ich bin mir nicht sicher, was "nahe genug" bedeutet, aber meine Tests ergaben, dass eine Iteration 2 Ziffern, zwei 7 Ziffern und drei 15 Ziffern ergaben, natürlich mit Doubles, so dass es aufgrund der Darstellung und der wahr Die Berechnung könnte genauer sein.

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)

Wie wäre es schließlich mit etwas Pi-Golf (800 Ziffern)? 160 Zeichen!

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

Answer 2

Dieses Programm gefällt mir sehr gut, weil es eine Annäherung durch Betrachtung des eigenen Bereichs vornimmt.

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()
{
            _-_-_-_
       _-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
        _-_-_-_-_-_-_-_
            _-_-_-_
}

Answer 3

Hier ist eine allgemeine Beschreibung einer Technik zur Berechnung von Pi, die ich in der High School gelernt habe.

Ich erzähle dies nur, weil ich denke, dass es einfach genug ist, dass sich jeder daran erinnern kann, und es lehrt Sie das Konzept der "Monte-Carlo"-Methoden - das sind statistische Methoden, um zu Antworten zu gelangen, die nicht sofort durch Zufallsprozesse ableitbar zu sein scheinen.

Zeichne ein Quadrat und beschrifte einen Quadranten (ein Viertel eines Halbkreises) innerhalb dieses Quadrats (ein Quadrant mit einem Radius, der der Seite des Quadrats entspricht, so dass er das Quadrat so weit wie möglich ausfüllt)

Werfen Sie nun einen Pfeil auf das Quadrat und notieren Sie, wo er landet, d. h. wählen Sie einen beliebigen Punkt innerhalb des Quadrats. Natürlich ist er innerhalb des Quadrats gelandet, aber ist er auch innerhalb des Halbkreises? Notiere diese Tatsache.

Wiederholen Sie diesen Vorgang viele Male - und Sie werden feststellen, dass es ein Verhältnis zwischen der Anzahl der Punkte innerhalb des Halbkreises und der Gesamtzahl der geworfenen Punkte gibt, nennen Sie dieses Verhältnis x.

Da der Flächeninhalt des Quadrats r mal r ist, kann man ableiten, dass der Flächeninhalt des Halbkreises x mal r mal r ist (also x mal r zum Quadrat). Somit ergibt x mal 4 pi.

Diese Methode lässt sich nicht schnell anwenden. Aber sie ist ein schönes Beispiel für eine Monte-Carlo-Methode. Und wenn Sie sich umsehen, werden Sie feststellen, dass viele Probleme, die sonst außerhalb Ihrer rechnerischen Fähigkeiten liegen, mit solchen Methoden gelöst werden können.

Answer 4

Der Vollständigkeit halber eine C++-Vorlage, die bei einem optimierten Build zur Kompilierzeit eine Annäherung an PI berechnet und zu einem einzigen Wert inlinet.

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

Bei I > 10 können optimierte Builds langsam sein, ebenso bei nicht optimierten Läufen. Bei 12 Iterationen gibt es meiner Meinung nach etwa 80k Aufrufe von value() (ohne Memoisierung).

Answer 5

Es gibt sogar ein ganzes Buch, das (unter anderem) folgenden Themen gewidmet ist schnell Methoden zur Berechnung von \pi Pi und die Hauptversammlung", von Jonathan und Peter Borwein ( erhältlich bei Amazon ).

Ich habe mich eingehend mit dem AGM und verwandten Algorithmen beschäftigt: Es ist sehr interessant (wenn auch manchmal nicht trivial).

Beachten Sie, dass die meisten modernen Algorithmen zur Berechnung von \pi benötigen Sie eine Bibliothek für Multipräzisionsarithmetik ( GMP ist eine gute Wahl, auch wenn es schon eine Weile her ist, dass ich es das letzte Mal benutzt habe).

Die Zeitkomplexität der besten Algorithmen liegt bei O(M(n)log(n)), wobei M(n) die Zeitkomplexität für die Multiplikation von zwei n-Bit-Ganzzahlen (M(n)=O(n log(n) log(log(n))) unter Verwendung von FFT-basierten Algorithmen ist, die normalerweise für die Berechnung von Ziffern von \pi und ein solcher Algorithmus ist in GMP implementiert).

Auch wenn die Mathematik hinter den Algorithmen nicht trivial ist, bestehen die Algorithmen selbst in der Regel aus ein paar Zeilen Pseudocode, und ihre Implementierung ist in der Regel sehr einfach (wenn Sie nicht Ihre eigene Multipräzisionsarithmetik schreiben wollen :-) ).