Algorytm precyzji pierwiastka kwadratowego?

Czy są znane algorytmy subkwadratowe do obliczania wartości pierwiastka kwadratowego z nliczby całkowitej?

Naiwny algorytm byłby podobny

def sqrt(x):
    r = 0
    i = x.bit_length() // 2
    while i >= 0:
        inc = (r << (i+1)) + (1 << (i*2))
        if inc <= x:
            x -= inc
            r += 1 << i
        i -= 1
    return r

Wymaga to O(n)iteracji, z których każda wiąże się z dodatkami, które są O(n)czasem, więc O(n^2)czas ogólnie. Czy jest coś szybciej? Wiem, że w przypadku mnożenia istnieją specjalne algorytmy, które działają lepiej niż czas kwadratowy, ale nie mogę znaleźć niczego dla pierwiastków kwadratowych.

Możesz użyć metody Newtona lub dowolnej z wielu innych metod do znalezienia aproksymacji do pierwiastków wielomianu .$p(x) = x^2 -c$

Szybkość zbieżności dla metody Newtona będzie kwadratowa, co oznacza, że ​​liczba poprawnych bitów podwaja się w każdej iteracji. Oznacza to, że wystarczą iteracje metody Newtona.$O(\lg n)$

Każda iteracja metody Newtona jest obliczana

$$x_{j+1} = x_j - (x_j^2 -c)/(2x_j) = 0.5 x_j + \frac{c}{2x_j}.$$

Złożoność bitowa mnożenia wynosi , aby pomnożyć dwa$\stackrel{~}{O}(b \lg b)$$b$-bit liczby całkowite (ignorowanie $\lg \lg b$czynniki). Złożoność bitowa dla podziału (na$b$bitów precyzji) jest taki sam. Dlatego każdą iterację można obliczyć w$\stackrel{~}{O}(n \lg n)$operacje. Mnożenie przez$O(\lg n)$ iteracje, okazuje się, że całkowity czas działania do obliczenia pierwiastka kwadratowego $n$ precyzja jest $\stackrel{~}{O}(n (\lg n)^2)$. To jest subkwadratowe.

Myślę, że dokładniejsza analiza pokazuje, że można to poprawić $\stackrel{~}{O}(n \lg n)$ czas działania (biorąc pod uwagę, że musimy tylko znać każdy z nich $x_j$ w ciągu około $j$ fragmenty precyzji, a nie $n$fragmenty precyzji). Jednak nawet bardziej podstawowa analiza już pokazuje czas działania, który jest wyraźnie subkwadratowy.

Jednym z problemów związanych z metodą Newtona jest to, że wymaga ona operacji dzielenia w każdej iteracji, która jest najwolniejszą podstawową operacją na liczbach całkowitych.

Jednak metoda Newtona dla odwrotnego pierwiastka kwadratowego nie. Gdyby$x$ to numer, dla którego chcesz znaleźć $\frac{1}{\sqrt x}$, iteruj:

$$r_{i+1} = \frac{1}{2} r_i (3 - x r_i^2)$$

Jest to często wyrażane jako:

$$w_i = r_i^2$$ $$d_i = 1 - w_i x$$ $$r_{i+1} = r_i + \frac{r_i d_i}{2}$$

To trzy operacje mnożenia. Podział na dwa można zastosować jako przesunięcie w prawo.

Problem w tym, że $r$nie jest liczbą całkowitą. Można jednak nim manipulować, wprowadzając ręcznie zmiennoprzecinkowe i wykonując kilka operacji przesunięcia, aby skompensować w razie potrzeby.

Najpierw przeskalujmy $x$:

$$x' = 2^{-2e} x$$

