Jacek Matulewski Zakład Mechaniki Kwantowej Instytut Fizyki, UMK 14 marca 2011

1. Obliczanie p Jacek Matulewski Zakład Mechaniki Kwantowej Instytut Fizyki, UMK 14 marca 2011

2. Obliczanie p Jacek Matulewski Zakład Mechaniki Kwantowej Instytut Fizyki, UMK 14 marca 2011

3. Liczby • Liczby wyrażają ilość(liczby naturalne) • Pitagoras i pitagorejczycy – kult i symbolika liczb

4. Liczby • Jeden • Monas • Bóg • Parmenidejski byt

5. Liczby • Dwa • Relacja „po” • Czas • Arytmetyka

6. Liczby • Trzy • Relacja „obok” • Przestrzeń • Geometria

7. Liczby • Pitagorejczycy użyli liczb także do mierzenia wielkości ciągłych(długości boków figur) • Długość przekątnejz twierdzenie Pitagorasa • Wielka tajemnica…

8. Liczby p

9. Liczby • Liczba p jest liczbą niewymiernątzn. nie można jej zapisać jako ułamka zwykłego (ilorazu dwóch liczb całkowitych) • Liczba p jest liczbą przestępnątzn. nie istnieje wielomian o współczynnikach całkowitych, którego p jest pierwiastkiem • Pitagorejczycy utrzymywali liczby niewymierne w tajemnicy!

10. Wartość liczby p • Wartość p: 3,14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165 …

11. Przybliżenia liczby p • Wartość: 3,14159 26535 89793 23846 26433 83279… • Przybliżenia:22/7 = 3,14285 71428 57142 85714…355/113 = 3,14159 29203 53982 30088…52163/16604 = 3,1415923873 76535 77451…103993/33102 = 3,14159 26530 11902 60407…3,14626 43699 41972 34232 91350…3,14159..(30 cyfr)

12. Rekordy

13. Rekordy

14. Metoda Archimedesa

15. Metoda Monte-Carlo • Okrąg: r = 1 • Pole koła: pr2 = p • Pole kwadratu:a2 = 4 • Stosunek pól: p/4

16. Metoda Monte-Carlo • Autor: Rudolf Wolf(szwajcarski astronom) • Dokładność (N – ilość prób) • To raczej słabo!

17. Metoda Monte-Carlo • Realizacje (kłopot z liczbami losowymi) • dziecko rzuca kamyki na plaży (KM i BM) • deszcz i donica w akwarium (M. Chmielarz) • mąka na stole

18. Metoda Monte-Carlo • Realizacje (kłopot z liczbami losowymi) • dziecko rzuca kamyki na plaży (KM i BM) • deszcz i donica w akwarium (M. Chmielarz) • mąka na stole • symulacja komputerowa • http://www.fizyka.umk.pl/~jacek/dydaktyka/pi/(C, C++, C#, Java, Pascal/Delphi, nawet ZX Basic)

19. Metoda Monte-Carlo • Implementacja ulong i=0, N=(ulong)1E8; double x, y; do { i++; x=random.NextDouble(); y=random.NextDouble(); if(x*x+y*y<1) trafienia++; } while (i<N); double pi=4*(double)trafienia/N;

20. Metoda Monte-Carlo Demo

21. Igła Buffona • Problem sformułowany w 1773 roku przez Georga-Luisa Leclerca • Nietrywialny przykład realizacji metody Monte-Carlo • Typowe zadanie rach. prawdopodobieństwa(prawdopodobieństwo geometryczne)

22. Igła Buffona • Rzucamy N razy igłą na dywan w pasy • Igła o długości l jest krótsza niż odległość pasów d (l < d) • Zliczamy ile razy igła przetnie którąś z linii

23. Igła Buffona y – rzut połowy igły na kierunek poziomy x – odległość środka igły od najbliższej linii Igła przetnie linię gdy y > x y W rzucie igłą mamy dwa losowe parametry o jednorodnym rozkładzie: – położenie środka igły, co się przekłada na wielkość x, przy czym – orientację mierzoną kątem ostrym q między igłą i pionem, Rzut połowy igły: Igła przetnie linię gdy q x

24. Igła Buffona Losujemy dwa parametry (x, q) – punkt w przestrzeni 2D o rozmiarach Prawdopodobieństwo przecięcia linii: obszar pod krzywą pole całego prostokąta Prawdopodobieństwo = częstość losowania

25. Metoda Eulera • Autor: Leonhard Euler(szwajcarski matematyk) • Dokładność 1/N(N – ilość prób) • Znacznie lepiej niż w metodzie Wolfa

26. Metoda Eulera

27. Metoda Eulera • Sumowanie szeregu • Suma to p2/6 Można z grubsza przyjąć, że ilość cyfr znaczących oszacowania wartości prówna jest rzędowi zsumowanych wyrazów. Zatem dla 105 wyrazów otrzymamy pięć cyfr znaczących.

28. Metoda Eulera • Implementacja ulong i=0, N=(ulong)1E8; double S=0; do { i++; S+=(1.0/(1.0*i*i)); } while (i<N); double pi=Math.Sqrt(6*S);

