Transcript pps

Obliczanie
Jacek Matulewski
Zakład Mechaniki Kwantowej
Instytut Fizyki, UMK
14 marca 2011
Obliczanie
Jacek Matulewski
Zakład Mechaniki Kwantowej
Instytut Fizyki, UMK
14 marca 2011
Liczby
• Liczby wyrażają ilość
(liczby naturalne)
• Pitagoras i
pitagorejczycy –
kult i symbolika liczb
Liczby
•
•
•
•
Jeden
Monas
Bóg
Parmenidejski byt
Liczby
•
•
•
•
Dwa
Relacja „po”
Czas
Arytmetyka
Liczby
•
•
•
•
Trzy
Relacja „obok”
Przestrzeń
Geometria
Liczby
• Pitagorejczycy użyli
liczb także do mierzenia
wielkości ciągłych
(długości boków figur)
• Długość przekątnej
z twierdzenie Pitagorasa
2
d  12  12  2
• Wielka tajemnica…
Liczby
2
p
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!
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 …
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,14159 23873 76535 77451…
103993/33102 = 3,14159 26530 11902 60407…
2 3
3,14626 43699 41972 34232 91350…
ln(6403203  744) / 163  3,14159..(30 cyfr)
Rekordy
Kto:
Kiedy:
Prec.
Wartość
Babilończycy
2000 p.n.e
1
3.125
Egipcjanie
2000 p.n.e
1
3.16045
Chińczycy
1200 p.n.e.
1
3
Biblia
550 p.n.e.
1
3 (1 Król. 7:23)
Archimedes
250 p.n.e.
3
3.141851±0.00100604
Ptolemeusz
150
3
3.14166
Tsu Ch’ung Chi
480
7
3.1415926
Fibonacci
1220
3
3.141818
Al-Kashi
1429
14
Viete
1593
9
Newton
1665
16
Sharp
1699
71
Rutheford, W.
1824
208 (152)
3.1415926536
W 1853 policzył 404
Rekordy
Kto:
Kiedy:
Prec.
Ferguson
1946, 1947
620, 808
Reitwiesner i in.
1949
2 037 (ENIAC)
Felton
1957, 1958
7 480, 10 021
Shanks i Wrench
1961
100 265
Guilloud i in
1966 – 1982
250 000 – 2 000 050
Kanada i inni
1981 – 1995
2 000 036 – 6 442 450 938
bracia Chudovsky
1989 – 1994
480 000 000 – 4 044 000 000
Y. Kanada
2002
1 241 100 000 000 (600 godzin)
Metoda Archimedesa
10
10
3 p 3
71
70
3,140845 p  3,142857
Metoda Monte-Carlo
• Okrąg: r = 1
• Pole koła:
pr2  p
• Pole kwadratu:
a2 = 4
• Stosunek pól: p/4
Metoda Monte-Carlo
• Autor: Rudolf Wolf
(szwajcarski astronom)
• Dokładność 1 / N
(N – ilość prób)
• To raczej słabo!
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
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)
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;
Metoda Monte-Carlo
Demo
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)
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
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
q
x
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 0  x  d 2
– orientację mierzoną kątem ostrym q
między igłą i pionem, 0  q  p 2
Rzut połowy igły: y 
l
sin(q )
2
Igła przetnie linię gdy x  y 
l
sin(q )
2
Igła Buffona
Losujemy dwa parametry (x, q) – punkt
w przestrzeni 2D o rozmiarach
[0, d 2] [0, p 2]
Prawdopodobieństwo przecięcia linii:
p2
l l
sin(q)dq

2 pod
obszar
22l krzywą
0
pole całego
dppddpprostokąta

2 42
Prawdopodobieństwo = częstość losowania
2l
n

pd N
2lN
p
nd
x y
l
sin(q )
2
Metoda Eulera
• Autor: Leonhard Euler
(szwajcarski matematyk)
1 p2


2
6
i 1 i