gdzie chcielibyśmy $x'$ być większym, ale blisko, $1$. Jeśli uruchomimy powyższy algorytm$x'$ zamiast $x$, znaleźliśmy $r = \frac{1}{\sqrt x'}$. Następnie,$\sqrt{x} = 2^e r x'$.

Teraz podzielmy się $r$ w mantysę i wykładnik potęgi:

$$r_i = 2^{-e_i} r'_i$$

gdzie $r'_i$jest liczbą całkowitą. Intuicyjnie,$e_i$ reprezentują precyzję odpowiedzi.

Wiemy, że metoda Newtona z grubsza podwaja liczbę dokładnych cyfr znaczących. Możemy więc wybrać:

$$e_{i+1} = 2e_i$$

Przy odrobinie manipulacji znajdujemy:

$$e_{i+1} = 2e_i$$ $$w_i = {r'_i}^2$$ $$x'_i = \frac{x}{2^{2e - e_{i+1}}}$$ $$d_i = 2^{e_{i+1}} - \frac{w_i' x'_i}{2^{e_{i+1}}}$$ $$r'_{i+1} = 2^{e_i} r'_i - \frac{r'_i d_i}{2^{e_i + 1}}$$

Przy każdej iteracji:

$$\sqrt{x} \approx \frac{r'_i x}{2^{e + e_i}}$$

Jako przykład spróbujmy obliczyć pierwiastek kwadratowy z $x = 2^{63}$. Wiemy, że odpowiedź brzmi$2^{31}\sqrt{2}$. Odwrotnym pierwiastkiem kwadratowym jest$\frac{1}{\sqrt{2}} 2^{-31}$, więc ustawimy $e = 31$ (to jest skala problemu) i na nasze początkowe przypuszczenia wybierzemy $r'_0 = 3$ i $e_0 = 2$. (To znaczy, wybieramy$\frac{3}{4}$ dla naszego wstępnego oszacowania na $\frac{1}{\sqrt{2}}$.)

Następnie:

$$e_1 = 4, r'_1 = 11$$ $$e_2 = 8, r'_2 = 180$$ $$e_3 = 16, r'_3 = 46338$$ $$e_4 = 32, r'_4 = 3037000481$$

Możemy ustalić, kiedy przerwać iterację, porównując $e_i$ do $e$; jeśli poprawnie obliczyłem,$e_i > 2e$powinno być wystarczająco dobre. Zatrzymamy się tutaj i znajdziemy:

$$\sqrt{2^{63}} \approx \frac{3037000481 \times 2^{63}}{2^{31+32}} = 3037000481$$

Prawidłowy pierwiastek kwadratowy z liczby całkowitej to $3037000499$, więc jesteśmy bardzo blisko. Możemy wykonać kolejną iterację lub zoptymalizować iterację końcową, która się nie podwaja$e_i$. Szczegóły pozostawia się jako ćwiczenie.

Aby przeanalizować złożoność tej metody, zwróć uwagę, że pomnożenie dwóch $b$-bit liczba całkowita bierze $O(b \log b)$operacje. Jednak tak to ułożyliśmy$r'_i < 2^{e_i}$. Więc mnożenie do obliczenia$w_i$ mnoży dwa $e_i$-bitowe liczby, aby wygenerować $e_{i+1}$-bit liczba, a pozostałe dwa mnożenia mnożą dwa $e_{i+1}$-bitowe liczby, aby wygenerować $2e_{i+1}$-bitowa liczba.

W każdym przypadku liczba operacji na iterację wynosi $O(e_i \log e_i)$, i tu są $O(\log e)$wymagane iteracje. Ostateczne pomnożenie jest rzędu$O(2e \log 2e)$operacje. Ogólna złożoność jest$O(e \log^2 e)$ operacje, które są podkwadratowe pod względem liczby bitów w $x$. To zaznacza wszystkie pola.

W tej analizie kryje się jednak ważna zasada, o której powinni pamiętać wszyscy pracujący z dużymi liczbami całkowitymi: ponieważ mnożenie jest liczbą superliniową pod względem liczby bitów, wszelkie operacje mnożenia należy wykonywać tylko na liczbach całkowitych o przybliżonej wielkości bieżącej precyzji (i , Mógłbym dodać, powinieneś spróbować pomnożyć razem liczby, które mają podobny rząd wielkości). Używanie liczb całkowitych większych niż to strata wysiłku. Czynniki stałe są ważne, a dla dużych liczb całkowitych mają one duże znaczenie.

W końcowej obserwacji dwa z mnożenia mają formę $\frac{ab}{2^c}$. Najwyraźniej nie ma sensu obliczać wszystkich bitów$ab$ tylko rzucić $c$z nich z przesunięciem w prawo. Wdrożenie inteligentnej metody mnożenia, która bierze to pod uwagę, również pozostawia się jako ćwiczenie.