Miller–Rabin primality test: Difference between revisions
Content deleted Content added
VincentARM (talk | contribs) add task to arm assembly raspberry pi |
|||
| (10 intermediate revisions by 4 users not shown) | |||
Line 63:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
</pre>
=={{header|AArch64 Assembly}}==
{{works with|as|Raspberry Pi 3B version Buster 64 bits <br> or android 64 bits with application Termux }}
<syntaxhighlight lang AArch64 Assembly>
/* ARM assembly AARCH64 Raspberry PI 3B */
/* program testmiller64B.s */
// optimisation : one routine
/************************************/
/* Constantes */
/************************************/
/* for this file see task include a file in language AArch64 assembly*/
.include "../includeConstantesARM64.inc"
.equ NBLOOP, 5 // loop number change this if necessary
// if modify, thinck to add test to little prime value
// in routine
//.include "../../ficmacros64.inc" // use for developper debugging
/*******************************************/
/* Initialized data */
/*******************************************/
.data
szMessStartPgm: .asciz "Program 64 bits start \n"
szMessEndPgm: .asciz "Program normal end.\n"
szMessPrime: .asciz " is prime !!!.\n"
szMessNotPrime: .asciz " is not prime !!!.\n"
szCarriageReturn: .asciz "\n"
/*******************************************/
/* UnInitialized data */
/*******************************************/
.bss
.align 4
sZoneConv: .skip 24
/*******************************************/
/* code section */
/*******************************************/
.text
.global main
main: // program start
ldr x0,qAdrszMessStartPgm // display start message
bl affichageMess
ldr x4,iStart // start number
ldr x5,iLimit // end number
tst x4,#1
cinc x4,x4,eq // start with odd number
1:
mov x0,x4
bl isPrimeMiller // test miller rabin
cmp x0,#0
beq 2f
mov x0,x4
ldr x1,qAdrsZoneConv
bl conversion10 // decimal conversion
ldr x0,qAdrsZoneConv
bl affichageMess
ldr x0,qAdrszMessPrime
bl affichageMess
b 3f
2:
ldr x0,qAdrszMessNotPrime
// bl affichageMess
3:
add x4,x4,#2
cmp x4,x5
blo 1b
ldr x0,qAdrszMessEndPgm // display end message
bl affichageMess
100: // standard end of the program
mov x0, #0 // return code
mov x8, #EXIT // request to exit program
svc 0 // perform system call
qAdrszMessStartPgm: .quad szMessStartPgm
qAdrszMessEndPgm: .quad szMessEndPgm
qAdrszCarriageReturn: .quad szCarriageReturn
qAdrsZoneConv: .quad sZoneConv
qAdrszMessPrime: .quad szMessPrime
qAdrszMessNotPrime: .quad szMessNotPrime
//iStart: .quad 0x0
//iLimit: .quad 0x100
iStart: .quad 0xFFFFFFFFFFFFFF00
iLimit: .quad 0xFFFFFFFFFFFFFFF0
//iStart: .quad 341550071728360
//iLimit: .quad 341550071728380
//359341550071728361
/***************************************************/
/* test miller rabin algorithme wikipedia */
/* unsigned */
/***************************************************/
/* x0 contains number */
/* x1 contains parameter */
/* x0 return 1 if prime 0 if composite */
isPrimeMiller:
stp x1,lr,[sp,-16]! // TODO: save à completer
stp x2,x3,[sp,-16]!
stp x4,x5,[sp,-16]!
stp x6,x7,[sp,-16]!
stp x8,x9,[sp,-16]!
cmp x0,#1 // control 0 or 1
csel x0,xzr,x0,ls
bls 100f
cmp x0,#2 // control = 2
mov x1,1
csel x0,x1,x0,eq
beq 100f
cmp x0,#3 // control = 3
csel x0,x1,x0,eq
beq 100f
cmp x0,#5 // control = 5
csel x0,x1,x0,eq
beq 100f
cmp x0,#7 // control = 7
csel x0,x1,x0,eq
beq 100f
cmp x0,#11 // control = 11
csel x0,x1,x0,eq
beq 100f
tst x0,#1
csel x0,xzr,x0,eq // even
beq 100f
mov x4,x0 // N
sub x3,x0,#1 // D
mov x2,#2
mov x6,#0 // S
1: // compute D * 2 power S
lsr x3,x3,#1 // D= D/2
add x6,x6,#1 // increment S
tst x3,#1 // D even ?
beq 1b
2:
mov x8,#0 // loop counter
sub x5,x0,#3
mov x7,3
3:
mov x0,x7
mov x1,x3 // exposant = D
mov x2,x4 // modulo N
bl moduloPur64
cmp x0,#1
beq 5f
sub x1,x4,#1 // n -1
cmp x0,x1
beq 5f
sub x9,x6,#1 // S - 1
4:
mov x2,x0
mul x0,x2,x2 // compute square lower
umulh x1,x2,x2 // compute square upper
mov x2,x4 // and compute modulo N
bl division64R2023
mov x0,x2
cmp x0,#1
csel x0,xzr,x0,eq // composite
beq 100f
sub x1,x4,#1 // n -1
cmp x0,x1
beq 5f
subs x9,x9,#1
bge 4b
mov x0,#0 // composite
b 100f
5:
add x7,x7,2
add x8,x8,#1
cmp x8,NBLOOP
blt 3b
mov x0,#1 // prime
100:
ldp x8,x9,[sp],16
ldp x6,x7,[sp],16
ldp x4,x5,[sp],16
ldp x2,x3,[sp],16
ldp x1,lr,[sp],16 // TODO: retaur à completer
ret
/********************************************************/
/* Calcul modulo de b puissance e modulo m */
/* Exemple 4 puissance 13 modulo 497 = 445 */
/********************************************************/
/* x0 nombre */
/* x1 exposant */
/* x2 modulo */
moduloPur64:
stp x1,lr,[sp,-16]! // save registres
stp x3,x4,[sp,-16]! // save registres
stp x5,x6,[sp,-16]! // save registres
stp x7,x8,[sp,-16]! // save registres
stp x9,x10,[sp,-16]! // save registres
cbz x0,100f
cbz x1,100f
mov x8,x0
mov x7,x1
mov x10,x2 // modulo
mov x6,1 // resultat
udiv x4,x8,x10
msub x9,x4,x10,x8 // contient le reste
1:
tst x7,1
beq 2f
mul x4,x9,x6
umulh x5,x9,x6
mov x6,x4
mov x0,x6
mov x1,x5
mov x2,x10
bl division64R2023
mov x6,x2
2:
mul x8,x9,x9
umulh x5,x9,x9
mov x0,x8
mov x1,x5
mov x2,x10
bl division64R2023
mov x9,x2
lsr x7,x7,1
cbnz x7,1b
mov x0,x6 // result
cmn x0,0 // clear carry not error
100:
ldp x9,x10,[sp],16 // restaur des 2 registres
ldp x7,x8,[sp],16 // restaur des 2 registres
ldp x5,x6,[sp],16 // restaur des 2 registres
ldp x3,x4,[sp],16 // restaur des 2 registres
ldp x1,lr,[sp],16 // restaur des 2 registres
ret // retour adresse lr x30
/***************************************************/
/* division number 128 bits in 2 registers by number 64 bits */
/* unsigned */
/***************************************************/
/* x0 contains lower part dividende */
/* x1 contains upper part dividende */
/* x2 contains divisor */
/* x0 return lower part quotient */
/* x1 return upper part quotient */
/* x2 return remainder */
division64R2023:
stp x3,lr,[sp,-16]!
stp x4,x5,[sp,-16]!
stp x6,x7,[sp,-16]!
mov x4,x2 // save divisor
mov x5,#0 // init upper part divisor
mov x2,x0 // save dividende
mov x3,x1
mov x0,#0 // init result
mov x1,#0
mov x6,#0 // init shift counter
1: // loop shift divisor
cmp x5,#0 // upper divisor <0
blt 2f
cmp x5,x3
bhi 2f
blo 11f
cmp x4,x2
bhi 2f // new divisor > dividende
11:
lsl x5,x5,#1 // shift left one bit upper divisor
tst x4,#0x8000000000000000
lsl x4,x4,#1 // shift left one bit lower divisor
orr x7,x5,#1
csel x5,x7,x5,ne // move bit 63 lower on upper
add x6,x6,#1 // increment shift counter
b 1b
2: // loop 2
lsl x1,x1,#1 // shift left one bit upper quotient
tst x0,#0x8000000000000000
lsl x0,x0,#1 // shift left one bit lower quotient
orr x7,x1,#1
csel x1,x7,x1,ne // move bit 63 lower on upper
cmp x5,x3 // compare divisor and dividende
bhi 3f
blo 21f
cmp x4,x2
bhi 3f
21:
subs x2,x2,x4 // < sub divisor from dividende lower
sbc x3,x3,x5 // and upper
orr x0,x0,#1 // move 1 on quotient
3:
lsr x4,x4,#1 // shift right one bit upper divisor
tst x5,1
lsr x5,x5,#1 // and lower
orr x7,x4,#0x8000000000000000 // move bit 0 upper to 31 bit lower
csel x4,x7,x4,ne // move bit 0 upper to 63 bit lower
subs x6,x6,#1 // decrement shift counter
bge 2b // if > 0 loop 2
100:
ldp x6,x7,[sp],16
ldp x4,x5,[sp],16
ldp x3,lr,[sp],16
ret
/***************************************************/
/* ROUTINES INCLUDE */
/***************************************************/
.include "../includeARM64.inc"
</syntaxhighlight>
{{Out}}
<pre>
Program 64 bits start
18446744073709551427 is prime !!!.
18446744073709551437 is prime !!!.
18446744073709551521 is prime !!!.
18446744073709551533 is prime !!!.
18446744073709551557 is prime !!!.
Program normal end.
</pre>
=={{header|Ada}}==
Line 5,367 ⟶ 5,680:
False
>>> </pre>
=={{header|Quackery}}==
Translated from the current (22 Sept 2023) pseudocode from [[wp:Miller-Rabin primality test#Miller–Rabin_test|Wikipedia]], which is:
'''Input #1''': ''n'' > 2, an odd integer to be tested for primality
'''Input #2''': ''k'', the number of rounds of testing to perform
'''Output''': “''composite''” if ''n'' is found to be composite, “''probably prime''” otherwise
'''let''' ''s'' > 0 and ''d'' odd > 0 such that ''n'' − 1 = 2<sup>''s''</sup>''d'' <u># by factoring out powers of 2 from ''n'' − 1</u>
'''repeat''' ''k'' '''times''':
''a'' ← random(2, ''n'' − 2) # ''n'' is always a probable prime to base 1 and ''n'' − 1
''x'' ← ''a''<sup>''d''</sup> mod ''n''
'''repeat''' ''s'' '''times''':
''y'' ← ''x''<sup>2</sup> mod ''n''
'''if''' ''y'' = 1 and ''x'' ≠ 1 and ''x'' ≠ ''n'' − 1 '''then''' <u># nontrivial square root of 1 modulo ''n''</u>
'''return''' “''composite''”
''x'' ← ''y''
'''if''' ''y'' ≠ 1 '''then'''
'''return''' “''composite''”
'''return''' “''probably prime''”
As Quackery does not have variables, I have included a "builder" (i.e. a compiler extension) to provide basic global variables, in order to facilitate comparison to the pseudocode.
<code>var myvar</code> will create a variable called <code>myvar</code> initialised with the value <code>0</code>, and additionally the words <code>->myvar</code> and <code>myvar-></code>, which assign a value to <code>myvar</code> and return the value of <code>myvar</code> respectively.
The variable <code>c</code> is set to <code>true</code> if the number being tested is composite. This is included as Quackery does not have a <code>break</code> that can exit nested iterative loops.
<code>eratosthenes</code> and <code>isprime</code> are defined at [[Sieve of Eratosthenes#Quackery]].
<code>**mod</code> is defined at [[Modular exponentiation#Quackery]].
<syntaxhighlight lang="Quackery"> [ nextword
dup $ "" = if
[ $ "var needs a name"
message put bail ]
$ " [ stack 0 ] is "
over join
$ " [ "
join over join
$ " replace ] is ->"
join over join
$ " [ "
join over join
$ " share ] is "
join swap join
$ "-> " join
swap join ] builds var ( [ $ --> [ $ )
10000 eratosthenes
[ 1 & not ] is even ( n --> b )
var n var d var s var a
var t var x var y var c
[ dup 10000 < iff isprime done
dup even iff
[ drop false ] done
dup ->n
[ 1 - 0 swap
[ dup even while
1 >> dip 1+ again ] ]
->d ->s
false ->c
20 times
[ n-> 2 - random 2 + ->a
s-> ->t
a-> d-> n-> **mod ->x
s-> times
[ x-> 2 n-> **mod ->y
y-> 1 =
x-> 1 !=
x-> n-> 1 - !=
and and iff
[ true ->c conclude ]
done
y-> ->x ]
c-> iff conclude done
y-> 1 != if
[ true ->c conclude ] ]
c-> not ] is prime ( n --> b )
say "Primes < 100:" cr
100 times [ i^ prime if [ i^ echo sp ] ]
cr cr
[] 1000 times
[ i^ 9223372036854774808 + dup
prime iff join else drop ]
say "There are "
dup size echo
say " primes between 9223372036854774808 and 9223372036854775808:"
cr cr
witheach [ echo i^ 1+ 4 mod iff sp else cr ]
cr cr
4547337172376300111955330758342147474062293202868155909489 dup echo
prime iff
[ say " is prime." ]
else
[ say " is not prime." ]
cr cr
4547337172376300111955330758342147474062293202868155909393 dup echo
prime iff
[ say "is prime." ]
else
[ say " is not prime." ]</syntaxhighlight>
{{out}}
<pre>Primes < 100:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
There are 23 primes between 9223372036854774808 and 9223372036854775808:
9223372036854774893 9223372036854774917 9223372036854774937 9223372036854774959
9223372036854775057 9223372036854775073 9223372036854775097 9223372036854775139
9223372036854775159 9223372036854775181 9223372036854775259 9223372036854775279
9223372036854775291 9223372036854775337 9223372036854775351 9223372036854775399
9223372036854775417 9223372036854775421 9223372036854775433 9223372036854775507
9223372036854775549 9223372036854775643 9223372036854775783
4547337172376300111955330758342147474062293202868155909489 is prime.
4547337172376300111955330758342147474062293202868155909393 is not prime.</pre>
=={{header|Racket}}==
Line 5,532 ⟶ 5,969:
for 1──►10000, K=10, Miller─Rabin primality test found 1229 primes {out of 1229}.
</pre>
Above program has 2 issues:
1) The REXX BIF '''random()''' is limited to the range 1 to 99999. Thus the statement '''?= random(2, nL)''' runs into trouble very soon.
2) In the statement '''x= ?**d // n''' (meaning ?^d mod n) the part '''?**d''' overflows already on modest values of ? and d. REXX does not have 'modular exponentation'.
These two issues prevent the testing of bigger primes.
Below version implements the full Miller-Rabin primality test included witnesses for deterministic test for x < 3317044064679887385961981. Improved functions '''Rand''' and '''Powermod''' (see elsewhere on RosettaCode) are added. Also some small functions '''Floor''' and '''Whole''' are included.
<syntaxhighlight lang="rexx">
parse version version
say version
glob. = ''
numeric digits 1000
say '25 numbers of the form 2^n-1, mostly Mersenne primes'
say 'Up to about 25 digits deterministic, above probabilistic'
say
mm = '2 3 5 7 11 13 17 19 23 31 61 89 97 107 113 127 131 521 607 1279 2203 2281 2293 3217 3221'
do nn = 1 to 25
a = Word(mm,nn); b = 2**a-1; l = Length(b)
call time('r'); p = Prime(b); e = time('e')
if l > 25 then
b = Left(b,12)'...'Right(b,13)
select
when p = 0 then
p = 'not'
when l < 26 then
p = 'for sure'
otherwise
p = 'probable'
end
say '2^('a'-1)' '=' b '('l' digits) is' p 'prime' '('format(e,,3) 'seconds)'
end
return
Prime:
/* Is a number prime? */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 2 then
return 'X'
/* Low primes also used as deterministic witnesses */
w1 = ' 2 3 5 7 11 13 17 19 23 29 31 37 41 '
/* Fast values */
w = x
if Pos(' 'w' ',w1) > 0 then
return 1
if x//2 = 0 then
return 0
if x//3 = 0 then
return 0
if Right(x,1) = 5 then
return 0
/* Miller-Rabin primality test */
numeric digits 2*Length(x)
d = x-1; e = d
/* Reduce n-1 by factors of 2 */
do s = -1 while d//2 = 0
d = d%2
end
/* Thresholds deterministic witnesses */
w2 = ' 2047 1373653 25326001 3215031751 2152302898747 3474749660383 341550071728321 0 ',
|| ' 3825123056546413051 0 0 318665857834031151167461 3317044064679887385961981 '
w = Words(w2)
/* Up to 13 deterministic trials */
if x < Word(w2,w) then do
do k = 1 to w
if x < Word(w2,k) then
leave
end
end
/* or 3 probabilistic trials */
else do
w1 = ' '
do k = 1 to 3
r = Rand(2,e)/1; w1 = w1||r||' '
end
k = k-1
end
/* Algorithm using generated witnesses */
do k = 1 to k
a = Word(w1,k); y = powermod(a,d,x)
if y = 1 then
iterate
if y = e then
iterate
do s
y = (y*y)//x
if y = 1 then
return 0
if y = e then
leave
end
if y <> e then
return 0
end
return 1
Floor:
/* Floor */
procedure expose glob.
arg x
/* Formula */
if Whole(x) then
return x
else
return Trunc(x)-(x<0)
Powermod:
/* Power modulus function x^y mod z */
procedure expose glob.
arg x,y,z
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(x) then
return 'X'
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
if y < 0 then
return 'X'
if z < 1 then
return 'X'
/* Special values */
if z = 1 then
return 0
/* Binary method */
numeric digits Max(Length(Format(x,,,0)),Length(Format(y,,,0)),2*Length(Format(z,,,0)))
b = x//z; r = 1
do while y > 0
if y//2 then
r = r*x//z
x = x*x//z; y = y%2
end
return r
Rand:
/* Random number 12 digits precision */
procedure expose glob.
arg x,y
/* Validity */
if x <> '' then do
if \ Whole(x) then
return 'X'
if \ Whole(x) then
return 'X'
if x >= y then
return 'X'
end
/* Fixed precision */
p = Digits(); numeric digits 30
/* Get and save first seed from system Date and Time */
if glob.rand = '' then do
a = Date('s'); b = Time('l')
glob.rand = Substr(b,10,3)||Substr(b,7,2)||Substr(b,4,2)||Substr(b,1,2)||Substr(a,7,2)||Substr(a,5,2)||Substr(a,3,2)
end
/* Calculate next random number method HP calculators */
glob.rand = Right((glob.rand*2851130928467),15)
/* Uniform deviate [0,1) */
z = '0.'Left(glob.rand,Length(glob.rand)-3)
numeric digits 12
if x = '' then
return z/1
else
return Floor(z*(y+1-x)+x)
Whole:
/* Is a number integer? */
procedure expose glob.
arg x
/* Formula */
return Datatype(x,'w')
</syntaxhighlight>
As in the previous version, k = 3 was hard coded and seems good enough for this test. A higher k gives more confidence on the probable primes, but also a longer running time.
This version will produce output (Windows 11, Intel i7, 4.5Ghz, 16G):
<pre>
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022
25 numbers of the form 2^n-1, mostly Mersenne primes
Up to about 25 digits deterministic, above probabilistic
2^2-1 = 3 (1 digits) is for sure prime (0.000 seconds)
2^3-1 = 7 (1 digits) is for sure prime (0.000 seconds)
2^5-1 = 31 (2 digits) is for sure prime (0.000 seconds)
2^7-1 = 127 (3 digits) is for sure prime (0.000 seconds)
2^11-1 = 2047 (4 digits) is not prime (0.000 seconds)
2^13-1 = 8191 (4 digits) is for sure prime (0.000 seconds)
2^17-1 = 131071 (6 digits) is for sure prime (0.000 seconds)
2^19-1 = 524287 (6 digits) is for sure prime (0.000 seconds)
2^23-1 = 8388607 (7 digits) is not prime (0.000 seconds)
2^31-1 = 2147483647 (10 digits) is for sure prime (0.000 seconds)
2^61-1 = 2305843009213693951 (19 digits) is for sure prime (0.000 seconds)
2^89-1 = 618970019642...0137449562111 (27 digits) is probable prime (0.016 seconds)
2^97-1 = 158456325028...5187087900671 (30 digits) is not prime (0.000 seconds)
2^107-1 = 162259276829...1578010288127 (33 digits) is probable prime (0.000 seconds)
2^113-1 = 103845937170...0992658440191 (35 digits) is not prime (0.000 seconds)
2^127-1 = 170141183460...3715884105727 (39 digits) is probable prime (0.016 seconds)
2^131-1 = 272225893536...9454145691647 (40 digits) is not prime (0.000 seconds)
2^521-1 = 686479766013...8291115057151 (157 digits) is probable prime (0.453 seconds)
2^607-1 = 531137992816...3219031728127 (183 digits) is probable prime (0.765 seconds)
2^1279-1 = 104079321946...5703168729087 (386 digits) is probable prime (9.407 seconds)
2^2203-1 = 147597991521...7686697771007 (664 digits) is probable prime (36.527 seconds)
2^2281-1 = 446087557183...2418132836351 (687 digits) is probable prime (37.575 seconds)
2^2293-1 = 182717463422...4672097697791 (691 digits) is not prime (16.324 seconds)
2^3217-1 = 259117086013...7362909315071 (969 digits) is probable prime (111.373 seconds)
2^3221-1 = 414587337621...7806549041151 (970 digits) is not prime (36.390 seconds)
</pre>
Above 1000 digits it becomes very slow.
=={{header|Ring}}==
Line 6,397 ⟶ 7,052:
I've therefore used this method to check the same numbers as in my Kotlin entry.
<syntaxhighlight lang="
var iters = 10
| |||