Schönhage-Straßen-Algorithmus - Schönhage–Strassen algorithm

Der Schönhage-Strassen-Algorithmus ist ein asymptotisch schneller Multiplikationsalgorithmus für große ganze Zahlen . Es wurde entwickelt von Arnold Schönhage und Volker Strassen im Jahr 1971. Die Laufzeit - Bit - Komplexität ist, in Big O - Notation , für zwei n - stellige Zahlen. Der Algorithmus verwendet rekursive Fast-Fourier-Transformationen in Ringen mit 2 ⁿ +1 Elementen, eine spezielle Art der zahlentheoretischen Transformation . ${\displaystyle O(n\cdot \log n\cdot \log \log n)}$

Der Schönhage-Strassen-Algorithmus war das asymptotisch schnellste bekannte Multiplikationsverfahren von 1971 bis 2007, als ein neues Verfahren, der Fürer-Algorithmus , mit geringerer asymptotischer Komplexität angekündigt wurde; der Algorithmus von Fürer erzielt jedoch derzeit nur bei astronomisch großen Werten einen Vorteil und wird nur in Basic Polynomial Algebra Subprograms (BPAS) verwendet (siehe Galaktische Algorithmen ).

In der Praxis beginnt der Schönhage-Strassen-Algorithmus ältere Methoden wie die Karatsuba- und die Toom-Cook-Multiplikation für Zahlen jenseits von 2 ^{2 15} bis 2 ^{2 17} (10.000 bis 40.000 Dezimalstellen) zu übertreffen . Die GNU Multi-Precision Library verwendet sie je nach Architektur für Werte von mindestens 1728 bis 7808 64-Bit-Wörtern (33.000 bis 150.000 Dezimalstellen). Es gibt eine Java-Implementierung von Schönhage–Strassen, die es über 74.000 Dezimalstellen verwendet.

Zu den Anwendungen des Schönhage-Strassen-Algorithmus gehören die mathematische Empirie wie die Great Internet Mersenne Prime Search und die Berechnung von Approximationen von π sowie praktische Anwendungen wie die Kronecker-Substitution , bei der die Multiplikation von Polynomen mit ganzzahligen Koeffizienten effizient auf große ganze Zahlen reduziert werden kann Multiplikation; dies wird in der Praxis von GMP-ECM für die elliptische Kurvenfaktorisierung von Lenstra verwendet .

Einzelheiten

In diesem Abschnitt wird ausführlich erläutert, wie Schönhage–Strassen umgesetzt wird. Es basiert hauptsächlich auf einem Überblick über die Methode von Crandall und Pomerance in ihren Primzahlen: A Computational Perspective . Diese Variante unterscheidet sich etwas von der ursprünglichen Methode von Schönhage darin, dass sie die diskrete gewichtete Transformation nutzt , um negazyklische Faltungen effizienter durchzuführen . Eine weitere Quelle für detaillierte Informationen Knuth ‚s The Art of Computer Programming . Der Abschnitt beginnt mit der Diskussion der Bausteine ​​und endet mit einer schrittweisen Beschreibung des Algorithmus selbst.

Windungen

Angenommen, wir multiplizieren zwei Zahlen wie 123 und 456 unter Verwendung einer langen Multiplikation mit Ziffern zur Basis B , jedoch ohne jegliches Tragen. Das Ergebnis könnte etwa so aussehen:

		1	2	3
	×	4	5	6
		6	12	18
	5	10	fünfzehn
4	8	12
4	13	28	27	18

Diese Folge (4, 13, 28, 27, 18) wird als azyklische oder lineare Faltung der beiden ursprünglichen Folgen (1,2,3) und (4,5,6) bezeichnet. Sobald wir die azyklische Faltung zweier Folgen haben, ist die Berechnung des Produkts der ursprünglichen Zahlen einfach: Wir führen einfach das Tragen durch (zum Beispiel behalten wir in der Spalte ganz rechts die 8 bei und fügen die 1 zu der Spalte mit 27 hinzu). Im Beispiel ergibt dies das richtige Produkt 56088.

Es gibt zwei andere Arten von Faltungen, die nützlich sein werden. Angenommen, die Eingabesequenzen haben n Elemente (hier 3). Dann hat die azyklische Faltung n + n -1 Elemente; wenn wir die ganz rechten n Elemente nehmen und die ganz linken n −1 Elemente hinzufügen , ergibt dies die zyklische Faltung :

	28	27	18
