Metody obliczania pierwiastka kwadratowego

Następująca wersja przejrzana tej strony, którą oznaczono 22 lis 2010, była oparta na tej wersji.

Oszacowanie

Wiele metod obliczania pierwiastka kwadratowego z dodatniej liczby rzeczywistej S wymaga wartości początkowej. Jeśli ta wartość jest zbyt odległa od aktualnej wartości pierwiastka, obliczenia będą znacznie wydłużone. w związku z tym jest wysoce pożądane aby mieć oszacowanie tej wielkości, które może być nawet bardzo niedokładne ale proste do wyznaczenia. Jeśli S ≥ 1, niech D będzie liczbą cyfr po lewej stronie przecinka dziesiętnego. Jeśli S < 1, niech D będzie ujemną liczbą zer bezpośrednio na prawo od przecinka dziesiętnego. Wtedy oszacowanie jest następujące:

Jeśli D jest nieparzyste, D = 2n + 1, to

{\sqrt {S}}\approx 2\cdot 10^{n}.

Jeśli D jest parzyste, D = 2n + 2, to

{\sqrt {S}}\approx 6\cdot 10^{n}.

Wybrane zostały dwa i sześć ponieważ ${\sqrt {\sqrt {1\cdot 10}}}={\sqrt[{4}]{10}}\approx 2\,$ i ${\sqrt {\sqrt {10\cdot 100}}}={\sqrt[{4}]{1000}}\approx 6\,$ .

Metoda babilońska

Prawdopodobnie pierwszy algorytm użyty do aproksymacji ${\sqrt {S}}$ jest znany jako metoda babilońska od babilońskich matematyków,^[1] lub metoda Herona od greckiego matematyka z pierwszego wieku Herona z Aleksandrii, który podał pierwszy jasny opis tej metody^[2]. Można ją uzyskać (ale jest znacznie starsza) z metody Newtona. Jest to algorytm zbieżny kwadratowo, co oznacza, że liczba prawidłowych cyfr aproksymacji średnio podwaja się z każdą iteracją. Kroki algorytmu są następujące:

Rozpocznij z dowolną dodatnią wartością początkową x₀ (im bliżej pierwiastka tym lepiej).
Niech x_n+1 będzie średnią z x_n i S / x_n (użycie średniej arytmetycznej do aproksymacji średniej geometrycznej).
Powtarzaj krok 2 aż zostanie osiągnięta pożądana dokładność.

Można go także przedstawić jako:

x_{0}\approx {\sqrt {S}}.

x_{n+1}={\frac {1}{2}}\left(x_{n}+{\frac {S}{x_{n}}}\right),

{\sqrt {S}}=\lim _{n\to \infty }x_{n}.

Przykład

Aby obliczyć ${\sqrt {S}}$ , gdzie S = 125348, do 6 cyfr znaczących, weźmy oszacowanie x₀ z przepisu wyżej. Liczba cyfr w S to D=6=2·2+2. Z tego n=2 i oszacowanie:

x_{0}=6\cdot 10^{2}=600.000\,.

x_{1}={\frac {1}{2}}\left(x_{0}+{\frac {S}{x_{0}}}\right)={\frac {1}{2}}\left(600.000+{\frac {125348}{600.000}}\right)=404.457\,.

x_{2}={\frac {1}{2}}\left(x_{1}+{\frac {S}{x_{1}}}\right)={\frac {1}{2}}\left(404.457+{\frac {125348}{404.457}}\right)=357.187\,.

x_{3}={\frac {1}{2}}\left(x_{2}+{\frac {S}{x_{2}}}\right)={\frac {1}{2}}\left(357.187+{\frac {125348}{357.187}}\right)=354.059\,.

x_{4}={\frac {1}{2}}\left(x_{3}+{\frac {S}{x_{3}}}\right)={\frac {1}{2}}\left(354.059+{\frac {125348}{354.059}}\right)=354.045\,.

x_{5}={\frac {1}{2}}\left(x_{4}+{\frac {S}{x_{4}}}\right)={\frac {1}{2}}\left(354.045+{\frac {125348}{354.045}}\right)=354.045\,.

