Pchełki Python: Test Millera-Rabina

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…

  1. import time
  2.  
  3. def try_composite(a, d, n, s):
  4.     if pow(a, d, n) == 1:
  5.         return False
  6.     for i in range(s):
  7.         if pow(a, 2 ** i * d, n) == n - 1:
  8.             return False
  9.     return True  # n  is definitely composite
  10.  
  11.  
  12. def is_prime(n, precision_for_huge_n=16):
  13.     if n in known_primes or n in (0, 1):
  14.         return True
  15.     if any((n % p) == 0 for p in known_primes):
  16.         return False
  17.     d, s = n - 1, 0
  18.     while not d % 2:
  19.         d, s = d >> 1, s + 1
  20.     if n < 1373653:
  21.         return not any(try_composite(a, d, n, s) for a in (2, 3))
  22.     if n < 25326001:
  23.         return not any(try_composite(a, d, n, s) for a in (2, 3, 5))
  24.     if n < 118670087467:
  25.         if n == 3215031751:
  26.             return False
  27.         return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7))
  28.     if n < 2152302898747:
  29.         return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11))
  30.     if n < 3474749660383:
  31.         return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13))
  32.     if n < 341550071728321:
  33.         return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17))
  34.     if n < 3825123056546413051:
  35.         return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17, 19, 23))
  36.     if n < pow(2, 64):
  37.         return not any(try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37))
  38.     # otherwise
  39.     return not any(try_composite(a, d, n, s) for a in known_primes[:precision_for_huge_n])
  40.  
  41.  
  42. known_primes = [2, 3]
  43. known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]
  44.  
  45. number = int("""
  46. 2814112013697373133393152975842584191818662382013600787892419349345515176682276313810715094745633257
  47. 0741987893085350715373424450164188818017893905487094143918572575715657587064784183567470706746334971
  48. 8805305087541682162432568055582607111069194660746087305696536083057159024277493422686618396630918543
  49. 3462514537484258655982386235046029227507801410907163348439547781093397260096909677091843944555754221
  50. 1154773437602069796500670878849934780129772778785328074322365540209315718023104299231675884324570361
  51. 0411085096043976903845036551402234962538366575120716966169735273223611192684645475170173452701137914
  52. 8175107820821297628946795631098960767492250494834254073334414121627833939461539212528932010726136689
  53. 2936888156654916713951747104526637091757536037741568557665153138276137272816966926335296663637872865
  54. 3976994160910777718359333600268012451763345149043959832482383645725121940639143263563922560455604239
  55. 6004307799361927379900586400420763092320813392262492942076312933268033818471555255820639308889948665
  56. 5702024038158563135789497797670462618453279567257672892052623117520147862478133318340150844753867605
  57. 2661221734057972123741448580372535546302200953630100814586752470460461886203909355520619532824095189
  58. 5107040793284825095462530151872823997171764140663315804309008611942578380931064748991594407476328437
  59. 7858488254239211706149382940294832571629792993889406958773754489480811083452933943278084527297898341
  60. 3514019391241966179948879521032823811274221870063454114974365728723284342636934880487899347196240339
  61. 3967857676150371600196650252168250117793178488012000505422821362550520509209724459895852366827477851
  62. 6191905032548531150294031321789890051957511943013402772827303906836511205878950601987531218821877886
  63. 5702400729178418651858997778851030674394589610864525876641569282566417447061615330514485227388454963
  64. 5059255410606458427323864109506687636314447514269094932953219924212594695157655009158521173420923275
  65. 8820633276254086179630329620335725635536040560978321115475359089884338169197476158171616066205573070
  66. 0037719473001343181556075015902784216490142254457122454693679323497089495466842543641234778537619431
  67. 0030139080568383420772628618722646109707506566928102800033961704343991962002059794565527774913883237
  68. 7567927200655437686407921774415592782723508230928436835343966791502296761018342437878204200872740286
  69. 1721268457638873360576949122410986659257736066624146728015898860552348634588088222785550570630927634
  70. 9415034547677180618296352866263005509222254318459768194126727603047460344175581029298320171226355234
  71. 4396768163099191275742063348077190218754138915808715290491878293084121334009104197563130215404784366
  72. 0417844675773899863208358620799223408516263437540677116970732321398828494377912217198595360589790229
  73. 1781768286548287878180415060635460047164104095483777201737468873324068550430695826210304316336385311
  74. 3840934900213323724634633739774274058966738275442031285748745819603352320056372293195923692881713752
  75. 7670226045091173506950402501666775521493207364365419948847701036390937200575789998958077577512662111
  76. 3057905717449417222016070530243916116705990451304256206318289297738303095152430549772239514964821601
  77. 8386288614463019360177105467775031092630309947473976185762073734477254414271353624283608636693271576
  78. 3598304544797181671880163986954752514630565557184371791687566914032072497856858671852758660243960233
  79. 5283513944980064327030278104224144971883680541689784796267391476087696392191""".replace("\n", ""))
  80.  
  81. start_time = time.time()
  82. print(is_prime(number))
  83. 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).

Dodaj komentarz

avatar
  Subscribe  
Powiadom o
%d bloggers like this: