Schönhage-Strassen-Algorithmus

aus Wikipedia, der freien Enzyklopädie

Du hast neue Nachrichten (Unterschied zur vorletzten Version).

Der Schönhage-Strassen-Algorithmus ist der effizienteste bislang bekannte Algorithmus zur Multiplikation zweier n-stelliger ganzer Zahlen. Er wurde 1971 von Arnold Schönhage und Volker Strassen entwickelt. Der Algorithmus basiert auf einer „superschnellen“ Variante der diskreten schnellen Fourier-Transformation sowie einem geschickten Wechsel zwischen der Restklassen- und der zyklischen Arithmetik in endlichen Zahlenringen.

Der Schönhage-Strassen-Algorithmus terminiert in $O \Big(n \cdot \log(n) \cdot \log \big(\log(n) \big) \Big)$ (s. Landau-Notation), wenn als Effizienzmaß die Bitkomplexität auf mehrbändigen Turingmaschinen, also die maximale Laufzeit des Algorithmus gemessen als benötigte Bitoperationen in Abhängigkeit von der Bitlänge $n$ der Eingabegrößen gewählt wird. Diese Komplexität stellt eine Verbesserung sowohl gegenüber dem „naiven“ aus der Schule bekannten Algorithmus der Laufzeit $O \left(n^2 \right)$ als auch gegenüber dem 1962 entwickelten Karatsuba-Algorithmus mit einer Laufzeit von $O \left(n^{\log_2 (3)} \right)$ sowie dessen verbesserter Variante, dem Toom-Cook-Algorithmus mit $O(n^{1+\varepsilon})$ Laufzeit dar.

[Bearbeiten] Bedeutung

Bis heute konnte kein effizienterer Algorithmus gefunden werden. Als untere Schranke gibt es für den allgemeinen Fall nur die (triviale) lineare Laufzeit, an die sich der Algorithmus mit wachsender Zahlenlänge annähert. Allerdings haben die Forscher Hinweise dafür gefunden, dass die Schranke $O(n \cdot \log(n))$ niemals unterboten werden kann. Selbst bei modernen Computern ist diese Methode der Berechnung erst bei Zahlen mit mehreren tausend Stellen effizienter als der Karatsuba-Algorithmus. Dies liegt wohl allerdings weniger am Overhead des Schönhage-Strassen-Algorithmus, sondern vielmehr an der seit Jahrzehnten typischen Designoptimierung der Computerprozessoren, die dem Erreichen schneller Gleitkommaoperationen den Vorzug vor der Arithmetik in endlichen Restklassenringen ganzer Zahlen gibt.

Für die Suche nach den Algorithmen mit der besten (Zeit-) Komplexität in der Computer-Algebra genießt der Schönhage-Strassen-Algorithmus zentrale Bedeutung.

[Bearbeiten] Algorithmus

[Bearbeiten] Grundidee und Terminologie

Um zwei ganze Zahlen $a$ und $b$ zu multiplizieren, wird im Groben folgendes Schema angewandt:

Aufspaltung der Zahlen (in Binärdarstellung) $a$ und $b$ in Stücke passender Länge
Schnelle diskrete Fourier-Transformation der beiden Stückfolgen
Komponentenweise Multiplikation der transformierten Stücke
Rücktransformation (inverse Fouriertransformation) der Ergebnisse
Zusammensetzen der Ergebnisstücke zur Ergebniszahl

Die im mittleren Schritt durchzuführenden „kleinen“ Multiplikation werden im rekursiven Sinne wiederum durch den Schönhage-Strassen-Algorithmus ausgeführt.

Bei der Umsetzung dieser Strategie sind allerdings viele „technische“ Schwierigkeiten zu meistern, unter anderem:

Die Multiplikation ist keine reine Faltung, sondern es kann auch zu Überträgen kommen; nach Durchführen der FT und iFT müssen diese passend behandelt werden.
Eine gewöhnliche FFT wäre nicht schnell genug, um die zu erzielende Komplexitätsschranke zu erreichen. Hier muss die schnellste aller nur irgend denkbaren FT-Varianten (sozusagen eine „superschnelle DFT“) benutzt werden, bei der die benötigten Multiplikationen als reine Shift-Operationen durchführbar sind.

Die Aufgabe der Multiplikation zweier ganzer Zahlen wird nun wie folgt konkretisiert:

Es seien die zwei zu multiplizierenden Zahlen $a, b\in\Bbb Z$ in Binärzifferdarstellung gegeben. Weiter sei $N$ die maximale Länge (also Binärziffernanzahl) der beiden Zahlen.

Nach passender Behandlung der Vorzeichen der beiden Zahlen sowie der trivialen Sonderfälle $a = 0$ und $b = 0$ (was mit linearem Aufwand $O (N)$ machbar ist) darf man davon ausgehen, dass $a,b\in\Bbb N$ natürliche Zahlen sind. Der Schönhage-Strassen-Algorithmus löst diese Aufgabe in $O(N\cdot\log N\cdot \log\log N)$ .

[Bearbeiten] Theoretische Vorbereitungen

[Bearbeiten] „Superschnelle DFT“

Die oben angesprochene „superschnelle DFT“, die das Kernstück des Algorithmus darstellt, muss etwas ausführlicher erläutert werden, da sie hier sehr speziell eingesetzt wird.

Es sei $R$ ein kommutativer unitärer Ring. In $R$ sei das Element $2$ eine Einheit; weiterhin sei $w\in R$ eine $2 n$ te Einheitswurzel (also $w^{2^n}=1$ ), die die Gleichheit $w^{2^{n-1}} = -1$ erfüllt. Dann lässt sich die Berechnung der diskreten Fouriertransformation (DFT) im Produktraum $R^{2^n}$ (dies ist eine Kurznotation für $R^{(2^n)}$ ; der Begriff Vektorraum ist hier nur für den Fall, dass $R$ ein Körper ist, üblich) wie folgt in einer schnellen Variante (als FFT) durchführen:

Zu berechnen ist für $a=(a_0,\ldots,a_{2^n-1})\in R^{2^n}$ die Transformierte $\hat a\in R^{2^n}$ mit

$\hat a_k = \sum_{j=0}^{2^n-1} a_j\cdot w^{j\cdot k}$ für $k=0,\ldots,2^n-1$ .

Indem wir die Indizes $k=\sum_{\nu=0}^{n-1} k_{\nu}\cdot 2^{\nu}$ und $j=\sum_{\nu=0}^{n-1}j_{\nu}\cdot 2^{n-1-\nu}$ in Binärdarstellung aufschreiben, wobei wir dies bei der Zahl $j$ in umgekehrter Reihenfolge tun, ist die Transformierte $\hat a$ wie folgt optimiert berechenbar:

Es seien

$A_0(j_0,\ldots j_{n-1}) = a_j$ für $j=0,\ldots,2^n-1$

und

$A_{r+1}(k_0,\ldots,k_r,j_{r+1},\ldots,j_{n-1}) =$

$= A_r(k_0,\ldots, k_{r-1},0,j_{r+1},\ldots,j_{n-1}) +$

$A_r(k_0,\ldots, k_{r-1},1,j_{r+1},\ldots,j_{n-1})\cdot w^{2^{n-1-r}\cdot (k_r 2^r+\dots +k_0 2^0)}$

$= \sum_{j_r=0}^1 A_r(k_0,\ldots, k_{r-1},j_r,\ldots,j_{n-1}) \cdot w^{j_r 2^{n-1-r}\cdot (k_r 2^r+\dots +k_0 2^0)}$ .

Die geschlossene Darstellung für diese Zwischenterme ist

$A_{r+1}(k_0,\ldots,k_r,j_{r+1},\ldots,j_{n-1}) =$

$= \sum_{j_r=0}^1 \ldots \sum_{j_0=0}^1 a_j \cdot w^{(j_0 2^{n-1}+\ldots + j_r 2^{n-1-r}) \cdot (k_r 2^r+\dots +k_0 2^0)}$ .

(Zum Nachrechnen dieser Darstellung beachte man $w^{(j_0 2^{n-1}+\ldots + j_r 2^{n-1-r})\cdot k_r 2^r} = 1$ ).

Diese Rekursion liefert die gewünschten Fourierkoeffizienten $\hat a_k = A_n(k_0,\ldots,k_{n-1})$ .

Aufgrund der Eigenschaft $w^{2^{n-1}}=-1$ können wir den Rekursionsschritt etwas berechnungsfreundlicher umformen zu

$A_{r+1}(k_0,\ldots,k_{r-1},0,j_{r+1},\ldots,j_{n-1}) =$

$= A_r(k_0,\ldots, k_{r-1},0,j_{r+1},\ldots,j_{n-1}) + A_r(k_0,\ldots, k_{r-1},1,j_{r+1},\ldots,j_{n-1})\cdot w^x$