+		4	13
	28	31	31

Wenn wir die zyklische Faltung weiterführen, ist das Ergebnis äquivalent zum Produkt der Eingaben mod B ⁿ  − 1. Im Beispiel 10 ³  − 1 = 999 ergibt die Weiterführung (28, 31, 31) 3141, und 3141 56088 (mod 999).

Umgekehrt, wenn wir die ganz rechten n Elemente nehmen und die ganz linken n −1 Elemente subtrahieren , ergibt dies die negazyklische Faltung :

	28	27	18
−		4	13
	28	23	5

Wenn wir die negazyklische Faltung weiterführen, ist das Ergebnis äquivalent zum Produkt der Eingaben mod B ⁿ  + 1. Im Beispiel 10 ³  + 1 = 1001, die Weiterführung (28, 23, 5) ergibt 3035, und 3035 ≡ 56088 (mod 1001). Die negazyklische Faltung kann negative Zahlen enthalten, die beim Tragen durch Borrowing eliminiert werden können, wie dies bei der langen Subtraktion der Fall ist.

Faltungstheorem

Wie andere auf der Fast Fourier-Transformation basierende Multiplikationsverfahren hängt Schönhage–Strassen grundlegend vom Faltungssatz ab , der eine effiziente Möglichkeit bietet, die zyklische Faltung zweier Folgen zu berechnen. Es sagt, dass:

Die zyklische Faltung von zwei Vektoren kann gefunden werden, indem man die diskrete Fourier-Transformation (DFT) von jedem von ihnen nimmt, die resultierenden Vektoren Element für Element multipliziert und dann die inverse diskrete Fourier-Transformation (IDFT) nimmt.

Oder in Symbolen:

Zyklische Faltung( X, Y ) = IDFT(DFT( X ) · DFT( Y ))

Wenn wir die DFT und IDFT mit einem schnellen Fourier-Transformationsalgorithmus berechnen und unseren Multiplikationsalgorithmus rekursiv aufrufen, um die Einträge der transformierten Vektoren DFT( X ) und DFT( Y ) zu multiplizieren , ergibt dies einen effizienten Algorithmus zum Berechnen der zyklischen Faltung.

Bei diesem Algorithmus ist es nützlicher, die negazyklische Faltung zu berechnen ; wie sich herausstellt, kann dies auch durch eine leicht modifizierte Version des Faltungstheorems (siehe diskrete gewichtete Transformation ) ermöglicht werden. Angenommen, die Vektoren X und Y haben die Länge n und a ist eine primitive Einheitswurzel der Ordnung 2 n (dh a ^{2 n} = 1 und a in alle kleineren Potenzen ist nicht 1). Dann können wir einen dritten Vektor A , den Gewichtsvektor , definieren als:

A = ( a ^j ), 0 ≤ j < n

A ⁻¹ = ( a ^{− j} ), 0 ≤ j < n

Nun können wir feststellen:

Negazyklische Faltung( X , Y ) = A ⁻¹ · IDFT(DFT( A · X ) · DFT( A · Y ))

Mit anderen Worten, es ist dasselbe wie zuvor, außer dass die Eingaben zuerst mit A multipliziert werden und das Ergebnis mit A ^-1 multipliziert wird .

Ringwahl

Die diskrete Fourier-Transformation ist eine abstrakte Operation, die in jedem algebraischen Ring ausgeführt werden kann ; normalerweise wird es in komplexen Zahlen durchgeführt, aber komplexe Arithmetik mit ausreichender Genauigkeit durchzuführen, um genaue Ergebnisse für die Multiplikation zu gewährleisten, ist langsam und fehleranfällig. Stattdessen verwenden wir den Ansatz der zahlentheoretischen Transformation , die darin besteht, die Transformation in den ganzen Zahlen mod N für eine ganze Zahl N durchzuführen .

Genauso wie es primitive Einheitswurzeln jeder Ordnung in der komplexen Ebene gibt, können wir für jede beliebige Ordnung n ein geeignetes N wählen, so dass b eine primitive Einheitswurzel der Ordnung n in den ganzen Zahlen mod N ist (mit anderen Worten, b ⁿ ≡ 1 (mod N ), und keine kleinere Potenz von b entspricht 1 mod N ).