W związku z tym ${\sqrt {125348}}\approx 354.045\,.$

Zbieżność

Zdefiniujmy błąd względny w x_n jako

\varepsilon _{n}={\frac {x_{n}}{\sqrt {S}}}-1

a tym samym

x_{n}={\sqrt {S}}\cdot (1+\varepsilon _{n}).

Można wykazać, że

\varepsilon _{n+1}={\frac {\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}

a tym samym, że

0\leq \varepsilon _{n+2}\leq \min \left\{{\frac {\varepsilon _{n+1}^{2}}{2}},{\frac {\varepsilon _{n+1}}{2}}\right\}

a w związku z tym zbieżność jest zapewniona pod warunkiem, że x₀ and S są obie dodatnie.

Najgorszy przypadek zbieżności

Przy zastosowaniu oszacowań z przepisu powyżej, najgorsze zbieżności są następujące:

{\begin{aligned}S&=1;&x_{0}&=2;&x_{1}&=1.250;&\varepsilon _{1}&=0.250.\\S&=10;&x_{0}&=2;&x_{1}&=3.500;&\varepsilon _{1}&<0.107.\\S&=10;&x_{0}&=6;&x_{1}&=3.833;&\varepsilon _{1}&<0.213.\\S&=100;&x_{0}&=6;&x_{1}&=11.333;&\varepsilon _{1}&<0.134.\end{aligned}}

Stąd w każdym przypadku,

\varepsilon _{1}\leq 2^{-2}.\,

\varepsilon _{2}<2^{-5}<10^{-1}.\,

\varepsilon _{3}<2^{-11}<10^{-3}.\,

\varepsilon _{4}<2^{-23}<10^{-6}.\,

\varepsilon _{5}<2^{-47}<10^{-14}.\,

\varepsilon _{6}<2^{-95}<10^{-28}.\,

\varepsilon _{7}<2^{-191}<10^{-57}.\,

\varepsilon _{8}<2^{-383}<10^{-115}.\,

Należy pamiętać, że błędy zaokrągleń mogą spowolnić zbieżność. Należy wyznaczać przynajmniej jedną cyfrę więcej niż pożądana dokładność x_n aby zminimalizować błędy zaokrągleń.

Tożsamość wykładnicza

Kalkulatory zwykle implementują dobre procedury na obliczanie funkcji ekspotencjalnej i logarytmu naturalnego, i z tego obliczają pierwiastek kwadratowy z S wykorzystując tożsamość

{\sqrt {S}}=e^{{\frac {1}{2}}\ln S}

.

Ta sama tożsamość jest używana do wyliczania pierwiastka kwadratowego przy użyciu tablic logarytmicznych lub suwaka logarytmicznego.

Metoda równego podziału

Innym prostym sposobem znalezienia pierwiastka kwadratowego jest metoda więcej/mniej, podobnie jak w metodzie równego podziału. Ta metoda polega na zgadywaniu liczby na podstawie znanego kwadratu, a następnie sprawdzaniu czy jej kwadrat jest zbyt wysoki lub zbyt niski i odpowiednim poprawianiu wyniku.

Załóżmy, że chcemy znaleźć pierwiastek kwadratowy z 20. Wiemy, że kwadrat 5 wynosi 25, a kwadrat 4 to 16, czyli wynik musi być pomiędzy 4 i 5. Na początku zgadujemy, że 4,5. Kwadrat 4,5 wynosi 20,25 a to jest za dużo, więc zgadujemy dla 4,4. Wynik równa się 19,36 a to za mało. Z tego jednak wynika, że pierwiastek jest pomiędzy 4,4 i 4,5. Kontynuując ten schemat możemy uzyskać taką dokładność jaką tylko chcemy. Aby uzyskać trzy cyfry po przecinku:

4.45² = 19.8025 (za mało)

4.47² = 19.9809 (za dużo, ale blisko)

4.48² = 20.0704 (za dużo)

4.475² = 20.025625 (za dużo)

4.473² = 20.007729 (za dużo, ale blisko)

4.472² = 19.998784 (za mało)

Teraz wiemy, że wynik jest pomiędzy 4,472 a 4,473. Przyjmujemy, że pierwiastek kwadratowy z 20 dla dokładności do trzech miejsc po przecinku wynosi 4,472.

Aproksymacja Bakhshali'ego

Ta metoda znajdowania aproksymacji pierwiastka kwadratowego została opisana w starożytnym manuskrypcie zwanego manuskrypt Bakhshali'ego. Jest ona równoważna dwóm iteracjom metody babilońskiej rozpoczynając od N. Oryginalna prezentacja jest następująca: Aby obliczyć ${\sqrt {S}}$ , niech N² będzie najlepszym przybliżeniem kwadratu S. Następnie oblicz:

d=S-N^{2}\,\!

P={\frac {d}{2N}}

A=N+P\,\!

{\sqrt {S}}\approx A-{\frac {P^{2}}{2A}}

To można także zapisać jako:

{\sqrt {S}}\approx N+{\frac {d}{2N}}-{\frac {d^{2}}{8N^{3}+4Nd}}={\frac {8N^{4}+8N^{2}d+d^{2}}{8N^{3}+4Nd}}

Przykład

Znajdźmy ${\sqrt {9.2345}}$ .

N=3\,\!

d=9.2345-3^{2}=0.2345\,\!

P={\frac {0.2345}{2\times 3}}=0.0391

A=3+0.0391=3.0391\,\!

{\sqrt {9.2345}}\approx 3.0391-{\frac {0.0391^{2}}{2\times 3.0391}}\approx 3.0388

Wyznaczanie cyfra po cyfrze

Jest to metoda sekwencyjnego znajdowania kolejnych cyfr pierwiastka kwadratowego. Jest ona wolniejsza niż metoda babilońska, lecz ma kilka zalet:

jest prosta dla ręcznych obliczeń;
każda cyfra wyznaczanego pierwiastka jest prawidłowa, to znaczy nie zmieni się w toku dalszych obliczeń;
jeśli pierwiastek kwadratowy ma skończone rozwinięcie, algorytm kończy się po znalezieniu ostatniej cyfry. Można to wykorzystać do sprawdzenia czy zadana liczba całkowita jest kwadratem.

Kostki Napiera zawierają wskazówki do wykonania tego algorytmu.

Dziesiątkowo

Zapisujemy oryginalną liczbę w postaci dziesiętnej. Liczby są zapisywane podobnie jak w dzieleniu pisemnym i podobnie wynik będzie zapisany nad linią. Następnie należy cyfry pogrupować w pary, rozpoczynając od przecinka dziesiętnego w obie strony. Przecinek dziesiętny będzie pozycją przecinka dziesiętnego wyniku nad linią. Jedna cyfra wyniku będzie się pojawiała nad każdą parą cyfr pod linią.

Rozpoczynamy od pierwszej pary cyfr od lewej i wykonujemy poniższą procedurę dla każdej kolejnej pary:

Rozpoczynamy od lewej przepisując najbardziej znaczącą parę cyfr jeszcze nie używaną (jeśli zostały wykorzystane wszystkie zapisujemy "00") i dopisujemy je z prawej strony reszty z poprzedniego etapu (w pierwszym kroku nie ma reszty). Innymi słowy, mnożymy resztę przez 100 i dodajemy 2 cyfry. Będzie to bieżąca wartość c.
Znajdujemy p, y i x następująco:
- Niech p będzie częścią pierwiastka aktualnie znalezioną, pomijamy przecinek dziesiętny. (W pierwszym kroku, p = 0).
- Określamy największą cyfrę x taką, że $y=(20\cdot p+x)\cdot x$ $y=(20\cdot p+x)\cdot x$ nie przekracza c.
  - Uwaga: 20p + x to zwyczajnie dwa razy p, z cyfrą x dołączoną po prawej).
  - Uwaga: Można znaleźć x zgadując, że to c/(20·p) i wykonując próbne obliczenie y, a następnie zmienić x w górę lub dół stosownie do potrzeb.
