28/07/2013

Ciekawa metoda obliczenia pierwiastka kwadratowego

Home

Istnieje wiele metod wyznaczania pierwiastka kwadratowego z liczb naturalnych. Ja napiszę o jednej, która mnie ostatnio urzekła i pokazuje piękno matematyki.

Otóż, wynikiem zastosowania tej metody jest ułamek, który stanowi przybliżenie pierwiastka z danej liczby. Dlaczego tylko przybliżenie? Pierwiastki kwadratowe z liczb naturalnych, nie będących kwadratem innej liczby naturalnej, są liczbami niewymiernymi. Liczb niewymiernych nie da się natomiast zapisać w postaci ułamka.

Metoda ta opiera się o tzw. ułamki łańcuchowe (ang. continued fraction). W skrócie chodzi o to, aby pierwiastek kwadratowy z danej liczby przedstawić w następujący sposób:
sqrt(S) = a0 + 1 / (a1 + 1 / (a2 + 1 / (a3 + 1 / (...))))
Nim dłuższa sekwencja tym otrzymamy dokładniejsze przybliżenie. Co bardzo ciekawe, dla każdej liczby niewymiernej sekwencja wartości a0, a1, a2... w pewnym momencie stworzy cykl. Na przykład dla sqrt(2) sekwencja ta ma postać:
[1, 1, 2, 1, 2, ...]
Natomiast dla sqrt(53):
[7, 3, 1, 1, 3, 14, 3, 1, 1, 3, 14, ...].
Jak wyznaczyć taką sekwencję? Algorytm nie jest bardzo skomplikowany i można go zapisać w kilkunastu liniach kodu. Nie będą go jednak tutaj przytaczał, ponieważ doskonały opis znajduje się już na Wikipedii. Zachęcam do wypróbowania swoich sił i jego zaimplementowania.

Skorzystać możemy również z silnika Wolframalpha. Komenda continued fraction sqrt(X) zwróci nam rozwinięcie pierwiastka z danej liczby X w postaci ułamka łańcuchowego. Natomiast komenda convergents[sqrt(X),N] obliczy N pierwszych przybliżeń pierwiastka z danej liczby X w oparciu o opisaną metodę.

Metoda ta z pewnością może przydać się w obliczeniach naukowych, kiedy bardzo istotna jest precyzja obliczeń i operacje na zmiennym przecinku są niedopuszczalne. W takim przypadku powinniśmy wyznaczyć odpowiednio dokładne przybliżenie sqrt(x) w postaci ułamka.

Na przykład ułamek sqrt(3) ~= 13623482 / 7865521 da nam tyle samo poprawnie obliczonych cyfr po przecinku co Math.Pow(3, -0.5). W tym przypadku sekwencja a0, a1, a2... miała 25 elementów. Każdy dodatkowy element sekwencji da nam większą precyzję. W skrajnym przypadku możemy skorzystać ze struktury BigInteger do reprezentacji licznika i mianownika.

Jeśli temat Was zainteresował to polecam rozwiązanie zadania 64, 65 lub 66 ze strony Project Euler.

0 comments:

Post a Comment