Quiche Lorraine [nieczynne]

Ponieważ było Pi dni niedawno, ja zauważyłem jak wiele z wyzwań , które proszą Cię do obliczenia PI.

Oczywiście, quiche Lorraine nie jest całkiem ciastem (możesz zdobyć wynik bonusowy¹ +1, jeśli odgadłeś wyzwanie z tytułu). W związku z tym Twoim zadaniem jest napisanie algorytmu lub metody, która na pierwszy rzut oka będzie wyglądała tak, jakby była zbliżona do Pi, ale z pewnością nie zbliży się do Pi.

Jest to trudne zadanie, więc upewnij się, że wyświetli 3.14 ... dla prostego przypadku testowego, np. Z 10 iteracjami twojego algorytmu. To także wyzwanie popularności, więc nie idź za oczywistością echo(pi)i powiedz, że zmiennoprzecinkowe IEEE 754 zaokrągla niektóre cyfry w górę lub w dół.

Zwycięzca dostaje quiche Lorraine².

¹ Ostrzeżenie: tak naprawdę nie jest to wynik bonusowy. Podejmując decyzję, zgadzasz się upiec mi ciasto przed Świętem Pi 2016

² Ostrzeżenie: quiche Lorraine służy jako metafora oznaczenia odpowiedzi jako „zaakceptowana”

Algorytm

Korzystając ze znanego wyniku:

definiujemy w Pythonie 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

Testowanie

>>> 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  # ?!!

Spojler

Borwein integralną jest idea Math praktycznego żartu. Chociaż powyższa tożsamość utrzymuje się na poziomie sinc (x / 13), następna wartość to w rzeczywistości:

Aby znaleźć pi, zintegrujemy to dobrze znane równanie różniczkowe:

Z warunkiem początkowym

Dobrze wiadomo, że ten problem z wartością początkową zbiega się do π, gdy t wzrasta bez ograniczeń. Tak więc wszystko, czego potrzebujemy, to zacząć od rozsądnego odgadnięcia czegoś między 0 a 2π, i możemy przeprowadzić integrację numeryczną. 3 jest bliskie π, więc wybierzemy y = 3, aby rozpocząć.

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

Oto kilka wyników dla każdej liczby kroków:

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

Jak to działa:

To równanie różniczkowe jest dobrze znane, ponieważ jego integracja jest niezwykle trudna. Podczas gdy dla małych wartości t integracja naiwna da akceptowalne wyniki, większość metod integracji wykazuje ekstremalną niestabilność, gdy t staje się bardzo duża.

Kod:

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

Zasadniczo odkryłem tę sekwencję przez przypadek. Zaczyna się jak 1, 1i każdy kolejny termin s(n)podany przez s(n) = 2*s(n - 1) - s(n - 2) * (1 + m*m). Wynik końcowy jest najmniejszy ni s(n) < 0pomnożony przez 2m. W miarę mzmniejszania się, powinno być coraz bardziej dokładne.

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

Jestem prawie pewien, że są to błędy zmiennoprzecinkowe (1 + m*m)zbliżające się do jednego, ale nie jestem pewien. Tak jak powiedziałem, natknąłem się na to przez przypadek. Nie jestem pewien jego oficjalnej nazwy. Nie próbuj tego ze związkiem mzbyt małe lub będzie działać zawsze (jeśli 1 + m*m == 1ze względu na mbycie tak małe).

Jeśli ktoś zna nazwę tej sekwencji lub dlaczego tak się zachowuje, byłbym wdzięczny.