Ostatnio zainteresowałem się metodami wyznaczania liczby π. W zasadzie, to ten temat chodził mi po głowie od pewnego czasu, ale dopiero teraz zagłębiłem się w to zagadnienie. Napisałem kilka programów które obliczają przybliżenie tej liczby, efekty są względnie dobre, na pewno można lepiej to zrobić, ale podzielę się kilkoma obserwacjami.
Liczba π jest liczbą nie wymierną (a przynajmniej na taką się wydaje), czyli nie można określić jej w pełni. Możemy jedynie poznać jej przybliżenie. Jak na razie udało się wykonać obliczenie 699 999 990 000 miejsc po przecinku, a dokonał tego Fabrice Bellard (według wikipedi). Mi udało się tylko 12, ale sprawdzałem na razie tylko proste metody.
Najłatwiejszą metodą jest wzór Leibnitza, który ma postać:
Tego typu szeregów dających przybliżenie π jest kilka, nie będę ich tutaj prezentował, można je znaleźć pod adresem: http://jknow.republika.pl/pi/pi.html
Przy założeniu że dx wynosiło 10^-6 w krótkim czasie obliczyłem 12-tą cyfrę liczby π. Do obliczania całki oznaczonej użyłem metody najmniejszych trapezów. I tutaj pojawił się problem, otóż zmienne typu double i float nie przedstawiają dokładnie tej wartości jaką się do nich wpiszę. Przy obliczeniach z taką dokładnością zaczęło to mieć wpływ na wynik. Próbowałem zastosować bibliotekę ttmath, jednak kompilator wyrzucał wile błędów i z braku czasu zaniechałem dalszy rozwój w tym kierunku. Generalnie ujmując to temat zmiennych zmiennoprzecinkowych jest to przebadania w przyszłym czasie.
Znalazłem jeszcze jedną metodą obliczania π i prawdopodobnie, że przetestuję ją w wolnym czasie. Metoda Monte-Carlo polega na generowaniu losowo punktów wewnątrz prostokąta o boku 1 i zliczaniu ile tych punktów znajdzie się w kole wpisanym w tym prostokącie. Szczegóły i obliczenia są przedstawione na stronie: http://www.mathematica.pl/?przyblizenie-liczby-pi-metoda-monte-carlo.,191 .
Wydaje mi się, że ta metoda pozwoli ominąć jak najdłużej liczby zmiennoprzecinkowe i dzięki temu nie będzie aż tyle błędów nakładających się po drodze.
Jak sprawdzę i tą metodę, to zdam prostą relację z prac. Na dzisiaj to tyle, jeżeli ktoś ma jakieś propozycje to proszę pisać w komentarzach. Dosyć zaciekawił mnie ten temat i chętnie zgłębię go jeszcze bardziej.
Ciekawa jest ta metoda średnich opisana na końcu tego artykułu http://jknow.republika.pl/pi/pi.html
Sprawdziłem z jaką dokładnością da się policzyć obliczając 400 elementów ciągu:
// calculating PI
// author: TS
// explanation: http://jknow.republika.pl/pi/pi.html
#include <ttmath/ttmath.h>
#include <iostream>
typedef ttmath::Big<1,15> Float;
const int max_loop = 400;
void calc_s(std::vector<Float> & tab)
{
Float pi, two, four, add;
int denominator;
bool plus;
pi = 4;
two = 2;
four = 4;
plus = false;
tab[0] = pi;
denominator = 3;
for(size_t i=1 ; i<tab.size() ; ++i)
{
add = four / Float(denominator);
pi += (plus) ? add : -add;
tab[i] = pi;
denominator += 2;
plus = !plus;
}
}
Float calc_avg(std::vector<Float> & tab)
{
Float two, pi;
std::vector<Float> tab2;
bool use_tab2;
size_t a, i;
tab2.resize(tab.size());
two = 2;
use_tab2 = false;
for(a=1 ; a<tab.size() ; ++a)
{
use_tab2 = !use_tab2; // first will be plus
for(i=0 ; i<tab.size()-a ; ++i)
{
if( use_tab2 )
tab2[i] = (tab[i] + tab[i+1]) / two;
else
tab[i] = (tab2[i] + tab2[i+1]) / two;
}
}
if( use_tab2 )
pi = tab2[0];
else
pi = tab[0];
return pi;
}
int main()
{
Float pi;
std::vector<Float> tab;
tab.resize(max_loop);
calc_s(tab);
pi = calc_avg(tab);
std::cout << pi << std::endl;
}
Gdyby się źle wkleiło to tutaj całość:
http://wklej.se/obliczaniepi
Wyszło prawidłowo około 120 cyfr po przecinku, użyłem biblioteki ttmath z typem zmiennoprzecinkowym. Tylko uwaga że przy liczeniu średnich to mamy złożoność O(n^2).
Wielkie dzięki za przesłanie tego programu. Trochę popracowałem w tym kierunku i udało mi się obliczyć 1006 cyfr liczby pi.
Po drodze natrafiłem na ciekawe zagadnienie. Chodzi o to, skąd mamy wiedzieć, ile wynosi faktycznie to pi. Skąd mamy wiedzieć czy to co policzyliśmy jest prawidłowe. Jedynym pewnym rozwiązaniem jakie przyszło mi do głowy było obliczenie z jeszcze większą dokładnością drugi raz. W tedy ta część rozszerzenia dziesiętnego, która jest w obu wynikach jest taka sama, jest prawidłowa.
W moim przypadku wyglądało to tak, że program był podzielony na dwa procesy i równolegle wykonywane były oba obliczania, a następnie zapisywałem je do pliku, gdzie już ręcznie dokonywałem porównania.
Dla informacji, obliczenia tych 1006 cyfr trwały 3,5h.