Recepty z programátorské kuchařky
Teorie čísel

Dnes si budeme povídat o různých užitečných vlastnostech celých čísel, především o dělitelnosti a kongruencích. Mohlo by se zdát, že to nemá s informatikou nic společného, ale překvapivě v informatice zakopáváme o teorii čísel takřka na každém kroku. Někdy se jedná o hledání velkých prvočísel, jindy o rychlé násobení čísel s miliony cifer nebo o všudypřítomnou asymetrickou šifru RSA.

Začneme vyjasněním základních pojmů, postupně se prokoušeme kongruencemi k hledání největšího společného dělitele a Bézoutových koeficientů, chvilku se zamyslíme nad prvočísly a nakonec si také ukážeme, jak to všechno souvisí s čínskou armádou.

Definice na úvod

Množinu celých čísel si označíme Z a každé její podmnožině {0, 1, … , n-1} budeme říkat Z_n.

Často nás bude zajímat dělitelnost: a\b (nebo a|b) budeme značit, že číslo a je dělitelem čísla b (nebude-li hrozit mýlka, čteme prostě „a dělí b“).

Pro největšího společného dělitele dvou čísel zavedeme symbol nsd(a,b). Pokud nsd(a,b) = 1, říkáme, že čísla a a b jsou nesoudělná, zkráceně a ⊥ b. Když budou naopak a a b soudělná, napíšeme a ∥ b. Podobně nejmenší společný násobek dvou čísel označíme nsn(a,b) a všimneme si, že je roven a · b / nsd(a,b).

Často nás může zajímat největší společný dělitel a nejmenší společný násobek více než dvou čísel. Není těžké si rozmyslet, že nsd(a, b, c) = nsd(a, nsd(b, c)), nsd(a, b, c, d) = nsd(a, nsd(b, nsd(c, d))) a obecně nsd(a₁,… ,a_k) = nsd(a₁,nsd(a₂,… ,a_k)). Obdobný vztah platí pro nejmenší společný násobek.

Není-li jedno číslo dělitelné druhým, znamená to, že při celočíselném dělení vznikne zbytek: například pokud vydělíme 23 / 8, dostaneme zbytek 7, protože 23 = 8  · 2 + 7. Obecně zbytkem po dělení a/b nazveme hodnotu z v rovnici a = b  · x + z, kde x je celé číslo a z je nezáporné celé číslo menší než b. Obvyklé programovací jazyky mívají takovouto operaci zabudovanou a říkají jí modulo. Programátoři počítají zbytky po dělení často nějakou konstrukcí podobnou z = a % b, zatímco matematici spíše píší z = a mod b.

Dodejme, že pro záporná čísla už není definice zbytku po dělení tak jednoznačná: v některých programovacích jazycích je (-7)%3 = -1, jiné se shodnou s naší definicí na tom, že vyjde 2. Přitom oba výsledky vycházejí z jedné rovnice a = b  · x + z. Aby měla jednoznačné řešení, požadovali jsme 0≤ z<b. Lze na to ovšem jít i jinak: řekneme, že x má být celočíselný podíl a/b. Ten jde ale definovat dvěma způsoby: buď se zaokrouhlením dolů (což se shodne s naší definicí), nebo se zaokrouhlením k nule (to dělá většina procesorů), což dá pro záporné a záporný zbytek. Zkuste zjistit (a vysvětlit), jak je to ve vašem oblíbeném jazyce a co se stane, když je záporné i číslo, kterým modulíme.

Kongruence

Když čísla p a q dávají stejný zbytek po dělení číslem m, píšeme

a čteme „p je kongruentní s q modulo m“. To platí právě tehdy, je-li rozdíl p-q dělitelný m.

Zápis kongruence tak trochu připomíná rovnici. To není náhoda – kongruence totiž můžeme upravovat podobně jako rovnice.

Součet dvou kongruencí: Pokud a ≡ A a b ≡ B, pak také platí a+b ≡ A+B (to vše modulo totéž m). Že je to pravda, nahlédneme snadno. Napišme si a jako n_a  · m + z_a a čísla A, b, B obdobně. Pak dostaneme:

