Jak komputer liczy pierwiastek kwadratowy?

01 Historia

Pierwiastek kwadratowy liczono na długo przed powstaniem matematyki symbolicznej — metody, które dziś programujemy, mają tysiące lat.

• ok. 1800–1600 p.n.e. — Babilonia Gliniana tabliczka YBC 7289 zawiera rysunek kwadratu z przekątną i zapisane w systemie sześćdziesiątkowym przybliżenie 1;24,51,10, czyli 1,41421296… — wartość \(\sqrt{2}\) z błędem rzędu \(10^{-6}\). To najstarszy znany ślad procedury iteracyjnego uśredniania, która dziś nosi nazwę „metody babilońskiej”.
• ok. 1650 p.n.e. — Egipt Papirus Rhinda zawiera zadania wymagające przybliżania pierwiastków przy pomocy ułamków jednostkowych — inne, mniej systematyczne podejście do tego samego problemu.
• I w. n.e. — Aleksandria Heron z Aleksandrii w dziele Metrica opisuje wprost algorytm: wziąć dowolne przybliżenie, uśrednić je z ilorazem liczby i tego przybliżenia, i powtarzać. Stąd druga nazwa metody — metoda Herona.
• III–VII w. — Indie Rękopis Bakhshali zawiera analogiczny wzór iteracyjny do przybliżania pierwiastków, rozwijany niezależnie w tradycji indyjskiej.
• Starożytne Chiny „Dziewięć rozdziałów sztuki matematycznej” opisuje algorytm ekstrakcji pierwiastka „cyfra po cyfrze”, wykonywany ręcznie na liczydle — metoda ucząca się do dziś na lekcjach jako „pisemne obliczanie pierwiastka”.
• 1669–1690 — Newton i Raphson Ogólna metoda znajdowania miejsc zerowych funkcji różniczkowalnych. Jak pokazujemy w sekcji Podstawy matematyczne, metoda babilońska jest jej szczególnym przypadkiem — tylko o dwa tysiące lat starszym.
• 1999 — „Fast inverse square root” W silniku gry Quake III Arena pojawia się słynna sztuczka bitowa z magiczną stałą 0x5f3759df, przybliżająca \(1/\sqrt{x}\) bez dzielenia — jedna iteracja metody Newtona wystarczała, by dociągnąć wynik do użytecznej dokładności (patrz Inne metody).
• Dziś Współczesne procesory mają sprzętową instrukcję pierwiastka (np. SQRTSD na x86, FSQRT na ARM), zaimplementowaną w krzemie algorytmami cyfra-po-cyfrze lub Goldschmidta — potomkami tych samych idei co metoda babilońska.

02 Podstawy matematyczne

Skąd właściwie bierze się wzór \(\tfrac12(x + S/x)\)? Wyprowadźmy go z ogólnej metody Newtona i sprawdźmy, dlaczego zbiega tak szybko.

Definicja

Pierwiastkiem kwadratowym z liczby nieujemnej \(S\) nazywamy liczbę \(x \ge 0\) spełniającą równanie \(x^2 = S\). Zagadnienie obliczeniowe sprowadza się więc do znalezienia miejsca zerowego funkcji:

Wyprowadzenie z metody Newtona

Ogólny wzór iteracyjny Newtona–Raphsona dla równania \(f(x) = 0\) ma postać:

