```  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.
84 #
85 #     This program is free software: you can redistribute it and/or modify
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()
```