DIPLOMOVÁ PRÁCE. Hana Straková Testování perfektních mocnin

Univerzita Karlova v Praze Matematicko-fyzikální fakulta

DIPLOMOVÁ PRÁCE

Hana Straková Testování perfektních mocnin Katedra algebry

Vedoucí diplomové práce: RNDr. David Stanovský, Ph D. Studijní program: Informatika Studijní obor: Softwarové systémy

2008

Děkuji RNDr. Davidu Stanovskému, Ph D. za odborné vedení práce a za cenné připomínky a podněty. Dále děkuji Braňovi Bošanskému za pomoc s úpravou textu.

Prohlašuji, že jsem svou diplomovou práci napsala samostatně a výhradně s použitím citovaných pramenů. Souhlasím se zapůjčováním práce a jejím zveřejňováním. V Praze dne 15. 4. 2008

Hana Straková

2

Obsah 1 Úvod 1.1 Perfektní mocniny . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Značení a definice . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Struktura práce . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Algoritmy Bacha & Sorensona 2.1 Algoritmus A - binární vyhledávání a metoda . . . . . . . . . . . . . . . . . . 2.2 Algoritmus B - sítové testy . . . . . . . 2.3 Algoritmus C - zkusmé dělení . . . . .

10 upravená Newtonova . . . . . . . . . . . . . 10 . . . . . . . . . . . . . 11 . . . . . . . . . . . . . 12

3 Algoritmus Daniela J. Bernsteina 3.1 Operace na číslech s plovoucí řádovou čárkou 3.1.1 Sčítání, odčítání, násobení . . . . . . 3.1.2 Dělení, zaokrouhlování, umocňování . 3.2 Bernsteinův algoritmus . . . . . . . . . . . . 3.2.1 Funkce Nrootb . . . . . . . . . . . . . 3.2.2 Funkce Ppower . . . . . . . . . . . . 3.2.3 Funkce Compare . . . . . . . . . . . 3.2.4 Hlavní cyklus Bernsteinova algoritmu 4 Knihovna GMP 4.1 O knihovně GMP . 4.2 Aritmetické operace 4.2.1 Násobení . . 4.2.2 Dělení . . . 4.2.3 Umocňování

. . . . . . . . . . v knihovně GMP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . .

6 6 8 9

. . . . .

. . . . .

. . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . .

15 15 15 16 17 17 19 20 21

. . . . .

22 22 23 23 24 24

5 Implementace algoritmů 25 5.1 Implementace algoritmů Bacha & Sorensona . . . . . . . . . . 27 5.2 Implementace Bernsteinova algoritmu . . . . . . . . . . . . . . 29 3

6 Složitost algoritmů a výsledky měření 6.1 Algoritmus A . . . . . . . . . . . . . . . . . . . . 6.1.1 Asymptotická složitost . . . . . . . . . . . 6.1.2 Výsledky měření . . . . . . . . . . . . . . 6.2 Algoritmus B . . . . . . . . . . . . . . . . . . . . 6.2.1 Čas běhu programu v průměrném případě 6.2.2 Výsledky měření . . . . . . . . . . . . . . 6.3 Algoritmus C . . . . . . . . . . . . . . . . . . . . 6.3.1 Čas běhu programu v průměrném případě 6.3.2 Výsledky měření . . . . . . . . . . . . . . 6.3.3 Algoritmy B a C a mocniny prvočísla . . . 6.4 Bernsteinův algoritmus . . . . . . . . . . . . . . . 6.4.1 Složitost . . . . . . . . . . . . . . . . . . . 6.4.2 Výsledky měření . . . . . . . . . . . . . . 6.5 Funkce mpz_perfect_power_p . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

31 32 32 34 35 35 39 40 40 42 43 45 45 46 47

7 Závěr

50

A Knihovna GMP

55

B Implementace algoritmů 59 B.1 Implementace algoritmů Bacha & Sorensona . . . . . . . . . . 60 B.2 Implementace Bernsteinova algoritmu . . . . . . . . . . . . . . 62

4

Název práce: Testování perfektních mocnin Autor: Hana Straková Katedra(ústav): Katedra algebry Vedoucí diplomové práce: RNDr. David Stanovský, Ph D. e-mail vedoucího: [email protected]ff.cuni.cz Abstrakt: Perfektní mocnina je takové přirozené číslo, které lze zapsat jako netriviální mocninu jiného přirozeného čísla. Přirozené číslo n tedy nazveme perfektní mocninou, pokud existují přirozená čísla x a k, obě větší než 2, taková, že n = xk . Testy, zda je číslo perfektní mocnina, jsou důležité pro faktorizaci čísel nebo testy prvočíselnosti, které samy neumějí rozlišit, zda je vstupní číslo mocnina jiného čísla. Tato práce pojednává o dvou algoritmech na testování perfektních mocnin, prvním od E.Bacha & J. Sorensona a druhém od Daniela J. Bernsteina. Cílem práce je implementace algoritmů v jazyce C s pomocí knihovny GMP, srovnání teoretických výsledků pro jednotlivé algoritmy s praktickými měřeními a porovnání algoritmů mezi sebou z hlediska teoretického přístupu a rychlosti. Klíčová slova: perfektní mocnina, výpočetní teorie čísel, knihovna GMP Title: Perfect power testing Author: Hana Straková Department: Department of algebra Supervisor: RNDr. David Stanovský, Ph D. Supervisor’s e-mail address: [email protected]ff.cuni.cz Abstract: A positive integer n is a perfect power if there exist integers x and k, both at least 2, such that n = xk . Perfect power testing is important as preprocessing for number factorization and prime number testing, because many algorithms for that are not able to distinguish between prime number and power of prime number, so it is necessary to test it by perfect power tests. This thesis includes comparison of two algorithms for perfect power testing, one by Daniel J. Bernstein and the other by E. Bach & J. Sorenson. The goal is to implement described algorithms in C language with GMP library for multiple-precision arithmetics, to compare the theoretical results and running times of implemented algorithms. Keywords: perfect power, computational number theory, GMP library

5

Kapitola 1 Úvod 1.1

Perfektní mocniny

Kryptografie a informační bezpečnost nabývají v současné době stále více na významu. Zatímco dříve byly slova jako šifra, klíč šifry, certifikáty vcelku neznámými pojmy, dneska se s nimi setkáváme prakticky denně, především ve spojení s informačními technologiemi. Méně známé už jsou detaily principů, na kterých jsou šifry založeny, informace o průběhu šifrování a dešifrování dat a zabezpečení informací. Jednu z velkých rolí hrají v kryptografii generátory náhodných čísel a testy prvočíselnosti. Na vlastnostech prvočísel stojí mnoho matematických teorií a prvočísla jsou základem pro mnoho šifer. Zmiňme alespoň nejpoužívanější a nejznámější asymetrickou šifru RSA. Důležitou část veřejného klíče této šifry tvoří číslo, které vznikne vynásobením dvou stejně velkých prvočísel (tak velkých, že výsledné číslo má zpravidla velikost 1024 nebo 2048 bitů). Znalost těchto dvou prvočísel tvoří pak soukromý klíč šifry. Prvočísla potřebná pro šifry se získávají tak, že se vygenerují libovolná čísla s požadovaným počtem cifer (pseudo)náhodným generátorem, a pak se na tato vygenerovaná čísla spustí test prvočíselnosti. Od testů prvočíselnosti očekáváme, že nám zjistí, zda vstupní číslo je nebo není prvočíslo. Problém nastává tehdy, pokud test prvočíselnosti není schopen rozlišit mezi prvočíslem a mocninou prvočísla, čili pozitivní odpověď je, že vstupní číslo je prvočíslo nebo mocnina prvočísla. Tak je tomu např. u moderního deterministického testu prvočíselnosti AKS (základní informace o testu AKS lze nalézt např. v [11]). V těchto případech se využije testování, zda vstupní číslo není mocninou nějakého jiného přirozeného čísla, tedy zda není perfektní mocninou. Přirozené číslo n nazýváme perfektní mocninou, pokud existují dvě přirozená čísla x a k, obě větší nebo rovna 2, taková, že platí vztah n = xk . Testy perfektních mocnin se používají nejen v testech prvočíselnosti, ale 6

i v algoritmech na faktorizaci čísel. Jeden z příkladů, kdy je potřeba před faktorizací čísla provést test, zda je číslo perfektní mocnina, je algoritmus faktorizující vstupní číslo n pomocí kongruence x2 ≡ y 2 (mod n). Po vyřešení kongruence x2 ≡ y 2 (mod n) a nalezení jejího řešení x a y, kde x je různé od ±y, je spočítán největší společný dělitel čísel n a x−y, čímž dostaneme faktor čísla n. Taková řešení kongruence x a y budou existovat pouze v případě, že n není mocnina prvočísla. Je tedy nutno otestovat n, zda není perfektní mocninou. Pokud ne, můžeme hledat řešení kongruence, a pokud ano, tak máme rovnou nalezenu faktorizaci vstupního čísla n. Tato práce zkoumá čtyři algoritmy na testování perfektních mocnin. První tři algoritmy (označované jako algoritmus A, algoritmus B a algoritmus C) jsou popsány v práci Sieve Algorithms for Perfect Power Testing od Erica Bacha a Jonathana Sorensona [1] a čtvrtý algoritmus (označovaný jako Bernsteinův algoritmus) je popsán v práci Detecting Perfect Powers in Essentially Linear Time od Daniela J. Bernsteina [2]. Testy perfektních mocnin lze obecně rozdělit na tři základní druhy: detekční, dekompoziční a klasifikační. Detekční algoritmus je takový, který pro vstupní přirozené číslo n zjistí, zda n je nebo není perfektní mocnina. Dekompoziční algoritmus kromě toho, že zjistí, zda dané n je nebo není perfektní mocnina, najde v případě pozitivní odpovědi takové x, že xk = n. Klasifikační algoritmus se liší od dekompozičního v tom, že exponent k nalezený algoritmem je maximální. Algoritmy, kterými se zabýváme v této práci, spadají do kategorie dekompozičních algoritmů. Cílem práce je srovnání dvou odlišných přístupů k testování perfektních mocnin, implementace algoritmů v jazyce C s pomocí knihovny GMP (GNU Multiple Precision Arithmetic Library), změření reálného času běhu programů a srovnání výsledků s teoretickými výpočty. Součástí práce je i zkoumání chování programů, pokud jsou vstupními čísly pouze mocniny prvočísel. V rámci práce bylo rovněž provedeno srovnání sledovaných algoritmů s funkcí pro testování perfektních mocnin implementovanou v knihovně GMP. Bach & Sorenson ve svém článku předpokládali použití klasického algoritmu pro násobení, který má kvadratickou složitost vzhledem k velikosti násobených čísel. Implementace algoritmů v této práci ovšem používají k násobení čísel Karatsubův algoritmus s asymptotickou složitostí O(log n1,585 ) pro vstupní číslo n (více informací o Karatsubově algoritmu např. na [12] nebo v [5]). Důvod je takový, že v dnešní kryptografii se nejvíce používají čísla do velikosti 1 000 cifer v desítkovém zápisu a pro čísla této velikosti je Karatsubův algoritmus nejefektivnější algoritmus pro násobení. V kapitole 6 jsou mimo praktických výsledků uvedeny i důkazy asymptotických složitostí algoritmů Bacha & Sorensona pro případ použití Karatsubova algoritmu.

7

1.2

Značení a definice

V této části uvedeme značení, definice a věty používané v celé práci. Označení log n je použito pro logaritmus o základu 2, v ostatních případech je základ explicitně napsán. gcd(x, y) označuje největšího společného dělitele čísel x a y. V algoritmech Bacha & Sorensona je použit symbol k, který má následující význam: pe k n znamená, že platí pe |n (pe dělí n beze zbytku), ale neplatí pe+1 |n. U některých vět a lemmat je uvedeno označení (ERH). Tím je zdůrazněn fakt, že pro správnost věty či lemmatu je předpokládána platnost Rozšířené Riemannovy Hypotézy (Extended Riemann Hypothesis) (více informací na [10]). Pojmem „klasická arimetika“ označujeme aritmetické operace s následující asymptotickou časovou složitostí. Předpokládáme, že a, b, m jsou přirozená čísla. Pak pro výpočet a + b, a − b potřebujeme čas O(log a + log b) a pro výpočet a · b, ⌊a/b⌋, a mod b potřebujeme čas O(log a log b). Ze složitostí základních operací plynou i složitosti dalších operací, jako umocňování, kde pro výpočet c = ab potřebujeme čas O(log2 c) a modulární umocňování ab mod m, pro které potřebujeme čas O(log a log m + log b log2 m). V důkazech složitostí algoritmů v kapitole 6 budeme používat následující věty: Věta o prvočíslech (Prime Number Theorem) se týká hustoty prvočísel, čili jak jsou prvočísla distribuována. To, co budeme v důkazech složitosti algoritmů potřebovat, je počet prvočísel menších než nějaké číslo x. Z věty o prvočíslech plyne: x π(x) ∼ (1.1) ln x kde π(x) je počet prvočísel p ≤ x. Z Dirichletovy věty, někdy také nazývané „Dirichletova věta o prvočíslech“ (Dirichlet’s Prime Number Theorem), použijeme skutečnost, že v aritmetické posloupnosti čísel ve tvaru an + b, pro n = 1, 2, ..., je nekonečně mnoho prvočísel, za předpokladu, že gcd(a, b) = 1. V důkazech v kapitole 6 použijeme Čebyševovu funkci (viz [13]) k odhadu X log p ∼ log n p≤log n

8

Pro pochopení sítových testů v algoritmech Bacha & Sorensona, je potřeba malá Fermatova věta (Fermat’s little theorem), která tvrdí, že platí ap−1 ≡ 1