Dla \(f(x) = x^2 - S\) pochodna wynosi \(f'(x) = 2x\), zatem:

Metoda babilońska to więc dokładnie metoda Newtona zastosowana do równania \(x^2 = S\) — z tą różnicą, że Babilończycy odkryli ją intuicyjnie, niemal 3000 lat przed Newtonem.

Interpretacja geometryczna: uśrednianie

Jeżeli \(x\) jest przybliżeniem \(\sqrt{S}\) z nadmiarem, to \(S/x\) jest przybliżeniem z niedomiarem (i odwrotnie) — bo ich iloczyn zawsze wynosi dokładnie \(S\). Nowe przybliżenie to po prostu średnia arytmetyczna tych dwóch liczb, która zawsze leży bliżej prawdziwego pierwiastka niż gorsze z dwóch wyjściowych przybliżeń. Z nierówności między średnią arytmetyczną a geometryczną (AM–GM) wynika, że dla \(n \ge 1\) zawsze \(x_n \ge \sqrt{S}\), a ciąg \((x_n)\) jest malejący i zbieżny do \(\sqrt{S}\).

Dlaczego zbiega tak szybko

Niech \(e_n = x_n - \sqrt{S}\) oznacza błąd po \(n\)-tej iteracji. Podstawiając do wzoru iteracyjnego i upraszczając, otrzymujemy:

Błąd w kolejnym kroku jest w przybliżeniu proporcjonalny do kwadratu poprzedniego błędu — to tzw. zbieżność kwadratowa. W praktyce oznacza to, że liczba poprawnych cyfr znaczących podwaja się z każdą iteracją. Poniższa tabela (dla \(S = 10\), \(x_0 = 10\)) pokazuje to zjawisko: już po 6 iteracjach osiągana jest pełna precyzja liczby zmiennoprzecinkowej typu double (ok. 15–16 cyfr znaczących).

03 Metoda babilońska

Praktyczna strona algorytmu: jak go zatrzymać, ile kosztuje, i jak wybrać przybliżenie startowe.

Wypróbuj sam — kalkulator na żywo

Wpisz dowolną liczbę nieujemną — poniższa tabela pokaże, krok po kroku, jak metoda babilońska (uruchomiona bezpośrednio w Twojej przeglądarce, w JavaScript) dochodzi do wyniku, i porówna go z wbudowaną funkcją Math.sqrt().

04 Inne metody (dla chętnych)

Metoda babilońska nie jest jedynym sposobem — oto krótki przegląd alternatyw, historycznych i współczesnych.

Wszystkie te metody łączy jedno: sprowadzają obliczenie pierwiastka do prostszych operacji (dodawanie, dzielenie, przesunięcia bitowe) wykonywanych iteracyjnie — różnią się jedynie tempem zbieżności i kosztem pojedynczego kroku.

05 Realizacja metody babilońskiej

Ten sam algorytm w czterech wcieleniach: pseudokod, C++, Python i arkusz kalkulacyjny.

Pseudokod

funkcja pierwiastek_babilonski(S, eps = 1e-15):
    jeśli S < 0:  zwróć błąd
    jeśli S = 0:  zwróć 0

    x ← S                      // przybliżenie startowe
    powtarzaj:
        poprzednie ← x
        x ← 0.5 * (x + S / x)  // krok Herona
    dopóki |x - poprzednie| > eps

    zwróć x

C++

double pierwiastekBabilonski(double S, double eps = 1e-15) {
    if (S < 0) return NAN;
    if (S == 0) return 0.0;

    double x = S;                 // przyblizenie poczatkowe
    double poprzednie;

    do {
        poprzednie = x;
        x = 0.5 * (x + S / x);    // wzor iteracyjny Herona
    } while (std::fabs(x - poprzednie) > eps);

    return x;
}

Przykładowy wynik uruchomienia (pełny kod w code/cpp/babilonska.cpp):

sqrt(2.000000000000000)
  metoda babilonska : 1.414213562373095
  std::sqrt          : 1.414213562373095
  roznica            : 0.000000000000000

sqrt(999999.000000000000000)
  metoda babilonska : 999.999499999875070
  std::sqrt          : 999.999499999874956
  roznica            : 0.000000000000114

Python

def pierwiastek_babilonski(s: float, eps: float = 1e-15) -> float:
    if s < 0:
        raise ValueError("Liczba pod pierwiastkiem nie moze byc ujemna")
    if s == 0:
        return 0.0

    x = s
    while True:
        poprzednie = x
        x = 0.5 * (x + s / x)
        if abs(x - poprzednie) <= eps:
            return x

Przykładowy wynik uruchomienia (pełny kod w code/python/babilonska.py):

sqrt(2)
  metoda babilonska : 1.414213562373095
  math.sqrt          : 1.414213562373095
  roznica            : 0.000000000000000

sqrt(999999)
  metoda babilonska : 999.999499999875070
  math.sqrt          : 999.999499999874956
  roznica            : 0.000000000000114

Arkusz kalkulacyjny (Excel)

Plik code/excel/metoda_babilonska.xlsx zawiera dwa arkusze. Pierwszy rozpisuje kolejne iteracje wiersz po wierszu — komórka S jest parametrem wejściowym, a każdy wiersz odwołuje się formułą do wiersza poprzedniego:

Drugi arkusz, „Porównanie kilku liczb”, zestawia dla kilku wartości S wynik metody babilońskiej (jako jedna, 24-krotnie zagnieżdżona formuła Herona) z wbudowaną funkcją =SQRT(S) oraz różnicę bezwzględną między nimi — obie kolumny liczone są w pełni formułami, bez wpisywania gotowych wyników.

06 Eksperyment

Czy „ręczna” metoda babilońska jest w praktyce wolniejsza od funkcji bibliotecznej — i o ile?

Metodologia

• Zakres danych: liczby całkowite 1…1 000 000.
• Dokładność: \(\varepsilon = 10^{-15}\) (15 miejsc po przecinku) dla metody babilońskiej.
• Mierzona wielkość: czas wykonania całej serii 1 000 000 obliczeń (nie pojedynczego wywołania).
• Powtórzenia: 5 pomiarów każdej metody, w jednym uruchomieniu programu (aby wyeliminować narzut startu procesu).
• Język / kompilacja: C++17, g++ -O2, pomiar std::chrono::high_resolution_clock.
• Środowisko: Apple M4, macOS 15 (Darwin 24.6), Apple Clang 17.
• Kontrola poprawności: obie metody sumują wyniki do zmiennej pomocniczej — suma kontrolna obu metod zgadza się do wyświetlanej precyzji we wszystkich próbach, co potwierdza, że porównujemy poprawne, równoważne wyniki (kod: code/cpp/eksperyment.cpp).

07 Wyniki

Pięć powtórzeń serii 1 000 000 obliczeń, C++, wynik w sekundach.

Dodatkowo: to samo w Pythonie

Dla porównania między językami przeprowadzono identyczny eksperyment w Pythonie (code/python/eksperyment.py) — nieobowiązkowe rozszerzenie, poza wymaganym zakresem zadania, ale dobrze ilustrujące wpływ narzutu interpretera na algorytmy iteracyjne.

08 Wnioski

1. Metoda babilońska jest poprawna — dla całego przebadanego zakresu (1…1 000 000) wyniki obu metod zgadzają się do 13.–15. miejsca po przecinku; jedyne różnice wynikają z akumulacji błędów zaokrągleń liczb zmiennoprzecinkowych, a nie z błędu algorytmu.
2. Jest jednak wyraźnie wolniejsza — w C++ około 11–12 razy, a w Pythonie prawie 19 razy wolniejsza od wbudowanej funkcji sqrt(). Różnica ta wynika z tego, że biblioteczny sqrt() korzysta bezpośrednio ze sprzętowej instrukcji procesora (jeden takt zegara), podczas gdy metoda babilońska wykonuje kilka pełnych iteracji z dzieleniem zmiennoprzecinkowym za każdym razem.
3. Narzut interpretera potęguje różnicę — w Pythonie krotność spowolnienia jest wyższa niż w C++ (≈19× wobec ≈11×), ponieważ każda iteracja pętli while w czystym Pythonie kosztuje dodatkowo czas na interpretację bajtkodu — koszt ten mnoży się przez liczbę iteracji potrzebnych do zbieżności.
4. Zbieżność kwadratowa działa zgodnie z teorią — mimo że algorytm wykonuje kilka pełnych „przebiegów” (nie jedną operację), liczba iteracji potrzebna do pełnej precyzji double pozostaje mała (6–7) niezależnie od wielkości liczby w badanym zakresie, co potwierdza wyprowadzenie z sekcji Podstawy matematyczne.
5. Praktyczny wniosek: w kodzie produkcyjnym należy zawsze używać funkcji bibliotecznej sqrt() — jest szybsza, dokładniejsza w skrajnych przypadkach (nieskończoności, liczby ujemne, podprzepełnienia) i przetestowana przez dekady użycia. Wartość metody babilońskiej jest dziś przede wszystkim edukacyjna — uczy, jak „pod maską” działają metody numeryczne, na których opierają się też znacznie bardziej zaawansowane algorytmy (np. metoda Newtona w uczeniu maszynowym czy optymalizacji).