Der Algorithmus wird die meiste Zeit damit verbringen, rekursive Multiplikationen kleinerer Zahlen durchzuführen; bei einem naiven Algorithmus treten diese an mehreren Stellen auf:

1. Innerhalb des Fast-Fourier-Transformationsalgorithmus, bei dem die primitive Einheitswurzel b wiederholt gepowert, quadriert und mit anderen Werten multipliziert wird.
2. Beim Ziehen von Potenzen der primitiven Einheitswurzel a zur Bildung des Gewichtsvektors A und beim Multiplizieren von A oder A ⁻¹ mit anderen Vektoren.
3. Beim Ausführen einer Element-für-Element-Multiplikation der transformierten Vektoren.

Die wichtigste Erkenntnis von Schönhage-Strassen besteht darin, N, den Modul, gleich 2 ⁿ  + 1 für eine ganze Zahl n zu wählen, die ein Vielfaches der Anzahl der zu transformierenden Teile P = 2 ^p ist. Dies hat in Standardsystemen, die große Ganzzahlen in binärer Form darstellen, eine Reihe von Vorteilen:

• Jeder Wert kann schnell modulo 2 ⁿ  + 1 reduziert werden, indem nur Verschiebungen und Additionen verwendet werden, wie im nächsten Abschnitt erklärt .
• Alle von der Transformation verwendeten Einheitswurzeln können in der Form 2 ^{k geschrieben werden} ; folglich können wir jede Zahl mit einer Einheitswurzel multiplizieren oder dividieren, indem wir eine Verschiebung verwenden, und eine Einheitswurzel potenzieren oder quadrieren, indem wir nur auf ihren Exponenten operieren.
• Die Element-für-Element-rekursiven Multiplikationen der transformierten Vektoren können mit einer negazyklischen Faltung durchgeführt werden, die schneller als eine azyklische Faltung ist und bereits "umsonst" den Effekt hat, ihr Ergebnis mod 2 ⁿ  + 1 zu reduzieren .

Um die rekursiven Multiplikationen bequem zu machen, werden wir Schönhage–Strassen als einen spezialisierten Multiplikationsalgorithmus formulieren, der nicht nur das Produkt zweier Zahlen berechnet, sondern das Produkt zweier Zahlen mod 2 ⁿ  + 1 für ein gegebenes n . Dies ist keine Einschränkung der Allgemeinheit, da man n immer so groß wählen kann, dass das Produkt mod 2 ⁿ  + 1 einfach das Produkt ist.

Schichtoptimierungen

Im Verlauf des Algorithmus gibt es viele Fälle, in denen die Multiplikation oder Division mit einer Zweierpotenz (einschließlich aller Einheitswurzeln) gewinnbringend durch wenige Verschiebungen und Additionen ersetzt werden kann. Dies macht sich die Beobachtung zunutze, dass:

(2 ⁿ ) ^k ≡ (−1) ^k mod (2 ⁿ + 1)

Beachten Sie, dass eine k- stellige Zahl zur Basis 2 ⁿ in Positionsnotation geschrieben werden kann als ( d _{k −1} ,..., d ₁ , d ₀ ). Es stellt die Zahl dar . Beachten Sie auch, dass für jedes d _i 0 ≤ d _i < 2 ^{n gilt} . ${\displaystyle \scriptstyle \sum_{i=0}^{k-1}d_{i}\cdot (2^{n})^{i}}$

Dies macht es einfach, eine Zahl, die in binär mod 2 ⁿ  + 1 dargestellt wird, zu reduzieren : Nimm die ganz rechten (niedrigstwertigen) n Bits, subtrahiere die nächsten n Bits, addiere die nächsten n Bits und so weiter, bis die Bits erschöpft sind. Wenn der resultierende Wert immer noch nicht zwischen 0 und 2 ^{n liegt} , normalisieren Sie ihn, indem Sie ein Vielfaches des Moduls 2 ⁿ  + 1 addieren oder subtrahieren . Wenn zum Beispiel n = 3 ist (also der Modul 2 ³ +1 = 9) und die reduzierte Zahl 656 ist, haben wir:

656 = 1010010000 ₂ 000 ₂ − 010 ₂ + 010 ₂ − 1 ₂ = 0 − 2 + 2 − 1 = −1 ≡ 8 (mod 2 ³ + 1).