Protože a ≡ A a b ≡ B, musí být z_a = z_A a z_b = z_B. Můžeme si tedy všimnout, že rozdíl mezi a + b a A + B je ((n_a + n_b) - (n_A + n_B))  · m. To je násobek m, takže a + b ≡ A + B.

Rozdíl dvou kongruencí: Nahlédneme obdobně.

Přičtení téhož čísla k oběma stranám: Pokud a ≡ A, pak platí a+k ≡ A+k pro libovolné k. Stačí totiž přičíst evidentně platnou kongruenci k ≡ k.

Přičtení násobku m k jedné straně: Z a ≡ A plyne a+k ≡ A pro libovolné k, které je násobkem modulu m. Přičítáme totiž kongruenci k ≡ 0.

Vynásobení dvou kongruencí: Z a ≡ A a b ≡ B plyne ab ≡ AB. Stejně jako u součtů a rozdílů, i zde stačí čísla rozepsat na součty násobků m a zbytků po dělení m:

Přitom opět víme, že z_a = z_A a z_b = z_B.

Vynásobení obou stran kongruence tímtéž číslem: Pokud a ≡ A, platí také ax ≡ Ax pro libovolné x. To plyne z násobení kongruencí x ≡ x.

Ekvivalentnost úprav: Běžné úpravy rovnic jsou takzvaně ekvivalentní – to znamená, že fungují oběma směry, takže řešení rovnic ani neubírají, ani nepřidávají. Jak je to s kongruencemi? Sčítání kongruencí ekvivalentní musí být, protože opačný směr odpovídá odečtení kongruencí, což víme, že je také korektní úprava.

U násobení kongruencí to už tak jasné není. Zkusme zjistit, jestli je pravda, že z kongruence ax ≡ Ax plyne a ≡ A. Pro x=0 to jistě neplatí, ale co když zvolíme jiné x?

Vyzkoušíme to třeba na následujícím příkladě:

Takováto kongruence má jednoduché řešení: a je každé celé číslo, které dostaneme sečtením 5 a nějakého násobku 14. To se dá zapsat třeba takhle:

tudíž a může být například 5, 19 nebo 33.

Vyzkoušíme nyní obě strany vynásobit… třeba trojkou:

Tato kongruence platí pro všechna a, která po vynásobení 3 dávají modulo 14 zbytek 1. Když máme nějaké a z první kongruence, je ve tvaru 5 + k  · 14. Vynásobíme-li ho 3, dostaneme 15 + 3  · k  · 14. To je určitě kongruentní s 1 modulo 14, takže žádné řešení jsme neztratili a po chvíli uvažování zjistíme, že jsme ani žádné nepřidali.

Další pokus: původní kongruenci vynásobíme místo trojky dvojkou. Dostaneme:

Řešení původní kongruence pořád sedí, ale nová kongruence platí například i pro a=12. Posuďte sami:

Ouha, najednou vynásobení obou stran konstantou není ekvivalentní úprava!

Proč nám násobení trojkou fungovalo, ale násobení dvojkou si vymýšlí kořeny navíc? Postupně se ukáže, že násobení k je ekvivalentní úprava právě tehdy, když k ⊥ 14 (či obecněji k⊥ m, počítáme-li modulo m).

Než k tomu dojdeme, nejdřív na chvíli odbočíme k největším společným dělitelům.

Euklidův algoritmus

Největšího společného dělitele dvou čísel můžeme vypočítat pomocí prvočíselného rozkladu, ale to je pro velká čísla velmi pomalé. Daleko lepší je použít prastarý Euklidův algoritmus. (Jmenuje se podle starověkého matematika Euklida, v jehož díle Základy se nachází první dochovaná verze. Jedná se zřejmě o nejstarší netriviální algoritmus, jaký se s drobnými úpravami používá dodnes. Dokonce je pravděpodobné, že Euklidés pouze sepsal dávno známý trik.)

Pojďme si odvodit, jak Euklidův algoritmus funguje. Nahlédneme, že pro libovolná čísla a,b (a>b) platí:

Proč je to pravda? Dokážeme, že dvojice (a,b) a (a-b,b) sdílejí dokonce všechny společné dělitele, takže i toho největšího:

