- Schönhage-Strassen-Algorithmus
-
Der Schönhage-Strassen-Algorithmus ist ein Algorithmus zur Multiplikation zweier n-stelliger ganzer Zahlen. Er wurde 1971 von Arnold Schönhage und Volker Strassen entwickelt.[1] 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 war von 1971 bis 2007 der effizienteste bekannte Algorithmus zur Multiplikation großer Zahlen; 2007 veröffentlichte Martin Fürer eine Weiterentwicklung des Algorithmus mit einer noch niedrigeren asymptotischen Komplexität.[2]
Der Schönhage-Strassen-Algorithmus terminiert in
(siehe 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
als auch gegenüber dem 1962 entwickelten Karatsuba-Algorithmus mit einer Laufzeit von
sowie dessen verbesserter Variante, dem Toom-Cook-Algorithmus mit O(n1 + ε) Laufzeit dar.
Inhaltsverzeichnis
Bedeutung
Bis 2007 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
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.
Algorithmus
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“ Multiplikationen werden im rekursiven Sinne wiederum durch den Schönhage-Strassen-Algorithmus ausgeführt.
Um zu verstehen, warum das Ergebnis das Produkt der Zahlen a und b ist, betrachtet man die Polynome
und
Setzt man X = 2 ein, so erhält man gerade die Binärdarstellung der Zahlen a und b. Zu berechnen ist c = C(2) für das Produktpolynom
Wir bestimmen die Fouriertransformierte der Koeffiziententupel von A und B:
für
für
Anders gesagt wertet man die beiden Polynome an den Stellen wk aus. Multipliziert man nun diese Funktionswerte, so ergeben sich die entsprechenden Funktionswerte des Produktpolynoms
.
Um das Polynom C selbst zu gewinnen, müssen wir die Transformation rückgängig machen:
für
für
für
Nach Definition der Einheitswurzeln gilt
. Diese genügt folgender Identität geometrischer Summen von Einheitswurzeln:
für l = 0,...,2n + 1 − 1
denn
für
Somit gilt:
für
Im Artikel Diskrete Fourier-Transformation sind die mathematische Grundlagen dieser Transformation weiter ausgeführt. Da bei der Transformation 2n Summen mit jeweils 2n Termen entstehen, haben wir bei einer klassischen Berechnung der Terme (etwa durch das Horner-Schema) nach wie vor eine quadratische Laufzeit. Mittels der schnellen Fourier-Transformation kann man diese Werte schneller berechnen. Diese Berechnung beruht auf folgendem Teile und herrsche Prinzip:
Man setzt Teillösungen mittels einfacher Operationen (Addition und einfache Multiplikation) zusammen. Damit können die Transformationen in Zeit
berechnet werden. Durch das Runden der komplexen Einheitswurzeln auf feste Stellenlänge ergeben sich jedoch Rechenfehler. Um diese auszugleichen, muss für ein resultierendes Bit mit mindestens O(log N) Bits gerechnet werden. Daraus ergibt sich eine Gesamtlaufzeit von
. Bei der Schönhage-Strassen-Variante rechnen wir stattdessen in einem Restklassenring und vermeiden damit die Rechenfehler der komplexen Zahlen.
Des Weiteren ist die Multiplikation keine reine Faltung, sondern es kann auch zu Überträgen kommen; nach Durchführen der FT und iFT müssen diese passend behandelt werden.
Die Aufgabe der Multiplikation zweier ganzer Zahlen wird nun wie folgt konkretisiert:
Es seien die zwei zu multiplizierenden Zahlen
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
natürliche Zahlen sind. Der Schönhage-Strassen-Algorithmus löst diese Aufgabe in
.
Theoretische Vorbereitungen
„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
eine 2nte Einheitswurzel (also
), die die Gleichheit
erfüllt. Dann lässt sich die Berechnung der diskreten Fouriertransformation (DFT) im Produktraum
(dies ist eine Kurznotation für
; 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
die Transformierte
mit
für
.
Indem wir die Indizes
und
in Binärdarstellung aufschreiben, wobei wir dies bei der Zahl j in umgekehrter Reihenfolge tun, ist die Transformierte
wie folgt optimiert berechenbar:
Es seien
für
und
.
Die geschlossene Darstellung für diese Zwischenterme ist
.
(Zum Nachrechnen dieser Darstellung beachte man
).
Diese Rekursion liefert die gewünschten Fourierkoeffizienten
.
Aufgrund der Eigenschaft
können wir den Rekursionsschritt etwas berechnungsfreundlicher umformen zu
und
mit dem gleichen Exponenten
.
Die Umkehrtransformation, also die inverse FFT, gelingt, da wir vorausgesetzt haben, dass 2 im Ring R invertierbar ist:
sowie
,
wobei wiederum
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
nur für k0 = 1 und schränken wir die weiteren Schritte der Rekursion ebenso auf k0 = 1 ein, so berechnen wir gerade alle
für ungerade Werte k. Will man umgekehrt aus diesen
für ungerade k (das sind 2n − 1 Stück) lediglich die Differenzen
der ursprünglichen aj 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
mit Fermatzahlen
benötigt.
Hinweis zur Notation: Für den Restklassenring
benutzen wir hier die kürzere Schreibweise
, die lediglich im Kontext der p-adischen Zahlen zu Verwechslungen führen könnte.
Als Einheitswurzel wird im Ring
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
; allerdings sind sie nicht als reine Shift-Operationen durchführbar, da das Reduzieren eines größeren Zwischenergebnisses modulo Fn 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 Shift-Operation durchführbar ist. Diesen Trick kann man in einem schönen Struktursatz über Restklassen- und zyklische Arithmetik in endlichen Zahlenringen zusammenfassen.
Struktursatz über zyklische Arithmetik
Der Struktursatz über zyklische Arithmetik lässt sich formal wie folgt fassen:
Für eine Zweierpotenz D = 2n mit einer natürlichen Zahl
gilt
.
Hierbei bezeichnet
die durch die Repräsentanten
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
bezeichnet die Restklassen modulo der Zahl D2, 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 = D2 − 1 nicht den Wert
, sondern den Wert
. 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
kann ebenso ersetzt werden durch das zyklische Rechnen im größeren Zahlenraum
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 D2 − 1) die Zahl 0 aus dem Restklassenring
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 = D2 − 1 die optimale Wahl für die Größe des zyklischen Einbettungsraumes darstellt.
Der klassische Restklassenring
wäre für die Einbettung dagegen nicht geeignet, denn in diesem Ring gilt 22n = D2 = 0, d. h. die Zahl 2 ist in diesem Ring ein Nullteiler.
Durchführung
Haben wir die zu multiplizierenden Zahlen a,b mit
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:
Rekursionsschritt für ungerades m
Diesen Schritt der Rückführung von m = 2n − 1 auf n führen wir mit der Komplexität
durch.
Es seien
mit m = 2n − 1 und der Fermatzahl
zu multiplizieren. Wir werden in diesem Schritt die Rückführung auf die Fermatzahl
vollziehen.
Für die zu den beiden Fermatzahlen gehörenden Zweierpotenzen führen wir die Abkürzungen
und
ein. Die halbierte Stellenzahl von D wird unsere Stückelungsgröße werden, d. h. wir entwickeln a und b nach Potenzen von
:
und
,
wobei für die „Einzelstücke“
gilt. In Binärdarstellung entspricht diese Zerlegung einer einfachen Gruppierung der Bitfolgen in Stücke der Länge 2n − 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
und
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
und
.
Vermöge des oben erwähnten Struktursatzes zur zyklischen Arithmetik wechseln wir nun vom Restklassenring
über zum Quotientenraum
mit der zyklischen Arithmetik. In diesem Raum errechnet sich für die Multiplikationsaufgabe
,
wobei wir im letzten Schritt die Eigenschaft
in unserem zyklischen „Rechenraum“ benutzt haben.
Zusammenfassend erhält die Multiplikation also die Form
mit den Ergebniskoeffizienten
.
Wir können ck < 2n + 1D nach oben abschätzen.
Nun folgt eine Umschreibung der Summenformel, damit wir uns bei der anzuwendenden FFT auf eine „halbierte FFT“ beschränken können.
Es gilt
, also ist
mit
in
. Durch passende Addition können wir den Wertebereich ins Positive verschieben, es ist nämlich
, und mit der Definition
gilt
.
Für die nichttrivialen zk (Indizes 0 bis 2n − 1) gilt die Abschätzung 0 < zk < 2n + 2D < 2n + 2Fn. Da die beiden Zahlen 2n + 2 und Fn teilerfremd ist, genügt zur Bestimmung der zk die Berechnung der Reste
und
.
Hat man nämlich die Reste
und
bestimmt, so kann man in Komplexität O(2n) wie folgt rechnen: Berechne erst
und dann
.
Bestimmung der Reste modulo 2n+2
Hier wenden wir einen für die Computeralgebra sehr typischen Trick an: Wir setzen die Stückfolgen ai und bi 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
und
in
. Wir bilden nun
und
und haben dabei
. Das Produkt
enthält dann in disjunkten Stücken der Bitlänge 3n + 5 die Summen
mit
, denn es ist
. Für die Terme ck unserer ursprünglichen Multiplikationsaufgabe
sehen wir
.
Für die zu bestimmenden Reste
erhalten wir
in
.
Der Komplexitätsaufwand für die Bildung aller αj,βj sowie der Extraktion der ηk ist O(22n); die Multiplikation
„kostet“ ψ(2n + 1(3n + 5)), insgesamt ist dies also O(22n).
Bestimmung der Reste modulo (D+1)
Hier kommt die DFT zum Einsatz. Wir unterziehen die „Vektoren“
und
mit
der DFT in
mit
und der Zahl 2 als 2n + 1-ter Einheitswurzel. Da wir nur die Differenzen
benötigen, genügt die „halbierte DFT“:
- DFT zur Bestimmung der
und
nur für die ungeraden k mit
- 2n Multiplikationen
für alle ungeraden k
- Inverse DFT zur Gewinnung aller Differenzen
aus den
für ungerade k
Der Komplexitätsaufwand hierfür besteht aus O(n2n) Schritten des Einzelaufwands O(2n) für die DFT (gesamt also O(n22n)); hinzu kommen das Aufaddieren von 2n + 1D sowie die Reduktionen modulo (D + 1) für die Gewinnung der
, was in O(22n) bewältigt werden kann.
Rekursionsschritt für gerades m
Auch für diesen Schritt der Rückführung von m = 2n − 2 auf n wird die Komplexität
erreicht.
Es seien
mit m = 2n − 2 und der Fermatzahl
zu multiplizieren. Wir werden auch in diesem Schritt die Rückführung auf die Fermatzahl
vollziehen.
Für die zu den beiden Fermatzahlen gehörenden Zweierpotenzen führen wir analog die Abkürzungen
und
ein. Wiederum wird die halbierte Stellenzahl von D unsere Stückelungsgröße werden, d. h. wir entwickeln a und b nach Potenzen von
:
und
,
wobei für die „Einzelstücke“
gilt.
Wie oben verlängern wir die Zahlendarstellung auf Zweierpotenzlänge zu
und analog für b.
Unter abermaliger Zuhilfenahme des Struktursatzes zur zyklischen Arithmetik wechseln wir nun vom Restklassenring
über zum Quotientenraum
mit der zyklischen Arithmetik.
Damit können wir wieder
mit den Ergebniskoeffizienten
darstellen. Dabei können wir ck < 2nD nach oben abschätzen.
Aus
können wir wieder
folgern, und mit
gilt
mit 0 < zk < 2n + 1D. Für die nichttrivialen zk (Indizes 0 bis 2n − 1 − 1) gilt die Abschätzung 0 < zk < 2n + 1D < 2n + 1Fn. Wegen der Teilerfremdheit der beiden Zahlen 2n + 1 und Fn genügt es wieder zur Bestimmung der zk, die Reste
und
zu berechnen.
Bestimmung der Reste modulo 2n+1
Wir wenden wieder den Trick der Einfügung von „Sicherheitsabständen“ an: Es seien also
und
in
. Wir bilden
und
und haben dabei
. Das Produkt
enthält dann in disjunkten Stücken der Bitlänge 3n + 2 die Summen
mit
. Für die gesuchten ck unserer ursprünglichen Multiplikationsaufgabe
sehen wir
.
Für die zu bestimmenden Reste
erhalten wir
in
.
Bestimmung der Reste modulo (D+1)
Mit
unterziehen wir wieder die „Vektoren“
und
mit
der DFT in
, wobei wir diesmal die Zahl 4 als 2n-te Einheitswurzel wählen. Da wir nur die Differenzen
benötigen, genügt hier wiederum die „halbierte DFT“:
- DFT zur Bestimmung der
und
nur für die ungeraden k mit
- 2n − 1 Multiplikationen
für alle ungeraden k
- Inverse DFT zur Gewinnung aller Differenzen
aus den
für ungerade k
Zusammenfassung
Startend mit a und b mit Ziffernlänge n wird durch die dargestellte Rekursion eine Komplexität von
erreicht.
Quellen
- ↑ Arnold Schönhage und Volker Strassen: Schnelle Multiplikation großer Zahlen, Computing 7, 1971, Springer Verlag, S. 281–292
- ↑ Martin Fürer: Faster integer multiplication, STOC 2007 Proceedings, S. 57–66.
Literatur
- Arnold Schönhage: Asymptotically fast algorithms for the numerical multiplication and division of polynomials with complex coefficients. In: J. Calmet. (Hrsg.): EUROCAM '82: European Computer Algebra Conference (Marseille, France, April 1982). Lect. Notes Comp. Sci. 144, Springer, 1982.
- Donald E. Knuth: The Art of Computer Programming, Vol. 2: Seminumerical Algorithms. 3. Auflage. Addison-Wesley, 1998, ISBN 0-201-89684-2.
- Chee Yap, Chen Li: QuickMul: Practical FFT-based Integer Multiplication, 2000. Vereinfachung des Schönhage-Strassen-Algorithmus für praktische Anwendungen.
- Michael T. Goodrich, Roberto Tamassia: Algorithm Design Foundations, Analysis, and Internet Examples. Eine Einführung zur FFT mit einer Java-Implementierung des QuickMul-Algorithmus.
- Daniel J. Bernstein: Multidigit multiplication for mathematicians. 11. August 2001 (Zusammenfassung verschiedener Techniken zur Polynom- und Langzahlmultiplikation, http://cr.yp.to/papers.html#m3).
Weblinks
- Homepage von Prof. Dr. Arnold Schönhage
- Weltrekord-Rechenmethode kommt zu späten Ehren (Presseinformation der Universität Bonn vom 21. Dezember 2004)
Wikimedia Foundation.