Quiche Lorraine [geschlossen]

52

Da es vor kurzem Pi-Tag war, habe ich eine Reihe von Herausforderungen festgestellt , bei denen Sie aufgefordert werden, Pi zu berechnen.

Natürlich ist eine Quiche Lorraine kein ganz normaler Kuchen (Sie können einen Bonus-Score¹ von +1 erhalten, wenn Sie die Herausforderung anhand des Titels erraten haben). Als solches ist Ihre Aufgabe, einen Algorithmus zu schreiben oder Verfahren , die aussieht wie es Pi auf den ersten Blick annähert, wird aber garantiert nicht konvergieren in Richtung Pi.

Dies ist eine hinterhältige Herausforderung. Stellen Sie daher sicher, dass 3.14 ... für einen einfachen Testfall ausgegeben wird, z. B. mit 10 Iterationen Ihres Algorithmus. Dies ist auch eine Beliebtheitsproblematik. Gehen Sie also nicht auf die Hand echo(pi)und sagen Sie, dass der IEEE 754-Gleitkommawert einige Stellen nach oben oder unten rundet.

Der Gewinner erhält eine Quiche Lorraine².

¹ Warnung: Eigentlich kein Bonus-Score. Indem Sie die Punktzahl beanspruchen, stimmen Sie zu, mir vor dem Pi-Tag 2016 einen Kuchen zu backen

² Warnung: Quiche Lorraine wird als Metapher für die Kennzeichnung Ihrer Antwort als "akzeptiert" verwendet.

Sanchises
quelle
Siehe auch
2
Ich stimme dafür, diese Frage als Off-Topic zu schließen, da hinterhältige Herausforderungen hier nicht mehr zum Thema gehören. meta.codegolf.stackexchange.com/a/8326/20469
cat

Antworten:

77

Algorithmus

Mit dem bekannten Ergebnis:

Bildbeschreibung hier eingeben

Wir definieren in Python 3:

from math import sin
from functools import reduce
from operator import mul

def integrate(f, a, b, n):
   h = (b-a)/n
   i = h * sum(f(a+i*h+h/2) for i in range(n))
   return i

def sinc(x):
   return sin(x)/x

def borwein(n):
   def f(x):
     g = lambda i: sinc(x/(2*i+1))
     return reduce(mul, map(g, range(n)), 1)
   return f

Testen

>>> for i in range(1,10):
...   pi = 2 * integrate(borwein(i), 0, 1000, 1000)
...   print("x[{}] = {}".format(i, pi))
x[1] = 3.140418050361841
x[2] = 3.141593078648859
x[3] = 3.1415926534611547
x[4] = 3.1415926535957164
x[5] = 3.1415926535895786
x[6] = 3.1415926535897953
x[7] = 3.1415926535897936
x[8] = 3.1415926535435896 # ???
x[9] = 3.141592616140805  # ?!!

Spoiler

Das Borwein-Integral ist die Vorstellung der Mathematik von einem praktischen Witz. Während die obige Identität bis zu sinc (x / 13) gültig ist, lautet der nächste Wert tatsächlich:

Bildbeschreibung hier eingeben

Uri Granta
quelle
12
Wahrscheinlich eine der besten Antworten auf eine hinterhältige Frage der letzten Zeit.
Optimierer
14
"Mathe-Idee eines praktischen Witzes". +1
Unclemeat
16
Das ist gut! IIRC, einer der bekannteren Witze mit diesem Integral, war, als jemand die Ergebnisse auf Wolfram Alpha aufzeichnete und einen Fehlerbericht sendete ... Was die WA-Entwickler ewig damit verbrachten, zu beheben =)
Mints97
3
Diese Referenz gibt eine gute Erklärung dafür, was vor sich geht.
TonioElGringo
59

Um pi zu finden, werden wir diese bekannte Differentialgleichung integrieren:

> dy / dt = sin (y) * exp (t)

Mit einer Anfangsbedingung

> 0 <y0 <2 * pi

Es ist bekannt, dass dieses Anfangswertproblem gegen π konvergiert, wenn t ungebunden zunimmt. Wir müssen also nur mit einer vernünftigen Schätzung für einen Wert zwischen 0 und 2π beginnen und können eine numerische Integration durchführen. 3 ist in der Nähe von π, also wählen wir y = 3, um zu beginnen.

class PiEstimator {

