Newton-Raphson-Division

Das Newton–Raphson-Divisions-Verfahren benutzt das Newton-Verfahren, um den Kehrwert eines Nenners $N\in \mathbb {N}$ zu finden und diesen mit einem Zähler $Z\in \mathbb {Z}$ zu multiplizieren für das Ergebnis des Quotienten $Q:=Z/N$.

Wegen der besonderen Bedeutung für die Computertechnik wird das Verfahren im Folgenden für das Dualsystem vorgestellt. Es lässt sich aber auch bei jeder anderen Basis anwenden.

Schritte

Die Schritte des Newton–Raphson-Divisionsverfahrens sind:

1. Finden einer ersten Näherung $X_{0}$ für den Kehrwert $1/N$ des Nenners $N$.
2. Berechnen immer besserer Näherungen $X_{1},X_{2},\ldots ,X_{S}$ des Kehrwerts. Hier wird vom Newton-Verfahren Gebrauch gemacht.
3. Berechnen des Quotienten durch Multiplikation des Zählers mit dem Kehrwert des Nenners: $Q=ZX_{S}$.

Newton-Verfahren

Die Anwendung des Newton-Verfahrens benötigt eine Funktion $f(X)$, die eine Nullstelle bei $X=1/N$ hat. Die naheliegende Funktion $f(X)=NX-1$ hat triviale für das Newton-Verfahren untaugliche Ableitungen. Eine brauchbare Funktion ist $f(X)=1/X-N$ mit $f'(X)=-1/X^{2}$. Wegen $f'(X)<0$ schneidet der Graph der Funktion die $X$-Achse transversal, d. h. nicht-berührend. Für die Newton–Iteration ist

$X_{i+1}=X_{i}-{f(X_{i}) \over f'(X_{i})}=X_{i}-{1/X_{i}-N \over -1/X_{i}^{2}}=X_{i}+X_{i}(1-NX_{i})=X_{i}(2-NX_{i})$,

was ausgehend von $X_{i}$ ausschließlich über Multiplikation und Subtraktion (oder zwei fused multiply-add-Operationen) berechnet werden kann.

Konvergenzgeschwindigkeit

Wenn der Fehler als $\epsilon _{i}=NX_{i}-1$ definiert ist, dann ist:

${\begin{aligned}\epsilon _{i+1}&=NX_{i+1}-1\\&=N(X_{i}(2-NX_{i}))-1\\&=2NX_{i}-N^{2}X_{i}^{2}-1\\&=-(N^{2}X_{i}^{2}-2NX_{i}+1)\\&=-(NX_{i}-1)^{2}\\&=-{\epsilon _{i}}^{2}.\\\end{aligned}}$

Diese Quadrierung des Fehlers bei jedem Schritt – die sog. quadratische Konvergenz des Newton-Verfahrens – sorgt dafür, dass die Anzahl der korrekten Ziffern sich bei jeder Iteration in etwa verdoppelt; eine Eigenschaft die beim Rechnen mit Langzahlen besonders wertvoll ist.

Da für diese Methode die Konvergenzgeschwindigkeit exakt quadratisch ist, folgt, dass

$S=\left\lceil \log _{2}{\frac {P+1}{\log _{2}17}}\right\rceil \,$

Schritte für eine Genauigkeit von $P$ Binärstellen ausreichen. Das sind 3 für single und 4 für sowohl double wie extended precision IEEE 754 Formate.

Finden einer ersten Näherung

Durch Bitverschiebungen kann der Nenner $N$ ins Intervall $0{,}5\leq N<1$ gebracht werden. Dieselbe Anzahl Shifts sollte der Zähler $Z$ erfahren, so dass der Quotient ungeändert bleibt. Danach kann man eine lineare Approximation der Form

${\frac {1}{N}}:\approx X_{0}:=T_{1}+T_{2}N$   mit   $T_{1}:={\frac {48}{17}}$   und   $T_{2}:=-{\frac {32}{17}}$

anwenden, um das Verfahren zu initialisieren.

Die Koeffizienten $T_{1}$ und $T_{2}$ dieser linearen (Polynomgrad 1) Approximation ergeben sich folgendermaßen. Der relative Fehler ist $|(T_{1}+T_{2}N-1/N)/(1/N)|=|N(T_{1}+T_{2}N)-1|$. Der Minimalwert des maximalen solchen Fehlers im Intervall $[0{,}5\,;\,1[$ wird gegeben durch den Alternantensatz von Tschebyscheff angewendet auf die Funktion $F(N)=N(T_{1}+T_{2}N)-1$. Das lokale Extremum von $F(N)$ findet statt, wenn $F'(N)=0$ ist, was die Lösung $N=-T_{1}/(2T_{2})$ hat. Nach dem genannten Alternantensatz muss diese Funktion am Extremum (im Inneren) das umgekehrte Vorzeichen als an den Rändern des Intervalls haben, also $F(1/2)=F(1)=-F(-T_{1}/(2T_{2}))$. Für die zwei Unbekannten in den zwei Gleichungen ergibt sich die Lösung $T_{1}=48/17$ und $T_{2}=-32/17$, und der maximale relative Fehler ist $F(1)=1/17$. Nach dieser Approximation ist der relative Fehler des Anfangswertes

$\vert \epsilon _{0}\vert \leq {1 \over 17}\approx 0{,}059.$

Pseudocode

Das Folgende berechnet den Quotienten von $Z$ und $N$ mit einer Genauigkeit von $P$ Binärstellen:

Drücke N aus als M × 2e mit 1 ≤ M < 2 (Standard-Gleitkomma-Darstellung)
N' := N / 2e+1             // Bitverschiebungen resp. Verkleinerung des Exponenten
Z' := Z / 2e+1
X := 48/17 − 32/17 × N'   // erste Näherung mit der gleichen Genauigkeit wie N
repeat $\left\lceil \log _{2}{\frac {P+1}{\log _{2}17}}\right\rceil \,$ times   // kann für fixes P vorausberechnet werden
    X := X + X × (1 - N' × X)
end
return Z' × X

Diese Methode benötigt bspw. für eine double-precision Gleitkomma-Division 10 Multiplikationen, 9 Additionen und 2 Shifts.

Literatur

• Mário P. Véstias and Horácio C. Neto Decimal Division Using the Newton–Raphson Method and Radix-1000 Arithmetic
• Thomas L. Rodeheffer Software Integer Division

Ähnliche Verfahren

• Goldschmidt-Division
• SRT-Division