Pchełki Python: Test Millera-Rabina

https://xpil.eu/kmO0x 📋

Dziś pokażę, w jaki sposób sprawdzić pierwszość naprawdę wielkich liczb.

Niedawno, w ramach odrdzewiania swojej (i tak raczej mizernej) wiedzy o Pythonie postanowiłem policzyć sobie wszystkie liczby pierwsze od dwójki do miliarda.

Uczyniłem swój kod odpornym na pady zasilania - co jakiś czas zapisywał on bieżący stan (czyli listę aktualnie wyliczonych liczb pierwszych) do pliku, który przy starcie sobie odczytywał. Zrobiłem mu nawet w miarę sympatyczny interfejs informujący użytkownika, ile liczb pierwszych już znaleziono, jaka jest ostatnia znaleziona i tak dalej.

Tak się zabawiając dotarłem niepostrzeżenie do sześćdziesięciu milionów liczb pierwszych. Zapisanie takiego zbioru danych na dysk zajmowało już około 4 sekund, podobnie odczyt. Stwierdziłem, że dalej to bez sensu, bo i tak mogę sobie teraz sprawdzić (naiwnym algorytmem) pierwszość dowolnej liczby od jeden aż do iluś-tam biliardów. Z tym, że te większe to jednak sprawdzały się naprawdę wolno. Nic dziwnego, wykonanie miliona dzieleń może trochę zająć.

Zainteresowałem się więc wreszcie "na poważnie" algorytmami bardziej zaawansowanymi. Nie, żebym spróbował zrozumieć stojącą za nimi teorię - o, co to, to nie, za głupi jestem. Musiałbym zrozumieć uogólnioną hipotezę Riemanna oraz próby jej udowodnienia, musiałbym biegle operować terminologią krzywych eliptycznych a także liczbami zespolonymi (te ostatnie trochę rozumiem, ale bez przesady). Rozumiem kongruencję, rozumiem dzielenie modulo, rozumiem nawet potęgowanie modulo, ale to tylko strzępki - a do zrozumienia większości "poważnych" algorytmów na sprawdzenie pierwszości trzeba dużo, dużo więcej.

Na szczęście pomimo wymienionych braków w wiedzy nadal mam działającą klawiaturę, a na niej klawisze Ctrl, C oraz V. Znalazłem więc kod w Pythonie, skopiowałem go sobie do PyCharm-a i sprawdziłem, czy działa.

Kurdę. Działa.

I nagle te wszytkie moje poprzednie wysiłki, żeby wyliczyć miliard liczb pierwszych, okazały się funta kłaków warte. Sito Eratostenesa jest fajne do pewnego momentu, a potem staje się po prostu niepraktyczne. Za to taki test Millera-Rabina... Przez chwilę poczułem się jak bokser wagi piórkowej obserwujący walkę Mechagodzilli z Optimusem Prime. Czy coś.

Żeby nie przedłużać, napiszę jeszcze, że za pomocą owego testu, na średniej klasy komputerze byłem w stanie potwierdzić pierwszość liczby o długości ponad tysiąca cyfr w czasie kilku-kilkunastu sekund, a sprawdzenie pierwszości liczby o czterech tysiącach cyfr zajęło około minuty.

Zadziwiające.

Kod jest dość prosty i przejrzysty. Jak to w Pythonie...

import time

def try_composite(a, d, n, s):
    if pow(a, d, n) == 1:
        return False
    for i in range(s):
        if pow(a, 2 ** i * d, n) == n - 1:
            return False
    return True  # n  is definitely composite


def is_prime(n, precision_for_huge_n=16):
    if n in known_primes or n in (0, 1):
        return True
    if any((n % p) == 0 for p in known_primes):
        return False
    d, s = n - 1, 0
    while not d % 2:
        d, s = d >> 1, s + 1
    if n < 1373653:
        return not any(try_composite(a, d, n, s) for a in (2, 3))
    if n < 25326001:
        return not any(try_composite(a, d, n, s) for a in (2, 3, 5))
    if n < 118670087467:
        if n == 3215031751:
            return False
        return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7))
    if n < 2152302898747:
        return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11))
    if n < 3474749660383:
        return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13))
    if n < 341550071728321:
        return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17))
    if n < 3825123056546413051:
        return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17, 19, 23))
    if n < pow(2, 64):
        return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37))
    # otherwise
    return not any(try_composite(a, d, n, s) for a in known_primes[:precision_for_huge_n])


