Algorytm we wcześniejszym poście wymagał znalezienia pierwiastka kwadratowego z liczby, liczonego dwukrotnie podczas inicjacji.
Publikowałem już sposób na liczenie pierwiastka rok temu, w listopadzie 2012 roku, lecz teraz znalazłem nieco szybszą technikę.
Ustalmy liczbę n oraz drugą liczbę p położoną między pierwiastkiem sześciennym oraz kwadratowym z n. Nic więcej o liczbie p nie potrzeba wiedzieć, jest ona pierwszym przybliżeniem pierwiastka. Resztę załatwią konwersje, a być może także wzór skróconego mnożenia.
Dzielimy n z resztą przez p (dwukrotnie), by uzyskać postać
n = (a*p+b)*p+c (1)
Jest to zapis liczby trójcyfrowej.
Możemy sprawdzić, czy a, b, c stanowią współczynniki trójmianu kwadratowego d^2+2de+e^2+f, jeśli tak, pierwiastek można przybliżyć przez
(dp+e)^2+f.
Jeśli jednak b>2de, też możemy posłużyć się wzorem, do którego należy dodać poprawkę (b-2de)*p. Sposób też prowadzi do celu, lecz jest znacznie dłuższy. Zbieżność do pierwiastka jest wolniejsza.
Zatem pierwsza faza, sprowadzamy a do postaci podzielnej przez 4, b do liczby parzystej za pomocą formuł:
a-1 AND b+=p; b-1 AND c+=p. (2)
Np. mając (a,b,c) = (5, 4, 34) dla p=125 uzyskujemy z pomocą formuł (2)
(a,b,c) = (5-1, 4+125, 34) = (4, 129, 34) = (4, 129-1, 34+125) = (4, 128, 159)
Korzystamy z konwersji podwajania systemu
(a,b,c)_p = (a/4, b/2, c)_(2p) (3)
czyli w przykładzie mamy (a,b,c) = (1, 64, 159)_250 .
Kiedy już a<4, oraz nie dopasujemy wzoru skróconego mnożenia, dążymy do zmniejszenia b. W odróżnieniu od zeszłorocznego algorytmu najlepiej to zrobić za pomocą konwersji o floor(a*b/3). Konwersja ta zmniejsza b o co najmniej połowę, może także zmniejszyć a. Branie połowy (a*b/2) może doprowadzić do zmniejszenie b do wartości ujemnej, zaś chcemy zachować a=1 jako niezmiennik, oraz wartość liczby trójcyfrowej powinna pozostać większa niż pierwiastek kwadratowy.
Przykład dla n=8934053
Przyjmijmy p=1000, wtedy (a,b,c) = (8, 934, 53)_1000
Pierwsza konwersja (3) sprowadza do (a,b,c) = (2, 467, 53)_2000
konwersja o floor(a*b/3) = 311
2 467 53
2 467-2*311
2-1 2156 53 // dodaję do b (2000+311) = 2311
1 2156-311 53-2156*311
1 1845-291 -670463 // dodaję do c 291*p
1 1554 2038
oraz (a,b,c) = (1, 1554, 2038)_2311
konwersja o floor(a*b/3) = 518, nowa podstawa p = 2311+518 = 2829
1 1554 2038
1 1554-518
1 1036 2038
1 518 2038-1036*518 // dodajemy do c 189*p
1 518-189 -534610
1 329 71
konwersja o floor(a*b/3) = 109, nowa podstawa p = 2829+109 = 2938
1 329 71
1 329-109
1 220 71
1 220-109 71-220*109
1 111-9 -23909 // dodajamy do c 9*p
1 102 2533
Sprawdzamy wzór skróconego mnożenia (1+51)^2 = 1+102+2601 oraz porównujemy. Niezgodność na pozycji c, jest za dużo o 68. Ale jest dobrze.
Zatem (a,b,c)_p = (1*p+51)^2-68 = (2938+51)^2-68 = 2989^2-68.
Zmniejszamy kwadrat
2989^2-68 = 2988^2+5976+1-68 = 2988^2+5909.
Znaleziony został pierwiastek całkowitoliczbowy 2988 oraz reszty: 5909 z niedomiarem, 68 z nadmiarem.
Dokładniejszą wartość pierwiastka mozna policzyć za pomocą interpolacji 2988 + 5909/(5909+68) = 2988,988623; dokładniejsza wartość 2988,9886249365352845723814
Błąd dopiero na szóstej pozycji po przecinku.
Brak komentarzy:
Prześlij komentarz