und

$A_{r+1}(k_0,\ldots,k_{r-1},1,j_{r+1},\ldots,j_{n-1}) =$

$= A_r(k_0,\ldots, k_{r-1},0,j_{r+1},\ldots,j_{n-1}) - A_r(k_0,\ldots, k_{r-1},1,j_{r+1},\ldots,j_{n-1})\cdot w^x$

mit dem gleichen Exponenten $x=2^{n-1-r}\cdot (k_{r-1} 2^{r-1}+\ldots + k_0 2^0)$ .

Die Umkehrtransformation, also die inverse FFT, gelingt, da wir vorausgesetzt haben, dass $2$ im Ring $R$ invertierbar ist:

$A_r(k_0,\ldots,k_{r-1},0,j_{r+1},\ldots,j_{n-1}) =$

$= 2^{-1}\big( A_{r+1}(k_0,\ldots,k_{r-1},0,j_{r+1},\ldots,j_{n-1}) + A_{r+1}(k_0,\ldots,k_{r-1},1,j_{r+1},\ldots,j_{n-1})$

sowie

$A_r(k_0,\ldots,k_{r-1},1,j_{r+1},\ldots,j_{n-1}) =$

$= 2^{-1} w^{-x}\big( A_{r+1}(k_0,\ldots,k_{r-1},0,j_{r+1},\ldots,j_{n-1}) - A_{r+1}(k_0,\ldots,k_{r-1},1,j_{r+1},\ldots,j_{n-1})$ ,

wobei wiederum $x = 2^{n-1-r}\cdot (k_{r-1} 2^{r-1}+\ldots + k_0 2^0)$ ist.

In der Anwendung im Schönhage-Strassen-Algorithmus wird tatsächlich nur eine „halbierte FFT“ benötigt; gemeint ist damit folgendes: Beginnen wir im 1. Schritt der Rekursion mit der Berechnung

$A_1(1,j_1,\ldots,j_{n-1}) = A_0(0,j_1,\ldots, j_{n-1}) - A_0(1,j_1,\ldots, j_{n-1}) = a_j - a_{j+2^{n-1}}$

nur für $k 0 = 1$ und schränken wir die weiteren Schritte der Rekursion ebenso auf $k 0 = 1$ ein, so berechnen wir gerade alle $\hat a_k$ für ungerade Werte $k$ . Will man umgekehrt aus diesen $\hat a_k$ für ungerade $k$ (das sind $2 n - 1$ Stück) lediglich die Differenzen $a_j - a_{j+2^{n-1}}$ der ursprünglichen $a j$ zurückgewinnen, so genügt auch in der Rückrichtung die „halbierte“ Rekursion.

Im Schönhage-Strassen-Algorithmus wird die geschilderte schnelle Fouriertransformation für endliche Zahlenringe $\Bbb Z_{F_n}$ mit Fermatzahlen $F_n = 2^{(2^n)}+1$ benötigt.

Hinweis zur Notation: Für den Restklassenring $\Bbb Z / k\Bbb Z$ benutzen wir hier die kürzere Schreibweise $\Bbb Z_k$ , die lediglich im Kontext der p-adischen Zahlen zu Verwechslungen führen könnte.

Als Einheitswurzel wird im Ring $\Bbb Z_{F_n}$ die Zahl $2$ (oder je nach Kontext auch eine geeignete Potenz von 2) zum Einsatz kommen. Die beim FFT-Algorithmus durchzuführenden Multiplikationen sind dann von der Form $x\cdot 2^k$ ; allerdings sind sie nicht als reine Schift-Operationen durchführbar, da das Reduzieren eines größeren Zwischenergebnisses modulo $F n$ noch nachgeschoben werden muss. Hier greift eine der brillanten Ideen von Schönhage und Strassen: Sie betten den Ring (ausgestattet mit der Restklassenarithmetik) passend in einen größeren, mit der zyklischen Arithmetik ausgestatteten Überring ein. Dieser Überring hat eine 2-Potenz als Ordnung, so dass in ihm die entsprechende Multiplikation tatsächlich als reine Schift-Operation durchführbar ist. Diesen Trick kann man in einem schönen Struktursatz über Restklassen- und zyklische Arithmetik in endlichen Zahlenringen zusammenfassen.

[Bearbeiten] Struktursatz über zyklische Arithmetik

Der Struktursatz über zyklische Arithmetik lässt sich formal wie folgt fassen:

Für eine Zweierpotenz $D = 2 n$ mit einer natürlichen Zahl $n\in\Bbb N$ gilt

$(\Bbb Z_{D+1},+,\cdot) \cong (\Bbb Z_{D^2}, \oplus,\otimes) / (D+1)\cdot \Bbb Z$ .

Hierbei bezeichnet $(\Bbb Z_{D+1},+,\cdot)$ die durch die Repäsentanten $0,\ldots,D$ darstellbaren Restklassen modulo $D + 1$ ausgestattet mit der Restklassenarithmetik, d. h. mit der Addition und Multiplikation modulo $D + 1$ . Die in diesem Restklassenring vorkommenden Zahlen können mit $n + 1$ Binärziffern dargestellt werden.

Die auf der rechten Seite vorkommende Struktur $(\Bbb Z_{D^2}, \oplus,\otimes)$ bezeichnet die Restklassen modulo der Zahl $D 2$ , die allerdings nicht mit der Restklassenarithmetik, sondern abweichend mit der zyklischen Arithmetik ausgestattet werden. Hierbei werden bei Zwischenergebnissen, die zu groß werden, Überträge aufgehoben und auf das Endergebnis additiv aufgeschlagen. Dies entspricht in Binärzifferdarstellung einer Verschiebung der „überständigen“ Binärziffern (rechtsbündig an die niedrigsten Zifferpositionen gestellt) mit nachfolgender Addition. Beispielsweise ergibt die Addition $a + 1$ mit $a = D 2 - 1$ nicht den Wert $D^2\equiv 0$ , sondern den Wert $D^2+1\equiv 1$ . Aus der so erhaltenen Zahlenstruktur mit zyklischer Arithmetik wird nun noch der Faktorring modulo $D + 1$ gebildet. Es werden also die Endergebnisse noch modulo $D + 1$ reduziert.

Damit besagt dieser Struktursatz folgendes: Das modulo-Rechnen in $\Bbb Z_{D+1}$ kann ebenso ersetzt werden durch das zyklische Rechnen im größeren Zahlenraum $\Bbb Z_{D^2}$ mit nachfolgendem Reduzieren modulo $D + 1$ .

Entscheidend für das Gelingen der in diesem Struktursatz vorgestellten Einbettung ist die Eigenschaft, dass die größte darstellbare Zahl $F$ im zyklischen Zahlenraum (hier ist dies die Zahl $D 2 - 1$ ) die Zahl $0$ aus dem Restklassenring $\Bbb Z_{D+1}$ repräsentiert. Hierfür ist die Bedingung $(D + 1) | F$ notwendig. Damit die zyklische Arithmetik aber überhaupt sinnvoll definiert werden kann, muss andererseits $F + 1$ eine Zweierpotenz sein. Zusammen ergibt sich, dass $F = D 2 - 1$ die optimale Wahl für die Größe des zyklischen Einbettungsraumes darstellt.

Der klassische Restklassenring $(\Bbb Z_{D^2},+,\cdot)$ wäre für die Einbettung dagegen nicht geeignet, denn in diesem Ring gilt $2 2 n = D 2 = 0$ , d. h. die Zahl $2$ ist in diesem Ring ein Nullteiler.

[Bearbeiten] Durchführung

Haben wir die zu multiplizierenden Zahlen $a, b$ mit $N\leq 2^m$ Binärziffern vorliegen, so führen wir je nachdem, ob $m$ gerade oder ungerade ist, unterschiedliche Rekursionsschritte aus, um die Stellenzahl in einem Einzelschritt zu logarithmieren:

[Bearbeiten] Rekursionsschritt für ungerades m

Diesen Schritt der Rückführung von $m = 2 n - 1$ auf $n$ führen wir mit der Komplexität $O\left(2^n \psi(2^n) + n\cdot 2^{2n}\right)$ durch.

Es seien $a,b\in\Z_{F_m}$ mit $m = 2 n - 1$ und der Fermatzahl $F_m = 2^{2^m}+1$ zu multiplizieren. Wir werden in diesem Schritt die Rückführung auf die Fermatzahl $F_n = 2^{2^n}+1$ vollziehen.

Für die zu den beiden Fermatzahlen gehörenden Zweierpotenzen führen wir die Abkürzungen

$D=F_n-1 = 2^{2^n}$

und

$E = F_m-1 = 2^{2^{2n-1}}=D^{2^{n-1}}$

ein. Die halbierte Stellenzahl von $D$ wird unsere Stückelungsgröße werden, d. h. wir entwickeln $a$ und $b$ nach Potenzen von $\sqrt D$ :

$a = \sum_{i=0}^{2^n} a_i\cdot \sqrt D^i\quad$ und $\quad b = \sum_{i=0}^{2^n} b_i\cdot \sqrt D^i$ ,

wobei für die „Einzelstücke“ $0 \leq a_i, b_i < \sqrt D$ gilt. In Binärdarstellung entspricht diese Zerlegung einer einfachen Gruppierung der Bitfolgen in Stücke der Länge $2 n - 1$ Bits.

Eine kleine Schwäche des Algorithmus (die allerdings der erreichten Komplexitätsschranke keinen Abbruch tut) offenbart sich jetzt. Um die „superschnelle DFT“ auf die Stückfolgen $(a_0,\ldots,a_{2^n})$ und $(b_0,\ldots,b_{2^n})$ anwenden zu können, müssen diese zur nächsten Zweierpotenzlänge mit Nullen aufgefüllt werden; die Zahlendarstellung wird also künstlich verlängert zu

$a = \sum_{i=0}^{2^{n+1}-1} a_i\cdot \sqrt D^i$

und analog für $b$ .

Vermöge des oben erwähnten Struktursatzes zur zyklischen Arithmetik wechseln wir nun vom Restklassenring $\Bbb Z_{E+1}$ über zum Quotientenraum $(\Bbb Z_{E^2},\oplus,\otimes) /(E+1)\cdot\Bbb Z$ mit der zyklischen Arithmetik.

In diesem Raum erhält die Multiplikation die Form

$a\cdot b = \sum_{k=0}^{2^{n+1}-1} c_k \sqrt D^k$

mit den Ergebniskoeffizienten

$c_k = \sum_{r,s=0\atop r+s\equiv k\mod 2^{n+1}}^{2^{n+1}-1} a_r\cdot b_s$ .

Wir können $c k < 2 n + 1 D$ nach oben abschätzen. Für die Darstellung der $c k$ haben wir die Eigenschaft $\sqrt D^{2^{n+1}} = D^{2^n} = E^2 = 1$ in unserem zyklischen „Rechenraum“ benutzt.

Nun folgt eine Umschreibung der Summenformel, damit wir uns bei der anzuwendenden FFT auf eine „halbierte FFT“ beschränken können.

Es gilt $c_{k+2^n} \sqrt D^{k+2^n} = c_{k+2^n} \sqrt D^k E = - c_{k+2^n}\sqrt D^k$ , also ist

$a\cdot b = \sum_{k=0}^{2^n-1} (c_k - c_{k+2^n}) \sqrt D^k$

mit $-2^{n+1}D < c_k - c_{k+2^n} < 2^{n+1}D$ in $\Bbb Z_{E+1}$ . Durch passende Addition können wir den Wertebereich ins Positive verschieben, es ist nämlich $0 < c_k - c_{k+2^n} + 2^{n+1}D < 2^{n+2}D$ , und mit der Definition

$z_k = \Big\lbrace {c_k - c_{k+2^n} + 2^{n+1}D, \quad 0\leq k < 2^n \atop 2^{n+1}D, \quad\quad 2^n\leq k < 2^{n+1} }$

gilt

$a\cdot b = \sum_{k=0}^{2^{n+1}-1} z_k \sqrt D^k$ .

Für die nichttrivialen $z k$ (Indizes $0$ bis $2 n - 1$ ) gilt die Abschätzung $0 < z k < 2 n + 2 D < 2 n + 2 F n$ . Da die beiden Zahlen $2 n + 2$ und $F n$ teilerfremd ist, genügt zur Bestimmung der $z k$ die Berechnung der Reste $z_k \mod 2^{n+2}$ und $z_k \mod F_n$ .

Hat man nämlich die Reste $z_k=\xi \mod F_n$ und $z_k=\eta \mod 2^{n+2}$ bestimmt, so kann man in Komplexität $O (2 n)$ wie folgt rechnen: Berechne erst $\delta = \eta - \xi \mod 2^{n+2}$ und dann $z_k = \xi + \delta\cdot (D+1) = \xi + \delta\cdot F_n$ .

[Bearbeiten] Bestimmung der Reste modulo 2ⁿ⁺²

Hier wenden wir einen für die Computeralgebra sehr typischen Trick an: Wir setzen die Stückfolgen $a i$ und $b i$ durch Einfügen genügend langer Nullsequenzen mit „Sicherheitsabständen“ so zusammen, dass nach Produktbildung die Einzelergebnisse ebenfalls noch ohne Überlappungen in Stücken aneinandergereiht sind. Es seien also $\alpha_j = a_j\mod 2^{n+2}$ und $\beta_j = b_j\mod 2^{n+2}$ in $\Bbb Z_{2^{n+2}}$ . Wir bilden nun

$u = \sum_{k=0}^{2^{n+1}-1}\alpha_k 2^{k(3n+5)}$ und $v = \sum_{k=0}^{2^{n+1}-1}\beta_k 2^{k(3n+5)}$

und haben dabei $0 \leq u, v < 2^{2^{n+1}(3n+5)}$ . Das Produkt $u\cdot v$ enthält dann in disjunkten Stücken der Bitlänge $3 n + 5$ die Summen

$\gamma_k = \sum_{r,s=0 \atop r+s=k}^{2^{n+1}-1}\alpha_r\cdot \beta_s$

mit $0\leq k < 2^{n+2}$ , denn es ist $0\leq \gamma_k < 2^{3n+5}$ . Für die Terme $c k$ unserer ursprünglichen Multiplikationsaufgabe $a\cdot b$ sehen wir

$c_k = \gamma_k + \gamma_{k+2^{n+1}}\mod 2^{n+2}$ .

Für die zu bestimmenden Reste $\eta_k = z_k \mod 2^{n+2}$ erhalten wir

$\eta_k = \gamma_k - \gamma_{k+2^n } + \gamma_{k+2\cdot 2^n} - \gamma_{k+3\cdot 2^n}$ in $\Bbb Z_{2^{n+2}}$ .

Der Komplexitätsaufwand für die Bildung aller $α j,β j$ sowie der Extraktion der $η k$ ist $O (2 2 n)$ ; die Multiplikation $u\cdot v$ „kostet“ $ψ(2 n + 1 (3 n + 5))$ , insgesamt ist dies also $O (2 2 n)$ .

[Bearbeiten] Bestimmung der Reste modulo (D+1)

Hier kommt die DFT zum Einsatz. Wir unterziehen die „Vektoren“ $(a_0,\ldots a_{2^{n+1}-1})$ und $(b_0,\ldots b_{2^{n+1}-1})$ mit $0 \leq a_k, b_k < \sqrt D$ der DFT in $R^{2^{n+1}}$ mit $R = \Bbb Z_{D+1}$ und der Zahl $2$ als $2 n + 1$ -ter Einheitswurzel. Da wir nur die Differenzen $c_k - c_{k+2^n}$ benötigen, genügt die „halbierte DFT“:

DFT zur Bestimmung der $\hat a_k$ und $\hat b_k$ nur für die ungeraden $k$ mit $0\leq k < 2^{n+1}$
$2 n$ Multiplikationen $\hat c_k = \hat a_k\cdot \hat b_k$ für alle ungeraden $k$
Inverse DFT zur Gewinnung aller Differenzen $c_k - c_{k+2^n}$ aus den $\hat c_k$ für ungerade $k$

Der Komplexitätsaufwand hierfür besteht aus $O (n 2 n)$ Schritten des Einzelaufwands $O (2 n)$ für die DFT (gesamt also $O (n 2 2 n)$ ; hinzu kommen das Aufaddieren von $2 n + 1 D$ sowie die Reduktionen modulo $(D + 1)$ für die Gewinnung der $z_k \mod (D+1)$ , was in $O (2 2 n)$ bewältigt werden kann.

[Bearbeiten] Rekursionsschritt für gerades m

Auch für diesen Schritt der Rückführung von $m = 2 n - 2$ auf $n$ wird die Komplexität $O\left(2^n \psi(2^n) + n\cdot 2^{2n}\right)$ erreicht.

Es seien $a,b\in\Z_{F_m}$ mit $m = 2 n - 2$ und der Fermatzahl $F_m = 2^{2^m}+1$ zu multiplizieren. Wir werden auch in diesem Schritt die Rückführung auf die Fermatzahl $F_n = 2^{2^n}+1$ vollziehen.

Für die zu den beiden Fermatzahlen gehörenden Zweierpotenzen führen wir analog die Abkürzungen

$D=F_n-1 = 2^{2^n}$

und

$E = F_m-1 = 2^{2^{2n-2}}=D^{2^{n-2}}$

ein. Wiederum wird die halbierte Stellenzahl von $D$ unsere Stückelungsgröße werden, d. h. wir entwickeln $a$ und $b$ nach Potenzen von $\sqrt D$ :

$a = \sum_{i=0}^{2^{n-1}} a_i\cdot \sqrt D^i\quad$ und $\quad b = \sum_{i=0}^{2^{n-1}} b_i\cdot \sqrt D^i$ ,

wobei für die „Einzelstücke“ $0 \leq a_i, b_i < \sqrt D$ gilt.

Wie oben verlängern wir die Zahlendarstellung auf Zweierpotenzlänge zu

$a = \sum_{i=0}^{2^n-1} a_i\cdot \sqrt D^i$

und analog für $b$ .

Unter abermaliger Zuhilfenahme des Struktursatzes zur zyklischen Arithmetik wechseln wir nun vom Restklassenring $\Bbb Z_{E+1}$ über zum Quotientenraum $(\Bbb Z_{E^2},\oplus,\otimes) /(E+1)\cdot\Bbb Z$ mit der zyklischen Arithmetik.

Damit können wir wieder

$a\cdot b = \sum_{k=0}^{2^n-1} c_k \sqrt D^k$

mit den Ergebniskoeffizienten

$c_k = \sum_{r,s=0\atop r+s\equiv k\mod 2^n}^{2^n-1} a_r\cdot b_s$

darstellen. Dabei können wir $c k < 2 n D$ nach oben abschätzen.

Aus $\sqrt D^{2^{n-1}} = E^2=1$ können wir wieder

$0 < c_k - c_{k+2^{n-1}} + 2^n D < 2^{n+1} D$

folgern, und mit

$z_k = \Big\lbrace {c_k - c_{k+2^{n-1}} + 2^n D, \quad 0\leq k < 2^{n-1} \atop 2^n D, \quad\quad 2^{n-1}\leq k < 2^n}$

gilt

$a\cdot b = \sum_{k=0}^{2^n-1} z_k \sqrt D^k$

mit $0 < z k < 2 n + 1 D$ . Für die nichttrivialen $z k$ (Indizes $0$ bis $2 n - 1 - 1$ ) gilt die Abschätzung $0 < z k < 2 n + 1 D < 2 n + 1 F n$ . Wegen der Teilerfremdheit der beiden Zahlen $2 n + 1$ und $F n$ genügt es wieder zur Bestimmung der $z k$ , die Reste $z_k \mod 2^{n+2}$ und $z_k \mod F_n$ zu berechnen.

[Bearbeiten] Bestimmung der Reste modulo 2ⁿ⁺¹

Wir wenden wieder den Trick der Einfügung von „Sicherheitsabständen“ an: Es seien also $\alpha_j = a_j\mod 2^{n+1}$ und $\beta_j = b_j\mod 2^{n+1}$ in $\Bbb Z_{2^{n+1}}$ . Wir bilden

$u = \sum_{k=0}^{2^n-1}\alpha_k 2^{k(3n+2)}$ und $v = \sum_{k=0}^{2^n-1}\beta_k 2^{k(3n+2)}$

und haben dabei $0 \leq u, v < 2^{2^n(3n+2)}$ . Das Produkt $u\cdot v$ enthält dann in disjunkten Stücken der Bitlänge $3 n + 2$ die Summen

$\gamma_k = \sum_{r,s=0 \atop r+s=k}^{2^n-1}\alpha_r\cdot \beta_s$

mit $0\leq k < 2^{n+1}$ . Für die gesuchten $c k$ unserer ursprünglichen Multiplikationsaufgabe $a\cdot b$ sehen wir

$c_k = \gamma_k + \gamma_{k+2^n}\mod 2^{n+1}$ .

Für die zu bestimmenden Reste $\eta_k = z_k \mod 2^{n+1}$ erhalten wir

$\eta_k = \gamma_k - \gamma_{k+2^{n-1}} + \gamma_{k+2\cdot 2^{n-1}} - \gamma_{k+3\cdot 2^{n-1}}$ in $\Bbb Z_{2^{n+1}}$ .

[Bearbeiten] Bestimmung der Reste modulo (D+1)

Mit $R = \Bbb Z_{D+1}$ unterziehen wir wieder die „Vektoren“ $(a_0,\ldots a_{2^n-1})$ und $(b_0,\ldots b_{2^n-1})$ mit $0 \leq a_k, b_k < \sqrt D$ der DFT in $R^{2^n}$ , wobei wir diesmal die Zahl $4$ als $2 n$ -te Einheitswurzel wählen. Da wir nur die Differenzen $c_k - c_{k+2^{n-1}}$ benötigen, genügt hier wiederum die „halbierte DFT“:

DFT zur Bestimmung der $\hat a_k$ und $\hat b_k$ nur für die ungeraden $k$ mit $0\leq k < 2^n$
$2 n - 1$ Multiplikationen $\hat c_k = \hat a_k\cdot \hat b_k$ für alle ungeraden $k$
Inverse DFT zur Gewinnung aller Differenzen $c_k - c_{k+2^{n-1}}$ aus den $\hat c_k$ für ungerade $k$

[Bearbeiten] Zusammenfassung

Startend von $a, b$ mit Ziffernlänge $N$ wird durch die dargestellte Rekursion eine Komplexität von $O (N log N loglog N)$ erreicht.

[Bearbeiten] Literatur

Schönhage, Arnold und Strassen, Volker: Schnelle Multiplikation ganzer Zahlen, Computing 7, 1971, Springer Verlag, S. 281-292
Arnold Schönhage: Asymptotically fast algorithms for the numerical multiplication and division of polynomials with complex coefficients. In: In EUROCAM '82: European Computer Algebra Conference (Marseille, France, April 1982), ed. J. Calmet.. Lect. Notes Comp. Sci. 144, 1982, S. 3-15
Knuth, Donald E.: The Art of Computer Programming, Vol. 2: Seminumerical Algorithms (3rd Edition), Addison-Wesley 1998, ISBN 0-201-89684-2 (engl.)

[Bearbeiten] Weblinks

Homepage von Prof. Dr. Arnold Schönhage
Weltrekord-Rechenmethode kommt zu späten Ehren (Presseinformation der Universität Bonn vom 21. Dezember 2004)

Von „http://de.wikipedia.org../../../s/c/h/Sch%C3%B6nhage-Strassen-Algorithmus_5925.html“

Kategorien: Algebra | Algorithmus | Arithmetik

Schönhage-Strassen-Algorithmus

aus Wikipedia, der freien Enzyklopädie

Inhaltsverzeichnis

[Bearbeiten] Bedeutung

[Bearbeiten] Algorithmus

[Bearbeiten] Grundidee und Terminologie

[Bearbeiten] Theoretische Vorbereitungen

[Bearbeiten] „Superschnelle DFT“

[Bearbeiten] Struktursatz über zyklische Arithmetik

[Bearbeiten] Durchführung

[Bearbeiten] Rekursionsschritt für ungerades m

[Bearbeiten] Bestimmung der Reste modulo 2ⁿ⁺²

[Bearbeiten] Bestimmung der Reste modulo (D+1)

[Bearbeiten] Rekursionsschritt für gerades m

[Bearbeiten] Bestimmung der Reste modulo 2ⁿ⁺¹

[Bearbeiten] Bestimmung der Reste modulo (D+1)

[Bearbeiten] Zusammenfassung

[Bearbeiten] Literatur

[Bearbeiten] Weblinks

Views

Navigation

Mitmachen

Suche

Andere Sprachen

Schönhage-Strassen-Algorithmus

aus Wikipedia, der freien Enzyklopädie

Inhaltsverzeichnis

[Bearbeiten] Bedeutung

[Bearbeiten] Algorithmus

[Bearbeiten] Grundidee und Terminologie

[Bearbeiten] Theoretische Vorbereitungen

[Bearbeiten] „Superschnelle DFT“

[Bearbeiten] Struktursatz über zyklische Arithmetik

[Bearbeiten] Durchführung

[Bearbeiten] Rekursionsschritt für ungerades m

[Bearbeiten] Bestimmung der Reste modulo 2n+2

[Bearbeiten] Bestimmung der Reste modulo (D+1)

[Bearbeiten] Rekursionsschritt für gerades m

[Bearbeiten] Bestimmung der Reste modulo 2n+1

[Bearbeiten] Bestimmung der Reste modulo (D+1)

[Bearbeiten] Zusammenfassung

[Bearbeiten] Literatur

[Bearbeiten] Weblinks

Views

Navigation

Mitmachen

Suche

Andere Sprachen

[Bearbeiten] Bestimmung der Reste modulo 2ⁿ⁺²

[Bearbeiten] Bestimmung der Reste modulo 2ⁿ⁺¹