Metody obliczania pierwiastka kwadratowego

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 faktycznej 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 ${\displaystyle {\sqrt {S}}\approx 2\cdot 10^{n}.}$

Jeśli D jest parzyste, D = 2n + 2, to ${\displaystyle {\sqrt {S}}\approx 6\cdot 10^{n}.}$

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

Metoda babilońska

(c) Crotalus horridus at the English Wikipedia, CC-BY-SA-3.0

Wykres z zastosowania metody babilońskiej do aproksymacji pierwiastka kwadratowego ze 100 (10) przy użyciu wartości początkowych x₀ = 50, x₀ = 1, i x₀ = −5. Ujemne wartości początkowe dają ujemny wynik.

Prawdopodobnie pierwszy algorytm użyty do aproksymacji ${\displaystyle {\sqrt {S}}}$ jest znany jako metoda babilońska od babilońskich matematyków^[a], lub metoda Herona od greckiego matematyka z pierwszego wieku Herona z Aleksandrii, który podał pierwszy jasny opis tej metody^[1]. 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:

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

Można go także przedstawić jako:

${\displaystyle x_{0}\approx {\sqrt {S}},}$

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

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

Przykład

Aby obliczyć ${\displaystyle {\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:

${\displaystyle x_{0}=6\cdot 10^{2}=600{,}000}$

${\displaystyle 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,}$

${\displaystyle 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,}$

${\displaystyle 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,}$

${\displaystyle 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,}$

${\displaystyle 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 ${\displaystyle {\sqrt {125348}}\approx 354,045.}$

Zbieżność

Zdefiniujmy błąd względny w x_n jako

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

a tym samym

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

Można wykazać, że

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

a tym samym, że

${\displaystyle 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₀ i S są obie dodatnie.

Najgorszy przypadek zbieżności

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

${\displaystyle {\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,

${\displaystyle \varepsilon _{1}\leq 2^{-2},}$

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

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

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

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

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

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

${\displaystyle \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 wykładniczej i logarytmu naturalnego, i z tego obliczają pierwiastek kwadratowy z S wykorzystując tożsamość

${\displaystyle {\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:

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

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

${\displaystyle d=S-N^{2}}$

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

${\displaystyle A=N+P}$

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

To można także zapisać jako:

${\displaystyle {\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 ${\displaystyle {\sqrt {9.2345}}.}$

${\displaystyle N=3}$

${\displaystyle d=9,2345-3^{2}=0,2345}$

${\displaystyle P={\frac {0,2345}{2\times 3}}=0,0391}$

${\displaystyle A=3+0,0391=3,0391}$

${\displaystyle {\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ę pierwiastkowaną w postaci dziesiętnej. Liczby są zapisywane podobnie jak w dzieleniu pisemnym i podobnie wynik pierwiastkowania będzie zapisany nad linią. Następnie należy cyfry liczby pierwiastkowanej 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:

1. 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 liczbę wyrażoną przez kolejne 2 cyfry. Będzie to bieżąca wartość c.
2. Znajdujemy p, y i x następująco:
   • Niech p będzie dotychczas obliczoną częścią wyniku pierwiastkowania zapisaną nad kreską, pomijamy przecinek dziesiętny. (W pierwszym kroku, p = 0).
   • Określamy największą cyfrę x taką, że ${\displaystyle y=(20\cdot p+x)\cdot x}$ nie jest większe od c (czyli może być mniejsze lub równe).
     • Uwaga: 20p + x to oczywiście w zapisie cyfrowym dwa razy p z cyfrą x dołączoną po prawej.
     • Uwaga: Można znaleźć x biorąc liczbę p, czyli dotychczas zapisany wynik nad kreską, i zgadując, że kolejna cyfra x to c/(20·p), a następnie wykonując obliczenie y z ewentualną zmianą x w górę lub dół stosownie do porównania obliczonego y z c.
   • Umieszczamy cyfrę ${\displaystyle x}$ jako następną cyfrę wyniku pierwiastkowania.
3. Odejmujemy y od c i otrzymujemy nową resztę.
4. 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ę, ${\displaystyle e,}$ która dodana z prawej bieżącego rozwiązania ${\displaystyle r,}$ takiego, że ${\displaystyle (r+e)\cdot (r+e)\leq x,}$ gdzie ${\displaystyle x}$ jest pożądaną wartością pierwiastka. Rozwijając, otrzymujemy ${\displaystyle r\cdot r+2re+e\cdot e\leq x.}$ Bieżąca wartość ${\displaystyle r\cdot r}$ – lub, zwykle, reszta – mogą być przyrostowo aktualizowane w systemie dwójkowym, gdzie wartość ${\displaystyle e}$ jest jednobitowa, a operacje wymagane do obliczenia ${\displaystyle 2\cdot r\cdot e}$ i ${\displaystyle e\cdot e}$ można zastąpić szybszymi operacjami przesunięć bitów. W wyniku otrzymujemy prosty algorytm komputerowy^[2]:

    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;
    }

Zobacz też

• szybka odwrotność pierwiastka kwadratowego

Uwagi

Przypisy