known_primes = [2, 3]
known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]

number = int("""
2814112013697373133393152975842584191818662382013600787892419349345515176682276313810715094745633257
0741987893085350715373424450164188818017893905487094143918572575715657587064784183567470706746334971
8805305087541682162432568055582607111069194660746087305696536083057159024277493422686618396630918543
3462514537484258655982386235046029227507801410907163348439547781093397260096909677091843944555754221
1154773437602069796500670878849934780129772778785328074322365540209315718023104299231675884324570361
0411085096043976903845036551402234962538366575120716966169735273223611192684645475170173452701137914
8175107820821297628946795631098960767492250494834254073334414121627833939461539212528932010726136689
2936888156654916713951747104526637091757536037741568557665153138276137272816966926335296663637872865
3976994160910777718359333600268012451763345149043959832482383645725121940639143263563922560455604239
6004307799361927379900586400420763092320813392262492942076312933268033818471555255820639308889948665
5702024038158563135789497797670462618453279567257672892052623117520147862478133318340150844753867605
2661221734057972123741448580372535546302200953630100814586752470460461886203909355520619532824095189
5107040793284825095462530151872823997171764140663315804309008611942578380931064748991594407476328437
7858488254239211706149382940294832571629792993889406958773754489480811083452933943278084527297898341
3514019391241966179948879521032823811274221870063454114974365728723284342636934880487899347196240339
3967857676150371600196650252168250117793178488012000505422821362550520509209724459895852366827477851
6191905032548531150294031321789890051957511943013402772827303906836511205878950601987531218821877886
5702400729178418651858997778851030674394589610864525876641569282566417447061615330514485227388454963
5059255410606458427323864109506687636314447514269094932953219924212594695157655009158521173420923275
8820633276254086179630329620335725635536040560978321115475359089884338169197476158171616066205573070
0037719473001343181556075015902784216490142254457122454693679323497089495466842543641234778537619431
0030139080568383420772628618722646109707506566928102800033961704343991962002059794565527774913883237
7567927200655437686407921774415592782723508230928436835343966791502296761018342437878204200872740286
1721268457638873360576949122410986659257736066624146728015898860552348634588088222785550570630927634
9415034547677180618296352866263005509222254318459768194126727603047460344175581029298320171226355234
4396768163099191275742063348077190218754138915808715290491878293084121334009104197563130215404784366
0417844675773899863208358620799223408516263437540677116970732321398828494377912217198595360589790229
1781768286548287878180415060635460047164104095483777201737468873324068550430695826210304316336385311
3840934900213323724634633739774274058966738275442031285748745819603352320056372293195923692881713752
7670226045091173506950402501666775521493207364365419948847701036390937200575789998958077577512662111
3057905717449417222016070530243916116705990451304256206318289297738303095152430549772239514964821601
8386288614463019360177105467775031092630309947473976185762073734477254414271353624283608636693271576
3598304544797181671880163986954752514630565557184371791687566914032072497856858671852758660243960233
5283513944980064327030278104224144971883680541689784796267391476087696392191""".replace("\n", ""))

start_time = time.time()
print(is_prime(number))
print("--- %s seconds ---" % (time.time() - start_time))

Jeżeli ktoś dotarł w lekturze aż tutaj, to albo mu się bardzo nudzi, albo jej. Tak czy siak, przypomnę jeszcze, że powyżej \(2^{64}\) powyższy test jest niedeterministyczny (można regulować prawdopodobieństwo błędu za pomocą drugiego argumentu funkcji is_prime() - czym większy, tym pewniej, ale również tym wolniej).

https://xpil.eu/kmO0x 📋