Nagyon hosszú tizedes tört reciprokának kiszámításához hol találok algoritmust?

A Newton-Raphson algoritmus általánosan azt csinálja, hogy megoldja ezt az egyenletet:

f(x) = 0

ahol f(x) egy deriválható függvény. Az algoritmus így vezethető le:

f'(x) = Δy/Δx = (f(x_(n+1)) - f(x_n))/(x_(n+1) - x_n)

x_(n+1)-nél n+1 az x alsó indexében van.

Na most ha n+1 közelít a végtelenhez, vagyis x_(n+1) közelít a keresett x-hez, akkor f(x_(n+1)) közelíteni fog a nullához, hisz f(x)=0 a keresett x-nél. Ezért így is folytathatjuk a fenti egyenletet:

f'(x) = (0 - f(x_n))/(x_(n+1) - x_n)

Ezt átrendezve ez jön ki:

x_(n+1) = x_n - f(x_n)/f'(x)

Ez már maga a Newton-Raphson képlet. Valamilyen x_0 értékkel kell elindulni, aztán a fenti képlettel iterálni az x_1, x_2, stb közelítéseket.

Normál esetben (szóval jól viselkedő függvényeknél) minden egyes iteráció duplázza az értékes jegyek számát, tehát az algoritmus gyorsan konvergál. A kezdőérték megválasztása azért hat a konvergenciára... Ezekbe a részletekbe most ne menjünk bele.

A reciprok számolásához ezt a függvényt lehet választani:

f(x) = 1/x - b

Ennek a függvénynek a gyöke (vagyis az f(x)=0 egyenlet megoldása) x = 1/b. Szóval ezt a függvényt választva éppen 1/b lesz az iteráció "végén" x értéke.

Hogyan megy ezzel az iteráció? A függvény deriváltja f'(x) = -1/x², tehát:

x_(n+1) = x_n - (1/x_n - b)/(-1/x_n²)

Beszorozva a nevező reciprokával ez jön ki:

x_(n+1) = x_n - (-x_n + b·x_n²)

x_(n+1) = 2·x_n - b·x_n²

x_(n+1) = x_n·(2 - b·x_n)

Ez a gyors reciprok algoritmus iterációjának a képlete.

Az érdekesség az, hogy most már csak egy kivonás meg két szorzás maradt az iterációs lépéseknél, szorozni meg gondolom van olyan algoritmusod, ami tud milliós tizedesjegyű számokat is gyorsan. (Bár gyors szorzó algoritmus sincs csak úgy...)

Kezdőérték:

Ha mondjuk kettes számrendszerben számolsz (gondolom úgy lesz), akkor a számot először 1/2 és 1 közé kell normálni. Vagyis shift-elgetni kell addig, amíg olyan nem lesz, és számolni, hogy hányszor shift-eltél. A végén ugyanannyival kell majd a reciprokot is shift-elni. Ugyanis k-szor jobbra shift-elve b reciproka helyett c=b/2^k reciprokát számolod így ki, az meg 1/c = 2^k/b lesz, amit k-szor kell jobbra shift-elni 1/b-hez.

1/c számolásához x_0 értékét így érdemes választani:

x_0 = 48/17 - 32/17·c

Ekkor a kezdő hiba:

c=1 esetén x_0 = 16/17 lesz, ami 1/17-del kisebb a pontos értéknél

c=1/2 esetén x_0 = 32/17 lesz, ami 2/17-del kisebb a pontos értéknél.

A hiba alatt ε=c·x-1 értékét szokták érteni, az mindkét esetben 1/17 (c·x pontosan 1 kellene legyen).

Más értékeknél is ε_0 ≤ 1/17.

Mivel a hiba minden iterációban négyzetelődik (ε_(n+1) = -ε_n²), ezért L lépés után

ε_L ≤ (1/17)^(2^L)

Ennyi iterációs lépés kell, hogy N biten pontos legyen az eredmény:

(1/17)^(2^L) ≈ (1/2)^(N+1)

17^(2^L) ≈ 2^(N+1)

(2^L)·log_2(17) ≈ N+1

L = [log_2((N+1)/log_2(17)]

Ahol [x] az x érteke felfelé kerekítve.