• Dokładność 1/N
(N – ilość prób)
• Znacznie lepiej niż
w metodzie Wolfa
Metoda Eulera
Metoda Eulera
• Sumowanie szeregu
n
1 1 1 1
1
S n   2  2  2  2  ... 2
1 2 3
n
i 1 i
• Suma S to p2/6
Można z grubsza przyjąć, że ilość cyfr znaczących oszacowania wartości p
równa jest rzędowi zsumowanych wyrazów.
Zatem dla 105 wyrazów otrzymamy pięć cyfr znaczących.
6  Sn  p  6  Sn  1 / n  6  Sn  1 / n
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);
Metoda Eulera
• Wartości częściowych sum
n
6  S n  6
i 1
1
i2
1
2,44948974278318
10
3,04936163598207
2
2,73861278752583
100
3,13207653180911
3
2,85773803324704
1 000
3,14063805620599
4
2,92261298612503
10 000
3,14149716394721
5
2,96338770103857
100 000
3,14158310432646
6
2,99137649474842
1 000 000
3,14159169866051
7
3,01177394784621
10 000 000
3,14159255809590
8
3,02729785665784
100 000 000
3,14159264498239
9
3,03950758956105
1 000 000 000
przekracza precyzję
liczb podwójnej
precyzji (double)!
Metoda Eulera
Demo
Metody obliczania p
• Inne szeregi:
http://mathworld.wolfram.com/PiFormulas.html
– szereg Madhava’ego, Gregory’ego, Reihe’a
i Leibniza (1671):
p  (1)i
1 1 1

 1     ...
4 i 0 2i  1
3 5 7
– szereg Abrahama Sharpa (1717):
2(1)i 31/ 2i
p 
2i  1
i 0

Metody obliczania p
• Inne szeregi:
http://mathworld.wolfram.com/PiFormulas.html

2
(1)i (1)i
1 1 1 1 1
p 

 1       ...
4
4i  1
3 5 7 9 11
i 0 4i  3
p 3
4
p2
(1)i
1
1
1




 ...
2  3 4 4  5 6 6  7 8
i 0 2i  22i  32i  4

1
1 1
p2 n
1
1 1

 1    ...
  2  1    ...
2
8
9 25
6
4 9
i 1 2i  1
i 1 i
n
Metody obliczania p
• Inne szeregi:
http://mathworld.wolfram.com/PiFormulas.html
- Szereg Wallisa
(i!) 2 2i 1 1 2 1 4 2  8
2 2
p 



 ...  2    ...
1
6 120
3 15
i 0 2i  1!

- Iloczyn Wallisa (1655)
p

4i 2
2 2 4 4 6 6 8 8
 2
         ...
