Schönhage – Strassenův algoritmus - Schönhage–Strassen algorithm

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

The Schönhage – Strassenův algoritmus je asymptoticky rychlý multiplikační algoritmus pro velké celá čísla. Byl vyvinut společností Arnold Schönhage a Volker Strassen v roce 1971.^[1] Runtime trochu složitost je v Velká O notace, ${ displaystyle O (n cdot log n cdot log log n)}$ pro dva n-místná čísla. Algoritmus používá rekurzivní Rychlé Fourierovy transformace v prsteny s 2ⁿ+1 prvky, konkrétní typ číselná teoretická transformace.

Algoritmus Schönhage – Strassen byl asymptoticky nejrychlejší multiplikační metoda známá od roku 1971 do roku 2007, kdy nová metoda, Fürerův algoritmus, bylo oznámeno s nižší asymptotickou složitostí;^[2] Fürerův algoritmus však v současné době dosahuje výhody pouze pro astronomicky velké hodnoty a v praxi se nepoužívá (viz Galaktické algoritmy ).

V praxi začíná algoritmus Schönhage – Strassen překonávat starší metody, jako např Karatsuba a Násobení Toom – Cook pro čísla nad 2²¹⁵ až 2²¹⁷ (10 000 až 40 000 desetinných míst).^[3][4][5] The GNU Multi-Precision Library používá jej pro hodnoty alespoň 1728 až 7808 64bitových slov (33 000 až 150 000 desetinných míst), v závislosti na architektuře.^[6] Existuje implementace Schönhage – Strassen v jazyce Java, která ji používá nad 74 000 desetinných míst.^[7]

Mezi aplikace patří Schönhage – Strassenův algoritmus matematický empirismus, tak jako Skvělé internetové vyhledávání Mersenne Prime a výpočetní technika aproximace π, stejně jako praktické aplikace jako Střídání v týmu Kronecker, ve kterém lze násobení polynomů s celočíselnými koeficienty efektivně snížit na velké celočíselné násobení; toto v praxi používá GMP-ECM pro Faktorizace eliptické křivky Lenstra.^[8]

Detaily

Tato část podrobně vysvětluje, jak je Schönhage – Strassen implementován. Je založen především na přehledu metody Crandalla a Pomerance v jejich Prvočísla: Výpočetní perspektiva.^[9] Tato varianta se od Schönhageovy původní metody poněkud liší tím, že využívá diskrétní vážená transformace vystupovat negacyklické konvoluce efektivněji. Dalším zdrojem podrobných informací je Knuth je Umění počítačového programování.^[10] Sekce začíná diskusí o stavebních blocích a končí podrobným popisem samotného algoritmu.

Závity

Předpokládejme, že vynásobíme dvě čísla jako 123 a 456 pomocí dlouhé násobení základnou B číslic, ale bez jakéhokoli přenášení. Výsledek může vypadat asi takto:

