Prime check : Miller-Rabin primality test (with bases)

Contributions to this software library are always welcome. Please ensure that you post program listings rather than .raw files. They give a reasonable idea of what your program does without having to load them into a DM42 and you can also include comments in your code. Check out the following link for a decoder/encoder: https://technical.swissmicros.com/decoders/dm42/

You can then copy/paste the listing and post it in "code" tags.
Boub65
Posts: 231
Joined: Tue Sep 12, 2017 4:34 pm
Location: Rabat, Morocco

Prime check : Miller-Rabin primality test (with bases)

Post by Boub65 »

Here is a DM42 program that checks if X register is a (big) prime number using the Miller-Rabin algorithm with a set of bases.
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
EDIT: BTW, I couldn't upload the .raw file due to some limitations 🤔
Sincèrement, Sincerely, 73,
Boubker

DM15L, DM41L, DM42 #00855 (domes upgraded), DM41X #00707
HP48SX (with dark screen), HP42s, HP32SII (1990 with fraction bug), HP41C/CV
TI-89 titanium, CASIO fx-cg50 and Numworks (to play with micropython)
User avatar
Walter
Posts: 3070
Joined: Tue May 02, 2017 11:13 am
Location: On a mission close to DRS, Germany

Re: Prime check : Miller-Rabin primality test (with bases)

Post by Walter »

This test is hard-wired in the WP43S. Feel free to check PRIME?. Just FYI.
WP43 SN00000, 34S, and 31S for obvious reasons; HP-35, 45, ..., 35S, 15CE, DM16L S/N# 00093, DM42β SN:00041
Boub65
Posts: 231
Joined: Tue Sep 12, 2017 4:34 pm
Location: Rabat, Morocco

Re: Prime check : Miller-Rabin primality test (with bases)

Post by Boub65 »

Hello Walter,
Walter wrote:
Wed Aug 11, 2021 11:32 pm
This test is hard-wired in the WP43S. Feel free to check PRIME?. Just FYI.
Up to what number is the test deterministic in WP43s?
or are you using random numbers as "base" to cover the full range of "integer"?
Sincèrement, Sincerely, 73,
Boubker

DM15L, DM41L, DM42 #00855 (domes upgraded), DM41X #00707
HP48SX (with dark screen), HP42s, HP32SII (1990 with fraction bug), HP41C/CV
TI-89 titanium, CASIO fx-cg50 and Numworks (to play with micropython)
User avatar
Over_score
Posts: 160
Joined: Fri May 05, 2017 9:37 pm
Location: France

Re: Prime check : Miller-Rabin primality test (with bases)

Post by Over_score »

We use gmplib-6.2.1, see here https://gmplib.org/manual/Number-Theoretic-Functions, and we use reps=25
DM42 SN00284 & SN03835 running C47, HP34C, HP41CV, HP42S, HP35s, WP34S, HP Prime
Boub65
Posts: 231
Joined: Tue Sep 12, 2017 4:34 pm
Location: Rabat, Morocco

Re: Prime check : Miller-Rabin primality test (with bases)

Post by Boub65 »

Merci Over_score,
Interesting... the Bailly-PSW test is my next RPN programming challenge 😜
Sincèrement, Sincerely, 73,
Boubker

DM15L, DM41L, DM42 #00855 (domes upgraded), DM41X #00707
HP48SX (with dark screen), HP42s, HP32SII (1990 with fraction bug), HP41C/CV
TI-89 titanium, CASIO fx-cg50 and Numworks (to play with micropython)
User avatar
Walter
Posts: 3070
Joined: Tue May 02, 2017 11:13 am
Location: On a mission close to DRS, Germany

Re: Prime check : Miller-Rabin primality test (with bases)

Post by Walter »

Boub65 wrote:
Thu Aug 12, 2021 1:58 am
Up to what number is the test deterministic in WP43s?
or are you using random numbers as "base" to cover the full range of "integer"?
As usual, the documentation is your friend. Please see p. 164 in the Reference Manual.
WP43 SN00000, 34S, and 31S for obvious reasons; HP-35, 45, ..., 35S, 15CE, DM16L S/N# 00093, DM42β SN:00041
Boub65
Posts: 231
Joined: Tue Sep 12, 2017 4:34 pm
Location: Rabat, Morocco

Re: Prime check : Miller-Rabin primality test (with bases)

Post by Boub65 »

@walter
@over_score

About GMP :
"GMP versions up to and including 6.1.2 did not use the Baillie-PSW primality test. In those older versions of GMP, this function performed reps Miller-Rabin tests."

So the older version of GMP where indeed limited to 3E24, like my algorithm...
... But, version 6.2.1 that is used in WP43S uses the Baillie-PSW test to check for primes...

