1 #!/usr/bin/env python3
  2 #
  3 #============================================================================
  4 #
  5 # NAME
  6 #
  7 #     millerRabin.py
  8 #
  9 # DESCRIPTION
 10 #
 11 #     Miller-Rabin test for primality.
 12 #
 13 # USAGE
 14 #
 15 #     To test n for primality using random number a, 
 16 #
 17 #     python millerRabin.py n a 1
 18 #
 19 # EXAMPLES
 20 #
 21 #     python millerRabin.py 104729 47340 1
 22 #         Testing n = 104729 for primality using random integer x = 47340
 23 #         Remove powers of 2
 24 #         k = 1 n1 = 52364
 25 #         k = 2 n1 = 26182
 26 #         k = 3 n1 = 13091
 27 #         q = 3 k = 13091
 28 #         Generate and test the sequence x^(q^k) mod n. 
 29 #         j = 0 y = 87711
 30 #         j = 1 y = 36639
 31 #         j = 2 y = 104728
 32 #         Probably prime.  Found y = n-1 or j = 0 and y = 1  y = 104728 j = 2
 33 #
 34 #     python millerRabin.py 9999612 32423 1
 35 #         Testing n = 9999612 for primality using random integer x = 32423
 36 #         Remove powers of 2
 37 #         q = 0 k = 9999611
 38 #         Generate and test the sequence x^(q^k) mod n. 
 39 #         Definitely composite.  Fell off the end of the loop.
 40 #
 41 #     ./millerRabin.py 8398 1884460498967805432001612672369307101507474835976431925948333387748670120353629453261347843140212808\
 42 #                           570505767386771290423087216156597588216186445958479269565424431335013281 122 1
 43 #         Testing n = 188446049896780543200161267236930710150747483597643192594833338774867012035362945326134784314021280857050576\
 44 #                     7386771290423087216156597588216186445958479269565424431335013281 for primality using random integer x = 122
 45 #         Remove powers of 2
 46 #         k = 1 n1 = 942230249483902716000806336184653550753737417988215962974166693874335060176814726630673921570106404285252883693\
 47 #                     385645211543608078298794108093222979239634782712215667506640
 48 #         k = 2 n1 = 47111512474195135800040316809232677537686870899410798148708334693716753008840736331533696078505320214262644184\
 49 #                     6692822605771804039149397054046611489619817391356107833753320
 50 #         k = 3 n1 = 235557562370975679000201584046163387688434354497053990743541673468583765044203681657668480392526601071313220923\
 51 #                     346411302885902019574698527023305744809908695678053916876660
 52 #         k = 4 n1 = 1177787811854878395001007920230816938442171772485269953717708367342918825221018408288342401962633005356566104616\
 53 #                     73205651442951009787349263511652872404954347839026958438330
 54 #         k = 5 n1 = 5888939059274391975005039601154084692210858862426349768588541836714594126105092041441712009813165026782830523083\
 55 #                     6602825721475504893674631755826436202477173919513479219165
 56 #         q = 5 k = 5888939059274391975005039601154084692210858862426349768588541836714594126105092041441712009813165026782830523083\
 57 #                     6602825721475504893674631755826436202477173919513479219165
 58 #         Generate and test the sequence x^(q^k) mod n. 
 59 #         j = 0 y = 4562627375418921997495128389472904369506849162624397308228643567229422125069828207111411946939507039364910257582\
 60 #                     53679056110147196100972654986996407565104255499753038631160
 61 #         j = 1 y = 181921118629873824198891977536932720607952192628256115040458970624355555142399332847486840509884746125806162856\
 62 #                     2013490792783220259040379029588830209459993575609245826016160
 63 #         j = 2 y = 1675291436909496660109561278257723410018131370324303168442668644288721609236515790057247572821758832994480351123\
 64 #                     772012853725200818096891590375783891432523321208439264772913
 65 #         j = 3 y = 1884460498967805432001612672369307101507474835976431925948333387748670120353629453261347843140212808570505767386\
 66 #                     771290423087216156597588216186445958479269565424431335013280
 67 #         Probably prime.  Found y = n-1 or j = 0 and y = 1  y = 18844604989678054320016126723693071015074748359764319259483333877486\
 68 #                     70120353629453261347843140212808570505767386771290423087216156597588216186445958479269565424431335013280 j = 3
 69 #
 70 # AUTHOR
 71 #
 72 #     Sean E. O'Connor        3 Apr 2015  Version 1.0 released.
 73 #
 74 # NOTES
 75 #
 76 #    Python interpreter:               https://www.python.org/
 77 #    Python tutorial and reference:    htttp://docs.python.org/lib/lib.html
 78 #    Python regular expression howto:  https://docs.python.org/3.7/howto/regex.html
 79 #
 80 # LEGAL
 81 #
 82 #     millerRabin.py Version 1.0 - A Python program to do the Miller-Rabin test for primality.
 83 #     Copyright (C) 2017-2021 by Sean Erik O'Connor.  All Rights Reserved.
 84 #
 85 #     This program is free software: you can redistribute it and/or modify
 86 #     it under the terms of the GNU General Public License as published by
 87 #     the Free Software Foundation, either version 3 of the License, or
 88 #     (at your option) any later version.
 89 #
 90 #     This program is distributed in the hope that it will be useful,
 91 #     but WITHOUT ANY WARRANTY; without even the implied warranty of
 92 #     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 93 #     GNU General Public License for more details.
 94 #
 95 #     You should have received a copy of the GNU General Public License
 96 #     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 97 #
 98 #     The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com
 99 #     with !DOT! replaced by . and the !AT! replaced by @
