Quiche Lorraine [nieczynne]

52

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”

Sanchises
źródło
Powiązane: link
Sp3000,
2
Głosuję za zamknięciem tego pytania jako nie na temat, ponieważ nieuczciwe wyzwania nie są już tutaj na temat. meta.codegolf.stackexchange.com/a/8326/20469
kot

Odpowiedzi:

77

Algorytm

Korzystając ze znanego wyniku:

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

Uri Granta
źródło
12
Prawdopodobnie jedna z najlepszych odpowiedzi na podstępne pytanie w ostatnim czasie.
Optymalizator
14
„pomysł matematyki na praktyczny żart”. +1
niestosowne
16
To jest dobre! IIRC, jeden z bardziej znanych żartów z tą całką miał miejsce, gdy ktoś zarejestrował wyniki aż do dziwnego na Wolfram Alpha i wysłał raport o błędzie ... Który
twórcy
3
To odniesienie daje dobre wyjaśnienie tego, co się dzieje.
TonioElGringo
59

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

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

Z warunkiem początkowym

> 0 <y0 <2 * pi

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.

AJMansfield
źródło
4
@UriZarfaty Ten artykuł na Wikipedii całkiem dobrze to wyjaśnia: en.wikipedia.org/wiki/Stiff_equation
AJMansfield
1
Co to jest n? ...
Cole Johnson
1
@AJMansfield Miałem na myśli: nigdzie nie jest zadeklarowany. Twoje forspowolnienie używa t, ale twoja pętla używa n.
Cole Johnson
1
@ColeJohnson Właśnie to naprawiłem.
AJMansfield
2
Myślę, że twoje równanie różniczkowe powinno brzmieć dy / dt = sin (y) * exp (t).
David Zhang
6

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.

soktinpk
źródło
Myślę, że jest to spowodowane anulowaniem, czyli utratą cyfr po odjęciu dwóch prawie równych liczb. S1 i s2 są prawie równe po iteracji.
Sanchises
1
Muszę się jeszcze dowiedzieć, jak to działa, ale przypomina mi coś, co kiedyś zrobiłem: wielokrotnie pobierałem sumę hałaśliwego sygnału i znormalizowałem go do wartości 0, maksymalna wartość 1. To zbiegałoby się w falę sinusoidalną, ponieważ jest to jedyny sygnał, który jest jego własną pochodną anty-pochodną (z przesunięciem fazowym).
Sanchises
Zapytałem o to na stronie math.SE i otrzymałem odpowiedź.
Sanchises