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

Der Chudnovsky-Algorithmus ist ziemlich schnell, wenn es Ihnen nichts ausmacht, eine Quadratwurzel und ein paar Inversionen durchzuführen. Er konvergiert bei doppelter Genauigkeit in nur 2 Iterationen.

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

Ergebnisse:

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

Answer 2

Ich denke, der Wert von Pi ist das Verhältnis zwischen dem Umfang und dem Radius des Kreises.

Es kann einfach durch eine normale mathematische Berechnung erreicht werden

Answer 3

Besserer Ansatz

Um die Ausgabe von Standardkonstanten wie pi oder die Standardkonzepte zu verwenden, sollten wir zunächst die in der von Ihnen verwendeten Sprache verfügbaren eingebauten Methoden verwenden. Sie geben einen Wert auf die schnellste und beste Weise zurück. Ich verwende Python, um den schnellsten Weg zu finden, den Wert von Pi zu ermitteln.

• pi-Variable der mathematischen Bibliothek . Die mathematische Bibliothek speichert die Variable pi als Konstante.

math_pi.py

import math
print math.pi

Führen Sie das Skript mit dem Zeitdienstprogramm von Linux aus /usr/bin/time -v python math_pi.py

Sortie :

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

• Verwendung der Arc-Cos-Methode in der Mathematik

acos_pi.py

import math
print math.acos(-1)

Führen Sie das Skript mit dem Zeitdienstprogramm von Linux aus /usr/bin/time -v python acos_pi.py

Sortie :

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

• utiliser BBP-Formel

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

Führen Sie das Skript mit dem Zeitdienstprogramm von Linux aus /usr/bin/time -v python bbp_pi.py

Sortie :

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

Daher ist es am besten, die von der Sprache bereitgestellten eingebauten Methoden zu verwenden, da sie die schnellste und beste Möglichkeit sind, die Ausgabe zu erhalten. In Python verwenden Sie math.pi