Tato posloupnost (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í, výpočet produktu původních čísel je snadný: pouze provedeme přenášení (například ve sloupci úplně vpravo jsme ponechali 8 a přidali 1 do sloupce obsahujícího 27). V tomto příkladu se získá správný produkt 56088.

Budou užitečné dva další typy konvolucí. Předpokládejme, že vstupní sekvence mají n prvky (zde 3). Pak má acyklická konvoluce n+n-1 prvky; vezmeme-li úplně doprava n prvky a přidat úplně vlevo n−1 prvků, výsledkem je cyklická konvoluce:

Pokud provádíme cyklickou konvoluci, výsledek je ekvivalentní součinu vstupů mod Bⁿ - 1. V příkladu 10³ - 1 = 999, při provádění (28, 31, 31) se získá 3141 a 3141 ≡ 56088 (mod 999).

Naopak, vezmeme-li úplně doprava n prvky a odčítat úplně vlevo n−1 prvků, výsledkem je negacyklická konvoluce:

Pokud provádíme nesení negacyklické konvoluce, výsledek je ekvivalentní k součinu vstupů mod Bⁿ + 1. V příkladu 10³ + 1 = 1001, při provádění pokračování (28, 23, 5) se získá 3035 a 3035 ≡ 56088 (mod 1001). Negacyklická konvoluce může obsahovat záporná čísla, která lze během přenášení pomocí výpůjčky eliminovat, jak je tomu u dlouhého odečítání.

Konvoluční věta

Jako ostatní multiplikační metody založené na rychlé Fourierově transformaci, Schönhage – Strassen zásadně závisí na konvoluční věta, který poskytuje efektivní způsob výpočtu cyklické konvoluce dvou sekvencí. Uvádí, že:

Cyklickou konvoluci dvou vektorů lze zjistit pomocí diskrétní Fourierova transformace (DFT) každého z nich, vynásobení výsledného vektorového prvku prvkem a následná inverzní diskrétní Fourierova transformace (IDFT).

Nebo v symbolech:

Cyklická konverze (X, Y) = IDFT (DFT (X) · DFT (Y))

Pokud vypočítáme DFT a IDFT pomocí a rychlá Fourierova transformace algoritmus a rekurzivně vyvolat náš multiplikační algoritmus, aby se vynásobily položky transformovaných vektorů DFT (X) a DFT (Y), čímž se získá efektivní algoritmus pro výpočet cyklické konvoluce.

V tomto algoritmu bude užitečnější vypočítat negacyclic konvoluce; jak se ukázalo, mírně upravená verze věty o konvoluci (viz diskrétní vážená transformace ) lze také povolit. Předpokládejme, že vektory X a Y mají délku n, a A je primitivní kořen jednoty z objednat 2n (to znamená, A²ⁿ = 1 a A pro všechny menší síly není 1). Pak můžeme definovat třetí vektor A, volal hmotnostní vektor, tak jako:

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

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

Nyní můžeme konstatovat:

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

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

Volba prstenu

Diskrétní Fourierova transformace je abstraktní operace, kterou lze provést v jakékoli algebraický prsten; obvykle se provádí v komplexních číslech, ale ve skutečnosti provádění komplexní aritmetiky s dostatečnou přesností, aby se zajistilo, že přesné výsledky pro násobení jsou pomalé a náchylné k chybám. Místo toho použijeme přístup číselná teoretická transformace, což je provedení transformace v režimu celých čísel N pro celé číslo N.

Stejně jako existují primitivní kořeny jednoty každého řádu v komplexní rovině, dané libovolným řádem n můžeme vybrat 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ší síla b odpovídá 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 mnoha místech:

1. Uvnitř algoritmu rychlé Fourierovy transformace, kde je primitivní kořen jednoty b je opakovaně napájen, čtvercový a vynásoben dalšími hodnotami.
2. Při převzetí pravomocí primitivního kořene jednoty A k vytvoření váhového vektoru A a při vynásobení A nebo A⁻¹ jinými vektory.
3. Při provádění násobení elementů po elementech transformovaných vektorů.

Klíčovým poznatkem Schönhage – Strassena je volba N, modulu, který se bude rovnat 2ⁿ + 1 pro celé číslo n to je násobek počtu kusů P=2^str být transformován. To má řadu výhod ve standardních systémech, které představují velká celá čísla v binární podobě:

• Libovolnou hodnotu lze rychle snížit modulo 2ⁿ + 1 pouze s použitím směn a přidání, jak je vysvětleno v další část.
• Všechny kořeny jednoty použité transformací lze zapsat ve formě 2^k; následně můžeme libovolné číslo vynásobit nebo rozdělit kořenem jednoty pomocí posunu a mocninu nebo druhou mocninu kořene jednoty pracovat pouze na jejím exponentu.
• Rekurzivní multiplikace transformovaných vektorů po prvcích lze provádět pomocí negacyklické konvoluce, která je rychlejší než acyklická konvoluce a již má „zdarma“ účinek snižování jeho výsledku mod 2ⁿ + 1.

Aby rekurzivní multiplikace byly pohodlné, vytvoříme 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. To není ztráta obecnosti, protože člověk si může vždy vybrat n dostatečně velký, aby produkt modⁿ + 1 je jednoduše produkt.

Optimalizace řazení

V průběhu algoritmu existuje mnoho případů, kdy lze násobení nebo dělení mocí 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 a k-místné číslo v základně 2ⁿ napsáno v poziční notace lze vyjádřit jako (d_k-1,...,d₁,d₀). Představuje číslo ${ displaystyle scriptstyle sum _ {i = 0} ^ {k-1} d_ {i} cdot (2 ^ {n}) ^ {i}}$. Všimněte si také, že pro každého d_i máme 0≤d_i < 2ⁿ.

Díky tomu je snadné snížit počet zastoupených v binárním režimu 2ⁿ + 1: vzít úplně vpravo (nejméně významný) n bity, odečtěte další n bity, přidejte další n bity atd., dokud se bity nevyčerpají. Pokud výsledná hodnota stále není mezi 0 a 2ⁿ, normalizujte jej přidáním nebo odečtením násobku modulu 2ⁿ + 1. Například pokud n= 3 (a modul je tedy 2³+1 = 9) a snížený počet je 656, máme:

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

Kromě toho je možné provádět velmi velké posuny, aniž by došlo ke konstrukci posunutého výsledku. Předpokládejme, že máme číslo A mezi 0 a 2ⁿ, a chcete jej vynásobit 2^k. Dělení k podle n shledáváme 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(Posun doleva r) (mod 2ⁿ + 1).

Protože A je ≤ 2ⁿ a r < n, Posun doleva r má maximálně 2n−1 bitů, a proto je potřeba pouze jeden posun a odčítání (následované normalizací).

Nakonec vydělte 2^k, pozorujte, že druhou mocninou první ekvivalence se získá:

2²ⁿ ≡ 1 (mod 2ⁿ + 1)

Proto,

A / 2^k = A (2^−k) ≡ A (2^{2n − k}) = Posun doleva (2n − k) (mod 2ⁿ + 1).

Algoritmus

Algoritmus sleduje rozdělení, vyhodnocení (vpřed FFT), bodové násobení, interpolaci (inverzní FFT) a kombinování fází podobných metodám Karatsuba a Toom-Cook.

Daná vstupní čísla X a ya celé číslo N, následující algoritmus spočítá produkt xy mod 2^N + 1. Pokud je N dostatečně velké, jedná se jednoduše o produkt.

1. Rozdělte každé vstupní číslo na vektory X a Y 2^k každý díl, kde 2^k rozděluje N. (např. 12345678 → (12, 34, 56, 78)).
2. Aby bylo možné dosáhnout pokroku, je nutné použít menší N pro rekurzivní násobení. Za tímto účelem vyberte n jako nejmenší celé číslo alespoň 2N/2^k + k a dělitelné 2^k.
3. Vypočítejte součin X a Y mod 2ⁿ +1 pomocí negacyklické konvoluce:
   1. Vynásobte X a Y váhovým vektorem A pomocí posunů (posuňte jth položka vlevo od jn/2^k).
   2. Vypočítejte DFT X a Y pomocí FFT s teoretickým počtem (proveďte všechny násobení pomocí posunů; pro 2^k-tý kořen jednoty, použijte 2^2n/2k).
   3. Rekurzivně použijte tento algoritmus k ná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šechny násobení pomocí posunů). To odpovídá interpolační fázi.
   5. Vynásobte výsledný vektor C číslem A.⁻¹ pomocí směn.
   6. Upravte značky: některé prvky výsledku mohou být negativní. Vypočítáme největší možnou kladnou hodnotu pro jth prvek C, (j + 1)2^2N/2k, a pokud to překročí, odečteme modul 2ⁿ + 1.
