Algoritmus Schönhage – Strassen - Schönhage–Strassen algorithm

Algoritmus Schönhage – Strassen je založen na metodě rychlé Fourierovy transformace (FFT) celočíselného násobení . Tento obrázek ukazuje násobení 1234 × 5678 = 7006652 pomocí jednoduché metody FFT. Jsou použity číslově teoretické transformace v celých číslech modulo 337, přičemž jako 8. kořen jednoty je vybráno 85. Základna 10 se používá místo báze 2 ^w pro ilustrativní účely. Schönhage – Strassen to vylepšuje pomocí negacyklických konvolucí.

Schönhage ( r ) a Strassen ( l ) hrají šachy v Oberwolfachu , 1979

Algoritmus Schönhage-Strassen je asymptoticky rychlý násobení algoritmus velkých čísel . To byl vyvinut Arnold Schönhage a Volker Strassen v roce 1971. run-time bit složitost je ve velké O notace , dva n -digit čísel. Algoritmus používá rekurzivní rychlé Fourierovy transformace v prstencích s prvky 2 ⁿ +1, což je specifický typ teoretické transformace čísel . ${\ Displaystyle O (n \ cdot \ log n \ cdot \ log \ log n)}$

Algoritmus Schönhage – Strassen byl asymptoticky nejrychlejší multiplikační metoda známá od roku 1971 do roku 2007, kdy byla oznámena nová metoda, Fürerův algoritmus , s nižší asymptotickou složitostí; Fürerův algoritmus však v současnosti dosahuje výhody pouze u astronomicky velkých hodnot a používá se pouze v podprogramech Basic Polynomial Algebra Subprograms (BPAS) (viz Galaktické algoritmy ).

V praxi začíná algoritmus Schönhage – Strassen překonávat starší metody, jako je násobení Karatsuba a Toom – Cook, pro čísla nad 2 ^{2 ¹⁵} až 2 ^{2 ¹⁷} (10 000 až 40 000 desetinných číslic). Knihovna GNU Multi-Precision Library ji používá pro hodnoty nejméně 1728 až 7808 64bitových slov (33 000 až 150 000 desetinných číslic), v závislosti na architektuře. Existuje implementace Schönhage – Strassen v jazyce Java, která ji používá nad 74 000 desetinných míst.

Aplikace algoritmu Schönhage – Strassen zahrnují matematický empirismus , jako je Great Internet Mersenne Prime Search a výpočetní aproximace π , stejně jako praktické aplikace, jako je substituce Kronecker , ve které lze účinně redukovat násobení polynomů s celočíselnými koeficienty na velká celá čísla násobení; toto v praxi používá GMP-ECM pro faktorizaci eliptické křivky Lenstra .

Podrobnosti

Tato část podrobně vysvětluje, jak je Schönhage – Strassen implementován. Je založen především na přehledu metod Crandallem a Pomerance v jejich Prime Numbers: A Computational Perspective . Tato varianta se poněkud liší od původní metody Schönhage v tom, že využívá diskrétní váženou transformaci k efektivnějšímu provádění negacyklických konvolucí . Dalším zdrojem pro detailní informace je Knuth ‚s Umění programování . Sekce začíná diskusí o stavebních blocích a končí podrobným popisem samotného algoritmu.

Konvoluce

Předpokládejme, že vynásobíme dvě čísla jako 123 a 456 pomocí dlouhého násobení se základními číslicemi B , ale bez provedení jakéhokoli přenosu. Výsledek může vypadat nějak takto:

		1	2	3
	×	4	5	6

		6	12	18
	5	10	15
4	8	12

4	13	28	27	18

Tato sekvence (4, 13, 28, 27, 18) se nazývá acyklická nebo lineární konvoluce dvou původních sekvencí (1,2,3) a (4,5,6). Jakmile máme acyklickou konvoluci dvou sekvencí, je výpočet součinu původních čísel snadný: pouze provedeme nošení (například v pravém sloupci ponecháme 8 a přidáme 1 do sloupce obsahujícího 27). V příkladu to poskytne správný produkt 56088.

Užitečné jsou dva další typy konvolucí. Předpokládejme, že vstupní sekvence mají n prvků (zde 3). Pak má acyklická konvoluce n + n −1 prvků; vezmeme -li úplně n prvků n a nejvíce vlevo n −1 prvků, vznikne cyklická konvoluce :

	28	27	18
+		4	13

	28	31	31

Provádíme -li cyklickou konvoluci, výsledek je ekvivalentní součinu vstupů mod B ⁿ - 1. V příkladu 10 ³ - 1 = 999, provádějící pokračování (28, 31, 31) poskytne 3141 a 3141 ≡ 56088 (mod 999).