2 i 0 4i  1 1 3 3 5 5 7 7 9
Metoda iteracyjna Salamina-Brenta
• Metody używane w XX wieku:
D.H. Bailey, J.M. Borwein, P.B. Borwein, S. Plouffe The Quest for Pi (1996)
- Algorytm iteracyjny Salamina i Brenta (1976)
Ustalamy: a0  1 , b0  1 2 , s0  1 2
Dla k=1,2,3,… obliczamy
ak 
ak 1  bk 1
2
bk  ak 1bk 1
ck  ak2  bk2
sk  sk 1  2k ck
Jeden z algorytmów używanych do bicia rekordów
ilości policzonych cyfr p (obliczenia komputerowe)
2ak2
pk 
sk
pk zbiega do p kwadratowo,
czyli każda iteracja podwaja
ilość prawidłowych cyfr.
Metody iteracyjne braci Borwein
• Algorytm braci Borwein (1985):
ciąg zbiegający do p kubicznie
• Algorytm dwukwadratowy braci Borwein
uż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)
Metody iteracyjne braci Borwein
• Yasumasa Kanada w swoim biurze
Metody iteracyjne braci Borwein
Źródło: D.H. Bailey, J.M. Borwein, P.B. Borwein, S. Plouffe The Quest for Pi (1996)
Wartość liczby p
• Wartość p:
3,243F6 A8885 A308D 31319 8A2E0 37073
44A40 93822 299F3 1D008 2EFA9 8EC4E
6C894 52821 E638D 1377B E5466 CF34E
90C6C C0AC2 9B7C9 7C50D D3F84
D5B5B 54709 17921 6D5D9 8979F B1BD1
310BA 698DF B5AC2 FFD72 DBD01
ADFB7 B8E1A FED6A 267E9 6BA7C
9045F 12C7F 9924A 19947 B3916 CF708 …
Algorytmy kurkowe(?)
• Niepewna nazwa: ang. spigot algoritm
(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, 2 , itp.
• Algorytmy te można zaimplementować bez
liczb rzeczywistych (tylko liczby całkowite)
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 to
algorytm Baileya-Borweina-Plouffe’a
(BBP; 1995, 1997, 2003) bazujący na wzorze:

1
p  i
i 0 16
2
1
1 
 4





 8i  1 8i  4 8i  5 8i  6 
• Pozwala na znalezienie konkretnej serii cyfr p!
Algorytm BBP
David Bailey
Peter Borwein
Simon Plouffe
Algorytm BBP
• Przygotowania:

1
p  i
i 0 16
2
1
1 
 4





 8i  1 8i  4 8i  5 8i  6 
p  4S (1)  2S (4)  S (5)  S (6)
1  1
S ( j )   i 
i 0 16  8i 



j
Algorytm BBP
• Aby obliczyć n-tą cyfrę liczby l
w reprezentacji dziesiętnej można użyć triku:
 
10 l  10 l
n
n
n
10 l mod 1
l*10^n – trunc(l*10^n)
10000*p – (int)(10000*p) = 0,9265358979
l*10^n % 1
• Aby zaokrąglić liczbę l do n cyfr po przecinku:
10 l 
n
trunc(l*10^n) / 10^n
n
(int)(10000*p) / 10000 = 3,1415
10
Algorytm BBP
• Aby obliczyć n-tą cyfrę liczby l
w reprezentacji heksydecymalnej należy:
 
16 l  16 l
n
n
n
16 l mod 1
l*16^n – trunc(l*16^n)
65536*p – (int)(65536*p) = 0,6A8885A308D31319
l*16^n % 1
• Aby zaokrąglić liczbę l do n cyfr po przecinku:
16 l 
n
trunc(l*16^n) / 16^n
n
(int)(65536*p) / 65536 = 3,243F
16
Algorytm BBP
• Aby obliczyć n-tą cyfrę liczby l
w reprezentacji heksydecymalnej należy:
 
10 l  10 l
n
n
n
10 l mod 1
l*10^n – trunc(l*10^n)
10000*p – (int)(10000*p) = 0,6A8885A308D31319
l*10^n % 1
• Aby zaokrąglić liczbę l do n cyfr po przecinku:
10 l 
n
trunc(l*10^n) / 10^n
n
(int)(10000*p) / 10000 = 3,243F
10
Algorytm BBP
• Cyfry sum S( j) począwszy od n-tej

1
S ( j)   i
i 0 16
 1

 8i 


j
n i
16
n
16
 S ( j)   8i  j 
i 0

część
ułamkowa
Obliczanie potęg za pomocą algorytmu
szybkiego potęgowania dla l. całkowit.
(exponentiation by squaring)

16ni mod8i  j
16ni


8i  j
i 0
i  n 18i  j
n
Wartość p od
n-tej cyfry HEX
p  4S (1)  2S (4)  S (5)  S (6)
Bardzo szybko
zbieżny szereg
(sprawdzamy tylko,
czy nie zmienia się
ostatnia cyfra serii)
Algorytm BBP
Demo
Motywacja
Po co tak dokładnie liczyć przybliżenia p?
Motywacja
Po co tak dokładnie liczyć przybliżenia p?
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 …