100 #
101 #============================================================================
102 
103 import sys
104 
105 def main():
106 
107     print( "Running Python Version {0:d}.{1:d}.{2:d}".format( sys.version_info[ 0 ], sys.version_info[ 1 ], sys.version_info[ 2 ] ) )
108 
109     if len( sys.argv ) != 4:
110         print( "Usage: python MillerRabin.py n x d" )
111         print( "    n   the number to be tested for primality" )
112         print( "    x   any random number" )
113         print( "    d   nonzero turns on debug prints" )
114         sys.exit()
115     else:
116         n = int( sys.argv[ 1 ] )
117         x = int( sys.argv[ 2 ] )
118         if int( sys.argv[ 3 ] ) == 0:
119             debug = False
120         else:
121             debug = True
122 
123         if debug:
124             print( "Testing n = {0:d} for primality using random integer x = {1:d}".format( n, x ) )
125 
126         if MillerRabin( n, x, debug ):
127             print( "Probably prime" )
128         else:
129             print( "Definitely composite" )
130 
131         sys.exit( 0 )
132 
133 def MillerRabin( n, x, debug ):
134     n1 = n - 1
135 
136     # Remove powers of 2. 
137     k = 0
138     even = True
139 
140     if debug:
141         print( "Remove powers of 2" )
142 
143     while n1 % 2 == 0:
144         n1 //= 2
145         k += 1
146         if debug:
147             print( "k = {0:d} n1 = {1:d}".format( k, n1 ))
148 
149     q = n1
150 
151     if debug:
152         print( "q = {0:d} k = {1:d}".format( k, n1 ) )
153 
154     y = pow( x, q, n )
155 
156     if debug:
157         if debug:
158             print( "Generate and test the sequence x^(q^k) mod n. " ) ;
159 
160     for j in range (0,k):
161         if debug:
162             print( "j = {0:d} y = {1:d}".format( j, y ))
163 
164         if (j == 0 and y == 1) or y == n-1:
165             if debug:
166                 print( "Probably prime.  Found y = n-1 or j = 0 and y = 1  y = {0:d} j = {1:d}".format( y, j ) )
167             return True
168 
169         if j > 0 and y == 1:
170             if debug:
171                 print( "Definitely composite.  Found j > 0 and y = 1" )
172             return False
173 
174         y = pow( y, 2, n )
175 
176     if debug:
177         print( "Definitely composite.  Fell off the end of the loop." )
178     return False
179 
180 # Call the main.
181 if __name__ == '__main__':
182     main()