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-2024 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 seanerikoconnor!AT!gmail!DOT!com
99# with !DOT! replaced by . and the !AT! replaced by @
100#
101#============================================================================
102
103import sys
104
105def 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
133def 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.
181if __name__ == '__main__':
182 main()