- Umieszczamy cyfrę $x$ jako następną cyfrę pierwiastka.
Odejmujemy y od c i otrzymujemy nową resztę.
Jeśli reszta wynosi zero i nie ma więcej cyfr do spisania, to algorytm jest zakończony. W innym przypadku wracamy do punktu 1 i kontynuujemy następną iterację.

Przykłady

Znajdź pierwiastek kwadratowy z 152.2756.

          1  2. 3  4 
       /
     \/  01 52.27 56                            

         01                   1*1 <= 1 < 2*2                 x = 1
         01                     y = x*x = 1*1 = 1
         00 52                22*2 <= 52 < 23*3              x = 2
         00 44                  y = (20+x)*x = 22*2 = 44                      
            08 27             243*3 <= 827 < 244*4           x = 3       
            07 29               y = (240+x)*x = 243*3 = 729
               98 56          2464*4 <= 9856 < 2465*5        x = 4       
               98 56            y = (2460+x)*x = 2464*4 = 9856
               00 00          Algorytm zakończony: Wynik to 12.34

Znajdź pierwiastek kwadratowy z 2.

          1. 4  1  4  2
       /
     \/  02.00 00 00 00

         02                  1*1 <= 2 < 2*2                 x = 1
         01                    y = x*x = 1*1 = 1
         01 00               24*4 <= 100 < 25*5             x = 4
         00 96                 y = (20+x)*x = 24*4 = 96                      
            04 00            281*1 <= 400 < 282*2           x = 1       
            02 81              y = (280+x)*x = 281*1 = 281
            01 19 00         2824*4 <= 11900 < 2825*5       x = 4       
            01 12 96           y = (2820+x)*x = 2824*4 = 11296
               06 04 00      28282*2 <= 60400 < 28283*3     x = 2
                             Osiągnięto pożądaną dokładność 
                             Pierwiastek kwadratowy z 2 wynosi około 1.4142

