Newton-Raphson-Division

Das Newton–Raphson-Divisions-Verfahren benutzt das Newton-Verfahren, um den Kehrwert eines Nenners ${\displaystyle N\in \mathbb {N} }$ zu finden und diesen mit einem Zähler ${\displaystyle Z\in \mathbb {Z} }$ zu multiplizieren für das Ergebnis des Quotienten ${\displaystyle 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.

Die Schritte des Newton–Raphson-Divisionsverfahrens sind:

1. Finden einer ersten Näherung ${\displaystyle X_{0}}$ für den Kehrwert ${\displaystyle 1/N}$ des Nenners ${\displaystyle N}$.
2. Berechnen immer besserer Näherungen ${\displaystyle 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: ${\displaystyle Q=ZX_{S}}$.

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

${\displaystyle 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 ${\displaystyle X_{i}}$ ausschließlich über Multiplikation und Subtraktion (oder zwei fused multiply-add-Operationen) berechnet werden kann.

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

${\displaystyle {\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

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

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

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

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

anwenden, um das Verfahren zu initialisieren.

Die Koeffizienten ${\displaystyle T_{1}}$ und ${\displaystyle T_{2}}$ dieser linearen (Polynomgrad 1) Approximation ergeben sich folgendermaßen. Der relative Fehler ist ${\displaystyle |(T_{1}+T_{2}N-1/N)/(1/N)|=|N(T_{1}+T_{2}N)-1|}$. Der Minimalwert des maximalen solchen Fehlers im Intervall ${\displaystyle [0{,}5\,;\,1[}$ wird gegeben durch den Alternantensatz von Tschebyscheff angewendet auf die Funktion ${\displaystyle F(N)=N(T_{1}+T_{2}N)-1}$. Das lokale Extremum von ${\displaystyle F(N)}$ findet statt, wenn ${\displaystyle F'(N)=0}$ ist, was die Lösung ${\displaystyle 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 ${\displaystyle F(1/2)=F(1)=-F(-T_{1}/(2T_{2}))}$. Für die zwei Unbekannten in den zwei Gleichungen ergibt sich die Lösung ${\displaystyle T_{1}=48/17}$ und ${\displaystyle T_{2}=-32/17}$, und der maximale relative Fehler ist ${\displaystyle F(1)=1/17}$. Nach dieser Approximation ist der relative Fehler des Anfangswertes

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

Das Folgende berechnet den Quotienten von ${\displaystyle Z}$ und ${\displaystyle N}$ mit einer Genauigkeit von ${\displaystyle 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 ${\displaystyle \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.

• 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

• Goldschmidt-Division
• SRT-Division