4. Nakonec proveďte režim nošení 2^N +1 k získání konečného výsledku.

Optimální počet kusů pro rozdělení vstupu je úměrný ${ displaystyle { sqrt {N}}}$, kde N je počet vstupních bitů a toto nastavení dosahuje doby běhu O (N log N log log N),^[1][9] tedy parametr k by měl být odpovídajícím způsobem nastaven. V praxi je empiricky nastaven na základě vstupních velikostí a architektury, obvykle na hodnotu mezi 4 a 16.^[8]

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

• Každý prvek vstupních vektorů má maximálně n/2^k bity;
• Produkt libovolných dvou vstupních vektorových prvků má maximálně 2n/2^k bity;
• Každý prvek konvoluce je součtem maximálně 2^k takové výrobky, a proto nesmí překročit 2n/2^k + k bity.
• n musí být dělitelné 2^k zajistit, aby v rekurzivních voláních byla podmínka „2^k rozděluje N"platí v kroku 1.

Optimalizace

Tato část vysvětluje řadu důležitých praktických optimalizací, které byly zohledněny při implementaci Schönhage – Strassen ve skutečných systémech. Je založen především na díle Gaudryho, Kruppy a Zimmermanna z roku 2007, které popisuje vylepšení GNU Multi-Precision Library.^[8]

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

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

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

Normalizace zahrnuje přidání nebo odečtení modulu 2ⁿ + 1; tato hodnota má nastaveny pouze dva bity, což znamená, že to lze provést průměrně v konstantním čase se specializovanou operací.

Iterativní algoritmy FFT, jako například Cooley – Tukey FFT algoritmus, i když se často používají pro FFT na vektorech komplexních čísel, mají tendenci vykazovat velmi špatnou mezipaměť lokalita s velkými vektorovými položkami použitými v Schönhage – Strassen. Přímočará rekurzivní implementace FFT, která není na místě, je úspěšnější, protože všechny operace zapadají do mezipaměti mimo určitý bod v hloubce volání, ale přesto umožňují suboptimální využití mezipaměti 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, jdou co nejdále do algoritmu na každém prvku vektoru, než se přesunou k dalšímu.

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

Nakonec jsou autoři opatrní při výběru správné hodnoty k pro různé rozsahy vstupních čísel s tím, že optimální hodnota k se může s rostoucí velikostí vstupu několikrát pohybovat tam a zpět mezi stejnými hodnotami.

Reference