Dwójkowo

Nieodłącznym krokiem w algorytmie cyfra po cyfrze jest szukaj i sprawdź: znajdź cyfrę, $\,e$ , która dodana z prawej bieżącego rozwiązania $\,r$ , takiego, że $\,(r+e)\cdot (r+e)\leq x$ , gdzie $\,x$ jest pożądaną wartością pierwiastka. Rozwijając, otrzymujemy $\,r\cdot r+2re+e\cdot e\leq x$ . Bieżąca wartość $\,r\cdot r$ —lub, zwykle, reszta—mogą być przyrostowo aktualizowane w systemie dwójowym, gdzie wartość $\,e$ jest jednobitowa, a operacje wymagane do obliczenia $\,2\cdot r\cdot e$ i $\,e\cdot e$ można zastąpić szybszymi operacjami przesunięć bitów. W wyniku otrzymujemy prosty algorytm komputerowy:^[3]

    short sqrt(short num) {
        short op = num;
        short res = 0;
        short one = 1 << 14; // The second-to-top bit is set: 1L<<30 for long
 
        // "one" starts at the highest power of four <= the argument.
        while (one > op)
            one >>= 2;
        
        while (one != 0) {
            if (op >= res + one) {
                op -= res + one;
                res = (res >> 1) + one;
            }
            else
              res >>= 1;
            one >>= 2;
        }
        return res;
    }

↑ Nie ma bezpośrednich dowodów wskazujących jak Babilończycy obliczali pierwiastki kwadratowe, mimo to są pewne przesłanki (pierwiastek kwadratowy z 2).
↑ Thomas Heath: A History of Greek Mathematics. T. 2. Oxford: Clarendon Press, 1921, s. 323-324.
↑ Fast integer square root by Mr. Woo's abacus algorithm

[1] Nie ma bezpośrednich dowodów wskazujących jak Babilończycy obliczali pierwiastki kwadratowe, mimo to są pewne przesłanki (pierwiastek kwadratowy z 2).

[2] Thomas Heath: A History of Greek Mathematics. T. 2. Oxford: Clarendon Press, 1921, s. 323-324.

[3] Fast integer square root by Mr. Woo's abacus algorithm

[1]

[2]

[3]