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

Diese Version (in Delphi) ist nichts Besonderes, aber sie ist zumindest schneller als die Version, die Nick Hodge in seinem Blog veröffentlicht hat :). Auf meinem Rechner dauert es etwa 16 Sekunden, eine Milliarde Iterationen durchzuführen, was einen Wert von 3.14159265 25879 (der richtige Teil ist fett gedruckt).

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.

Answer 2

Früher, als es noch kleine Wortgrößen und langsame oder nicht vorhandene Gleitkommaoperationen gab, haben wir so etwas gemacht:

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

Für Anwendungen, die keine große Präzision erfordern (z. B. Videospiele), ist dies sehr schnell und genau genug.

Answer 3

Wenn Sie möchten, dass berechnen Sie eine Annäherung an den Wert von (aus irgendeinem Grund), sollten Sie einen binären Extraktionsalgorithmus versuchen. Bellard's Verbesserung der BBP ergibt PI in O(N^2).

---

Wenn Sie möchten, dass erhalten eine Annäherung an den Wert von zu tun Berechnungen, dann:

PI = 3.141592654

Zugegeben, das ist nur ein Näherungswert und nicht ganz genau. Die Abweichung liegt bei etwas mehr als 0,00000000004102. (vier Zehntrillionstel, etwa 4 / _{10,000,000,000} ).

---

Wenn Sie Folgendes tun möchten Mathe mit , dann besorgen Sie sich einen Stift und Papier oder ein Computer-Algebra-Paket, und verwenden Sie den genauen Wert, .

Wenn Sie wirklich eine Formel wollen, dann ist diese ein Spaß:

\= - i ln(-1)

Answer 4

Im Grunde die C-Version von paperclip optimizer's Antwort, und viel einfacher:

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

Zur weiteren Vereinfachung verwendet dieser Algorithmus die Formel von Chudnovsky, die ich vollständig vereinfachen kann, wenn Sie den Code nicht wirklich verstehen.

Zusammenfassung: Wir nehmen eine Zahl von 1 bis 5 und fügen sie in eine Funktion ein, mit der wir PI erhalten. Dann werden dir 3 Zahlen gegeben: 545140134 (A), 13591409 (B), 640320 (D). Dann werden wir D als eine double sich 3 mal in ein anderes multiplizieren double (ID3). Wir werden dann die Quadratwurzel von ID3 in eine andere double (b) und vergeben 2 Nummern: 1 (p), den Wert von B (a). Beachten Sie, dass bei C die Groß- und Kleinschreibung keine Rolle spielt. Dann wird ein double (Summe) wird durch Multiplikation der Werte von p, a und b erstellt, die alle in double s. Dann beginnt eine Schleife bis zu der für die Funktion angegebenen Zahl und addiert den Wert von A zu a, der Wert von b wird mit ID3 multipliziert, der Wert von p wird mit mehreren Werten multipliziert, die Sie hoffentlich verstehen können, und wird auch durch mehrere Werte geteilt. Die Summe wird noch einmal mit p, a und b addiert und die Schleife wird so lange wiederholt, bis der Wert der Schleifennummer größer oder gleich 5 ist. Später wird die Summe mit 12 multipliziert und von der Funktion zurückgegeben, so dass wir das Ergebnis von PI erhalten.

Okay, das war lang, aber ich denke, Sie werden es verstehen...

Answer 5

Berechnen aus der Kreisfläche :-)

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