Naopak, vezmeme -li n nejpravších prvků n a odečteme n -1 prvků úplně vlevo , vznikne negacyklická konvoluce :

	28	27	18
-		4	13

	28	23	5

Pokud provádíme přenesení negacyklické konvoluce, výsledek je ekvivalentní součinu vstupů mod B ⁿ + 1. V příkladu 10 ³ + 1 = 1001, provádějící pokračování (28, 23, 5) poskytne 3035 a 3035 ≡ 56088 (mod 1001). Negacyklická konvoluce může obsahovat záporná čísla, která lze eliminovat při nošení pomocí půjčování, jak se to dělá při dlouhém odčítání.

Konvoluční věta

Stejně jako ostatní metody násobení založené na rychlé Fourierově transformaci , Schönhage – Strassen závisí v zásadě na konvoluční větě , která poskytuje účinný způsob výpočtu cyklické konvoluce dvou sekvencí. Uvádí, že:

Cyklickou konvoluci dvou vektorů lze nalézt tak, že z každého z nich vezmeme diskrétní Fourierovu transformaci (DFT), vynásobíme výsledné vektory prvek po prvku a poté provedeme inverzní diskrétní Fourierovu transformaci (IDFT).

Nebo v symbolech:

CyclicConvolution ( X, Y ) = IDFT (DFT ( X ) · DFT ( Y ))

Pokud vypočítáme DFT a IDFT pomocí rychlého Fourierova transformačního algoritmu a vyvoláme náš multiplikační algoritmus rekurzivně, abychom vynásobili záznamy transformovaných vektorů DFT ( X ) a DFT ( Y ), získá se účinný algoritmus pro výpočet cyklické konvoluce.

V tomto algoritmu bude užitečnější vypočítat negacyklickou konvoluci; jak se ukazuje, mírně upravená verze konvoluční věty (viz diskrétní vážená transformace ) to může také umožnit. Předpokládejme, že vektory X a Y mají délku n , a je primitivní kořen jednoty z řádu 2 n (to je, ²ⁿ = 1 a pro všechny menší sil není 1). Potom můžeme definovat třetí vektor A , nazývaný váhový vektor , jako:

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

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

Nyní můžeme konstatovat:

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

Jinými slovy, je to stejné jako dříve, kromě toho, že vstupy se nejprve vynásobí A a výsledek se vynásobí A ⁻¹ .

Volba prstenu

Diskrétní Fourierova transformace je abstraktní operace, kterou lze provést v jakémkoli algebraickém kruhu ; obvykle se provádí v komplexních číslech, ale ve skutečnosti se provádí komplexní aritmetika s dostatečnou přesností, aby se zajistilo, že přesné výsledky pro násobení jsou pomalé a náchylné k chybám. Namísto toho budeme používat přístup na čísla teoretické transformace , která je k provedení transformace v celých mod N pro nějaké celé číslo N .

Stejně jako existují primitivní kořeny jednoty každého řádu v komplexní rovině, při jakémkoli řádu n můžeme zvolit vhodné N takové, že b je primitivní kořen jednoty řádu n v celých číslech mod N (jinými slovy, b ⁿ ≡ 1 (mod N ) a žádná menší mocnina b je ekvivalentní 1 mod N ).

Algoritmus stráví většinu času prováděním rekurzivních násobení menších čísel; s naivním algoritmem se tyto vyskytují na několika místech:

Uvnitř algoritmu rychlé Fourierovy transformace, kde je primitivní kořen jednoty b opakovaně napájen, čtvercován a vynásoben jinými hodnotami.
Když vezmeme síly primitivního kořene jednoty a k vytvoření vektoru hmotnosti A a když vynásobíme A nebo A ⁻¹ jinými vektory.
Při provádění násobení transformovaných vektorů po jednotlivých prvcích.

Klíčovým vhledem do Schönhage – Strassen je zvolit N, modul, aby byl roven 2 ⁿ + 1 pro nějaké celé číslo n, které je násobkem počtu transformovaných kusů P = 2 ^p . To má ve standardních systémech řadu výhod, které představují velká celá čísla v binární formě:

Jakoukoli hodnotu lze rychle snížit modulo 2 ⁿ + 1 pouze pomocí posunů a sčítání, jak je vysvětleno v další části .
Všechny kořeny jednoty použité transformací lze zapsat ve tvaru 2 ^k ; v důsledku toho můžeme násobit nebo dělit libovolné číslo kořenem jednoty pomocí posunu a mocnit nebo odmocnit kořen jednoty tím, že budeme operovat pouze na jeho exponentu.
Rekurzivní násobení transformovaných vektorů po jednotlivých prvcích lze provádět pomocí negacyklické konvoluce, která je rychlejší než acyklická konvoluce a již má „zdarma“ účinek snížení jejího výsledku mod 2 ⁿ + 1.