• Nechť d je společným dělitelem a a b. Platí tedy a=a' · d, b=b' · d pro nějaká celá čísla a' a b'. Pak ovšem můžeme zapsat a-b jako (a'-b') · d, což je zase dělitelné číslem d.
• Nechť naopak d je společným dělitelem a-b a b. Opět zapíšeme a-b = c' · d, b = b' · d a získáme a = (a-b) + b = (c'+b') · d.

Euklidův algoritmus dostane na vstupu nějaká dvě čísla a a b a opakovaně odčítá menší z nich od většího. Jak už víme, tato operace zachovává největšího společného dělitele. Pokaždé se přitom součet a+b zmenší, takže po konečně mnoha krocích musíme jedno z čísel vynulovat. Pak víme, že největším společným dělitelem je druhé z nich (platí přeci nsd(0,x)=x).

Pojďme si takhle nějakého největšího společného dělitele spočítat. A abychom netroškařili, zkusme rovnou čísla 1518 a 945.

Zastavme se na chvíli. Teď bychom mohli pracně odečítat 30 od b, dokud bychom nenašli v b něco menšího než 30. Buďme trochu líní: když to budeme dělat dost dlouho, zbude nám v b zkrátka zbytek po dělení b číslem 30. Můžeme tedy místo odečítání modulit. Pokračujeme:

V a nám zbyla 0, takže nsd(1518, 945) = 3.

Pojďme si tento postup převést do návodu pro počítač. Při implementaci Euklidova algoritmu se hodí držet si v jedné proměnné pořád to větší z čísel a a b. Navíc můžeme využít toho, že každým krokem algoritmu se z většího čísla stane menší, takže je stačí prohodit a není potřeba znovu porovnávat. (Dokonce i porovnání před cyklem bychom si mohli ušetřit, kdyby nám nevadilo, že první průchod cyklem může projít „naprázdno“.)

def Euclid(a, b):
    # Prohodíme, je-li třeba
    if b > a:
        a, b = b, a

    # Zde je vždy a >= b
    while b > 0:
        # nsd(a % b, b) = nsd(a, b).
        a = a % b
        a, b = b, a

    # nsd(0, a) = a.
    return a

Jak rychle náš algoritmus běží? Podívejme se, co se stane, když pustíme dva kroky algoritmu na čísla a₁ a b₁ (a₁ ≥ b₁):

Dokážeme, že a₃ < a₁ / 2. Rozebereme přitom dva případy podle toho, jestli bylo b₁ menší, nebo větší než a₁ / 2:

• Pokud b₁ ≤ a₁ / 2, pak určitě platí a₃ < b₁, a tedy i a₃ < a₁ / 2. (Zde využíváme toho, že zbytek po dělení čehokoliv číslem b₁ musí být menší než b₁.)
• V opačném případě leží b₁ mezi a₁/2 a a₁, takže a₃ = a₁ mod b₁ = a₁ - b₁ < a₁ / 2. (Poslední rovnost platí, protože ⌊a₁ / b₁ ⌋= 1.)

Dokázali jsme tedy, že po dvou krocích algoritmu se větší z obou proměnných zmenší přinejmenším na polovinu a opět bude větší. Po O(log n) krocích tedy musí větší proměnná klesnout pod 1, čímž se algoritmus zastaví. Euklidův algoritmus proto provede O(log n) elementárních operací.

Jak dlouho ale trvá jedna elementární operace? Pokud počítáme s malými čísly, která se našemu počítači vejdou do celočíselné proměnné, zvládneme ji v konstantním čase. Jsou-li ovšem čísla větší, musíme ještě zohlednit složitost aritmetických operací: porovnání čísel a operace modulo. Když použijeme modulení pomocí školního dělení, které je kvadratické v počtu cifer, strávíme v každém z O(log n) kroků Euklidova algoritmu čas O(log ² n). Celková složitost algoritmu tedy vzroste na O(log ³ n).

Rozšířený Euklidův algoritmus

Právě jsme našli největšího společného dělitele d nějakých dvou obrovských čísel a a b. Jak ale přesvědčíme svého pochybovačného kolegu, že je náš výsledek správně? Snadno ověříme, že d dělí obě čísla. Ale jak ukážeme, že žádné větší číslo už a ani b nedělí? Překvapivě to jde jednoduše dosvědčit: pokud pro nějaká u a v platí