(mod p)

kde p je prvočíslo, a je celé číslo a p nedělí a.

1.3

Struktura práce

V kapitolách 2 a 3 jsou popsány studované algoritmy: kapitola 2 se zabývá třemi algoritmy Bacha & Sorensona a kapitola 3 je věnována Bernsteinovu algoritmu. K implementaci všech algoritmů byla použita knihovna GMP. Práce s touto knihovnou a její nastavení během sledování chování programů je popsáno v kapitole 4 a v dodatku A. V kapitole 4 lze nalézt informace o nastavení během měření a popis implementací aritmetických operací v této knihovně. V dodatku A jsou pak uvedeny podrobnosti k instalaci, způsob nastavování parametrů a zvyšování výkonu knihovny pomocí programu tuneup. Kapitola 5 a dodatek B jsou věnovány implementaci jednotlivých algoritmů. V kapitole 5 jsou zmíněny funkce použité v programech, princip generování čísel a prvočísel v programech, v dodatku B jsou podrobně popsány jednotlivé parametry funkcí včetně datových typů. Kapitola 6 analyzuje teoretickou složitost algoritmů a uvádí naměřené výsledky a jejich porovnání s teorií. V této kapitole jsou také dokázány složitosti algoritmů Bacha & Sorensona pro případ použití Karatsubova algoritmu pro násobení. Jsou tu uvedeny naměřené časy běhů jednotlivých programů se sledovanými algoritmy a pro některé algoritmy i výsledky pro případ, že na vstupu programu byly pouze mocniny prvočísel. Zde uvedeny výsledky testování funkce mpz_perfect_power_p, což je funkce pro testování perfektních mocnin implementovaná v knihovně GMP. Kapitola 7 ukazuje srovnání všech studovaných algoritmů z hlediska teoretických přístupů, teoretických časových náročností a naměřených rychlostí.

9

Kapitola 2 Algoritmy Bacha & Sorensona V článku Sieve Algorithms for Perfect Power Testing od Erica Bacha a Jonathana Sorensona [1] autoři vycházejí z nejjednoduššího algoritmu na testování perfektních mocnin, který postupně vylepšují přidáváním různých testů, čímž zlepšují asymptotickou složitost algoritmu. V této kapitole představíme jednotlivé algoritmy. Stanovení jejich složitostí je pak věnována kapitola 6.

2.1

Algoritmus A - binární vyhledávání a upravená Newtonova metoda

Nejjednodušší algoritmus na testování perfektních mocnin spočívá v tom, že jsou počítány všechny možné odmocniny vstupního čísla n a zkoumá se, zda odmocnina z čísla n není přesně rovna nějakému přirozenému číslu. Pokud je takové přirozené číslo nalezeno, je vstupní číslo perfektní mocnina. První jednoduché zlepšení spočívá v tom, že nemusíme kontrolovat všechny odmocniny vstupního čísla, ale stačí zkoumat pouze možné prvočíselné odmocniny. Podíváme se podrobněji na hledání odmocniny vstupního čísla (v algoritmu A označena jako x). Popíšeme dvě verze algoritmu A. Jednodušší verze algoritmu používá k nalezení hodnoty x binární vyhledávání. Binární vyhledávání, někdy také názýváno půlení intervalů, je jedna z vyhledávacích metod. V našem případě hledáme v intervalu přirozených čísel hodnotu x takovou, že x je p-tá odmocnina vstupního čísla n pro právě prověřovaný prvočíselný exponent p. Pro pevně zvolené p postupujeme následovně. Najdeme hodnotu ležící ve středu aktuálně prohledáváného intervalu, kterou umocníme na p-tou, a výsledek pak porovnáme se vstupním číslem n. Podle toho, zda je mocnina menší nebo větší než n, pokračujeme s polovinou intervalu ležící nalevo nebo napravo od středu intervalu. V každém kroku se tedy 10

délka prohledávaného intervalu zmenší na polovinu. AlgoritmusA(n) 1 for each p prvočíslo, p ≤ log n do 2 spočti x = ⌊n1/p ⌋; 3 if (n = xp ) then 4 n je perfektní p-tá mocnina; 5 stop; 6 n není perfektní mocnina;