Aby rekurzivní násobení bylo pohodlné, zformulujeme Schönhage – Strassen jako specializovaný multiplikační algoritmus pro výpočet nejen součinu dvou čísel, ale součinu dvou čísel mod 2 ⁿ + 1 pro některá daná n . Nejde o ztrátu obecnosti, protože vždy je možné zvolit dostatečně velké n , takže produktový mod 2 ⁿ + 1 je jednoduše produkt.

Optimalizace směn

V průběhu algoritmu existuje mnoho případů, kdy lze násobení nebo dělení mocninou dvou (včetně všech kořenů jednoty) výhodně nahradit malým počtem posunů a sčítání. To využívá pozorování, že:

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

Všimněte si, že k -číselné číslo v základu 2 ⁿ zapsané v pozičním zápisu lze vyjádřit jako ( d _{k −1} , ..., d ₁ , d ₀ ). Představuje číslo . Všimněte si také, že pro každé d _i máme 0 ≤ d _i <2 ⁿ . ${\ Displaystyle \ scriptstyle \ sum _ {i = 0}^{k-1} d_ {i} \ cdot (2^{n})^{i}}$

To zjednodušuje redukci čísla reprezentovaného v binárním módu $2 n + 1$ : vezměte nanejvýš vpravo (nejméně významný) n bitů, odečtěte dalších n bitů, přidejte dalších n bitů a tak dále, dokud se bity nevyčerpají. Pokud výsledná hodnota stále není mezi 0 a 2 ⁿ , normalizujte ji přičtením nebo odečtením násobku modulu $2 n + 1$ . Například pokud n = 3 (a modul je tedy 2 ³ +1 = 9) a zmenšené číslo je 656, máme:

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

Kromě toho je možné provádět velmi velké posuny, aniž byste museli konstruovat posunutý výsledek. Předpokládejme, že máme číslo A mezi 0 a 2 ⁿ , a chceme jej vynásobit 2 ^k . Dělení K o n najdeme K = qn + r s r < n . Z toho vyplývá, že:

A (2 ^k ) = A (2 ^{qn + r} ) = A [(2 ⁿ ) ^q (2 ^r )] ≡ (−1) ^q Shift Left (A, r ) (mod 2 ⁿ + 1).

Protože A je ≤ 2 ⁿ a r < n , má posun vlevo r maximálně 2 n −1 bitů, a proto je zapotřebí pouze jeden posun a odečtení (následované normalizací).

Nakonec dělíme 2 ^k , pozorujeme, že kvadratura první ekvivalence nad výnosy:

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

Proto,

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

Algoritmus

Algoritmus sleduje rozdělení, vyhodnocení (dopředný FFT), bodové násobení, interpolaci (inverzní FFT) a kombinuje fáze podobné metodám Karatsuba a Toom-Cook.

Vzhledem k vstupním číslům x a y a celému číslu N následující algoritmus vypočítá součin xy mod 2 ^N + 1. Za předpokladu, že N je dostatečně velké, je to jednoduše produkt.

Rozdělit každý číslo vstupu do vektorů X a Y 2 ^k- částí každý, kde 2 ^K dělení N . (např. 12345678 → (12, 34, 56, 78)).
Aby bylo možné dosáhnout pokroku, je nutné použít menší N pro rekurzivní násobení. Za tímto účelem zvolte N jako nejmenší celé číslo alespoň 2 N /2 ^k + k a dělitelné 2 ^k .
Vypočítejte součin X a Y mod 2 ^N + 1 pomocí negacyklické konvoluce:
1. Násobte X a Y váhovým vektorem A pomocí posunů (posuňte j. Vstup doleva o jn /2 ^k ).
2. Vypočtěte DFT X a Y Pomocí numerické-teoretické FFT (plnit všechny násobení pomocí směny, pro 2 ^K -tý kořen jednoty, použití 2 ^{2 n / 2 ^K} ).
3. Rekurzivně použijte tento algoritmus ke znásobení odpovídajících prvků transformovaných X a Y.
4. Vypočítejte IDFT výsledného vektoru, abyste získali výsledný vektor C (proveďte všechna násobení pomocí posunů). To odpovídá fázi interpolace.
5. Výsledný vektor C vynásobte A ⁻¹ pomocí posunů.
6. Upravte značky: některé prvky výsledku mohou být negativní. Vypočítáme největší možnou kladnou hodnotu pro j. Prvek C, $(j + 1) 2 2 N /2 k$ , a pokud ji překročí, odečteme modul 2 ⁿ + 1.
Nakonec proveďte nesoucí režim 2 ^N + 1, abyste získali konečný výsledek.