musí d být dělitelné každým společným dělitelem a a b, takže i číslem nsd(a,b). Nemůže tedy být menší než nsd(a,b).

Dobrá – kde taková u a v vzít? Kupodivu snadno: trochu upravíme Euklidův algoritmus. Nejprve ale prozradíme, že rovnici

se říká Bézoutova identita a číslům u a v Bézoutovy koeficienty.

Dokážeme, že spustíme-li Euklidův algoritmus na čísla a a b, v každém okamžiku se v proměnných a a b budou nacházet čísla tvaru α · a + β · b (kde α a β jsou nějaká celá čísla). Na začátku to triviálně platí, neboť a = a a b = b, a pokaždé, když se proměnné mění, buď se prohazují, nebo se jedna odčítá od druhé. Obě tyto operace z výrazů uvedeného tvaru dělají opět výrazy uvedeného tvaru. Takže i konečný výsledek algoritmu, tedy nsd(a,b), musí jít zapsat v takovém tvaru.

Algoritmus proto upravíme tak, aby si stále udržoval proměnné α_a, α_b, β_a a β_b a vždy platilo

Ke konci algoritmu je, jak víme, b = nsd(a,b), takže α_b a β_b jsou hledané Bézoutovy koeficienty. Opět si to vyzkoušejme na výpočtu nsd(1518,945):

Algoritmus tedy tvrdí, že hledaný nsd splňuje rovnost

Snadným výpočtem ověříme, že je to pravda. (Poznamenejme, že Bézoutova identita má nekonečně mnoho řešení. Jak byste našli ta další?)

Převeďme své myšlenky do zdrojového kódu. Do Aa, Ba, Ab, Bb budeme ukládat koeficienty α_a, β_a, α_b a β_b.

