Schönhage-Strassen-Algorithmus

Algorithmus zur Multiplikation

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 sehr schnellen 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 (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 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 Laufzeit dar.

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 der noch niedrigeren asymptotischen Komplexität , wobei der iterierte Logarithmus von n ist.[2] Durch Optimierungen des Algorithmus von Fürer erreichten David Harvey, Joris van der Hoeven und Grégoire Lecerf 2014 eine Verbesserung der asymptotischen Laufzeit auf .[3] Harvey und van der Hoeven stellten 2021 schließlich einen weiteren Algorithmus vor, der die von Schönhage und Strassen postulierte Laufzeit von erreicht.[4]

Bedeutung

Bearbeiten

Bis 2007 galt der Schönhage-Strassen-Algorithmus als effizientester bekannter Algorithmus für ganzzahlige Multiplikation. 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

Bearbeiten

Grundidee und Terminologie

Bearbeiten
 
Der Schönhage–Strassen-Algorithmus basiert auf der schnellen diskreten Fourier-Transformation (DFT). Dieses Beispiel zeigt die Berechnung von 1234 × 5678 = 7006652. Die Berechnung findet modulo 337 statt. Um die Anschaulichkeit zu verbessern, wird anstelle der Basis 2 mit Basis 10 gearbeitet.

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

  1. Aufspaltung der Zahlen (in Binärdarstellung)   und   in Stücke passender Länge
  2. Schnelle diskrete Fourier-Transformation (DFT) der beiden Stückfolgen
  3. Komponentenweise Multiplikation der transformierten Stücke
  4. Rücktransformation (inverse Fouriertransformation) der Ergebnisse
  5. 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   ein, so erhält man gerade die Binärdarstellung der Zahlen a und b. Zu berechnen ist   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   aus. Multipliziert man nun diese Funktionswerte, so ergeben sich die entsprechenden Funktionswerte des Produktpolynoms

 .

Um das Polynom   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  

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   Summen mit jeweils   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   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   die maximale Länge (also Binärziffernanzahl) der beiden Zahlen.

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

Theoretische Vorbereitungen

Bearbeiten

Superschnelle DFT

Bearbeiten

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   ein kommutativer unitärer Ring. In   sei das Element   eine Einheit; weiterhin sei   eine  te 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   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   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   im Ring   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   und schränken wir die weiteren Schritte der Rekursion ebenso auf   ein, so berechnen wir gerade alle   für ungerade Werte  . Will man umgekehrt aus diesen   für ungerade   (das sind   Stück) lediglich die Differenzen   der ursprünglichen   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   (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   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 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

Bearbeiten

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

Für eine Zweierpotenz   mit einer natürlichen Zahl   gilt

 .

Hierbei bezeichnet   die durch die Repräsentanten   darstellbaren Restklassen modulo   ausgestattet mit der Restklassenarithmetik, d. h. mit der Addition und Multiplikation modulo  . Die in diesem Restklassenring vorkommenden Zahlen können mit   Binärziffern dargestellt werden.

Die auf der rechten Seite vorkommende Struktur   bezeichnet die Restklassen modulo der Zahl  , 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   mit   nicht den Wert  , sondern den Wert  . Aus der so erhaltenen Zahlenstruktur mit zyklischer Arithmetik wird nun noch der Faktorring modulo   gebildet. Es werden also die Endergebnisse noch modulo   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  .

Entscheidend für das Gelingen der in diesem Struktursatz vorgestellten Einbettung ist die Eigenschaft, dass die größte darstellbare Zahl   im zyklischen Zahlenraum (hier ist dies die Zahl  ) die Zahl   aus dem Restklassenring   repräsentiert. Hierfür ist die Bedingung   notwendig. Damit die zyklische Arithmetik aber überhaupt sinnvoll definiert werden kann, muss andererseits   eine Zweierpotenz sein. Zusammen ergibt sich, dass   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  , d. h. die Zahl   ist in diesem Ring ein Nullteiler.

Durchführung

Bearbeiten

Haben wir die zu multiplizierenden Zahlen   mit   Binärziffern vorliegen, so führen wir je nachdem, ob   gerade oder ungerade ist, unterschiedliche Rekursionsschritte aus, um die Stellenzahl in einem Einzelschritt zu logarithmieren:

Rekursionsschritt für ungerades m

Bearbeiten

Diesen Schritt der Rückführung von   auf   führen wir mit der Komplexität   durch.

Es seien   mit   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   wird unsere Stückelungsgröße werden, d. h. wir entwickeln   und   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   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 diesem zyklischen Zahlenraum benutzt haben.

Zusammenfassend erhält die Multiplikation also die Form

 

mit den Ergebniskoeffizienten

 .

Wir können   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   (Indizes   bis  ) gilt die Abschätzung  . Da die beiden Zahlen   und   teilerfremd ist, genügt zur Bestimmung der   die Berechnung der Reste   und  .

Hat man nämlich die Reste   und   bestimmt, so kann man in Komplexität   wie folgt rechnen: Berechne erst   und dann  .

Bestimmung der Reste modulo 2n+2
Bearbeiten

Hier wenden wir einen für die Computeralgebra sehr typischen Trick an: Wir setzen die Stückfolgen   und   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   die Summen

 

mit  , denn es ist  . Für die Terme   unserer ursprünglichen Multiplikationsaufgabe   sehen wir

 .

Für die zu bestimmenden Reste   erhalten wir

  in  .

Der Komplexitätsaufwand für die Bildung aller   sowie der Extraktion der   ist  ; die Multiplikation   kostet  , insgesamt ist dies also  .

Bestimmung der Reste modulo (D+1)
Bearbeiten

Hier kommt die DFT zum Einsatz. Wir unterziehen die Vektoren   und   mit   der DFT in   mit   und der Zahl   als  -ter Einheitswurzel. Da wir nur die Differenzen   benötigen, genügt die halbierte DFT:

  • DFT zur Bestimmung der   und   nur für die ungeraden   mit  
  •   Multiplikationen   für alle ungeraden  
  • Inverse DFT zur Gewinnung aller Differenzen   aus den   für ungerade  

Der Komplexitätsaufwand hierfür besteht aus   Schritten des Einzelaufwands   für die DFT (gesamt also  ); hinzu kommen die Addition von   sowie die Reduktionen modulo   für die Gewinnung der  , was in   bewältigt werden kann.

Rekursionsschritt für gerades m

Bearbeiten

Auch für diesen Schritt der Rückführung von   auf   wird die Komplexität   erreicht.

Es seien   mit   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   unsere Stückelungsgröße werden, d. h. wir entwickeln   und   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  .

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   nach oben abschätzen.

Aus   können wir wieder

 

folgern, und mit

 

gilt

 

mit  . Für die nichttrivialen   (Indizes   bis  ) gilt die Abschätzung  . Wegen der Teilerfremdheit der beiden Zahlen   und   genügt es wieder zur Bestimmung der  , die Reste   und   zu berechnen.

Bestimmung der Reste modulo 2n+1
Bearbeiten

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   die Summen

 

mit  . Für die gesuchten   unserer ursprünglichen Multiplikationsaufgabe   sehen wir

 .

Für die zu bestimmenden Reste   erhalten wir

  in  .
Bestimmung der Reste modulo (D+1)
Bearbeiten

Mit   unterziehen wir wieder die Vektoren   und   mit   der DFT in  , wobei wir diesmal die Zahl   als  -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   mit  
  •   Multiplikationen   für alle ungeraden  
  • Inverse DFT zur Gewinnung aller Differenzen   aus den   für ungerade  

Zusammenfassung

Bearbeiten

Startend mit   und   mit Ziffernlänge   wird durch die dargestellte Rekursion eine Komplexität von   erreicht.

Abgewandelte Form

Bearbeiten

Zimmermann und Brent beschreiben eine Variante des Algorithmus, bei der die Laufzeit (in Abhängigkeit von der Länge der Eingabe) keine Sprünge macht, sondern stetiger verläuft. Dies wird erreicht, indem die DFT-Vektoren nicht aus  -stelligen Binärzahlen, sondern Zahlen der passenden Länge gebildet werden. Dadurch muss die Länge der zu transformierenden Vektoren keine Zweierpotenz sein.[5][6]

Literatur

Bearbeiten
Bearbeiten

Einzelnachweise

Bearbeiten
  1. Arnold Schönhage, Volker Strassen: Schnelle Multiplikation großer Zahlen. In: Computing, 7, 1971, S. 281–292, Springer Verlag
  2. Martin Fürer: Faster integer multiplication. STOC 2007 Proceedings, S. 57–66.
  3. David Harvey, Joris van der Hoeven, Grégoire Lecerf: Even faster integer multiplication. 2014, arxiv:1407.3360
  4. David Harvey, Joris van der Hoeven: Integer multiplication in time $O(n\mathrm{log}\, n)$. In: Annals of Mathematics. Band 193, Nr. 2, 1. März 2021, ISSN 0003-486X, doi:10.4007/annals.2021.193.2.4 (projecteuclid.org [abgerufen am 17. April 2024]).
  5. loria.fr (PDF; 1,9 MB) S. 56
  6. loria.fr (PDF)