Optimální počet kusů, na které se má vstup rozdělit, je úměrný , kde N je počet vstupních bitů, a tímto nastavením se dosáhne doby běhu O ( N log N log log N ), takže parametr k by měl být odpovídajícím způsobem nastaven. V praxi je nastaven empiricky na základě vstupních velikostí a architektury, obvykle na hodnotu mezi 4 a 16. ${\ displaystyle {\ sqrt {N}}}$

V kroku 2 se používá pozorování, které:

Každý prvek vstupních vektorů má nejvýše n /2 ^k bitů;
Součin jakýchkoli dvou vstupních vektorových prvků má nejvýše 2 n /2 ^k bitů;
Každý prvek konvoluce je součtem nejvýše 2 ^k takových produktů, a nesmí tedy překročit 2 n /2 ^k + k bitů.
n musí být dělitelné 2 ^k, aby bylo zajištěno, že v rekurzivních hovorech platí v kroku 1 podmínka „2 ^k dělí N “.

Optimalizace

Tato část vysvětluje řadu důležitých praktických optimalizací, které byly vzaty v úvahu při implementaci Schönhage – Strassen do reálných systémů. Je založen především na práci Gaudryho, Kruppy a Zimmermanna z roku 2007, která popisuje vylepšení vícepřesné knihovny GNU .

Pod určitým mezním bodem je efektivnější provádět rekurzivní násobení pomocí jiných algoritmů, například multiplikace Toom – Cook . Výsledky musí být sníženy mod 2 ⁿ + 1, což lze provést efektivně, jak je vysvětleno výše v optimalizaci směn s posuny a sčítání/odčítání.

Výpočet IDFT zahrnuje dělení každého záznamu primitivním kořenem jednoty 2 ^{2 n /2 ^k} , což je operace, která je často kombinována s vynásobením vektoru A ⁻¹ poté, protože oba zahrnují dělení mocninou dvou.

V systému, kde je velký počet reprezentován jako pole 2 ^w -bitových slov, je užitečné zajistit, aby velikost vektoru 2 ^k byla také násobkem bitů na slovo výběrem k ≥ w (např. Zvolit k ≥ 5 na 32bitový počítač a k ≥ 6 na 64bitovém počítači); to umožňuje rozdělit vstupy na kousky bez bitových posunů a poskytuje jednotnou reprezentaci pro hodnoty mod 2 ⁿ + 1, kde vysoké slovo může být pouze nula nebo jedna.

Normalizace zahrnuje sčítání nebo odčítání modulu 2 ⁿ + 1; tato hodnota má nastaveny pouze dva bity, což znamená, že to lze provádět v průměru ve stálém čase pomocí specializované operace.

Iterační algoritmy FFT, jako je algoritmus Cooley -Tukey FFT , přestože se často používají pro FFT na vektorech komplexních čísel, mají tendenci vykazovat velmi špatnou lokalitu mezipaměti s velkými vektorovými položkami používanými v Schönhage – Strassen. Přímočará rekurzivní, nemístná implementace FFT je úspěšnější, přičemž všechny operace zapadají do mezipaměti za určitý bod v hloubce volání, ale přesto suboptimálně využívá mezipaměť ve vyšších hloubkách volání. Gaudry, Kruppa a Zimmerman použili techniku kombinující Baileyho 4krokový algoritmus s vyššími radixovými transformacemi, které kombinují více rekurzivních kroků. Míchají také fáze a jdou co nejdále do algoritmu na každém prvku vektoru, než přejdou na další.

„Druhá odmocnina ze 2 triků“, kterou poprvé popsal Schönhage, je, že za předpokladu, že k ≥ 2, 2 ^{3 n /4} −2 ^{n /4} je druhá odmocnina ze 2 mod 2 ⁿ +1, a tak 4 n -tý kořen jednoty (od 2 ^{2 n} ≡ 1). To umožňuje prodloužit délku transformace z 2 ^k na 2 ^{k + 1} .

Nakonec si autoři dávají pozor na výběr správné hodnoty k pro různé rozsahy vstupních čísel s tím, že optimální hodnota k se může mezi stejnými hodnotami několikrát pohybovat tam a zpět, jak se velikost vstupu zvyšuje.

Languages

In other projects