Darüber hinaus ist es möglich, sehr große Verschiebungen zu bewirken, ohne jemals das verschobene Ergebnis zu konstruieren. Angenommen, wir haben eine Zahl A zwischen 0 und 2 ⁿ und möchten sie mit 2 ^k multiplizieren . Wenn wir k durch n teilen, erhalten wir k = qn + r mit r < n . Es folgt dem:

A(2 ^k ) = A(2 ^{qn + r} ) = A[(2 ⁿ ) ^q (2 ^r )] (−1) ^q Nach links verschieben (A, r ) (mod 2 ⁿ + 1).

Da A ≤ 2 ⁿ und r < n ist , hat eine Linksverschiebung r höchstens 2 n −1 Bits, und daher ist nur eine Verschiebung und Subtraktion (gefolgt von einer Normalisierung) erforderlich.

Um schließlich durch 2 ^k zu dividieren , beobachte, dass die Quadrierung der ersten obigen Äquivalenz ergibt:

2 ^{2 n} ≡ 1 (mod 2 ⁿ + 1)

Somit,

A/2 ^k = A(2 ^{− k} ) A(2 ^{2 n − k} ) = Shift_Left(A, 2 n − k ) (mod 2 ⁿ + 1).

Algorithmus

Der Algorithmus folgt einer Aufteilung, Auswertung (Vorwärts-FFT), punktweiser Multiplikation, Interpolation (inverse FFT) und kombiniert Phasen ähnlich den Methoden von Karatsuba und Toom-Cook.

Bei gegebenen Eingangszahlen x und y und einer ganzen Zahl N berechnet der folgende Algorithmus das Produkt xy mod 2 ^N  + 1. Vorausgesetzt, N ist ausreichend groß, ist dies einfach das Produkt.

1. Teilen Sie jede Eingabezahl in Vektoren X und Y von jeweils 2 ^k Teilen auf, wobei 2 ^k N teilt . (zB 12345678 → (12, 34, 56, 78)).
2. Um Fortschritte zu erzielen, ist es notwendig, ein kleineres N für rekursive Multiplikationen zu verwenden. Wählen Sie dazu N als kleinste ganze Zahl, die mindestens 2 N /2 ^k + k und durch 2 ^k teilbar ist .
3. Berechnen Sie das Produkt von X und Y mod 2 ^N  + 1 unter Verwendung der negazyklischen Faltung:
   1. Multiplizieren Sie X und Y jeweils mit dem Gewichtungsvektor A unter Verwendung von Verschiebungen (verschieben Sie den j- ten Eintrag nach links um jn / ^2k ).
   2. Berechnen Sie die DFT von X und Y mit der zahlentheoretischen FFT (führen Sie alle Multiplikationen mit Verschiebungen durch; für die ^2k- te Einheitswurzel verwenden Sie 2 ^{2 n /2 k} ).
   3. Wenden Sie diesen Algorithmus rekursiv an, um entsprechende Elemente der transformierten X und Y zu multiplizieren.
   4. Berechnen Sie die IDFT des resultierenden Vektors, um den Ergebnisvektor C zu erhalten (führen Sie alle Multiplikationen mit Verschiebungen durch). Dies entspricht der Interpolationsphase.
   5. Multiplizieren Sie den Ergebnisvektor C mit A ^{−1 unter} Verwendung von Verschiebungen.
   6. Vorzeichen anpassen: Einige Elemente des Ergebnisses können negativ sein. Wir berechnen den größtmöglichen positiven Wert für das j- te Element von C, ( j + 1)2 ^{2 N /2 k} , und wenn er diesen überschreitet, ziehen wir den Modulus 2 ⁿ  + 1 ab.
4. Führen Sie schließlich das Tragen von Mod 2 ^N  + 1 durch, um das Endergebnis zu erhalten.

Die optimale Stückzahl zum Aufteilen der Eingabe ist proportional zu , wobei N die Anzahl der Eingabebits ist und diese Einstellung die Laufzeit von O( N log N log log N ) erreicht, daher sollte der Parameter k entsprechend eingestellt werden. In der Praxis wird er empirisch basierend auf den Eingangsgrößen und der Architektur festgelegt, typischerweise auf einen Wert zwischen 4 und 16. ${\displaystyle {\sqrt {N}}}$

In Schritt 2 wird die Beobachtung verwendet, dass:

