The algorithm works up to 3 E 24 (well over 2 ^ 64).
More on the Miller-Rabin primality test here : https://en.wikipedia.org/wiki/Miller%E2 ... ality_test
The program uses the following 13 elements base : a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, and 41 the primality test is DETERMINISTIC up to 3.317.044.064.679.887.385.961.981 (approx 3 E 24)
3.317.044.064.679.887.385.961.981, is the biggest prime number that passes the test (found in 1 minute and 35 seconds @80MHz CPU on USB).
2^64 - 59 is also a prime number (found in 1 minute 13 seconds @80MHz CPU on USB).
2^80 - 65 is also a prime number (found in 2 minutes @80MHz CPU on USB).
Comparison with DM41X :
4.938.271.579 (DM41X limit as seen here viewtopic.php?f=30&t=3072) is found to be prime instantaneously!
if you try the composite number, 4.938.271.579 x 5.027.672.681 = 24.828.013.109.097.033.299 the program finds it is not prime in 4 seconds!
In the program, you can also find two interesting routines :
- POW that calculates (Z**Y) mod X for big numbers
- AXB%MOD that calculates (Z x Y) mod X for big numbers
Code: Select all
00 { 780-Byte Prgm }
01▸LBL "MILLRAB"
02 STO "N"
@ n = 4 then not prime
03 4
04 X<>Y
05 X=Y?
06 GTO "MRPRIM0"
@ n < 4 then prime (1, 2 and 3)
07 X<Y?
08 GTO "MRPRIM1"
@ we store the base numbers (a)
@ 13 base numbers 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, and 41
09 14.00001
10 STO 00
11 41
12 STO 01
13 37
14 STO 02
15 31
16 STO 03
17 29
18 STO 04
19 23
20 STO 05
21 19
22 STO 06
23 17
24 STO 07
25 13
26 STO 08
27 11
28 STO 09
29 7
30 STO 10
31 5
32 STO 11
33 3
34 STO 12
35 2
36 STO 13
@ d = n - 1
37 RCL "N"
38 1
39 -
40 STO "D"
@ we find r such that n = 2^d * r + 1 (r>=1)
@ while (d % 2 = 0)
41▸LBL 00
42 RCL "D"
43 2
44 MOD
45 X≠0?
46 GTO "MRMAIN"
@ d = d / 2
47 2
48 STO÷ "D"
49 GTO 00
@ loop for each base number
50▸LBL "MRMAIN"
51 DSE 00
52 GTO 05
@ are we at index = 0, we have finished, number is prime
53 GTO "MRPRIM1"
54▸LBL 05
@ base number
55 RCL IND 00
56 RCL "D"
57 RCL "N"
58 XEQ "MRTEST"
@ is the test composite
59 X=0?
@ then not prime
60 GTO "MRPRIM0"
@ if the test is not composite
@ we loop for next base nummber
61 GTO "MRMAIN"
62▸LBL "MRPRIM0"
63 "NOT PRIME"
64 0
65 BEEP
66 RTN
67▸LBL "MRPRIM1"
68 "PRIME"
69 1
70 BEEP
71 RTN
@ Miller-Rabin test
@ a is Z
@ d is Y
@ n is X
72▸LBL "MRTEST"
73 LSTO "N"
74 R↓
75 LSTO "D"
76 R↓
@ x = (a ** d) % n
77 LSTO "A"
78 RCL "D"
79 RCL "N"
80 XEQ "POW"
81 LSTO "X"
@ x = 1 then true
82 1
83 X=Y?
84 GTO "MRTEST1"
@ x = n - 1 then true
85 R↓
86 RCL "N"
87 1
88 -
89 X=Y?
90 GTO "MRTEST1"
@ while (d != n - 1) loop
91▸LBL 07
92 RCL "D"
93 RCL "N"
94 1
95 -
96 X=Y?
@ if we exit the loop, d = n - 1 then false
97 GTO "MRTEST0"
@ x = (x *x) % n
98 RCL "X"
99 X^2
100 1ᴇ34
101 X≤Y?
102 GTO 06
103 R↓
104 RCL "N"
105 MOD
106 GTO 08
107▸LBL 06
108 RCL "X"
109 RCL "X"
110 RCL "N"
111 XEQ "AXB%MOD"
112▸LBL 08
113 LSTO "X"
@ d *=2
114 2
115 STO× "D"
116 RCL "D"
117 1ᴇ34
118 X≤Y?
119 GTO "OVFLOW"
@ x = 1 then false
120 1
121 RCL "X"
122 X=Y?
123 GTO "MRTEST0"
@ x = n - 1 then true
124 RCL "N"
125 1
126 -
127 X=Y?
128 GTO "MRTEST1"
@ we loop
129 GTO 07
@ return composite
130▸LBL "MRTEST0"
131 0
132 RTN
133▸LBL "MRTEST1"
134 1
135 RTN
@ computes (x ** y) % mod for big numbers
@ input :
@ mod in X
@ Y in Y
@ X in Z
@ output :
@ axb%mod in X
136▸LBL "POW"
137 LSTO "MOD"
138 R↓
139 LSTO "Y"
140 R↓
141 LSTO "X"
@ res = 1
142 1
143 LSTO "RES"
@ x = x % mod
144 RCL "X"
145 RCL "MOD"
146 MOD
147 LSTO "X"
@ loop while y > 0
148▸LBL 10
149 RCL "Y"
150 X≤0?
@ if y <= 0 we exit with result = res
151 GTO 96
@ is y even ?
152 RCL "Y"
153 2
154 MOD
155 X=0?
@ if Y is even we skip next calculation
156 GTO 11
@ y is odd
@ res = (res * x) % mod
157 RCL "RES"
158 RCL "X"
159 ×
160 1ᴇ34
161 X≤Y?
162 GTO 12
163 R↓
164 RCL "MOD"
165 MOD
166 GTO 13
167▸LBL 12
168 RCL "RES"
169 RCL "X"
170 RCL "MOD"
171 XEQ "AXB%MOD"
172▸LBL 13
173 LSTO "RES"
@ Y is even si we substract 1 before dividing by 2
174 1
175 STO- "Y"
176▸LBL 11
@ y >>= 1
177 2
178 STO÷ "Y"
@ x = (x *x) % mod
179 RCL "X"
180 RCL "X"
181 ×
182 1ᴇ34
183 X≤Y?
184 GTO 14
185 R↓
186 RCL "MOD"
187 MOD
188 GTO 15
189▸LBL 14
190 RCL "X"
191 RCL "X"
192 RCL "MOD"
193 XEQ "AXB%MOD"
194▸LBL 15
195 LSTO "X"
196 GTO 10
197▸LBL 96
198 RCL "RES"
199 RTN
@ computes (a x b) % mod for big numbers
@ input :
@ mod in X
@ b in Y
@ a in Z
@ output :
@ axb%mod in X
200▸LBL "AXB%MOD"
201 LSTO "MOD"
202 R↓
203 LSTO "B"
204 R↓
205 LSTO "A"
206 0
207 LSTO "RES"
@ main loop until a = 0
208▸LBL 20
209 RCL "A"
210 X=0?
@ a = 0 we exit with res
211 GTO 97
@ is a even?
212 2
213 MOD
214 X=0?
215 GTO 21
@ a is odd, res = (res + b) % mod
216 RCL "RES"
217 RCL "B"
218 +
219 1ᴇ34
220 X≤Y?
221 GTO "OVFLOW"
222 R↓
223 RCL "MOD"
224 MOD
225 LSTO "RES"
@ a was odd, we retreive 1 so we can divide by 2
226 1
227 STO- "A"
228▸LBL 21
@ a >>= 1
@ knowing that a is even we just divide by 2
229 2
230 STO÷ "A"
@ b = (b * 2) % mod
231 RCL "B"
232 2
233 ×
234 1ᴇ34
235 X≤Y?
236 GTO "OVFLOW"
237 R↓
238 RCL "MOD"
239 MOD
240 LSTO "B"
@ we loop until a = 0
241 GTO 20
@ end, return res
242▸LBL 97
243 RCL "RES"
244 RTN
245▸LBL "OVFLOW"
246 "OVERFLOW"
247 X<>Y
248 BEEP
249 STOP
250 END