def ExtEuclid(a, b):
    Aa, Ba = 1, 0  # a = 1 * a + 0 * b
    Ab, Bb = 0, 1  # b = 0 * a + 1 * b
    # Prohodíme, je-li třeba
    if b > a:
        a, b = b, a
        Aa, Ab = Ab, Aa
        Ba, Bb = Bb, Ba
    # Zde je vždy a >= b
    while b > 0:
        # Odečteme od proměnné a proměnnou b
        # tolikrát, kolikrát se tam vejde.
        # ("//" značí celočíselné dělení)
        Aa = Aa - (a // b) * Ab
        Ba = Ba - (a // b) * Bb
        # nsd(a % b, b) = nsd(a, b).
        a = a % b
        # Prohodíme
        a, b = b, a
        Aa, Ab = Ab, Aa
        Ba, Bb = Bb, Ba
    # nsd(0, a) = a.
    # Vrátíme také Bézoutovy koeficienty.
    return (a, Aa, Ba)

Žádná operace, kterou děláme s koeficienty pro proměnné a a b, netrvá asymptoticky déle než operace modulo. Přidáním počítání Bézoutových koeficientů si tedy časovou složitost Euklidova algoritmu nezhoršíme.

Řešení lineárních kongruencí

Bézoutovy koeficienty jsou užitečné také k řešení kongruencí. Pojďme si to na jedné kongruenci vyzkoušet.

Máme nakoupena 4 vajíčka. V obchodě se vajíčka prodávají pouze v balíčcích po 6 kusech, zatímco my je skladujeme v platech po 20 kusech. Kolik si musíme koupit balíčků, abychom neměli v žádném platu volno?

Přepišme si tento příklad do formy kongruence:

čili

To je totéž, jako že pro x a nějaké další celé číslo y platí

Ejhle, to je rovnice podobná Bézoutově identitě. Kdyby na její pravé straně byl nsd(6,20)=2, byla by to přesně Bézoutova identita a rozšířený Euklidův algoritmus by nám prozradil, že platí

Tím bychom měli vyřešeno.

Jenže v našem případě je na pravé straně 8krát víc, než bychom potřebovali. Tak obě strany Bézoutovy identity vynásobíme 8:

Řešením naší rovnice tedy je x=-24, y=8.

Tak hurá do obchodu nakoupit -24 balíčků vajec. Cože? Že záporné nemají? Nevadí – stačí si vzpomenout, že jsme původně počítali modulo 20, takže k x můžeme přičíst libovolný násobek 20 a dostaneme další řešení. Můžeme tedy jít třeba pro 16 balíčků. (Mimochodem, není to nejmenší počet balíčků, který vyhovuje úloze: 6 balíčků by také fungovalo. Rozmyslete si, jak najít nejmenší řešení kongruence.)

Teď už můžeme zformulovat obecný návod na řešení kongruence

s neznámou x. Kongruenci přepíšeme do tvaru

a označíme d = nsd(a,n). Rozlišíme 3 případy:

• d=b … tehdy jsou hledaná x a y rovná Bézoutovým koeficientům a najdeme je rozšířeným Euklidovým algoritmem.
• d \b … pak najdeme řešení x' a y' rovnice s d na pravé straně a položíme x = x'  · b/d a y = y'  · b/d.
• b není násobkem d … v tomto případě kongruence nemůže mít žádné řešení, neboť levá strana rovnice je pro každé x a y dělitelná d, zatímco pravá strana dělitelná d nikdy není.

Inverzní prvky modulo m

Vraťme se teď zpátky z dlouhé odbočky a zkusme se znovu zamyslet nad tím, kdy je vynásobení obou stran kongruence ve tvaru

konstantou k ekvivalentní úprava. Tak říkáme úpravě, která neubírá ani nepřidává řešení. Už máme dokázáno, že když x ≡ z, tak pro každé k platí i k  · x ≡ k  · z. Takže zbývá zajistit, aby každé řešení kongruence k  · x ≡ k  · z bylo i řešením x ≡ z.

Nejprve ukážeme, že pokud k je soudělné s m, je naše snaha předem ztracená. Označme d = nsd(k,m) > 1. Vezměme libovolnou dvojici x a z splňující kongruenci x ≡ z, což je totéž jako x-z ≡ 0. Nyní vytvořme novou dvojici x'=x a z' = z + m/d. Pro tu dostaneme

Ovšem kongruenci vynásobenou k tato nová dvojice stále splňuje:

Dokázali jsme tedy, že pokud číslo k, kterým násobíme obě strany kongruence, je soudělné s modulem m, nejedná se o ekvivalentní úpravu. Teď naopak ukážeme, že jsou-li k a m nesoudělná, ekvivalentní to je.

Nahlédneme, že kdykoliv k⊥ m, existuje nějaké číslo k^-1 ∈Z_m takové, že k · k^-1 ≡ 1. Tomuto číslu se říká inverzní prvek ke k (nebo také multiplikativní invers čísla k), a pokud jím kongruenci kx ≡ kz vynásobíme, získáme

což je kýžená kongruence x ≡ z.

Kongruenci k · k^-1 ≡ 1 přitom už umíme vyřešit – předchozí kapitola nám říká, že takové k^-1 existuje právě tehdy, je-li k⊥ m, a že se dá najít Euklidovým algoritmem. Dodejme ještě, že prvkům, které mají multiplikativní invers, se říká invertibilní prvky modulo m.

Konečná tělesa

Když speciálně zvolíme za m nějaké prvočíslo, budou všechny prvky Z_m kromě nuly invertibilní. Tím pádem se Z_m bude chovat dost podobně racionálním nebo reálným číslům. Má s nimi například tyto společné vlastnosti (sčítáním a násobením v případě Z_m myslíme operace modulo m):

• Sčítání je asociativní a komutativní.
• Pro každé a platí a+0 = a.
• Pro každé a existuje (-a) takové, že a + (-a) = 0.
• Násobení je asociativní a komutativní.
• Pro každé a je a · 1 = a.
• Pro každé nenulové a existuje a^-1 takové, že a · a^-1 = 1.
• Násobení a sčítání jsou distributivní: a  · (b - c) = a  · b - a  · c.

Obecněji, máme-li libovolnou množinu, můžeme v ní označit „jedničku“ a „nulu“ a „přibalit“ operace sčítání, násobení, „dej mi (-a)“ a „dej mi a^-1“. Operací přitom myslíme libovolnou funkci, která prvkům množiny nebo jejich dvojicím přiřazuje prvky. Pokud navíc pro naši množinu s operacemi platí všechny vyjmenované vlastnosti, říká se jí komutativní těleso.

Racionální, reálná i komplexní čísla jsou příklady takových těles a my jsme k nim přidali konečná tělesa velikosti prvočísla. (Na okraj poznamenejme, že je známo, že všechna konečná tělesa mají velikost mocniny prvočísla, což ovšem neznamená, že se vždy chovají jako celá čísla modulo nějakým m.)

Malá Fermatova věta

S prvočísly úzce souvisí takzvaná Malá Fermatova věta. Říká, že pokud je p prvočíslo a a libovolné číslo od 1 do p-1, tak

Tato věta má mnoho různých použití (třeba ve známém šifrovacím algoritmu RSA nebo níže v algoritmu na testování prvočíselnosti), nám se bude především hodit jako další způsob invertování čísel modulo prvočíslo:

takže a^p-2 je inverzní prvek k a.

Pojďme nyní Malou Fermatovu větu dokázat. Indukcí podle a budeme dokazovat ekvivalentní tvrzení

(Jelikož a je určitě nesoudělné s modulem p, tak už víme, že násobení obou stran kongruence je ekvivalentní úprava.)

Pro a = 1 je snadné vidět, že věta platí:

Teď uděláme indukční krok. Řekněme, že máme dokázáno, že naše věta platí pro nějaké a, a chceme ji dokázat i pro a+1. K tomu se bude hodit známá Binomická věta, která říká, že pro každé reálné x a y a přirozené n platí:

mající v čitateli i jmenovali zlomku právě i členů.

Indukce po nás chce, abychom dokázali, že (a+1)^p ≡ a. Rozepíšeme tedy levou stranu kongruence pomocí Binomické věty:

Tím je indukce hotova.

Fermatův test prvočíselnosti

Jako malou odměnu za dlouhý důkaz předvedeme, jak Malou Fermatovu větu využívat ke zjištění, zda je nějaké obrovské číslo n prvočíslem. Jistě bychom mohli zkoušet všechny kandidáty na dělitele od 2 do n-1 (nebo chytřeji do ⌊√n⌋), ale to by trvalo příliš dlouho.

Raději zkusíme vybrat nějaké náhodné a∈{1,… ,n-1} a spočítat, kolik je a^n-1 mod n. Pro prvočíselné n musí vyjít jednička, takže pokud vyjde něco jiného, usvědčili jsme n z toho, že není prvočíslem (aniž jsme našli jediného dělitele – zvláštní, že?).

Pokud pro toto konkrétní a jednička vyjde, samozřejmě to neznamená, že n je určitě prvočíslo. Vyzkoušíme proto několik různých a, a pokud test pro žádné z nich neselže, drze prohlásíme, že n je pravděpodobně prvočíslo.

Jak moc velká drzost to je? Překvapivě ne moc velká. Pro skoro každé složené číslo n platí, že alespoň polovina a-ček dosvědčí, že se nejedná o prvočíslo. Takže jeden pokus selže s pravděpodobností nejvýše 1/2 a pokud uděláme t pokusů, pravděpodobnost chybného výsledku je nanejvýš 1/2^t.

Jedinou výjimku z našeho pravidla tvoří takzvaná Carmichaelova čísla (nejmenší z nich je číslo 561). To jsou čísla, jejichž složenost prokážeme jen tehdy, když se strefíme do a soudělného s n, a takových a je velmi málo. Naštěstí není Carmichaelových čísel moc (relativně k prvočíslům), takže Fermatův test funguje docela spolehlivě.

Existují i důmyslnější testy, které se Carmichaelovými čísly obalamutit nenechají. Jejich popis, jakož i důkaz našeho tvrzení o spolehlivosti Fermatova testu, najdete v literatuře zmíněné na konci kuchařky.

Rychlé mocnění

Ve Fermatově testu nebo při počítání inverzí pomocí Malé Fermatovy věty potřebujeme spočítat a^k mod m pro velké k. Pokud budeme mocninu a^k počítat přímo podle definice, tedy jako a  · a  · …  · a, budeme potřebovat O(k) násobení, což je příliš.

Jednoduchou fintou lze počet operací snížit na O(log k). Například a¹⁶ můžeme spočítat jako:

Pro obecný exponent bude rychlejší umocňování vypadat takto:

def FastExp(a, k):
    # Nejdříve ošetříme triviální případy.
    if k == 0: return 1
    if k == 1: return a

    # Když je x sudé, vrátíme a^(k/2) * a^(k/2).
    # Když je x liché, vrátíme a * a^(k - 1).

    if k % 2 == 0:
        i = FastExp(a, k / 2)
        return i * i
    else:
        return a * FastExp(a, k - 1)

Každé volání FastExp pro sudé k jednou zavolá FastExp s polovičním k a jednou vynásobí dvě čísla. Když je k liché, převede se na sudé a provede se jedno vynásobení. FastExp tedy provede O(log k) násobení.

Je ale důležité uvědomit si, že kdybychom si neuložili výsledek a^k/2 do pomocné proměnné i, ale rovnou vraceli FastExp(a, k / 2) * FastExp(a, k / 2), byl by náš kód stejně pomalý, jako kdybychom počítali mocninu podle definice!

Pro použití ve Fermatově testu (nebo obecně na spočítání mocniny a^k mod m) stačí po každém násobení výsledek vymodulit m.

Síto na prvočísla

Už jsme zjistili, že se nám hodí umět najít prvočísla. Kde je ale vezmeme? Můžeme určitě zkoušet jedno číslo po druhém a pokaždé otestovat, jestli držíme prvočíslo (třeba Fermatovým testem nebo zkoušením všech dělitelů). Už staří Řekové ale znali algoritmus, který najde všechna prvočísla menší než n efektivněji. Říká se mu Eratosthenovo síto.

Síto funguje na docela jednoduchém principu. Budeme si uchovávat v paměti pro každé číslo od 2 do n příznak, jestli je prvočíslo, nebo složené. Začneme u dvojky a označíme všechny násobky 2 ležící mezi 4 a n jako složená čísla. Další prvočíslo je 3. Označíme všechny násobky 3 ležící od 6 do n jako složená čísla. Další číslo na řadě je 4. Když jsme ale vyškrtávali násobky 2, vyškrtli jsme i 4. Nebudeme tedy provádět nic a rovnou přejdeme na 5. Takto najdeme všechna prvočísla od 2 do n a stihneme to rychle.

Ukažme si ještě zdrojový kód v Pythonu:

def Eratosthenes(n):
    prvocislo = [ True ] * n

    for i in range(2, n):
        if prvocislo[i]:
            print("%d je prvocislo." % i)

            # Násobky prvočísla jsou složené
            j = i * 2
            while j < n:
                prvocislo[j] = False
            j += i

Jak dlouho síto poběží? Dá se dokázat, že jeho časová složitost činí O(n log log n), ale není to snadné. My si zde předvedeme jenom slabší odhad O(n log n). Zájemce o těžší důkaz menší složitosti odkazujeme na vzorové řešení úlohy 24-3-5.

Síto tráví čas O(n) hledáním prvočísel a mezitím škrtá jejich násobky. Když škrtáme násobky dvojky, vyškrtneme nejvýše n/2 čísel, když škrtáme násobky trojky, vyškrtneme jich nejvýše n/3, atd. Složitost Eratosthenova síta tedy bude shora omezena součtem

Sumě

se říká n-té harmonické číslo a dokážeme o něm, že leží v O(log n).

Uvažujme, o co se zvětší H_2n oproti H_n:

To je součet n členů, z nichž každý je menší než 1/n. Celý součet je tedy menší než 1.

Zjistili jsme tedy, že H_2n < H_n + 1. Funkce, které rostou takhle pomalu, jdou shora omezit nějakým logaritmem n, takže H_n = O(log n).

Proto složitost celého Eratosthenova síta činí O(n log n).

Čínská zbytková věta

Následující věta dostala své jméno po staročínském způsobu počítání vojáků. Čínská armáda je velká, a kdybychom chtěli počítat vojáky jednoho po druhém, trvalo by to dlouho. Pomáhalo prý armádu rozřadit do řad o velikostech m₁, m₂, …, m_n (součin všech m_i si označíme jako M bez indexu). Někdy zbyli nezařazení dva, někdy třicet, někdy se seřadili všichni. Tyto zbytky si označíme z₁, z₂, …, z_n.

A co Čínská zbytková věta říká? Tvrdí, že když jsou všechna m_i navzájem nesoudělná a počet vojáků je menší než M, lze ho ze zbytků z_i jednoznačně určit. Když například rozdělujeme vojáky do řad velikostí 2, 3, 5, 7, 11 a 13, můžeme zbytky z řad jednoznačně vyjádřit každý počet vojáků menší než 2  · 3  · 5  · 7  · 11  · 13 = 30 030.

Formálněji řečeno: Jsou-li dána navzájem nesoudělná přirozená čísla m₁, …, m_n (jejichž součin označíme M) a zbytky z₁, …, z_n, pak existuje právě jedno číslo x∈Z_M takové, že pro všechna i je

A jak dokážeme, že něco takového platí? Mějme 2 čísla a, b ∈Z_M taková, že mají stejné zbytky po dělení všemi m_i. Ukážeme, že musí nutně být stejná.

Víme, že pro všechna i platí

To podle definice kongruence znamená, že rozdíl a-b je dělitelný všemi m_i. Proto je dělitelný i nejmenším společným násobkem všech m_i, což ovšem díky nesoudělnosti musí být jejich součin M.

Máme tedy dvě čísla ze Z_M, jejichž rozdíl je dělitelný M. To nutně znamená, že jsou stejná.

Dokázali jsme tedy, že jedna sada zbytků z₁, …, z_n odpovídá jednoznačně určenému číslu x ∈Z_M, ale ještě nevíme, jak bez zkoušení všech možností toto x najít. Půjdeme na to od lesa.

Nejprve se hodí všimnout si toho, že když sečteme dvě čísla, sečtou se i jejich zbytky modulo všemi m_i.

Co kdybychom nyní dokázali sehnat čísla Q₁, …, Q_n taková, že Q_j je dělitelné všemi m_i kromě m_j a že Q_j ≡ 1 (mod m_j) ?

To by potom stačilo položit

Vskutku: počítáme-li x mod m_i, všechny členy z_j  · Q_j pro j≠ i vyjdou nulové a člen z_i  · Q_i bude roven z_i. To, že celý výsledek nakonec vymodulíme M, na věci nic nemění, protože přičtení či odečtení libovolného násobku M zbytek po dělení žádným m_i neovlivní.

Jak se ale k číslům Q_i dostaneme? Číslo Q_i má být dělitelné všemi m_j kromě m_i. Uvažujme tedy součin

Ten modulo každé m_j (j≠ i) dá nulu, zatímco modulo m_i nějaké číslo r_i nesoudělné s m_i (nesoudělné musí být, protože jinak by m_i bylo soudělné s některým m_j). Speciálně to znamená, že r_i není 0.

Toto Q_i už má požadované vlastnosti: Q_i mod m_j pro j≠ i vyjde nulové, protože Q_i je násobkem S_i, které bylo dělitelné m_j. A modulo m_i získáme

„Kouzelná“ čísla Q_i tedy dokážeme sestrojit a jejich zkombinováním i hledané x.

Pojďme si to teď zkusit v praxi. Chceme najít nejmenší x takové, že platí následující kongruence:

Spočítáme si nejdříve M a všechna S_i:

Teď zjistíme, kolik vychází každé S_i modulo m_i a určíme příslušné multiplikativní inverze (například pomocí rozšířeného Euklidova algoritmu):

Nakonec sečteme příslušné násobky Q_i a zjistíme x:

Výsledek opravdu vypadá správně:

Pár slov na závěr

Doufáme, že se vám naše povídání o teorii čísel líbilo a že jste poznali, že i tak základní objekty, jako jsou celá čísla, mají spousty zajímavých vlastností.

Přejete-li si dozvědět se více o prvočíselných testech nebo o RSA, můžeme navrhnout ke studiu textík Algoritmy okolo teorie čísel od jednoho z autorů kuchařky. Důkladný rozbor Eratosthenova síta a jiné zajímavosti o prvočíslech najdete v článku Tři věty o prvočíslech od téhož autora.

S teorií čísel také souvisí algebra, která zobecňuje různé poznatky na libovolné množiny opatřené nějakými operacemi (například tělesa). Máte-li o ni zájem, mohla by vám pomoci například skripta Základy algebry od Davida Stanovského.

Michal Pokorný a Martin Mareš