You can read also here : https://mathworld.wolfram.com/Baillie-P ... yTest.html
"No examples of composite numbers passing the test are known, and as of June 13, 2009, Jeff Gilchrist has confirmed that there are no Baillie-PSW pseudoprimes up to 10^(17). However, the elliptic curve primality proving program PRIMO checks all intermediate probable primes with this test, and if any were composite, the certification would necessarily have failed. Based on the fact that this has not occurred in three years of usage, PRIMO author M. Martin estimates that there is no composite less than about 10000 digits that can fool this test."

I repeat "PRIMO author M. Martin estimates that there is no composite less than about 10000 digits that can fool this test."

So the WP43S "PRIME?" instruction should be ok up to 10000 digits according to literature (with no mathematical proof)...

@walter... may be you should update "the manual" 😜
Sincèrement, Sincerely, 73,
Boubker

DM15L, DM41L, DM42 #00855 (domes upgraded), DM41X #00707
HP48SX (with dark screen), HP42s, HP32SII (1990 with fraction bug), HP41C/CV
TI-89 titanium, CASIO fx-cg50 and Numworks (to play with micropython)
User avatar
Walter
Posts: 3070
Joined: Tue May 02, 2017 11:13 am
Location: On a mission close to DRS, Germany

Re: Prime check : Miller-Rabin primality test (with bases)

Post by Walter »

@Boub65
Boub65 wrote:
Thu Aug 12, 2021 12:42 pm
"GMP versions up to and including 6.1.2 did not use the Baillie-PSW primality test. In those older versions of GMP, this function performed reps Miller-Rabin tests."

So the older version of GMP where indeed limited to 3E24, like my algorithm...
... But, version 6.2.1 that is used in WP43S uses the Baillie-PSW test to check for primes...

You can read also here : https://mathworld.wolfram.com/Baillie-P ... yTest.html
"No examples of composite numbers passing the test are known, and as of June 13, 2009, Jeff Gilchrist has confirmed that there are no Baillie-PSW pseudoprimes up to 10^(17). However, the elliptic curve primality proving program PRIMO checks all intermediate probable primes with this test, and if any were composite, the certification would necessarily have failed. Based on the fact that this has not occurred in three years of usage, PRIMO author M. Martin estimates that there is no composite less than about 10000 digits that can fool this test."
...
So the WP43S "PRIME?" instruction should be ok up to 10000 digits according to literature (with no mathematical proof)...

@walter... may be you should update "the manual" 😜
Merci pour votre gentile remarque.

I think we agree that 10^10000 >> 3.3x10^24 >> 10^17. And please allow me repeating "PRIMO author M. Martin estimates that there is no composite less than about 10000 digits that can fool this test." Hmmh. The test author claims a working range up to 10^10000 on the base that his test works up to 10^17 - that's a bit far from a mathematical proof, isn't it? Pretty keen IMHO.

I won't update this section of the manual based on "should be ok". And BTW, looking for primes >296 digits is of limited value since you can't output such integers from the calculator so far.
WP43 SN00000, 34S, and 31S for obvious reasons; HP-35, 45, ..., 35S, 15CE, DM16L S/N# 00093, DM42β SN:00041
mcc
Posts: 277
Joined: Fri Jun 23, 2017 5:10 am

Re: Prime check : Miller-Rabin primality test (with bases)

Post by mcc »

Hi,

I may have found a little bug in above program:

4 and 5 are both marked as not being prime which is about 50% correct. ;)

May be the reason for that is that lable MRPRIM0 and MRPRIM1 are reached via GTO and
the according parts of the program are left with RTN?

I used https://technical.swissmicros.com/decoders/dm42/ to translate the program
into the binary RAW format. The byte count was correct.

Cheers
mcc
DM 42 - SN: 00373, Firmware release v.:3.22. / DMCP 3.24. as compiled by SwissMicros
Boub65
Posts: 231
Joined: Tue Sep 12, 2017 4:34 pm
Location: Rabat, Morocco

Re: Prime check : Miller-Rabin primality test (with bases)

Post by Boub65 »

Hi mcc
You are right, for every element of the base over 4 (first test goes to 4) a = 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, and 41 it says not prime...
43 is ok (first prime after the base)

Will check program and correct...
Sincèrement, Sincerely, 73,
Boubker

DM15L, DM41L, DM42 #00855 (domes upgraded), DM41X #00707
HP48SX (with dark screen), HP42s, HP32SII (1990 with fraction bug), HP41C/CV
TI-89 titanium, CASIO fx-cg50 and Numworks (to play with micropython)
Post Reply