• Jedes Element der Eingangsvektoren hat höchstens n /2 ^k Bits;
• Das Produkt von zwei beliebigen Eingangsvektorelementen hat höchstens 2 n /2 ^k Bits;
• Jedes Element der Faltung ist die Summe von höchstens 2 ^k solchen Produkten, und so kann nicht mehr als 2 n / 2 ^k + k - Bits.
• n muss durch 2 ^k teilbar sein, um sicherzustellen, dass in den rekursiven Aufrufen die Bedingung "2 ^k teilt N " in Schritt 1 gilt.

Optimierungen

In diesem Abschnitt werden einige wichtige praktische Optimierungen erläutert, die bei der Implementierung von Schönhage–Straßen in realen Systemen berücksichtigt wurden. Es basiert hauptsächlich auf einer Arbeit von Gaudry, Kruppa und Zimmermann aus dem Jahr 2007, die Verbesserungen der GNU Multi-Precision Library beschreibt .

Unterhalb eines bestimmten Grenzwerts ist es effizienter, die rekursiven Multiplikationen mit anderen Algorithmen wie der Toom-Cook-Multiplikation durchzuführen . Die Ergebnisse müssen um mod 2 ⁿ  + 1 reduziert werden, was effizient wie oben in Verschiebungsoptimierungen mit Verschiebungen und Additionen/Subtraktionen erläutert durchgeführt werden kann .

Das Berechnen der IDFT beinhaltet das Dividieren jedes Eintrags durch die primitive Einheitswurzel 2 ^{2 n /2 k} , eine Operation, die häufig mit der anschließenden Multiplikation des Vektors mit A ^–1 kombiniert wird , da beide eine Division durch eine Zweierpotenz beinhalten.

In einem System , in dem eine große Zahl als ein Array von 2 dargestellt ist , ^w -Bit Worten, es ist nützlich , daß die Vektorgröße 2 sicherzustellen ^k ist auch ein Vielfaches der Bits pro Wort durch die Wahl k ≥ w (zB wählen k ≥ 5 auf ein 32-Bit-Computer und k ≥ 6 auf einem 64-Bit-Computer); dies ermöglicht es, die Eingänge ohne Bitverschiebungen in Stücke aufzubrechen, und liefert eine einheitliche Darstellung für Werte mod 2 ⁿ  + 1, bei denen das High-Wort nur null oder eins sein kann.

Die Normalisierung beinhaltet das Addieren oder Subtrahieren des Moduls 2 ⁿ  + 1; Dieser Wert hat nur zwei Bits gesetzt, was bedeutet, dass dies mit einer spezialisierten Operation im Durchschnitt in konstanter Zeit durchgeführt werden kann.

Iterative FFT-Algorithmen wie der Cooley-Tukey-FFT-Algorithmus , obwohl sie häufig für FFTs auf Vektoren mit komplexen Zahlen verwendet werden, neigen dazu, bei den großen Vektoreinträgen, die in Schönhage-Strassen verwendet werden , eine sehr schlechte Cache- Lokalität aufzuweisen. Die einfache rekursive, nicht-in-place-Implementierung von FFT ist erfolgreicher, wobei alle Operationen über einen bestimmten Punkt in der Aufruftiefe hinaus in den Cache passen, nutzt den Cache jedoch in höheren Aufruftiefen immer noch suboptimal. Gaudry, Kruppa und Zimmerman verwendeten eine Technik, die den 4-Schritt-Algorithmus von Bailey mit Transformationen mit höherem Radix kombiniert, die mehrere rekursive Schritte kombinieren. Sie mischen auch Phasen und gehen für jedes Element des Vektors so weit wie möglich in den Algorithmus ein, bevor sie zum nächsten übergehen.

Der von Schönhage erstmals beschriebene "Quadratwurzel-von-2-Trick" besteht darin, dass, vorausgesetzt, k 2, 2 ^{3 n /4} −2 ^{n /4} eine Quadratwurzel von 2 mod 2 ⁿ +1 ist, also a 4 n -te Einheitswurzel (da 2 ^{2 n} 1). Dadurch kann die Transformationslänge von 2 ^k auf 2 ^{k + 1 erweitert werden} .

Schließlich achten die Autoren darauf, den richtigen Wert von k für verschiedene Bereiche von Eingabezahlen zu wählen , da der optimale Wert von k mit zunehmender Eingabegröße mehrmals zwischen denselben Werten hin und her gehen kann.

Verweise