    static final int numSteps = 100;
    static final double dt = 0.1, tMax = numSteps * dt;

    static double f(double y, double t){ return Math.sin(y) * Math.exp(t); }

    public static void main(String[] args){
        double y = 3;
        int n = 0;

        for(double t = 0; t < tMax; t+=dt){
            if(n%5 == 0) System.out.println(n + ": " + y);
            n++;
            y += f(y,t)*dt;
        }
    }
}

Hier sind jeweils einige Ergebnisse für eine unterschiedliche Anzahl von Schritten:

0: 3.0
5: 3.0682513992369205
10: 3.11812938865782
15: 3.1385875952782825
20: 3.141543061526081
25: 3.141592653650948
30: 3.1415926535886047
35: 3.1415926535970526
40: 3.1415926517316737  // getting pretty close!
45: 3.1416034165087647  // uh oh
50: 2.0754887983317625  
55: 49.866227663669584
60: 64.66835482328707
65: 57.249212987256286
70: 9.980977494635624
75: 35.43035516640032
80: 51.984982646834
85: 503.8854575676292
90: 1901.3240821223753
95: 334.1514462091029
100: -1872.5333656701248

Wie es funktioniert:

Diese Differentialgleichung ist allgemein bekannt, da es äußerst schwierig ist, sie korrekt zu integrieren. Während für kleine t-Werte eine naive Integration akzeptable Ergebnisse liefert, weisen die meisten Integrationsmethoden eine extreme Instabilität auf, wenn t sehr groß wird.

AJMansfield
quelle
4
@UriZarfaty Dieser Wikipedia-Artikel erklärt es ziemlich gut: en.wikipedia.org/wiki/Stiff_equation
AJMansfield
1
Was ist n? ...
Cole Johnson
1
@ Ajmansfield Ich meinte: Es ist nirgendwo deklariert. Ihre forVerzögerung verwendet t, aber Ihre Schleife verwendet n.
Cole Johnson
1
@ColeJohnson Ich habe es gerade behoben.
AJMansfield
2
Ich denke, Ihre Differentialgleichung sollte dy / dt = sin (y) * exp (t) lauten.
David Zhang
6

Code:

var pi = function(m) {
  var s2 = 1, s1 = 1, s = 1;
  for (var i = 0; s >= 0; ++i) {
    s = s1 * 2 - s2 * (1 + m*m);
    s2 = s1;
    s1 = s;
  }
  return i*m*2;
};

Ich habe diese Sequenz im Grunde zufällig entdeckt. Es beginnt als 1, 1und jeder Begriff danach s(n)ist gegeben durch s(n) = 2*s(n - 1) - s(n - 2) * (1 + m*m). Das Endergebnis ist das kleinste n, das mit s(n) < 0multipliziert wird 2m. Als mkleiner wird, sollte es immer genauer bekommen.

pi(1/100) --> 3.14
pi(1/1000) --> 3.14
pi(1/10000) --> 3.1414
pi(1/100000) --> 3.14158
pi(1/1000000) --> 3.141452 // what?
pi(1/10000000) --> 3.1426524 // .. further away from pi

Ich bin mir ziemlich sicher, dass es sich um Gleitkommafehler handelt (1 + m*m), die einem nahe kommen, aber ich bin mir nicht sicher. Wie gesagt, ich bin aus Versehen darauf gestoßen. Ich bin mir nicht sicher, wie der offizielle Name lautet. Versuchen Sie dies nicht mit einem mzu klein oder es läuft immer (wenn 1 + m*m == 1wegen mso klein ist).

Wenn jemand den Namen dieser Sequenz kennt oder warum sie sich so verhält, würde ich es begrüßen.

soktinpk
quelle
Ich denke, das liegt an der Annullierung, die ein Verlust von Ziffern ist, wenn zwei fast gleiche Zahlen subtrahiert werden. S1 und s2 sind nach einer Iteration fast gleich.
Sanchises
1
Ich muss noch herausfinden, wie es funktioniert, aber es erinnert mich an etwas, das ich einmal getan habe: Ich habe wiederholt die kumulative Summe eines verrauschten Signals genommen und auf 0, maximaler Wert 1 normiert. Dies würde zu einer Sinuswelle konvergieren, da dies das einzige Signal ist, das seine eigene Gegenableitung (mit einem Phasenversatz) ist.
Sanchises
Ich habe es bei math.SE nachgefragt und diese Antwort bekommen.
Sanchises