Druhá verze algoritmu A používá k hledání odmocniny tzv. upravenou Newtonovu metodu. Newtonova metoda slouží k nalezení kořenu polynomu iterativním způsobem. Počáteční odhad kořene je v každé iteraci zpřesňován a pokračuje se v iterování, dokud není dosaženo požadované přesnosti odhadu kořene. V našem případě potřebujeme spočítat kořen polynomu f (y) = y p −n. i) V každé iteraci algoritmus spočítá přesnější odhad kořene yi+1 = yi − ff′(y , (yi ) kde yi je hodnota vypočítaná v předchozí iteraci. y0 je počáteční odhad kořene polynomu f (y). Upravenou Newtonovou metodou budeme nazývat kombinaci binárního vyhledávání a Newtonovy metody. Prvních log p bitů hledané odmocniny x je spočítáno binárním vyhledáváním a výsledek je použit jako vstupní odhad pro Newtonovu metodu, pomocí které je číslo x dopočítáno na požadovanou přesnost.

2.2

Algoritmus B - sítové testy

Algoritmus B vznikne „vylepšením“ algoritmu A. Hlavní myšlenka úpravy algoritmu A vychází ze skutečnosti, že většina čísel nejsou perfektní mocniny a tedy algoritmus poběží rychleji, pokud se nám podaří vyloučit hned z počátku ta čísla, která perfektními mocninami nejsou. K tomu nám poslouží tzv. sítové testy (sieve tests). Vstupní číslo n projde pro daný exponent p testem, pokud je možné, že n je perfektní p-tá mocnina. Pokud číslo n testem neprojde, můžeme si být jisti, že n perfektní p-tá mocnina není. Taková čísla, která určitě nejsou perfektní p-té mocniny, budeme dále označovat jako „neperfektní“ p-té mocniny. Očividným důsledkem malé Fermatovy věty je následující lemma:

11

Lemma 2.2.1. Nechť p a q jsou prvočísla, pro která platí q ≡ 1 (mod p). Dále předpokládejme, že n je perfektní p-tá mocnina a gcd(n, q) = 1. Potom n(q−1)/p ≡ 1 (mod q). Výpočet n(q−1)/p budeme označovat jako sítový test pro číslo n a prvočíslo p. Číslo q nazveme modulem pro sítový test. Číslo n projde testem, pokud n(q−1)/p ≡ 0, 1 (mod q). V opačném případě n perfektní p-tá mocnina být nemůže a lze přejít na testování dalšího exponentu. Sítových  testů provedeme 2 log log n pro každý exponent p a vstupní číslo n nejvýše . Pokud n prošlo log p všemi sítovými testy pro exponent p, pak je možné, že je to perfektní p 2 log log n  tá mocnina. Číslo je zvoleno jako kompromis mezi tím, aby bylo log p vyloučeno co nejvíce neperfektních p-tých mocnin, ale aby sítových testů nebylo provedeno příliš mnoho. Jeho bližší zdůvodnění lze nalézt v [1]. AlgoritmusB(n) 1 for each p prvočíslo, p ≤log n do  2 log log n sítových testů na n; 2 proveď log p 3 if (n prošlo všemi testy) then 4 spočti x = ⌊n1/p ⌋; 5 if (n = xp ) then 6 n je perfektní p-tá mocnina; 7 stop; 8 n není perfektní mocnina; V algoritmu B se pro nalezení odmocniny používá upravená Newtonova metoda.

2.3

Algoritmus C - zkusmé dělení

Algoritmus B popsaný v předešlé části jsme dostali doplněním základního algoritmu A na hledání perfektních mocnin o sítové testy, které vyloučily na začátku některá čísla, která nejsou perfektní mocniny pro nějaký exponent p. Algoritmus C získáme zlepšením algoritmu B přidáním zkusmého dělení (trial division). Zkusmým dělením je myšlena procedura, kde je testováno, zda r | n, pro prvočísla r ≤ b. b je číslo zvoleno tak, aby bylo o dost menší než n, řádově log n. Číslo b je zvoleno jako nejmenší číslo takové, pro které platí b log2 b ≥ log n 12

. Po provedení zkusmého dělení nastane jedna ze dvou situací: (a) Najdeme prvočíslo r takové, že r | n. V tomto případě spočítáme exponent e, pro který platí re k n. Je-li vstupní číslo n perfektní p-tá mocnina, pak musí platit, že p | e. Protože e je relativně malé číslo, takových p, že p dělí e, je málo. Tedy se oproti algoritmu B značně zmenší počet takových prvočíselných exponentů p, které musí být dále prověřeny. (b) Není nalezen žádný prvočíselný dělitel čísla n menší než b. Potom víme, že pokud je n perfektní p-tá mocnina, pak p-tá odmocnina z n (označme ji x), musí být větší než b. A tedy bp ≤ xp = n ⇒ logb bp ≤ logb n ⇒ p ≤ n logb n ⇒ p ≤ log , čímž se též zmenší počet exponentů k prověřování. log b V algoritmu C je pro hledání odmocniny opět použita upravená Newtonova metoda. q v algoritmu značí prvočíselný modulus pro sítové testy.

13

AlgoritmusC(n) 1 Spočti nejmenší přirozené číslo b : b log2 b ≥ log n; n 2 S ← {p : p ≤ log }; log b 3 for each r prvočíslo, r ≤ b do 4 if (r | n) then 5 najdi e : re k n; 6 S ← {p : p | e}; 7 ukonči dělení; 8 while (S 6= ∅) do 9 p = nejmenší  2 log logprvek  S; n 10 proveď sítových testů na n; log p 11 if (q | n) then 12 najdi e : q e k n ; 13 S ← S ∩ {p : p | e}; 14 if ((n prošlo všemi testy) & ( p ∈ S )) then 15 x = ⌊n1/p ⌋; 16 if (n = xp ) then 17 n je perfektní p-tá mocnina; 18 stop; 19 S ← S − {p}; 20 číslo n není perfektní mocnina;

14

Kapitola 3 Algoritmus Daniela J. Bernsteina Dalším algoritmem pro testování perfektních mocnin je algoritmus Daniela J. Bernsteina popsaný v článku Detecting Perfect Powers in Essentially Linear Time [2]. Bernstein díky upravení některých základních funkcí docílil v podstatě lineární složitosti algoritmu (essencially linear time). Protože v algoritmu nejsou použita pouze přirozená čísla, ale také čísla s plovoucí řádovou čárkou (floating-point numbers), první část této kapitoly je zaměřena na tato čísla a na operace s nimi. Potom následuje popis jednotlivých funkcí Bernsteinova algoritmu. Při popisu algoritmů a funkcí v celé kapitole platí, není-li uvedeno jinak: • b, k, n, n′ , x jsou přirozená čísla • a, a′ jsou celá čísla • p je prvočíslo • y, r, r′ jsou čísla s plovoucí řádovou čárkou

3.1

Operace na číslech s plovoucí řádovou čárkou

Číslo s plovoucí řádovou čárkou r je reprezentováno dvojicí (a, n), a ∈ Z a n ∈ N, tak, že platí r = 2a n. Reprezentaci čísla s plovoucí řádovou čárkou dostaneme vydělením mocninou čísla 2.

3.1.1

Sčítání, odčítání, násobení

Jsou-li r a r′ dvě čísla s plovoucí řádovou čárkou reprezentovaná dvojicemi ′ (a, n) a (a′ , n′ ), čili r = 2a n a r′ = 2a n′ , pak aritmetické operace sčítání, 15

odčítání a násobení čísel r a r′ vypadají následovně. ′

r + r′ = 2f (2a−f n + 2a −f n′ ) ′ r − r′ = 2f (2a−f n − 2a −f n′ ) ′ r · r′ = 2a+a nn′ f = min(a, a′ ) a předpokládáme, že r ≥ r′ .

3.1.2

Dělení, zaokrouhlování, umocňování

Funkce pro dělení divb (r, k) vrací takové číslo s plovoucí řádovou čárkou, které aproximuje podíl kr tak, že platí k divrb (r,k) ∈ (1, 1 + 21−b ). Funkci definujme takto: a

a+f −⌈log k⌉−b

divb (2 n, k) = 2



n 2f −⌈log k⌉−b k



f je číslo takové, že pro vstupní číslo r = 2a n platí 2f −1 ≤ n < 2f . Funkci pro zaokrouhlování čísla na b bitů označíme truncb (r). Zaokrouhlení na b bitů znamená, že truncr b (r) ∈ (1, 1 + 21−b ). Funkci definujme: truncb (2 n) = divb (2 n, 1) = 2 a

a

a+f −b



n 2f −b



f je číslo takové, že pro vstupní číslo r = 2a n platí 2f −1 ≤ n < 2f . Funkci pro umocňování značíme powb (r, k). Jedná se o b-bitovou aproximaci čísla rk , čili powb (r, k) ≤ rk < powb (r, k)(1 + 21−b )2k−1 . Funkce je definována následujícím způsobem: powb (r, 1) = truncb r powb (r, 2k) = truncb (powb (r, k)2 ) powb (r, 2k + 1) = truncb (powb (r, 2k)powb (r, 1))

16

3.2

Bernsteinův algoritmus

Bernsteinův algoritmus vychází ze stejné myšlenky jako algoritmus A Bacha & Sorensona. Jedná se o základní algoritmus na počítání perfektních mocnin, který má na vstupu přirozené číslo n:

1 for each p prvočíslo, p < ⌊log 2n⌋ do 2 zkontroluj, zda n je p-tá mocnina;

Názvy funkcí Bernsteinova algoritmu jsou změněny oproti původním názvům v článku tak, aby se funkce nepletly s názvy algoritmů Bacha & Sorensona. Pro lehčí orientaci v původním článku však u každé funkce uvedeme i název použitý v článku D. Bernsteina. K pochopení celého algoritmu je zapotřebí podrobně popsat tři funkce. • Nrootb (y, k) spočítá aproximaci čísla y −1/k • Ppower(n, k, y) zjistí, zda n je k-tá mocnina; v y je předem vypočtená aproximace hodnoty n−1 . (Původní označení funkce je Algoritmus K.) • Compare(n, x, k) zjistí, zda n = xk ; funkce Compare je použita ve funkci Ppower. (Původní označení funkce je Algoritmus C.)

3.2.1

Funkce Nrootb

Vstupními argumrnty funkce Nrootb jsou číslo s plovoucí řádovou čárkou y, k ∈ N a b ∈ N. Funkce spočítá s přesností b aproximaci čísla y −1/k . Zadanou přesností se podobně jako u ostatních Bernsteinových funkcí myslí, že výstup funkce splňuje podmínku: (1 − 2−b )Nrootb (y, k) < y −1/k < (1 + 2−b )Nrootb (y, k) . Pro počítání odmocniny zvolil Bernstein, stejně jako Bach & Sorenson, kombinaci binárního vyhledávání a Newtonovy metody. Výsledná funkce je pojmenována Nrootb a počítá aproximaci čísla y −1/k . Binární vyhledávání označíme BinarySearch (v článku označeno jako Algoritmus B) a Newtonovu metodu NewtonsMethod (v článku označena jako Algoritmus N). 17

V každém kroku binárního vyhledávání je prohledávaný interval, ve kterém hledaný kořen leží, zmenšen na polovinu. Nejprve je aproximována hodnota výrazu z k y a vybrána levá nebo pravá polovina intervalu podle toho, zda je výraz větší či menší než 1. Bernstein ještě přidává podmínku, že pokud je hodnota výrazu z k y příliš blízko 1 pro zkoumané z, nevybere pravou nebo levou polovinu intervalu, ale použije „střední“ část intervalu. Čili část intervalu okolo středu, jejíž velikost je opět polovina původního intervalu. Hodnota výrazu z k y není vyčíslována přesně, ale pouze aproximována, čímž je dosaženo velkého zrychlení výpočtu. BinarySearch(b , k , y ) 1 najdi g : 2g−1 < n ≤ 2g ; 2 a = ⌊ −g ⌋; k 3 B = ⌈log(66(2k + 1))⌉; 4 z = 2a + 2a−1 ; 5 j = 1; 6 while (j < b) do 7 r = truncB (powB (z, k)truncB (y)); 993 8 if (r ≤ 1024 ) then 9 z = z − 2a−j−1 ; 10 else if (r > 1) then 11 z = z + 2a−j−1 ; 12 j = j + 1; 13 return z;

Bernstein používá ve svém algoritmu mnoho různých konstant. Ne vždy je přesně vysvětleno, proč je použita právě daná konstanta. Většina konstant byla zvolena kvůli odhadům v důkazech složitostí a správností jednotlivých 993 funkcí. Jednou z takových konstant je např. i 1024 použitá v binárním vyhle993 dávání. Bernstein její použití zdůvodňuje tím, že je číslo 1024 je nejjednodušší −1 32 číslo s plovoucí řádovou čárkou v rozmezí mezi 33 a e 33 , což pak používá v důkazech rychlosti a správnosti funkce BinarySearch, které lze nalézt v [2]. Cílem Newtonovy metody je najít kořen polynomu h(z) = z −k y −1 − 1, čili číslo y −1/k . Prvotní odhad kořene z je v každém cyklu metody nahrazován k+1 přesnější hodnotou znew = z − hh(z) ). Algoritmus je iterován ′ (z) = z + (z − yz tak dlouho, dokud z nemá požadovanou přesnost. Níže popsaná funkce NewtonsMethod je tedy ta část Newtonovy metody, která je stále iterována, až je z dostatečně blízko správnému řešení. Každá iterace Newtonovy metody 18

NewtonsMethod(b , k , y , z ) 1 r2 = mul(truncb (z), k + 1); 2 r3 = truncb (powb (z, k + 1)truncb (y)); 3 r4 = divb (r2 − r3 , k); 4 return r4 ;

přibližně zdvojnásobí počet správných číslic výsledku. Funkce mul(y, k) značí násobení čísla s plovoucí řádovou čárkou y přirozeným číslem k. Binární vyhledávání a Newtonova metoda je spojena do výsledné funkce Nrootb obdobně, jako tomu bylo u upravené Newtonovy metody. V Bernsteinově algoritmu je prvních 3 + ⌈log k⌉ bitů kořene z odhadnuto pomocí binárního vyhledávání a zbytek upřesněn Newtonovou metodou.

N rootb (k , y ) 1 if b < 4 + ⌈log k⌉ then 2 z = BinarySearchb (y, k); 3 else k⌉) ⌉; 4 b′ = ⌈log 2k⌉ + ⌈ (b−⌈log 2 5 B = 2b′ + 4 − ⌈log k⌉; 6 if (b′ < 4 + ⌈log k⌉) then 7 z = BinarySearchb′ (y, k); 8 else 9 z = Nrootb′ (y, k); 10 z = NewtonsMethodB (y, k, z); 11 return z;

3.2.2

Funkce Ppower

Funkce Ppower, v článku původně označovaná jako Algoritmus K, zjišťuje, zda n = xk . Vstupními parametry funkce jsou přirozená čísla n a k a číslo s plovoucí řádovou čárkou y. y je pomocná proměnná, ve které je uložena předpočítaná aproximace čísla n−1 , která se pak využije během výpočtu. Funkce Ppower nejprve spočítá aproximaci r čísla n1/k . Pokud existuje přirozené číslo

19

x takové, že |x − r| < 41 , zkusí, zda xk = n.

Ppower(n, k , y ) 1 f = ⌊log 2n⌋; 2 b = 3 + ⌈ fk ⌉; 3 r = Nrootb (y, k); 4 najdi x ∈ N : |x − r| ≤ 85 ; 5 if ((x = 0) or (|x − r| ≥ 14 )) then return 0; 6 else 7 s = Compare(n, k, x); 8 if (s = 0) then return x; 9 else return 0;

Funkce Compare, podrobně popsaná v další části, vrací znaménko výrazu n − xk . Tedy pokud s je rovno 0, znamená to, že n = xk a tedy, že n je k-tá mocnina. Funkce Ppower vrací číslo n1/k , pokud je n k-tá mocnina, v ostatních případech vrací 0. Kompletní důkaz správnosti funkce lze najít v [2].

3.2.3

Funkce Compare

Compare(n, k , x ) 1 f = ⌊log 2n⌋; 2 b = 1; 3 while (b < f ) do 4 r = powb+⌈log 8k⌉ (x, k); 5 if (n < r) then return -1; 6 if (r(1 + 2−b ) ≤ n) then return 1; 7 b = min{2b, f }; Vstupními argumenty funkce Compare, která je v článku původně označená jako Algoritmus C, jsou přirozená čísla n, x a k. Funkce testuje rovnost čísel n a xk . Návratová hodnota funkce je znaménko výrazu n − xk . Od klasických funkcí na porovnávání čísel se funkce Compare liší v tom, že 20

výraz xk nevyčísluje přesně, ale používá z něj pouze tolik prvních číslic, kolik potřebuje. Čím rozdílnější čísla jsou, tím rychleji je to rozpoznáno, a mocnina nemusí být dopočítána přesně, čímž se algoritmus urychlí. Důkaz správnosti funkce lze opět najít v [2].

3.2.4

Hlavní cyklus Bernsteinova algoritmu

Jelikož byly již všechny potřebné funkce popsány, můžeme nyní napsat celý Bernsteinův algoritmus. Hlavní cyklus Bernsteinova algoritmu je v článku označován jako Algoritmus X, my ho budeme nazývat BernsteinsAlgorithm. BernsteinsAlgorithm(n) 1 f = ⌊log 2n⌋; 2 y = nroot3+⌈f /2⌉ (n, 1); 3 for each p prvočíslo, p < f do 4 x = Ppower(n, p, y); 5 if (x > 0) then 6 n je perfektní p-tá mocnina (n = xp ); 7 n není perfektní mocnina;

21

Kapitola 4 Knihovna GMP Všechny algoritmy popsané v předchozích kapitolách jsme implementovali v jazyce C s pomocí funkcí knihovny GMP (GNU Multi-Precision Library) pro práci s velkými čísly. Cílem této kapitoly je seznámení s knihovnou GMP, popis jejích základních charakteristik a způsobu použití při implementaci algoritmů. V této kapitole je také popsána implementace některých aritmetických operací v knihovně GMP imlementovány, jak je to důležité pro naše testování a jak vypadalo nastavení knihovny v průběhu měření času běhu programů. Některé dodatečné informace, jako např. poznámky k instalaci knihovny, informace o nastavování jednotlivých parametrů knihovny apod. jsou podrobněji popsány v dodatku A.

4.1

O knihovně GMP

GNU Multi-Precision Library je knihovna určená k práci s libovolně velkými čísly nebo spíše čísly s libovolnou přesností (arbitrary precision arithmetic). Velikost čísel, se kterými jsou knihovní funkce schopny pracovat, je omezena typicky jen používaným hardwarem, především velikostí dostupné paměti. V knihovně GMP jsou implementovány funkce pro různé typy čísel. Nejpoužívanější jsou přirozená čísla, racionální čísla a čísla s plovoucí řádovou čárkou. Velkou předností knihovny je efektivita implementovaných operací a tedy i rychlost. Kvůli rychlosti operací a téměř neomezené velikosti čísel se knihovna GMP hodí především pro kryptografické aplikace, pro aplikace sloužící k zabezpečení komunikace přes Internet, pro implementaci algoritmů z výpočetní teorie čísel apod. Z projektů používajících knihovnu GMP uveďme např. knihovnu MPFR pro práci s čísly s plovoucí řádovou čárkou, knihovnu NTL pro teorii čísel nebo projekt Kaffe Java Virtual Machine. Seznam projektů používajících knihovnu GMP lze nalézt na oficiálních stránkách této 22

knihovny [6].

4.2

Aritmetické operace v knihovně GMP

Rychlost knihovny GMP je založena především na skutečnosti, že se použité algoritmy pro aritmetické operace mění v závislosti na velikosti operandů nebo na vlastnostech operandů. Protože naším cílem je zkoumat chování algoritmů, je třeba vysvětlit, jak jsou jednotlivé operace, jako násobení, dělení nebo umocňování v knihovně GMP implementovány a jaké algoritmy jsme během sledování chování algoritmů použili.

4.2.1

Násobení

Pro násobení jsou v knihovně GMP implementovány celkem 4 algoritmy: klasický „školský“ algoritmus, Karatsubův algoritmus, algoritmus Toom3 a FFT. Popis jednotlivých algoritmů lze nalézt v [5] nebo v [7]. Použitý algoritmus pro násobení dvou čísel je vybrán v závislosti na velikosti operandů. Hranice použití jednotlivých algoritmů jsou dané nastavením parametrů: MUL_KARATSUBA_THRESHOLD MUL_TOOM3_THRESHOLD MUL_FFT_THRESHOLD Velikost čísel v parametrech se udává v tzv. „limbech“. Limbem se v knihovně GMP rozumí taková část čísla, která se vejde do jednoho strojového slova, tedy obvykle 32 bitů. Velikost limbu se dá nastavit parametry BITS_PER_LIMB a BYTES_PER_LIMB. Knihovna má implicitně nastaveny hodnoty konstant pro různé třídy procesorů a jsou vybrány ty konstanty, které nejlépe vyhovují danému hardwaru. Konstanty je však možno změnit. Způsob nastavení konstant a případného zlepšení výkonu knihovny pomocí programu tuneup jsou popsány v dodatku A. Za účelem zkoumání chování algoritmů jsme nastavili konstanty tak, aby pro čísla ve zkoumaném rozmezí byl pro násobení použit pouze Karatsubův algoritmus. Tento algoritmus má asymptotickou složitost O(loglog2 3 n) = O(log1,585 n), kde n je vstupní číslo.

23

4.2.2

Dělení

Pro dělení čísel je také v knihovně GMP implementováno několik algoritmů, ale změna použitého algoritmu není dána vzrůstající velikostí operandů, jak tomu bylo u násobení. Různé algoritmy jsou většinou použity pro speciální případy, např. pokud je dělenec dvakrát větší než dělitel nebo pokud předem víme, že dělení vyjde beze zbytku. V našem případě nemůžeme zaručit, že při dělení nastane pouze nějaký speciální případ, a tedy počítáme s tím, že asymptotická složitost dělení je kvadratická, jako u klasické aritmetiky.

4.2.3

Umocňování

Použitý algoritmus násobení ovlivňuje složitost dalších operací, především umocňování. Knihovna GMP používá pro umocňování algoritmus „umocni a vynásob“ (square and multiply) variantu „zleva doprava“ (left to right). Vstupními argumenty jsou přirozená čísla k a e. Výstupem funkce je číslo n = k e . Binární zápis exponentu e označíme (et , et−1 , ..., e0 )2 , čili ei značí jednotlivé bity v binárním zápisu exponentu e. Square_and_multiply(k , e) 1 res = 1; 2 for (i = t downto 0) do 3 if (ei = 1) then 4 A = k · res · res; 5 else res = res · res; 6 return res;

V algoritmu je provedeno celkem t + 1 cyklů a v každém z nich se násobí čísla nejvýše rovná výsledku n. t + 1 = log2 e je počet bitů exponentu e. Umocňování „umocni a vynásob“ má tedy asymptotickou složitost O(log e logα n), kde logα n je složitost násobení. V našem případě počítáme α = 1, 585, jelikož počítáme s Karatsubovým algoritmem.

24

Kapitola 5 Implementace algoritmů V této kapitole je popsána vlastní implementace jednotlivých algoritmů řešení funkcí, použité struktury a datové typy, názvy funkcí v programech apod. Další podrobnosti týkající se implementace, jako přesné názvy funkcí, parametry funkcí, datové typy parametrů apod. jsou popsány v dodatku B. Datové typy a proměnné V implementaci algoritmů popsaných v článku Bacha & Sorensona byly použity klasické datové typy jazyka C: int unsigned long signed long a datový typ knihovny GMP: mpz_t Datový typ mpz_t slouží k ukládání velkých přirozených čísel. Tento typ byl použit především pro uložení vstupních čísel. Vstupní čísla jsou uložena v poli input_arr zmíněného typu, tedy mpz_t. V Bernsteinově algoritmu byly kromě výše zmíněných typů použity ještě dva další datové typy. mpf_t je datový typ knihovny GMP, který slouží pro ukládání velkých čísel s plovoucí řádovou čárkou. 25

Dále byla použita struktura struct { long a; mpz_t n; } fp; pro ukládání čísel y = 2a n. V úvodu bylo zmíněno, že popisované algoritmy jsou dekompoziční, tedy že algoritmus vrací přirozená čísla x a p v případě, že vstupní číslo n = xp . Pro toto slouží v programech pole base a exponent. Pokud je vstupní číslo input_arr[i] rozpoznáno jako perfektní mocnina, pak je do base[i] uloženo číslo x a do exponent[i] číslo p. Generování prvočísel Ve všech algoritmech potřebujeme mít k dispozici seznam prvočísel. Prvočísla potřebná v programech jsou generována pomocí funkce knihovny GMP nextprime. Všechna prvočísla vygenerovaná touto funkcí a používaná v programech jsou uložena v poli primes, kde primes[0] = 2, primes[1] = 3, primes[2] = 5 atd. V algoritmu A a v Bernsteinově algoritmu jsou zapotřebí pouze prvočísla menší než log n, kde n je vstupní číslo. U algoritmů B a C jsou potřeba ještě další prvočísla pro vytváření tabulky modulů pro sítové testy. Generování mocnin prvočísel Jedním ze sledovaných výstupů je chování algoritmů B a C v situaci, kdy vstupními čísly jsou pouze mocniny prvočísel. Pro generování mocnin prvočísel o určité velikosti jsme použili dvě konstanty: LOWER_NUMBER a UPPER_NUMBER, které vymezují požadovanou velikost generovaných mocnin prvočísel. Ze seznamu prvočísel je náhodně vybráno jedno, označme ho p, a je nalezen nejmenší exponent k takový, že 2U P P ER_N U M BER < pk . Pak je náhodně vybrán exponent velikosti nejvýše k, a pokud 2U P P ER_N U M BER < pk , pak mocnina pk spadá do požadovaného rozmezí. Vygenerovali jsme tedy další vstupní číslo, které uložíme do pole input_number. Generování vstupních čísel Parametry pro generování vstupních čísel (např. kolik čísel má být vygenerováno a kolik cifer mají vygenerovaná čísla mít) se v programu nastavují konstantami: 26

NUMBERS počet generovaných čísel BASE_NUM číselná soustava, ve které jsou čísla generována DIGITS počet cifer generovaných čísel MAX konstanta pro statickou alokaci polí FIRST_PRIME první prvočíslo ROUNDS pro testovací účely; nastavení, kolikrát se mají výpočty opakovat Konstanty byly při testování nastaveny následujícím způsobem: BASE_NUM 10 MAX 1 000 000 FIRST_PRIME 2 ROUNDS 1 Konstanta MAX byla v programech použita pro omezení velikosti pole vstupních čísel a pole prvočísel. Její velikost byla odvozena od počtu prvočísel odpovídajícímu maximálnímu modulu pro sítové testy a počtu generovaných vstupních čísel. Vstupní čísla pro programy jsou generována po jednotlivých cifrách v soustavě určené hodnotou BASE_NUM. Cifry vstupních čísel jsou generovány pomocí funkce knihovny GMP mpz_urandomm a vstupní čísla jsou ukládána do pole input_arr. Měření času Čas programů byl měřen pomocí knihovny time.h a funkce clock(). U každého implementovaného algoritmu jsme sledovali závislost doby běhu programu na délce vstupních čísel. Výsledky měření jsou uvedeny v následující kapitole.

5.1

Implementace algoritmů Bacha & Sorensona

Binární vyhledávání, upravená Newtonova metoda Ve všech algoritmech Bacha & Sorensona je použito binární vyhledávání a téměř ve všech upravená Newtonova metoda, proto tu popíšeme některé detaily týkající se těchto procedur.

27

Pro binární vyhledávání je potřeba nastavit na začátku spodní a horní mez intervalu, kde se hodnota má hledat. V našem případě je spodní mez stanovena na hodnotu 2, horní na 2⌊log n/p⌋+1 . V algoritmu, kde není použita Newtonova metoda, je binární vyhledávání iterováno, dokud se spodní a horní mez nerovnají, nebo dokud není pro daný prvočíselný exponent p nalezeno y, že y p = n, kde n je vstupní číslo a p je prověřovaný exponent. V dalších algoritmech s upravenou Newtonovou metodou je pro dané prvočíslo p počítáno pouze log p iterací binárního vyhledávání a zbytek je dopočítán upravenou Newtonovou metodou. Pomocí upravené Newtonovy metody počítáme kořen rovnice f (y) = y p − n, kde y je hledaná odmocnina, p je prověřovaný prvočíselný exponent a n je vstupní číslo. První odhad y získáme z binárního vyhledávání, následující zpřesňované odhady jsou počítány Newtonovou metodou. V každé iteraci Newtonovy metody je spočten nový přesnější odhad ynext = y −

f (y) (p − 1)xp − n = f ′ (y) pxp−1

Pokud je počítáno s racionálními čísly, Newtonova metoda bývá implementována tak, že je iterace prováděna tak dlouho, dokud hodnota f (y) = y p − n není dostatečně blízko nuly. Protože náš algoritmus počítá pouze s přirozenými čísly, bylo třeba koncové podmínky upravit. V každé iteraci Newtonovy metody je spočítán zpřesněný odhad kořene y, který je umocněn na právě prověřovaný prvočíselný exponent p. Výsledná mocnina je porovnána se vstupním číslem. Pokud se rovnají, tak je vstupní číslo n perfektní mocnina, neboť je rovno y p , a výpočet končí. Výpočet také končí, pokud je dvakrát za sebou spočten stejný odhad kořene y, protože se odhad nedá v přirozených číslech lépe zpřesnit, a nebo pokud je y p < n. Zdůvodnění je takové, že funkce f (y) = y p − n je napravo od nuly konkávní a jelikož původní odhad z binárního vyhledávání je hledán v intervalu h2, 2⌊log n/p⌋+1 i, pohybujeme se napravo od nuly. Ke kořeni se tedy „blížíme zprava“ v tom smyslu, že odhad y stále zmenšujeme. Pokud platí y p < n, jsme už nalevo od kořene, přesný kořen jsme tedy minuli, a tedy n není perfektní p-tá mocnina. Funkce vrací 1, pokud je číslo perfektní mocnina, 0 v ostatních případech. Sítové testy V algoritmu B jsou oproti algoritmu A přidány sítové testy. Je tedy třeba generovat moduly pro tyto sítové testy, čili taková prvočísla q, že q ≡ 1 (mod p), a také reprezentovat tabulku takovýchto modulů. K uložení hodnot pro každé prvočíslo p slouží struktura 28

struct sieve_moduli { unsigned long number_sm; unsigned long * sm; }   log n Pro každé prvočíslo p ≤ log n potřebujeme celkem 2 log modulů, pro log p každý sítový test jeden. Tabulku uložíme do pole struktur sieve_moduli. sieve_moduli[i] slouží pro uložení modulů pro prvočíslo, které je uloženo  log n a v sm je uložen v p = primes[i]. V number_sm je uložena hodnota 2 log log p jejich seznam. Pro vytvoření takovéto tabulky pro sítové testy slouží funkce fill_sieve_table a funkce sieve_tests slouží pro provádění vlastních sítových testů. Zkusmé dělení Funkce týkající se zkusmého dělení přidané do algoritmu C jsou funkce trial_division, která zkouší nalézt prvočíselného dělitele vstupního čísla, a dále funkce find_e pro nalezení exponentu e pro nějaké prvočíslo r tak, že re k n, kde n je vstupní číslo. Další přidané funkce jsou určeny pro práci s polem S, ve kterém jsou uložena taková prvočísla p, pro která je možné, že platí n = y p , pro nějaké přirozené číslo y.

5.2

Implementace Bernsteinova algoritmu

V Bernsteinově algoritmu byly pro nalezení p-té odmocniny z čísla n zkombinovány Newtonova metoda a binární vyhledávání do funkce Nrootb , která počítá aproximaci čísla n−1/p . Pro b ≤ 3+⌈log p⌉ použije binární vyhledávání, pro vyšší b pak Newtonovu metodu. Na začátku binárního vyhledávání je třeba stanovit interval, ve kterém se číslo bude hledat. V Bernsteinově algoritmu se na začátku binárního vyhledávání spočítá funkcí find_g exponent g takový, že 2g−1 < n ≤ 2g , a položí se a = ⌊− gp ⌋, takže platí, že 2a ≤ n−1/p < 2a+1 , a z tohoto intevalu se pak nadále v binárním vyhledávání vychází. V Newtonově metodě je v každé iteraci nahrazen odhad kořene z novým, k+1 ). Hodnota je počítána přesnějším odhadem znew = z − hh(z) ′ (z) = z + (z − yz pomocí funkcí trunc_b, pow_b a div_b, čili výrazy nejsou vyčíslovány přesně, ale vždy jen podle parametru b. Ten udává, na kolik bitů se má výsledek spočítat, což urychluje výpočty. Velká část implementace Bernsteinova algoritmu jsou funkce pro práci s vlastním definovaným typem fp, který slouží pro ukládání čísel y = 2a n. 29

Bylo potřeba implementovat aritmetické operace na číslech takto reprezentovaných, porovnávání čísel a také konverzi z datového typu mpz_t na tuto strukturu. Funkce pro aritmetické operace s tímto datovým typem jsou následující: fp_add pro sčítání, fp_sub pro odčítání, fp_mul pro násobení čísel a dále div_b pro dělení, trunc_b pro zaokrouhlování a pow_b pro umocňování. Všechny tyto funkce byly do detailů popsány v kapitole 3. Další potřebné funkce jsou funkce cmp_fp_fp a cmp_fp_int pro porovnávání čísel reprezentovaných strukturou fp jak mezi sebou, tak i s čísly typu mpz_t. První slouží pro porovnávání dvou čísel r1 a r2 reprezentovaných strukturou fp, druhá funkce slouží pro porovnávání čísla r1 ve struktuře fp s číslem r2 typu mpz_t. V obou případech je návratová hodnota kladné číslo, pokud je r1 větší než r2, nula, pokud jsou si čísla rovná, a záporné číslo, pokud r1 je menší než r2. Ve funkcích pro dělení a zaokrouhlování nebo v binárním vyhledávání jsou v Bernsteinově algoritmu zapotřebí takové exponenty f a g, že 2f ≤ n < 2f , resp. 2g−1 < y ≤ 2g . V obou případech je exponent počítán pro číslo y = 2a n. K nalezení takovýchto exponentů slouží funkce find_f a find_g. Tyto funkce mají k dispozici předpočítaná pole pow2 a pow2_zap, kde jsou uloženy mocniny čísla 2 a pomocí kterých se dané exponenty hledají. Funkce find_x nalezne ke vstupnímu číslu r s plovoucí řádovou čárkou takové přirozené číslo x, že |x − r| ≤ 5/8.

30

Kapitola 6 Složitost algoritmů a výsledky měření Cílem této kapitoly je prezentovat výsledky testování algoritmů a porovnat naměřené hodnoty s teoretickými výpočty. U každého algoritmu z článku Bacha & Sorensona jsou dokázány teoretické výpočty pro případ Karatsubova násobení. Teoretické výsledky pro případ „klasického“ násobení jsou pouze uvedeny, dokázány jsou v [1]. V důkazech používáme značení O(log α n) pro asymptotickou složitost násobení, případně O(log α n log p) pro asymptotickou složitost umocňování čísla np . V našem případě je α = 1, 585, protože počítáme s Karatsubovým algoritmem. Nejdůležitější z hlediska kryptografie jsou čísla do cca 2 000 decimálních cifer. Dnes je např. pro šifru RSA používán klíč velikosti 1 024 nebo 2 048 bitů, což jsou čísla o přibližně 300 a 600 decimálních cifrách, a předpokládá se, že v blízké budoucnosti budou stačit klíče o velikost 4 096 bitů, což je okolo 1 500 decimálních cifer. Tato čísla jsou však pro sledování asymptotického chování algoritmů přiliš malá, proto jsou výsledky každého algoritmu prezentovány na číslech od 1 000 do 10 000 cifer. U algoritmů B a C se během testování ukázala i velikost 10 000 decimálních cifer nedostatečná pro stanovení asymptotického chování algoritmu, proto jsou u těchto dvou algoritmů prezentovány výsledky až do 15 000 decimálních cifer. U každého programu je ještě pro zajímavost uvedeno, jak dlouho trvalo zpracování 500 ciferných čísel, která jsou pro dnešní kryptografii podstatná. Důvodem, proč nemohla být čísla do 1 000 cifer zahrnuta do zkoumání chování algoritmů, je, že nejnižší hodnota parametru MUL_KARATSUBA_THRESHOLD, jehož význam je popsán v kapitole 4, je 10 limbů, což odpovídá přibližně 100 ciferným číslům. Protože Karatsubův algoritmus při násobení dělí čísla na poloviny, při násobení „malých“ čísel se ještě znatelně projeví použití „škol31

ského“ násobení čísel pod 100 cifer. V každém grafu jsou vyneseny časy (v sekundách) v závislosti na velikosti čísel (počet decimálních cifer). Pro každý algoritmus byl měřen jiný počet vstupních čísel; konkrétní počet je uveden u jednotlivých algoritmů. Pro každý algoritmus byly naměřené hodnoty proloženy funkcí vystihující jejich průběh. K tomu byla použita regresní analýza (metoda nejmenších čtverců) implementovaná v programu MS Excel. Pro algoritmy Bacha & Sorensona byly použity mocninné funkce, pro Bernsteinův algoritmus funkce lineární. Programy byly spuštěny na procesoru AMD Athlon 64 X2 Dual Core Processor 3800+ s 32-bitovým jádrem (konstanta BITS_PER_LIMB byla nastavena na hodnotu 32). Verze použité knihovny GMP byla 4.2.1.

6.1 6.1.1

Algoritmus A Asymptotická složitost

Pro důkaz složitosti algoritmu A s binárním vyhledáváním využijeme rovnost: X1 p≤b

p

= log log b + O(1)

(6.1)

kde p jsou prvočísla (viz [4]).

K důkazu složitosti algoritmu A s upravenou Newtonovou metodou využijeme následující lemma. Lemma 6.1.1. Buď f (x) = xp − n, r > 0 splňující f (r) = 0 a x0 splňující 0 ≤ (x0 /r)−1 ≤ p1 . Potom pro spočítání odhadu kořene s absolutní chybou menší než 1 potřebuje Newtonova metoda O(log( logp n )) iterací, za předpokladu, že jako počáteční odhad kořene bylo použito x0 . Důkaz lemmatu je uveden v [1]. Věta 6.1.1. (Složitost algoritmu A s binárním vyhledáváním a klasickou aritmetikou) Buď n přirozené číslo. Pak algoritmus A, je-li použito pouze binární vyhledávání, rozhodne, zda n je nebo není perfektní mocnina v čase O(log3 n log log log n) je-li použita klasická aritmetika. 32

Věta 6.1.2. (Složitost algoritmu A s binárním vyhledáváním a Karatsubovým algoritmem) Buď n přirozené číslo. Pak algoritmus A, je-li použito pouze binární vyhledávání, rozhodne, zda n je nebo není perfektní mocnina v čase O(log2,585 n log log n log log log n) je-li použit pro násobení Karatsubův algoritmus. Důkaz. V algoritmu A je pro spočítání p-té odmocniny z n zapotřebí O( logp n ) iterací binárního vyhledávání. V každé iteraci je počítána aproximace p-té odmocniny, která je v zápětí umocněna na p-tou. Umocňování má složitost O(log p logα n). Pro pevný exponent p je tedy na spočítání ⌊n1/p ⌋ třeba α n log p čas O( log n log ). Součet přes všechny prvočíselné exponenty p ≤ log n p získáme pomocí rovnosti 6.1: X log p log p = log1+α n = p p p≤log n p≤log n ! X 1 = O log1+α n log log n = O(log1+α n log log n log log log n) p p≤log n X

log1+α n

V našem případě předpokládáme použití Karatsubova algoritmu, tedy výsledná složitost algoritmu A pouze s binárním vyhledáváním je O(log2,585 n log log n log log log n)

Věta 6.1.3. (Složitost algoritmu A s upravenou Newtonovou metodou a klasickou aritmetikou) Buď n přirozené číslo. Pak algoritmus A, je-li použita upravená Newtonova metoda, rozhodne, zda n je nebo není perfektní mocnina v čase O(log3 n) pokud je použita klasická aritmetika.

33

Věta 6.1.4. (Složitost algoritmu A s upravenou Newtonovou metodou a Karatsubovým algoritmem) Buď n přirozené číslo. Pak algoritmus A, je-li použita upravená Newtonova metoda, rozhodne, zda n je nebo není perfektní mocnina v čase O(log2,585 n log log n) je-li použit pro násobení Karatsubův algoritmus. Důkaz. V algoritmu A s upravenou Newtonovou metodou je provedeno pro každé prvočíslo p ≤ log n pouze log p iterací binárního vyhledávání a O(log( logp n )) = O(log log n − log p) iterací Newtonovy metody, což plyne z lemmatu 6.1.1. Celkem je tedy pro jedno prvočíslo p ≤ log n provedeno O(log log n) iterací. Protože na každou iteraci je třeba čas O(log p logα n), z rovnice 6.1 a z Čebyševovy věty plyne výsledná složitost: X X logα n log log n log p = logα n log log n log p = p≤log n

p≤log n

= O(logα n log log n log n) = O(log1+α n log log n) Dosadíme-li opět za α exponent Karatsubova násobení, vyjde nám: O(log2,585 n log log n)

6.1.2

Výsledky měření

Pro obě varianty algoritmu A byla náhodně generována čísla o velikosti 1 000, 2 000, ... , 10 000 decimálních cifer. Pro každou velikost bylo vygenerováno 10 čísel. Vidíme, že algoritmy nejsou příliš efektivní, protože zpracovat 10 čísel jim trvalo poměrně dlouho (u 10 000 ciferných čísel algoritmu A s upravenou Newtonovou metodou 125,1 s, u algoritmu A s binárním vyhledáváním dokonce 245,7 s). „Neefektivnost“ algoritmů A je ještě markantnější při porovnání s časy implementovaných algoritmů B a C. Pro informaci ještě uveďme, že zpracování 1 000 čísel o 500 cifrách trvalo oběma algoritmům A mezi 10 a 15 s. Algoritmus A s upravenou Newtonovou metodou byl o něco rychlejší (kolem 11 s) než algoritmus A s binárním vyhledáváním (kolem 13 s). Na grafech 6.1 a 6.2 tedy vidíme naměřené časy běhu programů, jimi proloženou funkci a rovnici této funkce. Pro oba algoritmy vidíme, že výsledný 34

Graf 6.1: algoritmus A s binárním vyhledáváním exponent odpovídá teoretickým výpočtům. Pro algoritmus A s binárním vyhledáváním i s upravenou Newtonovou metodou je podle teoretických výpočtů exponent hlavního členu v asymptotické složitosti roven 2,585. Z grafů vidíme, že exponent pro algoritmus A s binárním vyhledáváním vychází 2,6151 a pro algoritmus A s upravenou Newtonovou metodou 2,5299. U obou variant algoritmů A křivka velmi pěkně vystihuje chování naměřených dat. Při prokládání křivky je spočtena hodnota koeficientu determinace (R-squared value), která ukazuje, jak přesně vystihuje proložená křivka naměřená data. Čím blíže je hodnota R2 hodnotě 1, tím lépe jsou data proloženou křivkou aproximována. U algoritmu A s binárním vyhledáváním je R2 = 0, 9995 a u algoritmu A s upravenou Newtonovou metodou R2 = 0, 9998. Závěrem můžeme říci, že chování obou variant algoritmu A naměřené v praxi odpovídá teoretickým výsledkům.

6.2 6.2.1

Algoritmus B Čas běhu programu v průměrném případě

Před tím, než přistoupíme k důkazům vět o složitosti algoritmu B se ještě na chvíli zastavíme u sítových testů. Pro každé vstupní číslo n a prvočíslo 35

Graf 6.2: algoritmus A s upravenou Newtonovou metodou   log n sítových testů. Pro každý takový sítový p ≤ log n je provedeno 2 log log p test s prvočíslem p je třeba mít spočten modulus q, tedy prvočíslo takové, že platí: q ≡ 1 (mod p). V souvislosti s moduly pro sítové testy q se nabízejí dvě otázky: (a) Existuje dostatek takových prvočísel q a jak jsou tato prvočísla veliká? (b) Jaká je pravděpodobnost, že n projde všemi sítovými testy pro dané prvočíslo p? Existence dostatečného počtu prvočíselných modulů plyne z Dirichletovy věty, kterou ovšem nemůžeme použít k zodpovězení druhé části bodu (a). Odhad velikosti maximálního modulu pro sítové testy je založen na předpokladu platnosti rozšířené Riemannovy hypotézy(ERH)(viz [10]). Bach & Sorenson ve svém článku [1] dokázali, že za předpokladu platnosti ERH má největší modulus, potřebný pro sítové testy, velikost O(log2 n log4 log n). Následně ještě dodávají, že je pravděpodobně dokázaný odhad je nadsazen. Podle nich je přesnější odhad na velikost největšího modulu pro sítové testy log n ln2 (log n), což ale v článku už není dokázáno. Na lemma o velikosti největšího modulu se budeme dál v textu odkazovat, proto ho vyslovíme v celém znění: 36

Lemma 6.2.1. (o velikosti největšího modulu pro sítové testy) Pokud je na vstupu algoritmu přirozené číslo menší než n, pak má největší modulus potřebný pro sítové testy velikost nejvýše O(log2 n log4 log n). Tím byly zodpovězeny otázky z bodu (a). Ještě zbývá odpovědět na otázku týkající se pravděpodobnosti, že prvočíslo p projde všemi sítovými testy. Důkaz lze nalézt opět v [1], ve kterém je opět předpokládána platnost ERH k omezení velikosti největšího modulu potřebného pro sítové testy. Lemma 6.2.2. (ERH) Buď n přirozené číslo vybrané rovnoměrně z intervalu délky L a přepokládejme pro každé takové n, L ≥ (log n)3 log log log n . Potom log n pravděpodobnost, že n projde všemi 2 log sítovými testy pro pevně zvo log p  log n . lený exponent p, je shora omezena O log log2 n

Při výpočtu složitostí algoritmů B se předpokládá, že algoritmy mají k dispozici předpočítanou tabulku modulů pro sítové testy. Takovou  tabulkou 2 log log n je myšlena struktura, kde je pro každé prvočíslo p seznam prvočísellog p ných modulů q, q ≡ 1 (mod p). Z uvedených lemmat a předpokladů, především za předpokladu platnosti rozšířené Riemannovy hypotézy pro odhad velikosti modulů pro sítové testy, plyne věta o časové náročnosti algoritmu B: Věta 6.2.1. (Složitost algoritmu B v průměrném případě) (ERH) Buď n přirozené číslo a L jako v lemmatu 6.2.2 a přepokládejme, že tabulka modulů pro sítové testy je již vytvořena. Pak Algoritmus B rozhodne, zda n je nebo není perfektní mocnina a časová náročnost algoritmu je O(log2 n) počítáno jako průměr přes všechny vstupy z intervalu L. U algoritmu B se na výsledné časové náročnosti neprojeví použití asymptoticky lepšího Karatsubova násobení, proto jsme uvedli pouze jednu větu o časové náročnosti algoritmu B. Uvedeme zde důkaz pro Karatsubův algoritmus, původní důkaz lze opět nalézt v [1]. V důkazu časové náročnosti algoritmu B použijeme rovnost:   X 1 B (6.2) =O 2 log p log B p≤B kde p jsou prvočísla. 37

Důkaz. (důkaz rovnosti 6.2) X 1 X 1 X = + log p √ log p √ p≤B p≤ B

B≤p≤B

√ = O( B) + O



√ 1 = O( B) + O log p

B 1 log B log B



=O

Druhá rovnost plyne z věty o prvočíselnosti.





B 1 √ log B log B

B log2 B



=



Důkaz. Abychom dostali horní odhad, budeme předpokládat, že všechny možné sítové testy jsou provedeny. Z lemmatu 6.2.1 o velikosti největšího modulu pro sítové testy q plyne log q = O(log log n). Vzhledem k předpokladu, že tabulka je předpočítána, je čas k nalezení příslušného modulu pro sítový test O(1). Ke spočtení výrazu n(q−1)/p (mod q) je třeba jedno dělení a modulární umocňování v čase log n log q + log3 q = O(log n log log n + log3 log n) = log n O(log n log log n). Tedy pokud je celkem provedeno 2 log sítových testů, log p celkový čas pro každý prvočíselný exponent p je   log n log2 log n O log p

.

Při počítání asymptotické složitosti algoritmu A s upravenou Newtonovou metodou jsme ukázali, že čas potřebný k aproximaci p-té odmocniny z n a následné spočítání její p-té mocniny je O(logα n log log n log p). Nicméně z lemmatu 6.2.2 je známa pravděpodobnost, že toto počítání vůbec bude   log log n prováděno, která je O log2 n . Tedy průměrný čas běhu programu potřebný pro tuto část je:  2  log log n log p O log0,415 n

Spočteme-li složitosti přes všechny exponenty p s použitím 6.2 a Čebyševovy věty, vyjde nám v prvním případě O(log2 n) a v druhém případě ! X log2 log n log p = O(log2 log n log1−0,415 n) = O(log2 log n log0,585 n) O 0,415 log n p≤log n

Čili výsledný průměrný čas vychází O(log2 n).

38

6.2.2

Výsledky měření

Stejně jako pro algoritmy A, i pro algoritmus B jsme sestavili graf s výsledky měření a nechali jsme naměřenými daty proložit mocninnou funkci. Opět byla měřena čísla od 1 000 do 10 000 cifer, tentokrát však bylo měřeno pro každou velikost zpracování 5 000 vygenerovaných čísel, což je oproti předchozím algoritmům A znatelné zlepšení. Algoritmus B byl schopen zpracovat pětsetkrát více čísel než algoritmus A s upravenou Newtonovou metodou za téměř stejný čas. Pro informaci opět uvedeme, že zpracování 50 000 čísel s 500 ciframi algoritmem B trvalo kolem 14 s.

Graf 6.3: algoritmus B Z grafu 6.3 je vidět, že mocninná funkce nevystihuje průběh naměřených dat tak přesně, jako tomu bylo u algoritmů A. Hlavně u 10 000 ciferných čísel se proložená křivka a naměřená data dost rozcházejí a u více ciferných čísel by rozdíl byl ještě viditelnější. Možné vysvětlení je takové, že námi zkoumaná čísla jsou ještě příliš malá na správnou aproximaci, protože se v těchto „malých“ číslech projevuje vliv jiného, např. lineárního členu. Tuto hypotézu podporuje i graf 6.4, na kterém jsou výsledky zpracování 5 000 čísel od 7 000 do 15 000 cifer. Proložená křivka už daleko lépe odpovídá naměřeným datům a i navržený exponent daleko více teoretickým výpočtům.

39

Graf 6.4: algoritmus B pro větší čísla

6.3 6.3.1

Algoritmus C Čas běhu programu v průměrném případě

Stejně jako u Algoritmu B i zde budeme předpokládat následující: (a) tabulka modulů pro sítové testy je k dispozici (b) platí ERH (c) platí vlastnosti modulů pro sítové testy uvedené u algoritmu B V důkazu věty o časové náročnosti algoritmu C, použijeme následující skutečnosti, které jsou dokázány v [1]: n ), z čehož plyne, že log b je O(log log n). 1. b je O( loglog 2 log n 2. Pravděpodobnost, že v algoritmu není nalezen žádný prvočíselný dělitel menší než b, je O( log1 b ), kde b je číslo počítané na začátku algoritmu a používané jako horní mez pro zkusmé dělení. Věta 6.3.1. (Složitost Algoritmu C v průměrném případě) Buď n přirozené číslo vybrané rovnoměrně z intervalu popsaného v lemmatu 6.2.2. Potom algoritmus C rozhodne, zda n je nebo není perfektní mocnina,

40

a průměrný čas k tomu potřebný je   log2 n O log2 log n

Opět, jako u algoritmu B, se na výsledném času neprojeví použití Karatsubova násobení. Původní důkaz lze najít v [1], zde opět dokážeme pouze případ Karatsubova násobení. Důkaz. Algoritmus C lze rozdělit na 4 hlavní části: zkusmé dělení, hledání exponentu e pro různá prvočísla, počítání aproximací p-tých odmocnin a jejich umocňování na p-tou a sítové testy. Výpočet výsledného času algoritmu C rozdělíme právě do těchto 4 částí. P Čas potřebný pro zkusmé dělení O( p≤b log n log p) v kombinaci s odhadem čísla b a větou o prvočíslech dává celkový čas na zkusmé dělení.   X b O(log n log p) = O log n log b( ) = log b p≤b =O .



log n log log n log n log2 log n log log n



=O



log2 n log2 log n



Podle lemmatu 6.2.1 žádný modulus pro sítové testy nebo dělitel ve zkusmém dělení není větší než O(log2+ǫ n), tedy čas potřebný pro nalezení maximálního exponentu e pro taková čísla je O(e log n log log n). Podle lemmatu uvedeného v [1] je střední hodnota e O(1), tedy předpokládaná práce pro nalezení maximálního exponentu je zanedbatelná. K odhadu času potřebného pro třetí část, tedy aproximování p-tých odmocnin a počítání jejich p-tých mocnin, lze použít výsledky spočtené u algoritmu B. U algoritmu C nebude potřebný čas větší než u algoritmu B, tedy 2 n log p z výpočtu složitosti algoritmu B použijeme odhad O logloglog . Tedy 0,415 n celkový čas je O(log2 log n log0,585 ). Poslední zbývající částí jsou sítové testy. U této části záleží, zda je během zkusmého dělení nalezen dělitel, či není. Pokud je nalezeneno číslo r ≤ b, které dělí n, potom maximální exponent e je nejvýše log n. Tedy algoritmu zbývá ověřit nejvýše log log n exponentů, protože e má nejvýše log log n prlog2 log n vočíselných dělitelů. Maximální čas pro jeden exponent p je O( log n log ) p podle složitosti algoritmu B a tedy celkový průměrný čas pro sítové testy v případě, že byl nalezen dělitel, je   X log n log2 log n O = O(log log n log n log2 log n) = O(log n log3 log n) log p p≤log n 41

 2  n . V případě, že žádný dělitel nalezen nebyl, je potřebný čas roven O log log b Požadovaný výsledek dostaneme vynásobením tohoto času pravděpodobností, že nebude nalezen žádný dělitel.

6.3.2

Výsledky měření

Graf 6.5: algoritmus C Stejným způsobem jako pro předchozí algoritmy jsme i pro algoritmus C změřili dobu běhu programu pro různé velikosti čísel. Implementovaný algoritmus C je nejrychlejší ze všech algoritmů Bacha & Sorensona, časy v grafu 6.5 odpovídají zpracování 50 000 čísel. Zpracování 100 000 čísel s 500 ciframi algoritmem C trvalo kolem 3 s. Do měření času opět nebyly zahrnuty „přípravné“ výpočty, jako předpočítání tabulky modulů pro sítové testy a generování. Pro 10 000 ciferná čísla trvalo předpočítání modulů pro sítové testy kolem 8 s, generování prvočísel zabralo asi 5 s. Z grafu 6.5 vidíme, že nastala podobná situace, jako u algoritmu B. Mocninná funkce nevystihuje průběh naměřených dat úplně přesně. Stejně jako u algorimu B, i zde jsme sledovali chování programu pro čísla do 15 000 cifer. Výsledky jsou znázorněny na grafu 6.6 a vidíme, že křivka proložená daty vystihuje průběh daleko lépe, než v grafu 6.5. Vidíme také, že navržený exponent je menší, než vyšlo v teoretických důkazech. 42

Graf 6.6: algoritmus C pro větší čísla

6.3.3

Algoritmy B a C a mocniny prvočísla

Teoretické výpočty u algoritmů B a C se zabývají složitostí v průměrném případě. „Průměrným“ vstupem budeme nadále rozumět čísla generovaná náhodně bez jakékoliv specifické vlastnosti; vstupy, na jakých byly programy testovány v předchozích případech. Protože se testy prvočíselnosti používají především na zjištění, zda je vstupní číslo mocninou prvočísla či nikoliv, zkusili jsme, jak se budou algoritmy B a C chovat, pokud na vstupu budou pouze mocniny prvočísel. Z podstaty algoritmů B i C vyplývá, že budou mocniny prvočísel zpracovávat pomaleji, než „průměrný“ vstup, protože hlavní součást algoritmů, sítové testy, jsou zařazeny kvůli vyloučení neperfektních mocnin. Z grafu 6.7 vidíme, že doba zpracování vstupních čísel je asi desekrát delší, než u „průměrného“ vstupu. Zpracování 500 mocnin prvočísel trvalo kolem 130 s, u „průměrného“ vstupu se za přibližně 160 s zpracovalo 5 000 čísel. Zpracování mocnin prvočísla trvá déle, protože více dvojic vstupní číslo n a prvočíselný exponent p projde sítovými testy, a tedy je provedeno více binárních vyhledávání a upravených Newtonových metod. Pro srovnání jsme vyzkoušeli, kolik takových dvojic vstupní číslo n a prvočíselný exponent p projde sítovými testy pro „průměrný“ vstup. Vygenerovali jsme náhodně 10 000 čísel o 1 000, 5 000 a 10 000 decimálních cifrách a sledovali jsme, kolik dvojic 43

Graf 6.7: algoritmus B s mocninami prvočísel vstupní číslo n a prvočíselný exponent p projde sítovými testy, čili kolikrát musí být v programu provedena upravená Newtonova metoda a binární vyhledávání. Z deseti takových měření nám vyšel průměr 11,8 krát pro 1 000 ciferná čísla, 1,6 krát pro 5 000 ciferná čísla a 0,9 krát pro 10 000 ciferná čísla. Pokud je na vstupu 10 000 mocnin prvočísel, projde sítovými testy alespoň jednou každé vstupní číslo n s nějakým exponentem p. Čili vidíme, že pro mocniny prvočísel je provedeno daleko více upravených Newtonových metod a binárních vyhledávání než pro „průměrný“ vstup. I přes znatelné zpomalení je algoritmus B stále vhodnější pro zpracování mocnin prvočísel, než algoritmus A. Pro ty prvočíselné exponenty p, které jsou vyloučeny sítovými testy, nemusí být prováděna upravená Newtonova metoda a binární vyhledávání. Pro srovnání rychlosti uvedeme, že zpracování 50 mocnin prvočísel o 10 000 decimálních cifrách algoritmem A trvalo kolem 540 s, což je asi padesátkrát pomalejší, než zpracování stejného počtu mocnin prvočísel algoritmem B. V grafu 6.8 můžeme vidět výsledky měření pro algoritmus C. Porovnámeli exponenty u proložených funkcí, vychází algoritmus C o něco lépe, než algoritmus B. Do 5 000 cifer jsou časy zpracování pro algoritmus B a C téměř stejné, algoritmus B vychází o něco lépe. Podle naměřených časů se od 5 000 ciferných čísel vyplatí použít algoritmus C, máme-li na vstupu mocniny prvočísel. 44

Graf 6.8: algoritmus C s mocninami prvočísel

6.4 6.4.1

Bernsteinův algoritmus Složitost

Bernsteinův článek [2] je celý koncipován jako důkaz věty o tom, že existuje dekompoziční algoritmus pro testování perfektních mocnin, který běží nejvýše v čase (log n)1+o(1) , čili algoritmus běží v „v podstatě lineárním čase“ (in essentially linear time). Tento algoritmus popsaný v kapitole 3 jsme taktéž naimplementovali pomocí knihovny GMP. Výsledky měření jsou prezentovány v další části. Bernstein upravil základní funkce takovým způsobem, že výsledná časová náročnost algoritmu je v podstatě lineární. Cílem následující části je tedy podrobně popsat funkce, které jsou v Bernsteinově algoritmu použity. Bernstein ve svém článku uvedl, že předpokládá použití rychlého násobení. Neuvedl konkrétně algoritmus, který používá, ale s velkou pravděpodobností předpokládal použití Schönhage-Strassenova algoritmu založeného na rychlé Fourierově transformaci (více informací v [3] nebo v [5]), který má asymptotickou složitost O(k log k log k), kde k je velikost vstupního čísla. Kompletní rozsáhlý důkaz linearity Bernsteinova algoritmu lze najít v [2], my zde uvedeme základní myšlenku důkazu.

45

Bernstein nejprve zavede funkci F (n), kterou pak používá k celému důkazu. X F (n) = (log p)max{1, log n − d(n, (round n1/p )p )} 2≤p≤log n

kde p je prvočíslo a round t vyjadřuje takové přirozené číslo, které je ve vzdálenosti nejvýše 12 od t. Funkce F (n) má, dle věty o prvočíslech, zhruba log n členů, za každé prvočíslo v rozmezí je počítán jeden člen. V každém ln log n členu část max{1, log n − d(n, (round n1/p )p )}, označme ji M AX, vyjadřuje, kolik bitů z n je shodných s nějakou blízkou p-tou mocninou. Pokud je n velmi blízko k p-té mocnině, pak hodnota části M AX je blízko k log n. Faktor log p vyjadřuje úsilí vynaložené na počítání p-tých mocnin. Bernstein ve svém článku analyzuje chování funkce F(n) a dokázal, že je omezená funkcí (log n)1+ǫ(n) , pro jistou funkci ǫ ∈ o(1). Časovou náročnost hlavního cyklu Bernsteinova algoritmu lze omezit funkcí lineární v F (n) + log n. Skutečnost, že F (n) je v podstatě lineární v log n, což Bernstein ve svém článku dokázal, dává výslednou linearitu časové náročnosti Bernsteinova algoritmu.

6.4.2

Výsledky měření

Stejným způsobem, jako u předchozích algoritmů Bacha & Sorensona, i u naimplementovaného Bersteinova algoritmu bylo změřeno, jak dlouho program zpracovává různě velká čísla na vstupu. Výsledky testování můžeme vidět na grafu 6.9. Implementovaný Bernsteinův algoritmus jsme testovali na 100 vygenerovaných čísel velikosti od 1 000 do 10 000 cifer. Porovnáme-li rychlost s naměřenými hodnotami u algoritmů Bacha & Sorensona, vyjde nám, že je Bernsteinův algoritmus zhruba osmsetkrát pomalejší, než algoritmus C, dvacetkrát pomalejší, než algoritmus B a zhruba dvacetkrát rychlejší než algoritmus A s upravenou Newtonovou metodou, mluvíme-li o našich implementacích jednotlivých algoritmů. Zpracování 1 000 čísel o 500 decimálních cifrách trvalo kolem 35 s. Musíme si však uvědomit, že samotným algoritmům Bacha & Sorensona, jejichž čas byl měřen, předcházely výpočty modulů pro sítové testy, které nebyly do časů započítávany. Naměřenými daty byla proložena lineární funkce, která data dobře aproximovala, soudě z hodnoty R2 = 0, 9992. Zajímavé je, že závislost doby běhu programu na velikosti dat vychází lineární navzdory tomu, že pro násobení nebyl použit Schönhage-Strassenův algoritmus, ale Karatsubův algoritmus, jehož asymptotická složitost je horší. Možným vysvětlením je, že v implementaci algoritmu pro většinu násobení byla použita funkce mpf_mul_2exp resp. mpz_mul_2exp, která pro zadaná čísla a typu unsigned long a n typu 46

mpf_t resp. mpz_t spočítá 2a n. Je možné, že jsou tyto funkce realizovány jinak, než s pomocí zvoleného, v našem případě Karatsubova násobení (např. v případě přirozených čísel jako bitový posun).

Graf 6.9: Bernsteinův algoritmus

6.5

Funkce mpz_perfect_power_p

V knihovně GMP je také implementovaná funkce pro testování perfektních mocnin, nazvaná mpz_perfect_power_p(number), kde number je typu mpz_t. Funkce vrací nenulové číslo, pokud je number perfektní mocnina, tedy number = xk . Na rozdíl od námi implementovaných algoritmů funkce v knihovně GMP nevrací ani základ x ani exponent k. Zjišťovali jsme, jak se funkce chová v závislosti na velikosti vstupních čísel. Na grafu 6.10 můžeme vidět časy zpracování 1 000 čísel velikosti od 1 000 do 10 000 decimálních cifer. Naměřeným datům nejlépe odpovídala mocninná funkce (R2 = 0, 999). Exponent u proložené funkce by nejlépe odpovídal algoritmům A. Exponent u algoritmu A s Newtonovou metodou vycházel experimentálně 2,5299, exponent u funkce mpz_perfect_power_p vyšel 2,4268. Funkce v knihovně GMP je ovšem značn rychlejší, než implementované algoritmy A. I přesto, že do časů našich programů nebyly započítávány přípravné

47

Graf 6.10: funkce implementovaná v knihovně GMP práce, a nevíme, co přesně funkce knihovny GMP potřebuje, zda generuje moduly pro sítové testy či nikoliv, jak získává prvočísla apod. V manuálu ke knihovně [7] je k algoritmu pro testování perfektních mocnin pouze zmíněno, že se k nalezení odmocniny používá Newtonova metoda a že algoritmus používá redukci počtu možných prvočíselných exponentů. Čili je velmi pravděpodobné, že funkce je implementována jako algoritmus A s upravenou Newtonovou metodou s tím rozdílem, že před tím je provedeno zkusmé dělení k redukci takových prvočíselných exponenů p, které musí být v programu proveřeny (podobně jako u algoritmu C Bacha & Sorensona). Tím by bylo vysvětleno zrychlení oproti implementacím algoritmů A. Na grafu 6.11 jsou výsledky pro běh funkce mpz_perfect_power_p, pokud na vstupu měla pouze mocniny prvočísel. Zajímavé je, že exponent u proložené křivky vychází téměř stejný jako u „průměrného“ vstupu a také časy zpracování jsou skoro stejné. To by potvrzovalo domněnku, že ve funkci nejsou použity sítové testy, které nejvíce ovlivňovaly rozdíl mezi časy programu na „průměrném“ vstupu a pokud vstup byl tvořen pouze mocninami prvočísel. „Průměrným“ vstupem zde rozumíme čísla generovaná náhodně bez jakékoliv specifické vlastnosti.

48

Graf 6.11: funkce implementovaná v knihovně GMP s mocninami prvočísel

49

Kapitola 7 Závěr Cílem práce bylo porovnat dva typy algoritmů na testování perfektních mocnin. Oba přístupy, tedy algoritmy Bacha & Sorensona a Bernsteinův algoritmus, byly srovnány z hlediska teoretického přístupu, teoretické časové náročnosti a naměřené rychlosti. Do srovnání zahrnujeme i funkci mpz_perfect_power_p pro testování perfektních mocnin implementovanou v knihovně GMP. Teoretický přístup Oba přístupy, v práci Bacha & Sorensona i v Bernsteinově práci, vycházely ze základního jednoduchého algoritmu na testování perfektních mocnin. Bach & Sorenson použili sítové testy založené na malé Fermatově větě, aby co nejvíce neperfektních p-tých mocnin bylo vyloučeno už na začátku a složitější testy se prováděly pouze s „kandidáty“ na perfektní mocniny. Dále přidali zkusmé dělení, aby se co nejvíce zmenšil počet prvočísel p, které musí být prověřeny, aby se zjistilo, zda vstupní číslo je nebo není p-tá mocnina . Bernstein se naopak nesnažil vyloučit čísla podle specifických vlastností, ale upravil každou operaci, např. porovnávání nebo umocňování, tak, aby trvala co nejkratší dobu. Zaměřil se na to, že např. pro porovnání dvou čísel není třeba znát všechny bity obou čísel, ale jen pár prvních, podle kterých lze v mnoha případech poznat, že se čísla nerovnají. Takovýmto zjednodušováním operací a nevyčíslováním hodnot hned na všechny bity dokázal ušetřit mnoho času. Z hlediska teoretického přístupu nemůžeme o funkci mpz_perfect_power_p mnoho říci, protože nemáme dostatek informací. Z popisu algoritmu v manuálu ke knihovně GMP [7] a podle výsledků je pravděpodobné, že algoritmus používá upravenou Newtonovu metodu a zkusmé dělení, ale jak jsou jednotlivé funkce implementovány, nevíme. 50

Teoretická časová náročnost a naměřená rychlost Srovnáme-li algoritmy z hlediska teoretické časové náročnosti, vychází Bernsteinův algoritmus lépe, než všechny algoritmy Bacha & Sorensona. Porovnáme-li však naměřené časy, zjistíme, že pomalejší než Bernsteinův algoritmus jsou jen oba algoritmy A. Tato skutečnost může být zapříčiněna tím, že se v Bernsteinově algoritmu pracuje s dynamicky alokovanými strukturami a ukazateli, tedy režie s pamětí je mnohem větší, a také práce se strukturami fp, převod z přirozených čísel a také počítání hodnoty reprezentované proměnnými typu fp zpomalují běh programu. Na druhou stranu u algoritmů Bacha & Sorensona nebylo do měřeného času počítáno vytváření tabulky modulů pro sítové testy. Pro informaci pro vstupy o velikosti 10 000 decimálních cifer trval výpočet modulů okolo 8,5 s. Chování algoritmu C bylo počítáno na „průměrný“ vstup (čísla generovaná náhodně bez jakékoliv specifické vlastnosti). Čili pokud na vstupu programu byly jen samé perfektní mocniny, strategie, že na začátku algoritmus vyloučí ta čísla, která nejsou perfektní mocninou, nepřinášela žádné výhody, a celý program běžel pomaleji. Výsledky pro běh programu, pokud na vstupu byly pouze mocniny prvočísla, jsou prezentovány v kapitole 6. Z hlediska rychlosti srovnáme ještě se všemi programy funkci mpz_perfect_power_p implementovanou v knihovně GMP. Čas zpracování 1 000 čísel byl kolem 130 s, což je asi čtyřikrát pomalejší, než algoritmus B. Algoritmy A i Bernsteinův algoritmus vycházejí ze srovnání o dost pomaleji. Bernsteinův algoritmus je asi šestkrát pomalejší, algoritmus A s upravenou Newtonovou metodou je asi stokrát pomalejší. Opět platí, že srovnání může být zavádějící, jelikož nemáme o funkci mpz_perfect_power_p dostatek informací na to, abychom věděli, co vše funkce ke svým výpočtům používá, jaká část změřeného času tvoří „přípravné“ práce jako generování prvočísel nebo v případě algoritmů B a C příprava tabulky modulů pro sítové testy. Další možnosti zkoumání Závěrem se dá říci, že implementované algoritmy splňují teoretické chování, i když zvláště u algoritmů B a C bylo testování 10 000 ciferných čísel stále málo a asymptotické chování se projevilo lépe až na vyšších číslech. Dalo by se tedy pracovat především na zrychlení algoritmů. V případě Bensteinova algoritmu by bylo možno zvážit, zda některé parametry nepřidat jako součást struktury fp, případně neimplementovat speciální funkci nroot1, vycházející z funkce nroot_b ale upravenou pro počítání aproximací n1 , případně nepřidat do výpočtů ještě nějaké testy na vyloučení neperfektních mocnin. Protože přístupy k vylepšování algoritmů na testování perfektních moc51

nin byly velmi rozdílné, další možnost výzkumu se nabízí ve skloubení Bernsteinova přístupu a přístupu Bacha & Sorensona. Další výzkumy se naskýtají v oblasti mocnin prvočísel, které jsou také u testů prvočíselnosti a faktorizačních algoritmů nejvíce potřeba.

52

Literatura [1] Eric Bach, Jonathan Sorenson (1993): Sieve algorithms for perfect power testing, Algorithmica 9, 313-328 [2] Daniel J. Bernstein (1998): Detecting perfect powers in essetially linear time, Mathematics of Computation, Vol. 67, No. 223, 1253-1283 [3] Donald E. Knuth (1997): The Art of Computer Programming, Volume 2: Seminumerical Algorithms(3rd Edition), Addison-Wesley Professional [4] J. B. Rosser, L. Shoenfeld (1967): Approximate formulas for some functions of prime numbers, Illinois J. Math, 6, 64-94 [5] D. Stanovský: Počítačová algebra 2007/08 http://www.karlin.mff.cuni.cz/ stanovsk/vyuka/skripta_palg.pdf [6] GMP web pages http://gmplib.org/ [7] GNU MP http://gmplib.org/manual/ [8] GMP Install Instruction for Windows Platform http://cs.nyu.edu/exact/core/gmp/ [9] GNU Lesser General Public License http://www.gnu.org/copyleft/lesser.html [10] Extended Riemann Hypothesis http://mathworld.wolfram.com/ExtendedRiemannHypothesis.html [11] AKS Primality Test http://mathworld.wolfram.com/AKSPrimalityTest.html [12] Karatsuba Multiplication http://mathworld.wolfram.com/KaratsubaMultiplication.html 53

[13] Chebyshev Functions http://mathworld.wolfram.com/ChebyshevFunctions.html

54

Dodatek A Knihovna GMP Cílem tohoto dodatku je shrnout všechny poznatky o knihovně GMP nabyté během práce, jako instalace knihovny, práce s knihovnou, nastavení parametrů případně zlepšování výkonu knihovny. Jak bylo napsáno na začátku kapitoly 4, GNU Multi-Precision Library je knihovna určena k práci s čísly s libovolnou přesností (arbitrary precision arithmetic) resp. přesností omezenou používaným hardwarem. Knihovna GMP je součástí projektu GNU a je distribuována pod licencí LGPL (viz [9]). Knihovna je volně šiřitelná, spadá pod tzv. free software, což znamená, že jsou k dispozici i zdrojové kódy knihovny, které je povoleno libovolně měnit, upravovat a distribuovat. Zdrojové kódy knihovny lze stáhnout z oficiální stránky [6]. Po stáhnutí zdrojových kódů je třeba knihovnu zkompilovat a nainstalovat. V manuálu ke knihovně je popsána instalace pro systémy UN*X, ale na Internetu lze nalézt i postupy pro kompilaci pro operační systémy Windows, např. na [8]. Automaticky je knihovna nainstalována do adresáře /usr/local/, ale použitím přepínače –prefix u říkazu ./configure lze umístění knihovny změnit. Označujme dále adresář, do kterého máme knihovnu GMP nainstalovánu, jako $gmp_install. V kapitole 4 v části o aritmetických operacích byly zmíněny konstanty, jejichž nastavení mění použití algoritmů, např. pro násobení. Parametry pro algoritmy, jako např. MUL_KARATSUBA_THRESHOLD nebo MUL_TOOM3_THRESHOLD, se nastavují v souboru gmp-mparam.h. Podíváme-li se do adresářové struktury knihovny, nalezneme více souborů pojmenovaných gmp-mparam.h. Během instalace knihovny se podle používaného počítače vybere takový adresář, který odpovídá právě používanému hardware. Např. pro procesory AMD Athlon je použit soubor $gmp_install/mpn/x86/k7/gmp-mparam.h, pro procesory Intel P6 (PentiumPro, Pentium II, Pentium III) je použitý soubor s parametry v adresáři 55

$gmp_install/mpn/x86/p6 apod. V příslušném souboru gmp-mparam.h jsou nastaveny výchozí hodnoty, které se používají k výpočtům na daném stroji, ale chceme-li výkon knihovny ještě zlepšit, lze použít program tuneup, jehož zdrojové kódy se nacházejí v adresáři $gmp_install/tune. Program navrhne optimální hodnoty parametrů pro používaný procesor. Je třeba podotknout, že pokud změníme parametry v souboru gmp-mparam.h, je třeba celou knihovnu opět překompilovat a nainstalovat od začátku, aby se změny opravdu projevily. Uvedeme si příklad parametrů nastavených v knihovně a parametrů upravených programem tuneup. Pro procesory AMD Athlon jsou výchozí hodnoty parametrů: #define MUL_KARATSUBA_THRESHOLD #define MUL_TOOM3_THRESHOLD

28 98

#define SQR_BASECASE_THRESHOLD #define SQR_KARATSUBA_THRESHOLD #define SQR_TOOM3_THRESHOLD

0 50 113

#define MULLOW_BASECASE_THRESHOLD #define MULLOW_DC_THRESHOLD #define MULLOW_MUL_N_THRESHOLD

4 88 357

#define DIV_SB_PREINV_THRESHOLD #define DIV_DC_THRESHOLD #define POWM_THRESHOLD

0 97 128

#define GCD_ACCEL_THRESHOLD #define GCDEXT_THRESHOLD #define JACOBI_BASE_METHOD #define #define #define #define #define

/* always (native) */

/* always */

3 44 1

USE_PREINV_DIVREM_1 USE_PREINV_MOD_1 DIVREM_2_THRESHOLD DIVEXACT_1_THRESHOLD MODEXACT_1_ODD_THRESHOLD

#define GET_STR_DC_THRESHOLD

1 1 0 0 0 22

56

/* /* /* /* /*

native native always always always

*/ */ */ (native) */ (native) */

#define GET_STR_PRECOMPUTE_THRESHOLD #define SET_STR_THRESHOLD #define 393216, #define #define

35 3987

MUL_FFT_TABLE { 592, 1184, 2432, 5632, 14336, 73728, 163840, 1572864, 6291456, 0 } MUL_FFT_MODF_THRESHOLD 584 MUL_FFT_THRESHOLD 3840

#define SQR_FFT_TABLE { 656, 1312, 2688, 5632, 14336, 57344, 163840, 393216, 1572864, 6291456, 0 } #define SQR_FFT_MODF_THRESHOLD 672 #define SQR_FFT_THRESHOLD 3840 Hodnoty parametrů navrhnuté programem tuneup:

#define MUL_KARATSUBA_THRESHOLD #define MUL_TOOM3_THRESHOLD

24 89

#define SQR_BASECASE_THRESHOLD #define SQR_KARATSUBA_THRESHOLD #define SQR_TOOM3_THRESHOLD

0 46 137

#define MULLOW_BASECASE_THRESHOLD #define MULLOW_DC_THRESHOLD #define MULLOW_MUL_N_THRESHOLD

5 60 286

#define DIV_SB_PREINV_THRESHOLD #define DIV_DC_THRESHOLD #define POWM_THRESHOLD

0 64 104

#define GCD_ACCEL_THRESHOLD #define GCDEXT_THRESHOLD #define JACOBI_BASE_METHOD #define #define #define #define #define

/* always (native) */

/* always */

3 27 1

USE_PREINV_DIVREM_1 USE_PREINV_MOD_1 DIVREM_2_THRESHOLD DIVEXACT_1_THRESHOLD MODEXACT_1_ODD_THRESHOLD 57

1 1 0 0 0

/* /* /* /* /*

native native always always always

*/ */ */ (native) */ (native) */

#define GET_STR_DC_THRESHOLD #define GET_STR_PRECOMPUTE_THRESHOLD #define SET_STR_THRESHOLD

21 35 3987

#define MUL_FFT_TABLE { 624, 928, 1920, 5632, 14336, 57344, 0 } #define MUL_FFT_MODF_THRESHOLD 640 #define MUL_FFT_THRESHOLD 3840 #define SQR_FFT_TABLE { 688, 928, 1920, 5632, 14336, 57344, 0 } #define SQR_FFT_MODF_THRESHOLD 704 #define SQR_FFT_THRESHOLD 3840

58

Dodatek B Implementace algoritmů Tento dodatek slouží k doplnění informací z kapitoly 5, týkajících se především funkcí a jejich parametrů. Generování prvočísel Ve všech třech algoritmech potřebujeme mít k dispozici seznam prvočísel. V algoritmu A a v Bernsteinově algoritmu jsou zapotřebí pouze prvočísla p ≤ log n pro cyklus testování, zda vstupní číslo n není p-tou mocninou. U algoritmů B a C jsou potřeba ještě další prvočísla pro vytváření tabulky modulů pro sítové testy. Prvočísla potřebná v programech jsou generována pomocí funkce knihovny GMP nextprime(next_prime, prime) kde jsou oba parametry typu mpz_t. Funkce parametru next_prime uloží první prvočíslo větší než zadaná hodnota prime. Prvočísla jsou generována do pole primes_mpz a následně konvertována na typ unsigned long, se kterým se pracuje v celém výpočtu. Výhodou konverze prvočísel na typ unsigned long je především to, že je potom možné používat knihovní funkce na umocňování mpz_pow_ui a mpz_ui_pow_ui, kde je exponent typu unsigned long. Prvočísla používaná v algoritmech jsou konvertována z typu mpz_t na typ unsigned long. Zdůvodníme, že tuto konverzi lze provést. Víme, že pro hlavní cyklus budeme potřebovat prvočíslo nejvýše log n, kde n je vstupní číslo. Tedy máme-li vstupní číslo n o velikosti nejvýše 10 000 cifer v 59

desítkovém zápise, pak log n je přibližně 10 000, a tedy pro uložení prvočísel můžeme použít typ unsigned long. V algoritmech B a C jsou prvočísla pooužívána nejen jako možné exponenty, ale také pro vytvoření tabulky modulů pro sítové testy. Moduly jsou uloženy opět v poli hodnot typu unsigned long. V tomto případě je největší potřebné prvočíslo v programu rovno maximálnímu potřebnému modulu pro sítové testy a pro jeho odhad použijeme lemma 6.2.1 a na ně navazující poznámku autorů, že praxe ukázala, že odhad z lemmatu je možno ještě zlepšit. V praxi vychází největší modulus potřebný pro sítové testy pro číslo s 10 000 ciframi 2 839 301. Praxe ukázala, že pro nějaké algoritmy bylo 10 000 cifer stále ještě málo, a proto byla provedena měření pro čísla až do 15 000 cifer. Nejvyšší modulus je pro takto velká čísla 5 677 081, takže výše uvedená zdůvodnění lze použít i pro čísla s 15 000 decimálními ciframi.

B.1

Implementace algoritmů Bacha & Sorensona

Binární vyhledávání, upravená Newtonova metoda Funkce sloužící pro binární vyhledávání v algoritmech Bacha & Sorensona je BS(lower_limit, upper_limit, i, input_number) kde lower_limit a upper_limit jsou hodnoty typu mpz_t a ostatní dva parametry jsou typu mpz_t. lower_limit a upper_limit jsou počáteční horní a spodní meze pro binární vyhledávání, i je pořadové číslo právě prověřovaného exponentu, kterým je primes[i], ale v textu ho budeme většinou označovat pouze jako p. input_number je pořadové číslo právě zkoumaného vstupu, kterým je input_arr[input_number]. Na začátku binárního vyhledávání je spodní mez stanovena na hodnotu 2, horní mez na hodnotu 2⌊log n/p⌋+1 . Funkce pro Newtonovu metodu je: NM (i, input_number) 60

kde význam parametrů i a input_number je stejný jako u binárního vyhledávání. Princip upravené Newtonovy metody byl podrobněji popsán v kapitole 5. Sítové testy Jak již víme z kapitoly 5, k uložení hodnot modulů pro sítové testy pro každé p slouží struktura struct sieve_moduli { unsigned long number_sm; unsigned long * sm; } Funkce pracující s moduly pro sítové testy jsou: fill_sieve_table(i, log_log_n)  2 logVytváří  tabulku modulů pro sítové testy. Pro prvočíslo primes[i] najde log n potřebných modulů, které uloží do struktury sieve_moduli. log p sieve_test(i, input_number) Provede pro exponent primes[i] a vstupní číslo input_arr[input_number] patřičný počet sítových testů. Zkusmé dělení V algoritmu C jsou funkce pro binární vyhledávání, upravenou Newtonovu metodu i práci s moduly pro sítové testy implementovány stejně jako u algoritmů A a B. Funkce přidané oproti předešlým algoritmům se týkají především zkusmého dělení, nalezení maximálního exponentu e pro nějaké prvočíslo r, že re k n, kde n je vstupní číslo, a pak funkce pro práci s polem S. Pole S je pole těch prvočísel p, pro které je možné, že platí n = y p , pro nějaké přirozené číslo y.

61

find_e(r) Najde přirozené číslo e takové, že re k n. r je typu unsigned long. trial_division(void) Vyzkouší pro všechna prvočísla r ≤ b, jestli r nedělí vstupní číslo n. b je nejmenší přirozené číslo takové, že pro něj platí b log2 b ≥ log n. Pokud je nalezen r dělitel n, pak je pomocí funkce find_e nalezen exponent e, že re k n a příslušným způsobem aktualizováno pole S funkcí update_S2. init_S(void), update_S1(void), update_S2(void), update_S3(void) Funkce slouží k práci s polem S. Pole S je pole hodnot 0 a 1 takových, že S[i] = 1 když i ∈ S. Pokud S[i] = 1, znamená to, že pro prvočíslo primes[i] = p může platit, že vstupní číslo je p-tá mocnina, a je třeba ho prověřit. Inicializace funkcí init_S spočívá v nastavení všech polí na 0 a nastavení proměnné size_S na 0 (v proměnné size_S si pamatujeme aktuální počet prvků pole S). Zbylé tři funkce slouží k aktualizaci prvků patřících do pole S podle toho, jestli byl nebo nebyl nalezen nějaký dělitel vstupního čísla n. update_S1 přidá do pole S všechny prvočísla p ≤ b, kde b je mez popsaná u funkce trial_division. update_S2 dá do pole S všechna prvočísla p, která dělí e. e je přirozené číslo takové, že re k n, kde r je nalezený prvočíselný dělitel vstupního čísla n. update_S3 vyhodí z pole S ty prvky, které nedělí e, kde e je exponent popsaný u funkce update_S3.

B.2

Implementace Bernsteinova algoritmu

Jednou z hlavních částí implementace Bernsteinova algoritmu jsou funkce pro práci s nově definovaným typem fp, jako je inicializace, aritmetické operace na číslech reprezentovaných touto strukturou a porovnávání čísel. fp *init_struct(void)

62

Slouží k inicializaci struktury fp. Vytvoří se nová struktura a alokuje se pro ni prostor v paměti. Funkce vrací ukazatel na nově vytvořenou strukturu. free_struct(fp_p) „uklidí“ strukturu danou ukazatelem f p_p. Uvolní paměť, kterou struktura již nebude potřebovat. Základní funkce pro aritmetické operace jsou: fp_add(res_p, a_p, b_p) fp_sub(res_p, a_p, b_p) fp_mul(res_p, a_p, b_p) jsou funkce pro sčítání, odčítání a násobení dvou čísel a a b reprezentovaných strukturami fp. Výsledek je uložen do proměnné res která je také typu fp. Všechny parametry funkcí jsou ukazatele na dané struktury. div_b ( res_p, b, r_p, k) trunc_b ( res_p, b, r_p) pow_b (res_p, b, r_p, k ) jsou funkce pro další operace s čísly: dělení, zaokrouhlování a umocňování. r_p je ukazatel na strukturu typu fp, b a k jsou typu unsigned int a výsledek se ukládá do struktury typu fp, která je určena ukazatelem res_p. int cmp_fp_fp(r1_p, r2_p) int cmp_fp_int(r_p, n) První funkce slouží pro porovnávání dvou čísel r1 a r2 (které jsou uloženy ve strukturách typu fp určené ukazateli r1_p a r2_p). Druhá slouží pro porovnávání čísla r reprezentovaního strukturou fp a čísla n typu mpz_t. Při porovnávání je spočítána hodnota čísel ve struktuře fp, která je následně porovnávána s druhou hodnotou. Funkce vrací kladné číslo, pokud r1 > r2, resp. r > n, nulu, pokud r1 = r2, resp. r = n a záporné číslo, pokud r1 < r2, resp. r < n.

63

mpz_to_fp(n_fp_p, n) slouží pro převod čísla n uloženého v proměnné typu mpz_t na reprezentaci ve struktuře fp. Vstupní číslo n se dělí 2 tak dlouho, dokud je sudé, tím se získá exponent f p → a, a zbytek se uloží do f p → n, Pokud je číslo liché, do f p → a je uložena 0 a celé číslo je uloženo do f p → n. Mnoho funkcí potřebuje znát exponent f takový, že 2f ≤ n < 2f , pro číslo n typu mpz_t, resp. g, že 2g−1 < n ≤ 2g , pro číslo y = 2a n. K nalezení takových exponentů slouží funkce unsigned long find_f(n) signed long find_g(y_p) Rozdíl mezi těmito funkcemi spočívá v odlišných vstupních datových typech. Vstupní parametr funkce find_f je typu mpz_t a tedy víme, že hledaný exponent f je kladný, protože vstupní parametr n > 1. Vstupní parametr druhé funkce find_g je ukazatel na strukturu fp a hledané g může být i záporné. Funkce mají k dispozici předpočítaná pole mocnin dvojky pow2 a pow2_zap, která používají pro nalezení správného exponentu. find_x(r_p, input_number) Funkce nalezne takové přirozené číslo x, že |x − r| ≤ 5/8. binary_search(res_p, b, y_p, k) newtons_method(res_p, b, y_p, uk, est_p) nroot_b(z_p, b, y_p, k) Algoritmus_C(input_number, k) Algoritmus_K(input_number, i, y_p) Funkce jsou implementovány tak, jak byly popsány v kapitole 3, s tím rozdílem, že většina z funkcí má název shodný s názvem v původním článku, ne s názvy uvedenými v této práci. Jediné, co se občas odlišuje, jsou datové typy vstupních čísel. Samotné algoritmy byly popsány pro vstupní čísla určitých datových typů, ale v kontextu pak s ostatními funkcemi bylo třeba datové typy změnit. Jako příklad uvedeme funkci N rootb ). V popisu samotné 64

funkce se předpokládá, že vstupní číslo je typu floating-point, tedy číslo s plovoucí řádovou čárkou. Ovšem pokud je pak používána v hlavním programu BernsteinsAlgorithm, jako vstupní číslo je jí předáváno přirozené číslo n. Za zmínku k programu stojí ještě jedna věc týkající se knihovny GMP. U typu mpf_t se dá nastavit jeho přesnost(precision). K nastavení přesnosti jednotlivých proměnných typu mpf_t lze při inicializaci použít funkci mpf_init2(number, prec), kde number je právě inicializované číslo typu mpf_t, prec je typu unsigned long udávající přesnost čísla number. Přesnost typu mpf_t lze v programech nastavit pomocí konstanty PREC.

65