29. Metoda Eulera • Wartości częściowych sum

30. Metoda Eulera Demo

31. Metody obliczania p • Inne szeregi: http://mathworld.wolfram.com/PiFormulas.html • szereg Madhava’ego, Gregory’ego, Reihe’a i Leibniza (1671): • szereg Abrahama Sharpa (1717):

32. Metody obliczania p • Inne szeregi: http://mathworld.wolfram.com/PiFormulas.html

33. Metody obliczania p • Inne szeregi: http://mathworld.wolfram.com/PiFormulas.html • Szereg Wallisa • Iloczyn Wallisa (1655)

34. Metoda iteracyjna Salamina-Brenta • Metody używane w XX wieku: D.H. Bailey, J.M. Borwein, P.B. Borwein, S. PlouffeThe Quest for Pi (1996) • Algorytm iteracyjny Salamina i Brenta (1976)Ustalamy: , , Dla k=1,2,3,… obliczamy pk zbiega do p kwadratowo,czyli każda iteracja podwajailość prawidłowych cyfr. Jeden z algorytmów używanych do bicia rekordów ilości policzonych cyfr p (obliczenia komputerowe)

35. Metody iteracyjne braci Borwein • Algorytm braci Borwein (1985):ciąg zbiegający do pkubicznie • Algorytm dwukwadratowy braci Borweinużywany przez Yasumasa Kanadę (Tokio) • Znaleźli wzór na ciągi dowolnego rzędu(co nie znaczy, że są szybsze od 4-go rzędu) • Dodatkowo mnożenia w wysokiej precyzji przyspieszane są dzięki użyciu algorytmu FFT (1965)

36. Metody iteracyjne braci Borwein • Yasumasa Kanada w swoim biurze

37. Metody iteracyjne braci Borwein Źródło: D.H. Bailey, J.M. Borwein, P.B. Borwein, S. PlouffeThe Quest for Pi (1996)

38. Wartość liczby p • Wartość p: 3,243F6 A8885 A308D 31319 8A2E0 37073 44A40 93822 299F3 1D008 2EFA9 8EC4E 6C894 52821 E638D 1377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF708 …

39. Algorytmy kurkowe(?) • Niepewna nazwa: ang. spigotalgoritm(spigot = kran przy beczce, szpunt) • Algorytm, który pozwala na liczenie kolejnych cyfr stałej matematycznej nie używając ich już po policzeniu (w przeciwieństwie do wcześniej przedstawionych algorytmów rekurencyjnych) • Takie algorytmy istnieją dla: p, e, , itp. • Algorytmy te można zaimplementować bez liczb rzeczywistych (tylko liczby całkowite)

40. Algorytmy kurkowe(?) • Do policzenia 5000 cyfr liczby p potrzeba tylko 600 000 000 operacji na liczbach całkowitych! • Koszt obliczenia n-tej cyfry nie zależy od n. • Algorytm kurkowy dla liczby p toalgorytm Baileya-Borweina-Plouffe’a(BBP; 1995, 1997, 2003) bazujący na wzorze: • Pozwala na znalezienie konkretnej serii cyfr p!

41. Algorytm BBP David Bailey Peter Borwein Simon Plouffe

42. Algorytm BBP • Przygotowania:

43. Algorytm BBP • Aby obliczyć n-tą cyfrę liczby lw reprezentacji dziesiętnej można użyć triku: • Aby zaokrąglić liczbę l do n cyfr po przecinku: l*10^n – trunc(l*10^n) 10000*p – (int)(10000*p) = 0,9265358979 l*10^n % 1 trunc(l*10^n) / 10^n (int)(10000*p) / 10000 = 3,1415

44. Algorytm BBP • Aby obliczyć n-tą cyfrę liczby lw reprezentacji heksydecymalnej należy: • Aby zaokrąglić liczbę l do n cyfr po przecinku: l*16^n – trunc(l*16^n) 65536*p – (int)(65536*p) = 0,6A8885A308D31319 l*16^n % 1 trunc(l*16^n) / 16^n (int)(65536*p) / 65536 = 3,243F

45. Algorytm BBP • Aby obliczyć n-tą cyfrę liczby lw reprezentacji heksydecymalnej należy: • Aby zaokrąglić liczbę l do n cyfr po przecinku: l*10^n – trunc(l*10^n) 10000*p – (int)(10000*p) = 0,6A8885A308D31319 l*10^n % 1 trunc(l*10^n) / 10^n (int)(10000*p) / 10000 = 3,243F

46. Algorytm BBP • Cyfry sum S(j) począwszy od n-tej Bardzo szybko zbieżny szereg (sprawdzamy tylko, czy nie zmienia się ostatnia cyfra serii) część ułamkowa Wartość p od n-tej cyfry HEX Obliczanie potęg za pomocą algorytmu szybkiego potęgowania dla l. całkowit.(exponentiation by squaring)

47. Algorytm BBP Demo

48. Motywacja Po co tak dokładnie liczyć przybliżenia p?

49. Motywacja Po co tak dokładnie liczyć przybliżenia p?

50. Koniec • Wartość p: 3,14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 44288 10975 66593 34461 28475 64823 37867 83165 …