Miller–Rabin primality test: Difference between revisions
m
→{{header|Wren}}: Minor tidy
Alextretyak (talk | contribs) (Added 11l) |
m (→{{header|Wren}}: Minor tidy) |
||
| (18 intermediate revisions by 8 users not shown) | |||
Line 15:
''x'' ← ''a''<sup>''d''</sup> mod ''n''
'''if''' ''x'' = 1 or ''x'' = ''n'' − 1 '''then''' '''do''' '''next''' LOOP
'''
''x'' ← ''x''<sup>2</sup> mod ''n''
'''if''' ''x'' = 1 '''then''' '''return''' ''composite''
Line 29:
{{trans|D}}
<
I n < 2 | n % 2 == 0
R n == 2
Line 57:
R 1B
print((2..29).filter(x -> isProbablePrime(x)))</
{{out}}
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 77 ⟶ 390:
such as for the [[Carmichael 3 strong pseudoprimes]] the [[Extensible prime generator]], and the [[Emirp primes]].
<
type Number is range <>;
package Miller_Rabin is
Line 85 ⟶ 398:
function Is_Prime (N : Number; K : Positive := 10) return Result_Type;
end Miller_Rabin;</
The implementation of that package is as follows:
<
package body Miller_Rabin is
Line 143 ⟶ 456:
end Is_Prime;
end Miller_Rabin;</
Finally, the program itself:
<
procedure Mr_Tst is
Line 171 ⟶ 484:
Ada.Text_IO.Put ("Enter the count of loops: "); Pos_IO.Get (K);
Ada.Text_IO.Put_Line ("What is it? " & Result_Type'Image (Is_Prime(N, K)));
end MR_Tst;</
{{out}}
Line 183 ⟶ 496:
Using the big integer implementation from a cryptographic library [https://github.com/cforler/Ada-Crypto-Library/].
<
procedure Miller_Rabin is
Line 268 ⟶ 581:
Ada.Text_IO.Put_Line("Prime(" & S & ")=" & Boolean'Image(Is_Prime(+S, K)));
Ada.Text_IO.Put_Line("Prime(" & T & ")=" & Boolean'Image(Is_Prime(+T, K)));
end Miller_Rabin;</
{{out}}
Line 276 ⟶ 589:
Using the built-in Miller-Rabin test from the same library:
<
procedure Miller_Rabin is
Line 301 ⟶ 614:
Ada.Text_IO.Put_Line("Prime(" & T & ")="
& Boolean'Image (LN.Mod_Utils.Passed_Miller_Rabin_Test(+T, K)));
end Miller_Rabin;</
The output is the same.
Line 310 ⟶ 623:
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}}
<!-- {{does not work with|ELLA ALGOL 68|Any (with appropriate job cards AND formatted transput statements removed) - tested with release 1.8.8d.fc9.i386 - ELLA has no FORMATted transput, also it generates a call to undefined C LONG externals }} -->
<
MODE LOOPINT = INT;
Line 347 ⟶ 660:
print((" ",whole(i,0)))
FI
OD</
{{out}}
<pre>
937 941 947 953 967 971 977 983 991 997
</pre>
=={{header|ARM Assembly}}==
{{works with|as|Raspberry Pi <br> or android 32 bits with application Termux}}
<syntaxhighlight lang ARM Assembly>
/* ARM assembly Raspberry PI */
/* program testmiller.s */
/* for constantes see task include a file in arm assembly */
/************************************/
/* Constantes */
/************************************/
.include "../constantes.inc"
.equ NBDIVISORS, 2000
//.include "../../ficmacros32.inc" @ use for developper debugging
/*******************************************/
/* Initialized data */
/*******************************************/
.data
szMessStartPgm: .asciz "Program 32 bits start \n"
szMessEndPgm: .asciz "Program normal end.\n"
szMessErrorArea: .asciz "\033[31mError : area divisors too small.\n"
szMessPrime: .asciz " is prime !!!.\n"
szMessNotPrime: .asciz " is not prime !!!.\n"
szCarriageReturn: .asciz "\n"
.align 4
iGraine: .int 123456
/*******************************************/
/* UnInitialized data */
/*******************************************/
.bss
.align 4
sZoneConv: .skip 24
/*******************************************/
/* code section */
/*******************************************/
.text
.global main
main: @ program start
ldr r0,iAdrszMessStartPgm @ display start message
bl affichageMess
ldr r4,iStart @ start number
ldr r5,iLimit @ end number
tst r4,#1
addeq r4,#1 @ start with odd number
1:
mov r0,r4
ldr r1,iAdrsZoneConv
bl conversion10 @ decimal conversion
ldr r0,iAdrsZoneConv
bl affichageMess
mov r0,r4
bl isPrimeMiller @ test miller rabin
cmp r0,#0
beq 2f
ldr r0,iAdrszMessPrime
bl affichageMess
b 3f
2:
ldr r0,iAdrszMessNotPrime
bl affichageMess
3:
add r4,r4,#2
cmp r4,r5
ble 1b
ldr r0,iAdrszMessEndPgm @ display end message
bl affichageMess
100: @ standard end of the program
mov r0, #0 @ return code
mov r7, #EXIT @ request to exit program
svc 0 @ perform system call
iAdrszMessStartPgm: .int szMessStartPgm
iAdrszMessEndPgm: .int szMessEndPgm
iAdrszCarriageReturn: .int szCarriageReturn
iAdrsZoneConv: .int sZoneConv
iAdrszMessPrime: .int szMessPrime
iAdrszMessNotPrime: .int szMessNotPrime
iStart: .int 4294967270
iLimit: .int 4294967295
/***************************************************/
/* check if a number is prime test miller rabin */
/***************************************************/
/* r0 contains the number */
/* r0 return 1 if prime 0 else */
@2147483647
@4294967297
@131071
isPrimeMiller:
push {r1-r6,lr} @ save registers
cmp r0,#1 @ control 0 or 1
movls r0,#0
bls 100f
cmp r0,#2 @ control = 2
moveq r0,#1
beq 100f
tst r0,#1
moveq r0,#0 @ even
beq 100f
mov r1,#5 @ loop number
bl testMiller
100:
pop {r1-r6,pc}
/***************************************************/
/* test miller rabin algorithme wikipedia */
/* unsigned */
/***************************************************/
/* r0 contains number */
/* r1 contains parameter */
/* r0 return 1 if prime 0 if composite */
testMiller:
push {r1-r9,lr} @ save registers
mov r4,r0 @ N
mov r7,r1 @ loop maxi
sub r3,r0,#1 @ D
mov r2,#2
mov r6,#0 @ S
1: @ compute D * 2 power S
lsr r3,#1 @ D= D/2
add r6,r6,#1 @ increment S
tst r3,#1 @ D even ?
beq 1b
2:
mov r8,#0 @ loop counter
sub r5,r0,#3
3:
mov r0,r5
bl genereraleas
add r0,r0,#2 @ alea (entre 2 et n-1)
mov r1,r3 @ exposant = D
mov r2,r4 @ modulo N
bl moduloPuR32
cmp r0,#1
beq 5f
sub r1,r4,#1 @ n -1
cmp r0,r1
beq 5f
sub r9,r6,#1 @ S - 1
4:
mov r2,r0
umull r0,r1,r2,r0 @ compute square
mov r2,r4 @ and compute modulo N
bl division32R2023
mov r0,r2
cmp r0,#1
moveq r0,#0 @ composite
beq 100f
sub r1,r4,#1 @ n -1
cmp r0,r1
beq 5f
subs r9,r9,#1
bge 4b
mov r0,#0 @ composite
b 100f
5:
add r8,r8,#1
cmp r8,r7
blt 3b
mov r0,#1 @ prime
100:
pop {r1-r9,pc}
/********************************************************/
/* Calcul modulo de b puissance e modulo m */
/* Exemple 4 puissance 13 modulo 497 = 445 */
/* */
/********************************************************/
/* r0 nombre */
/* r1 exposant */
/* r2 modulo */
/* r0 return result */
moduloPuR32:
push {r1-r6,lr} @ save registers
cmp r0,#0 @ control <> zero
beq 100f
cmp r2,#0 @ control <> zero
beq 100f
1:
mov r4,r2 @ save modulo
mov r5,r1 @ save exposant
mov r6,r0 @ save base
mov r3,#1 @ start result
mov r1,#0 @ division r0,r1 by r2
bl division32R2023
mov r6,r2 @ base <- remainder
2:
tst r5,#1 @ exposant even or odd
beq 3f
umull r0,r1,r6,r3 @ multiplication base
mov r2,r4 @ and compute modulo N
bl division32R2023
mov r3,r2 @ result <- remainder
3:
umull r0,r1,r6,r6 @ compute base square
mov r2,r4 @ and compute modulo N
bl division32R2023
mov r6,r2 @ base <- remainder
lsr r5,#1 @ left shift 1 bit
cmp r5,#0 @ end ?
bne 2b
mov r0,r3
100:
pop {r1-r6,pc}
/***************************************************/
/* division number 64 bits in 2 registers by number 32 bits */
/* unsigned */
/***************************************************/
/* r0 contains lower part dividende */
/* r1 contains upper part dividende */
/* r2 contains divisor */
/* r0 return lower part quotient */
/* r1 return upper part quotient */
/* r2 return remainder */
division32R2023:
push {r3-r6,lr} @ save registers
mov r4,r2 @ save divisor
mov r5,#0 @ init upper part divisor
mov r2,r0 @ save dividende
mov r3,r1
mov r0,#0 @ init result
mov r1,#0
mov r6,#0 @ init shift counter
1: @ loop shift divisor
cmp r5,#0 @ upper divisor <0
blt 2f
cmp r5,r3
cmpeq r4,r2
bhs 2f @ new divisor > dividende
lsl r5,#1 @ shift left one bit upper divisor
lsls r4,#1 @ shift left one bit lower divisor
orrcs r5,r5,#1 @ move bit 31 lower on upper
add r6,r6,#1 @ increment shift counter
b 1b
2: @ loop 2
lsl r1,#1 @ shift left one bit upper quotient
lsls r0,#1 @ shift left one bit lower quotient
orrcs r1,#1 @ move bit 31 lower on upper
cmp r5,r3 @ compare divisor and dividende
cmpeq r4,r2
bhi 3f
subs r2,r2,r4 @ < sub divisor from dividende lower
sbc r3,r3,r5 @ and upper
orr r0,r0,#1 @ move 1 on quotient
3:
lsr r4,r4,#1 @ shift right one bit upper divisor
lsrs r5,#1 @ and lower
orrcs r4,#0x80000000 @ move bit 0 upper to 31 bit lower
subs r6,#1 @ decrement shift counter
bge 2b @ if > 0 loop 2
100:
pop {r3-r6,pc}
/***************************************************/
/* Generation random number */
/***************************************************/
/* r0 contains limit */
genereraleas:
push {r1-r4,lr} @ save registers
ldr r4,iAdriGraine
ldr r2,[r4]
ldr r3,iNbDep1
mul r2,r3,r2
ldr r3,iNbDep1
add r2,r2,r3
str r2,[r4] @ save seed for next call
cmp r0,#0
beq 100f
mov r1,r0 @ divisor
mov r0,r2 @ dividende
bl division
mov r0,r3 @ résult = remainder
100: @ end function
pop {r1-r4,pc} @ restaur registers
iAdriGraine: .int iGraine
iNbDep1: .int 0x343FD
iNbDep2: .int 0x269EC3
/***************************************************/
/* ROUTINES INCLUDE */
/***************************************************/
.include "../affichage.inc"
</syntaxhighlight>
{{Out}}
<pre>
Program 32 bits start
4294967271 is not prime !!!.
4294967273 is not prime !!!.
4294967275 is not prime !!!.
4294967277 is not prime !!!.
4294967279 is prime !!!.
4294967281 is not prime !!!.
4294967283 is not prime !!!.
4294967285 is not prime !!!.
4294967287 is not prime !!!.
4294967289 is not prime !!!.
4294967291 is prime !!!.
4294967293 is not prime !!!.
4294967295 is not prime !!!.
Program normal end.
</pre>
=={{header|AutoHotkey}}==
ahk forum: [http://www.autohotkey.com/forum/post-276712.html#276712 discussion]
<
MsgBox % MillerRabin(999809,10) ; 1
MsgBox % MillerRabin(999727,10) ; 1
Line 395 ⟶ 1,016:
y := i&1 ? mod(y*z,m) : y, z := mod(z*z,m), i >>= 1
Return y
}</
=={{header|bc}}==
Line 401 ⟶ 1,022:
{{works with|OpenBSD bc}}
(A previous version worked with [[GNU bc]].)
<
scale = 0
Line 468 ⟶ 1,089:
}
}
quit</
=={{header|BQN}}==
The function <code>IsPrime</code> in bqn-libs [https://github.com/mlochbaum/bqn-libs/blob/master/primes.bqn primes.bqn] uses deterministic Miller-Rabin to test primality when trial division fails. The following function, derived from that library, selects witnesses at random. The left argument is the number of witnesses to test, with default 10.
<syntaxhighlight lang="bqn">_modMul ← { n _𝕣: n|× }
MillerRabin ← { 𝕊n: 10𝕊n ; iter 𝕊 n: !2|n
# n = 1 + d×2⋆s
s ← 0 {𝕨 2⊸|◶⟨+⟜1𝕊2⌊∘÷˜⊢,⊣⟩ 𝕩} n-1
d ← (n-1) ÷ 2⋆s
# Arithmetic mod n
Mul ← n _modMul
Pow ← Mul{𝔽´𝔽˜⍟(/2|⌊∘÷⟜2⍟(↕1+·⌊2⋆⁼⊢)𝕩)𝕨}
# Miller-Rabin test
MR ← {
1 =𝕩 ? 𝕨≠s ;
(n-1)=𝕩 ? 0 ;
𝕨≤1 ? 1 ;
(𝕨-1) 𝕊 Mul˜𝕩
}
C ← { 𝕊a: s MR a Pow d } # Is composite
{0:1; C •rand.Range⌾(-⟜2) n ? 0; 𝕊𝕩-1} iter
}</syntaxhighlight>
The simple definition of <code>_modMul</code> is inaccurate when intermediate results fall outside the exact integer range (this can happen for inputs around <code>2⋆26</code>). When replaced with the definition below, <code>MillerRabin</code> remains accurate for all inputs, as floating point can't represent odd numbers outside of integer range.
<syntaxhighlight lang="bqn"># Compute n|𝕨×𝕩 in high precision
_modMul ← { n _𝕣:
# Split each argument into two 26-bit numbers, with the remaining
# mantissa bit encoded in the sign of the lower-order part.
q←1+2⋆27
Split ← { h←(q×𝕩)(⊣--)𝕩 ⋄ ⟨𝕩-h,h⟩ }
# The product, and an error relative to precise split multiplication.
Mul ← × (⊣ ⋈ -⊸(+´)) ·⥊×⌜○Split
((n×<⟜0)⊸+ -⟜n+⊢)´ n | Mul
}</syntaxhighlight>
{{out}}
<syntaxhighlight lang="bqn"> MillerRabin 15485867
1
MillerRabin¨⊸/ 101+2×↕10
⟨ 101 103 107 109 113 ⟩</syntaxhighlight>
=={{header|Bracmat}}==
{{trans|bc}}
<
& ( rand
=
Line 541 ⟶ 1,207:
& !primes:? [-11 ?last
& out$!last
);</
{{out}}
<pre>937 941 947 953 967 971 977 983 991 997</pre>
Line 548 ⟶ 1,214:
{{libheader|GMP}}
'''miller-rabin.h'''
<
#define _MILLER_RABIN_H
#include <gmp.h>
bool miller_rabin_test(mpz_t n, int j);
#endif</
'''miller-rabin.c'''
{{trans|Fortran}}
For <code>decompose</code> (and header <tt>primedecompose.h</tt>),
see [[Prime decomposition#C|Prime decomposition]].
<
#include <gmp.h>
#include "primedecompose.h"
Line 619 ⟶ 1,285:
gmp_randclear(rs);
return res;
}</
'''Testing'''
<
#include <stdlib.h>
#include <stdbool.h>
Line 649 ⟶ 1,315:
mpz_clear(num);
return EXIT_SUCCESS;
}</
===Deterministic up to 341,550,071,728,321===
<
size_t power(size_t a, size_t n, size_t mod)
{
Line 724 ⟶ 1,390:
return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13);
return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13) && witness(n, s, d, 17);
}</
Inspiration from http://stackoverflow.com/questions/4424374/determining-if-a-number-is-prime
===Other version===
It should be a 64-bit deterministic version of the Miller-Rabin primality test.
<syntaxhighlight lang="c">
typedef unsigned long long int ulong;
ulong mul_mod(ulong a, ulong b, const ulong mod) {
ulong res = 0, c; // return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (c = b) >= mod - b ? c -= mod : 0, b += c);
return res % mod;
}
ulong pow_mod(ulong n, ulong exp, const ulong mod) {
ulong res = 1; // return (n ^ exp) % mod
for (n %= mod; exp; exp & 1 ? res = mul_mod(res, n, mod) : 0, n = mul_mod(n, n, mod), exp >>= 1);
return res;
}
int is_prime(ulong N) {
// Perform a Miller-Rabin test, it should be a deterministic version.
const ulong n_primes = 9, primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
for (ulong i = 0; i < n_primes; ++i)
if (N % primes[i] == 0) return N == primes[i];
if (N < primes[n_primes - 1]) return 0;
int res = 1, s = 0;
ulong t;
for (t = N - 1; ~t & 1; t >>= 1, ++s);
for (ulong i = 0; i < n_primes && res; ++i) {
ulong B = pow_mod(primes[i], t, N);
if (B != 1) {
for (int b = s; b-- && (res = B + 1 != N);)
B = mul_mod(B, B, N);
res = !res;
}
}
return res;
}
int main(void){
return is_prime(8193145868754512737);
}
</syntaxhighlight>
=={{header|C sharp|C#}}==
<
{
public static bool IsPrime(int n, int k)
Line 754 ⟶ 1,463:
return true;
}
}</
[https://stackoverflow.com/questions/7860802/miller-rabin-primality-test] Corrections made 6/21/2017
<br><br>
<
// Based on the Ruby implementation on this page.
public static class BigIntegerExtensions
Line 814 ⟶ 1,523:
return true;
}
}</
=={{header|C++}}==
This test is deterministic for n < 3'474'749'660'383
<syntaxhighlight lang="c++">
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>
std::vector<uint32_t> small_primes{ 2, 3 };
uint64_t add_modulus(const uint64_t& a, const uint64_t& b, const uint64_t& modulus) {
uint64_t am = ( a < modulus ) ? a : a % modulus;
if ( b == 0 ) {
return am;
}
uint64_t bm = ( b < modulus ) ? b : b % modulus;
uint64_t b_from_m = modulus - bm;
if ( am >= b_from_m ) {
return am - b_from_m;
}
return am + bm;
}
uint64_t multiply_modulus(const uint64_t& a, const uint64_t& b, const uint64_t& modulus) {
uint64_t am = ( a < modulus ) ? a : a % modulus;
uint64_t bm = ( b < modulus ) ? b : b % modulus;
if ( bm > am ) {
std::swap(am, bm);
}
uint64_t result;
while ( bm > 0 ) {
if ( ( bm & 1) == 1 ) {
result = add_modulus(result, am, modulus);
}
am = ( am << 1 ) - ( am >= ( modulus - am ) ? modulus : 0 );
bm >>= 1;
}
return result;
}
uint64_t exponentiation_modulus(const uint64_t& base, const uint64_t& exponent, const uint64_t& modulus) {
uint64_t b = base;
uint64_t e = exponent;
uint64_t result = 1;
while ( e > 0 ) {
if ( ( e & 1 ) == 1 ) {
result = multiply_modulus(result, b, modulus);
}
e >>= 1;
b = multiply_modulus(b, b, modulus);
}
return result;
}
bool is_composite(const uint32_t& a, const uint64_t& d, const uint64_t& n, const uint32_t& s) {
if ( exponentiation_modulus(a, d, n) == 1 ) {
return false;
}
for ( uint64_t i = 0; i < s; ++i ) {
if ( exponentiation_modulus(a, pow(2, i) * d, n) == n - 1 ) {
return false;
}
}
return true;
}
bool composite_test(const std::vector<uint32_t>& primes, const uint64_t& d, const uint64_t& n, const uint32_t& s) {
for ( const uint32_t& prime : primes ) {
if ( is_composite(prime, d, n, s) ) {
return true;
}
}
return false;
}
bool is_prime(const uint64_t& n) {
if ( n == 0 || n == 1 ) {
return false;
}
if ( std::find(small_primes.begin(), small_primes.end(), n) != small_primes.end() ) {
return true;
}
if ( std::any_of(small_primes.begin(), small_primes.end(), [n](uint32_t p) { return n % p == 0; }) ) {
return false;
}
uint64_t d = n - 1;
uint32_t s = 0;
while ( ! d % 2 ) {
d >>= 1;
s++;
}
if ( n < 1'373'653 ) {
return composite_test({ 2, 3 }, d, n, s);
}
if ( n < 25'326'001 ) {
return composite_test({ 2, 3, 5 }, d, n, s);
}
if ( n < 118'670'087'467 ) {
if ( n == 3'215'031'751 ) {
return false;
}
return composite_test({ 2, 3, 5, 7 }, d, n, s);
}
if ( n < 2'152'302'898'747 ) {
return composite_test({ 2, 3, 5, 7, 11 }, d, n, s);
}
if ( n < 3'474'749'660'383 ) {
return composite_test({ 2, 3, 5, 7, 11, 13 }, d, n, s);
}
if ( n < 341'550'071'728'321 ) {
return composite_test({ 2, 3, 5, 7, 11, 13, 17 }, d, n, s);
}
const std::vector<uint32_t> test_primes(small_primes.begin(), small_primes.begin() + 16);
return composite_test(test_primes, d, n, s);
}
void create_small_primes() {
for ( uint32_t i = 5; i < 1'000; i += 2 ) {
if ( is_prime(i) ) {
small_primes.emplace_back(i);
}
}
}
int main() {
create_small_primes();
for ( const uint64_t number : { 1'234'567'890'123'456'733, 1'234'567'890'123'456'737 } ) {
std::cout << "is_prime(" << number << ") = " << std::boolalpha << is_prime(number) << std::endl;
}
}
</syntaxhighlight>
{{ out }}
<pre>
is_prime(1234567890123456733) = false
is_prime(1234567890123456737) = true
</pre>
=={{header|Clojure}}==
===Random Approach===
<
(:require [clojure.math.numeric-tower :as math])
(:require [clojure.set :as set]))
Line 892 ⟶ 1,743:
(println "Is Prime?" 743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153
(random-test 743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153))
</syntaxhighlight>
{{Output}}
<pre>
Line 903 ⟶ 1,754:
</pre>
===Deterministic Approach===
<
(:require [clojure.math.numeric-tower :as math]))
Line 985 ⟶ 1,836:
(println "Is Prime?" 743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153
(deterministic-test 743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153))
</syntaxhighlight>
{{Output}}
<pre>
Line 998 ⟶ 1,849:
(deterministic-test 743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153))
</pre>
=={{header|Commodore BASIC}}==
This displays a minimum probability of primality = 1-1/4<sup>k</sup>, as the fraction of "strong liars" approaches 1/4 in the limit.
<syntaxhighlight lang="basic">100 PRINT CHR$(147); CHR$(18); "**** MILLER-RABIN PRIMALITY TEST ****": PRINT
110 INPUT "NUMBER TO TEST"; N$
120 N = VAL(N$): IF N < 2 THEN 110
130 IF 0 = (N AND 1) THEN PRINT "(EVEN)": GOTO 380
140 INPUT "ITERATIONS"; K$
150 K = VAL(K$): IF K < 1 THEN 140
160 PRINT
170 DEF FNMD(M) = M - N * INT(M / N)
180 D = N - 1
190 S = 0
200 D = D / 2: S = S + 1
210 IF 0 = (D AND 1) THEN 200
220 P = 1
230 FOR I = 1 TO K
240 : A = INT(RND(.) * (N - 2)) + 2
250 : X = 1
260 : FOR J = 1 TO D
270 : X = FNMD(X * A)
280 : NEXT J
290 : IF (X = 1) OR (X = (N - 1)) THEN 360
300 : FOR J = 1 TO S - 1
310 : X = FNMD(X * X)
320 : IF X = 1 THEN P = 0: GOTO 370
330 : IF X = N - 1 THEN 360
340 : NEXT J
350 : P = 0: GOTO 370
360 NEXT I
370 P = P * (1 - 1 / (4 * K))
380 IF P THEN PRINT "PROBABLY PRIME ( P >="; P; ")": END
390 PRINT "COMPOSITE."</syntaxhighlight>
{{Out}}
Sample runs.
<pre>**** MILLER-RABIN PRIMALITY TEST ****
NUMBER TO TEST? 1747
ITERATIONS? 2
PROBABLY PRIME ( P >= .875 )
READY.</pre>
<pre>**** MILLER-RABIN PRIMALITY TEST ****
NUMBER TO TEST? 32329
ITERATIONS? 2
COMPOSITE.
READY.</pre>
=={{header|Common Lisp}}==
<
"Return two values R and E such that NUMBER = DIVISOR^E * R,
and R is not divisible by DIVISOR."
Line 1,042 ⟶ 1,947:
thereis (= y (- n 1)))))))
(loop repeat k
always (strong-liar? (random-in-range 2 (- n 2)))))))))</
<pre>
CL-USER> (last (loop for i from 1 to 1000
Line 1,054 ⟶ 1,959:
=== Standard non-deterministic M-R test ===
<
module Primes
Line 1,093 ⟶ 1,998:
puts 341521.prime?(20) # => true
puts 341531.prime? # => false</
=== Deterministic M-R test ===
Line 1,099 ⟶ 2,004:
It is a direct translation of the Ruby version for arbitrary sized integers.
It is deterministic for all integers < 3_317_044_064_679_887_385_961_981.
<
# crystal build miller-rabin.cr -Ddisable_overflow --release
# crystal build -Ddisable_overflow miller-rabin.cr --release
Line 1,242 ⟶ 2,147:
n = "94366396730334173383107353049414959521528815310548187030165936229578960209523421808912459795329035203510284576187160076386643700441216547732914250578934261891510827140267043592007225160798348913639472564715055445201512461359359488795427875530231001298552452230535485049737222714000227878890892901228389026881".to_big_i
print "\n number = #{n} is prime? "; print " in ", tm{ print n.primemr? }, " secs"
puts</
=={{header|D}}==
{{trans|Ruby}}
<
bool isProbablePrime(in ulong n, in uint k=10) /*nothrow*/ @safe /*@nogc*/ {
Line 1,295 ⟶ 2,200:
iota(2, 30).filter!isProbablePrime.writeln;
}</
{{out}}
<pre>[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]</pre>
=={{header|E}}==
<
if (n <=> 2 || n <=> 3) { return true }
if (n <=> 1 || n %% 2 <=> 0) { return false }
Line 1,322 ⟶ 2,227:
}
return true
}</
<
print(i, " ")
}
println()</
=={{header|EchoLisp}}==
EchoLisp natively implement the '''prime?''' function = Miller-Rabin tests for big integers. The definition is as follows :
<
(lib 'bigint)
Line 1,370 ⟶ 2,275:
(prime? (1+ (factorial 100))) ;; native
→ #f
</syntaxhighlight>
=={{header|Elixir}}==
<
defmodule Prime do
use Application
Line 1,423 ⟶ 2,328:
end
end
</syntaxhighlight>
{{out}}
Line 1,440 ⟶ 2,345:
</pre>
The following larger examples all produce true:
<
miller_rabin?( 94366396730334173383107353049414959521528815310548187030165936229578960209523421808912459795329035203510284576187160076386643700441216547732914250578934261891510827140267043592007225160798348913639472564715055445201512461359359488795427875530231001298552452230535485049737222714000227878890892901228389026881, 1000 )
miller_rabin?( 138028649176899647846076023812164793645371887571371559091892986639999096471811910222267538577825033963552683101137782650479906670021895135954212738694784814783986671046107023185842481502719762055887490765764329237651328922972514308635045190654896041748716218441926626988737664133219271115413563418353821396401, 1000 )
Line 1,448 ⟶ 2,353:
miller_rabin?( 153410708946188157980279532372610756837706984448408515364579602515073276538040155990230789600191915021209039203172105094957316552912585741177975853552299222501069267567888742458519569317286299134843250075228359900070009684517875782331709619287588451883575354340318132216817231993558066067063143257425853927599, 1000 )
miller_rabin?( 103130593592068072608023213244858971741946977638988649427937324034014356815504971087381663169829571046157738503075005527471064224791270584831779395959349442093395294980019731027051356344056416276026592333932610954020105156667883269888206386119513058400355612571198438511950152690467372712488391425876725831041, 1000 )
</syntaxhighlight>
=={{header|Erlang}}==
Line 1,455 ⟶ 2,360:
to permit use of integers of arbitrary precision.
<
-module(miller_rabin).
Line 1,545 ⟶ 2,450:
Acc;
power(B, E, Acc) ->
power(B, E - 1, B * Acc).</
The above code optimised as follows:
Line 1,553 ⟶ 2,458:
53s to 11s on a quad-core 17 with 16 GB ram. The performance
gain from the improved exponentiation was not evaluated.
<
%%% @author Tony Wallace <tony@resurrection>
%%% @copyright (C) 2021, Tony Wallace
Line 1,737 ⟶ 2,642:
.
</syntaxhighlight>
=={{header|F_Sharp|F#}}==
<
// Miller primality test for n<3317044064679887385961981. Nigel Galloway: April 1st., 2021
let a=[(2047I,[2I]);(1373653I,[2I;3I]);(9080191I,[31I;73I]);(25326001I,[2I;3I;5I]);(3215031751I,[2I;3I;5I;7I]);(4759123141I,[2I;7I;61I]);(1122004669633I,[2I;13I;23I;1662803I]);
Line 1,752 ⟶ 2,657:
printfn "%A %A" (mrP 2147483647I)(mrP 844674407370955389I)
</syntaxhighlight>
{{out}}
<pre>
Line 1,760 ⟶ 2,665:
Forth only supports native ints (e.g. 64 bits on most modern machines), so this version uses a set of bases that is known to be deterministic for 64 bit integers (and possibly greater). Prior to the Miller Rabin check, the "prime?" word checks for divisibility by some small primes.
<syntaxhighlight lang="forth">
\ modular multiplication and exponentiation
\
Line 1,849 ⟶ 2,754:
then
then ;
</syntaxhighlight>
{{Out}}
Test on some Fermat numbers and some Mersenne numbers
Line 1,864 ⟶ 2,769:
{{works with|Fortran|95}}
For the module ''PrimeDecompose'', see [[Prime decomposition#Fortran|Prime decomposition]].
<
module Miller_Rabin
use PrimeDecompose
Line 1,920 ⟶ 2,825:
end function miller_rabin_test
end module Miller_Rabin</
'''Testing'''
<
use Miller_Rabin
implicit none
Line 1,945 ⟶ 2,850:
end subroutine do_test
end program TestMiller</
''Possible improvements'': create bindings to the [[:Category:GMP|GMP library]], change <code>integer(huge)</code> into something like <code>type(huge_integer)</code>, write a lots of interfaces to allow to use big nums naturally (so that the code will be unchanged, except for the changes said above)
===With some avoidance of overflow===
Integer overflow is a severe risk, and even 64-bit integers won't get you far when the formulae are translated as <code>MOD(A**D,N)</code> - what is needed is a method for raising to a power that incorporates the modulus along the way. There is no library routine for that, so... <
CONTAINS !Working only with in-built integers.
LOGICAL FUNCTION MRPRIME(N,TRIALS) !Could N be a prime number?
Line 2,025 ⟶ 2,930:
END DO
END</
Output:
<pre>
Line 2,143 ⟶ 3,048:
Using the task pseudo code
===Up to 2^63-1===
<
' compile with: fbc -s console
Line 2,243 ⟶ 3,148:
Print : Print "hit any key to end program"
Sleep
End</
{{out}}
<pre>9223372036854774893 9223372036854774917 9223372036854774937
Line 2,257 ⟶ 3,162:
===Using Big Integer library===
{{libheader|GMP}}
<
' compile with: fbc -s console
Line 2,371 ⟶ 3,276:
Print : Print "hit any key to end program"
Sleep
End</
{{out}}
<pre> 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
Line 2,397 ⟶ 3,302:
Direct implementation of the task algorithm.
<
def isProbablyPrimeMillerRabin( n, k ) =
Line 2,426 ⟶ 3,331:
for i <- 3..100
if isProbablyPrimeMillerRabin( i, 5 )
println( i )</
{{out}}
Line 2,464 ⟶ 3,369:
The main difference between this algorithm and the pseudocode in the task description is that k numbers are not chosen randomly, but instead are the three numbers 2, 7, and 61. These numbers provide a deterministic primality test up to 2^32.
<
import "log"
Line 2,529 ⟶ 3,434:
}
return true
}</
=={{header|Haskell}}==
Line 2,538 ⟶ 3,443:
* Original Rosetta code has been simplified to be easier to follow
Another Miller Rabin test can be found in D. Amos's Haskell for Math module [http://www.polyomino.f2s.com/david/haskell/numbertheory.html Primes.hs]
<
import System.Random
Line 2,582 ⟶ 3,487:
g i b y | even i = g (i `quot` 2) (b*b `rem` m) y
| otherwise = f (i-1) b (b*y `rem` m)
</syntaxhighlight>
{{out|Sample output}}
Line 2,599 ⟶ 3,504:
* The code above likely has better complexity.
<syntaxhighlight lang="haskell">
import Control.Monad (liftM)
import Data.Bits (Bits, testBit, shiftR)
Line 2,655 ⟶ 3,560:
[n,k] <- liftM (map (\x -> read x :: Integer) . words) getLine
print $ isPrime n k
</syntaxhighlight>
Line 2,676 ⟶ 3,581:
The following runs in both languages:
<
every n := !A do write(n," is ",(mrp(n,5),"probably prime")|"composite")
end
Line 2,708 ⟶ 3,613:
}
return [s,d]
end</
Sample run:
Line 2,728 ⟶ 3,633:
=={{header|Java}}==
The Miller-Rabin primality test is part of the standard library (java.math.BigInteger)
<
public class MillerRabinPrimalityTest {
Line 2,736 ⟶ 3,641:
System.out.println(n.toString() + " is " + (n.isProbablePrime(certainty) ? "probably prime" : "composite"));
}
}</
{{out|Sample output}}
<pre>java MillerRabinPrimalityTest 123456791234567891234567 1000000
Line 2,742 ⟶ 3,647:
This is a translation of the [http://rosettacode.org/wiki/Miller-Rabin_primality_test#Python:_Proved_correct_up_to_large_N Python solution] for a deterministic test for n < 341,550,071,728,321:
<
public class Prime {
Line 2,821 ⟶ 3,726:
}
}
</syntaxhighlight>
=={{header|JavaScript}}==
For the return values of this function, <code>true</code> means "probably prime" and <code>false</code> means "definitely composite."
<
if (n === 2 || n === 3) return true
if (n % 2 === 0 || n < 2) return false
Line 2,849 ⟶ 3,754:
}
return false
}</
=={{header|Julia}}==
The built-in <code>isprime</code> function uses the Miller-Rabin primality test. Here is the implementation of <code>isprime</code> from the Julia standard library (Julia version 0.2):
<
witnesses(n::Union(Uint8,Int8,Uint16,Int16)) = (2,3)
witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? (2,3) : (2,7,61)
Line 2,881 ⟶ 3,786:
return true
end
</syntaxhighlight>
=={{header|Kotlin}}==
Translating the pseudo-code directly rather than using the Java library method BigInteger.isProbablePrime(certainty):
<
import java.math.BigInteger
Line 2,937 ⟶ 3,842:
for (bi in bia)
println("$bi is ${if (isProbablyPrime(bi, k)) "probably prime" else "composite"}")
}</
{{out}}
Line 2,949 ⟶ 3,854:
=={{header|Liberty BASIC}}==
<syntaxhighlight lang="lb">
DIM mersenne(11)
mersenne(1)=7
Line 3,240 ⟶ 4,145:
End Function
</syntaxhighlight>
=={{header|Lua}}==
This implementation of the Miller-Rabin probabilistic primality test is based on the treatment in Chapter 10 of "A Computational Introduction to Number Theory and Algebra" by Victor Shoup.
<syntaxhighlight lang="lua"> function MRIsPrime(n, k)
-- If n is prime, returns true (without fail).
-- If n is not prime, then returns false with probability ≥ 4^(-k), true otherwise.
-- First, deal with 1 and multiples of 2.
if n == 1 then
return false
elseif n == 2 then
return true
elseif n%2 == 0 then
return false
end
-- Now n is odd and greater than 1.
-- Find the unique expression n = 1 + t*2^h for n, where t is odd and h≥1.
t = (n-1)/2
h = 1
while t%2 == 0 do
t = t/2
h = h + 1
end
for i = 1, k do
-- Generate a random integer between 1 and n-1 inclusive.
a = math.random(n-1)
-- Test whether a is an element of the set L, and return false if not.
if not IsInL(n, a, t, h) then
return false
end
end
-- All generated a were in the set L; thus with high probability, n is prime.
return true
end
function IsInL(n, a, t, h)
local b = PowerMod(a, t, n)
if b == 1 then
return true
end
for j = 0, h-1 do
if b == n-1 then
return true
elseif b == 1 then
return false
end
b = (b^2)%n
end
return false
end
function PowerMod(x, y, m)
-- Computes x^y mod m.
local z = 1
while y > 0 do
if y%2 == 0 then
x, y, z = (x^2)%m, y//2, z
else
x, y, z = (x^2)%m, y//2, (x*z)%m
end
end
return z
end </syntaxhighlight>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
Do[
a=RandomInteger[{2,n-1}]; x=PowerMod[a,d,n];
Line 3,251 ⟶ 4,220:
];
,{k}];
Print[test] ]</
{{out|Example output (not using the PrimeQ builtin)}}
<
->False</
=={{header|Maxima}}==
<
Line 3,304 ⟶ 4,273:
)
)
)$</
=={{header|Mercury}}==
Line 3,325 ⟶ 4,294:
found with instructions for use in Github.
<syntaxhighlight lang="mercury">
%----------------------------------------------------------------------%
:- module primality.
Line 3,545 ⟶ 4,514:
:- end_module test_is_prime.
</syntaxhighlight>
=={{header|Nim}}==
Line 3,551 ⟶ 4,520:
===Deterministic approach limited to uint32 values.===
<
## Nim currently doesn't have a BigInt standard library
## so we translate the version from Go which uses a
Line 3,625 ⟶ 4,594:
assert isPrime(492366587u32)
assert isPrime(1645333507u32)
</syntaxhighlight>
=== Correct M-R test implementation for using bases > input, deterministic for all integers < 2^64.===
<
# Compile as: $ nim c -d:release mrtest.nim
Line 3,771 ⟶ 4,740:
echo("\nnumber of primes < ",num, " are ", primes.len)
echo (epochTime()-te).formatFloat(ffDecimal, 6)
</syntaxhighlight>
=={{Header|OCaml}}==
Line 3,777 ⟶ 4,746:
A direct translation of the wikipedia pseudocode (with <tt>get_rd</tt> helper function translated from <tt>split</tt> in the scheme implementation). This code uses the Zarith and Bigint (bignum) libraries.
<syntaxhighlight lang="ocaml">
(* Translated from the wikipedia pseudo-code *)
let miller_rabin n ~iter:k =
Line 3,819 ⟶ 4,788:
in
loop 0 true
</syntaxhighlight>
=={{header|Oz}}==
Line 3,826 ⟶ 4,795:
the Mercury and Prolog versions on this page.
<syntaxhighlight lang="oz">
%--------------------------------------------------------------------------%
% module: Primality
Line 3,932 ⟶ 4,901:
% end_module Primality
</syntaxhighlight>
=={{header|PARI/GP}}==
===Built-in===
<
===Custom===
<
my(s = valuation(n-1, 2), d = Mod(b, n)^(n >> s));
if (d == 1, return(1));
Line 3,953 ⟶ 4,922:
);
1
};</
===Deterministic version===
A basic deterministic test can be obtained by an appeal to the ERH (as proposed by Gary Miller) and a result of Eric Bach (improving on Joseph Oesterlé). Calculations of Jan Feitsma can be used to speed calculations below 2<sup>64</sup> (by a factor of about 250).
<
Miller(n)={
if (n%2 == 0, return(n == 2)); \\ Handle even numbers
Line 3,975 ⟶ 4,944:
1
)
};</
=={{header|Perl}}==
===Custom===
<
sub is_prime {
Line 4,011 ⟶ 4,980:
}
print join ", ", grep { is_prime $_, 10 } (1 .. 1000);</
===Modules===
{{libheader|ntheory}}
While normally one would use <tt>is_prob_prime</tt>, <tt>is_prime</tt>, or <tt>is_provable_prime</tt>, which will do a [[wp:Baillie--PSW_primality_test|BPSW test]] and possibly more, we can use just the Miller-Rabin test if desired. For large values we can use an object (e.g. bigint, Math::GMP, Math::Pari, etc.) or just a numeric string.
<
sub is_prime_mr {
my $n = shift;
Line 4,025 ⟶ 4,994:
# Otherwise, perform a number of random base tests, and the result is a probable prime test.
return miller_rabin_random($n, 20);
}</
Math::Primality also has this functionality, though its function takes only one base and requires the input number to be less than the base.
<
sub is_prime_mr {
my $n = shift;
Line 4,036 ⟶ 5,005:
1;
}
for (1..100) { say if is_prime_mr($_) }</
Math::Pari can be used in a fashion similar to the Pari/GP custom function. The builtin accessed using a second argument to <tt>ispseudoprime</tt> was added to a later version of Pari (the Perl module uses version 2.1.7) so is not accessible directly from Perl.
Line 4,044 ⟶ 5,013:
Native-types deterministic version, fails (false negative) at 94,910,107 on 32-bit [fully tested, ie from 1],
and at 4,295,041,217 on 64-bit [only tested from 4,290,000,000] - those limits have now been hard-coded below.
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">powermod</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
Line 4,112 ⟶ 5,081:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%d is %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],{</span><span style="color: #008000;">"composite"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"prime"</span><span style="color: #0000FF;">}[</span><span style="color: #000000;">is_prime_mr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</
{{out}}
<pre>
Line 4,129 ⟶ 5,098:
{{trans|Ruby}}
While desktop/Phix uses a thin wrapper to the builtin gmp routine, the following is also available and is used (after transpilation) in mpfr.js, renamed as mpz_prime:
<!--<
<span style="color: #000080;font-style:italic;">-- this is transpiled (then manually copied) to mpz_prime() in mpfr.js:</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">modp47</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">NULL</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">w</span>
Line 4,232 ⟶ 5,201:
<span style="color: #008080;">return</span> <span style="color: #004600;">true</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<!--</
Either the standard shim or the above can be used as follows
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
Line 4,263 ⟶ 5,232:
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</
{{out}}
<pre>
Line 4,295 ⟶ 5,264:
=={{header|PHP}}==
<
function is_prime($n, $k) {
if ($n == 2)
Line 4,333 ⟶ 5,302:
echo "$i, ";
echo "\n";
?></
=={{header|PicoLisp}}==
<
(use (R D)
(while (=0 (setq R (abs (rand)))))
Line 4,382 ⟶ 5,351:
(do K
(NIL (_prim? N D S))
T ) ) ) )</
{{out}}
<pre>: (filter '((I) (prime? I)) (range 937 1000))
Line 4,394 ⟶ 5,363:
=={{header|Pike}}==
<syntaxhighlight lang="pike">
Line 4,486 ⟶ 5,455:
36261430139487433507414165833468680972181038593593271409697364115931523786727274410257181186996611100786935727 PRIME
15579763548573297857414066649875054392128789371879472432457450095645164702139048181789700140949438093329334293 PRIME
</syntaxhighlight>
=={{header|Prolog}}==
Line 4,493 ⟶ 5,462:
from the Erlang version on this page.
<
% is_prime/2 returns false if N is composite, true if N probably prime
Line 4,575 ⟶ 5,544:
; Next_Loop =:= S -> Result = false
; inner_loop(Next_Base, N, Next_Loop, S, Result)
).</
=={{header|PureBasic}}==
<
#Composite
#Probably_prime
Line 4,606 ⟶ 5,575:
Wend
ProcedureReturn #Probably_prime
EndProcedure</
=={{header|Python}}==
Line 4,614 ⟶ 5,583:
This versions will give answers with a very small probability of being false. That probability being dependent number of trials (automatically set to 8).
<
import random
Line 4,654 ⟶ 5,623:
return False
return True </
===Python: Proved correct up to large N===
Line 4,660 ⟶ 5,629:
<br>For 341550071728321 and beyond, I have followed the pattern in choosing <code>a</code> from the set of prime numbers.<br>While this uses the best sets known in 1993, there are [http://miller-rabin.appspot.com/ better sets known], and at most 7 are needed for 64-bit numbers.
<
if pow(a, d, n) == 1:
return False
Line 4,696 ⟶ 5,665:
_known_primes = [2, 3]
_known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]</
;Testing:
Line 4,711 ⟶ 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}}==
<
(define (miller-rabin-expmod base exp m)
(define (squaremod-with-check x)
Line 4,745 ⟶ 5,838:
(prime? 4547337172376300111955330758342147474062293202868155909489) ;-> outputs true
</syntaxhighlight>
=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2015-09-22}}
<syntaxhighlight lang="raku"
sub expmod(Int $a is copy, Int $b is copy, $n) {
my $c = 1;
Line 4,791 ⟶ 5,884:
}
say (1..1000).grep({ is_prime($_, 10) }).join(", "); </
=={{header|REXX}}==
Line 4,804 ⟶ 5,897:
<br>To make the program small, the ''true prime generator'' ('''GenP''') was coded to be small, but not particularly fast.
<
parse arg limit times seed . /*obtain optional arguments from the CL*/
if limit=='' | limit=="," then limit= 1000 /*Not specified? Then use the default.*/
Line 4,851 ⟶ 5,944:
if x\==nL then return 0 /*nope, it ain't prime nohows, noway. */
end /*k*/ /*maybe it's prime, maybe it ain't ··· */
return 1 /*coulda/woulda/shoulda be prime; yup.*/</
{{out|output|text= when using the input of: <tt> 10000 10 </tt>}}
<pre>
Line 4,878 ⟶ 5,971:
=={{header|Ring}}==
<
# Project : Miller–Rabin primality test
Line 4,927 ⟶ 6,020:
ok
end
</syntaxhighlight>
Output:
<pre>
Line 4,938 ⟶ 6,031:
===Standard Probabilistic===
From 2.5 Ruby has fast modular exponentiation built in. For alternatives prior to 2.5 please see [[Modular_exponentiation#Ruby]]
<
d = n - 1
s = 0
Line 4,959 ⟶ 6,052:
end
p primes = (3..1000).step(2).find_all {|i| miller_rabin_prime?(i,10)}</
{{out}}
<pre>[3, 5, 7, 11, 13, 17, ..., 971, 977, 983, 991, 997]</pre>
The following larger examples all produce true:
<
puts miller_rabin_prime?(138028649176899647846076023812164793645371887571371559091892986639999096471811910222267538577825033963552683101137782650479906670021895135954212738694784814783986671046107023185842481502719762055887490765764329237651328922972514308635045190654896041748716218441926626988737664133219271115413563418353821396401,1000)
puts miller_rabin_prime?(123301261697053560451930527879636974557474268923771832437126939266601921428796348203611050423256894847735769138870460373141723679005090549101566289920247264982095246187318303659027201708559916949810035265951104246512008259674244307851578647894027803356820480862664695522389066327012330793517771435385653616841,1000)
Line 4,969 ⟶ 6,062:
puts miller_rabin_prime?(132082885240291678440073580124226578272473600569147812319294626601995619845059779715619475871419551319029519794232989255381829366374647864619189704922722431776563860747714706040922215308646535910589305924065089149684429555813953571007126408164577035854428632242206880193165045777949624510896312005014225526731,1000)
puts miller_rabin_prime?(153410708946188157980279532372610756837706984448408515364579602515073276538040155990230789600191915021209039203172105094957316552912585741177975853552299222501069267567888742458519569317286299134843250075228359900070009684517875782331709619287588451883575354340318132216817231993558066067063143257425853927599,1000)
puts miller_rabin_prime?(103130593592068072608023213244858971741946977638988649427937324034014356815504971087381663169829571046157738503075005527471064224791270584831779395959349442093395294980019731027051356344056416276026592333932610954020105156667883269888206386119513058400355612571198438511950152690467372712488391425876725831041,1000)</
===Deterministic for integers < 3,317,044,064,679,887,385,961,981===
It extends '''class Integer''' to make it simpler to use.
<
# Returns true if +self+ is a prime number, else returns false.
def primemr?(k = 10)
Line 5,096 ⟶ 6,189:
n = 94366396730334173383107353049414959521528815310548187030165936229578960209523421808912459795329035203510284576187160076386643700441216547732914250578934261891510827140267043592007225160798348913639472564715055445201512461359359488795427875530231001298552452230535485049737222714000227878890892901228389026881
print "\n number = #{n} is prime? "; print " in ", tm{ print n.primemr? }, " secs"
puts</
=={{header|Run BASIC}}==
Line 5,102 ⟶ 6,195:
''This code has not been fully tested. Remove this comment after review.''
<
input "Input test:";k
Line 5,156 ⟶ 6,249:
wend
[funEnd]
END FUNCTION</
=={{header|Rust}}==
<
num = "0.2.0"
rand = "0.6.5"
Line 5,322 ⟶ 6,415:
// that n really is a prime number, so return true:
true
}</
'''Test code:'''
<
let n = 1234687;
let result = is_prime(&n);
Line 5,350 ⟶ 6,443:
let result = is_prime(&n);
println!("Q: Is {} prime? A: {}", n, result);
}</
{{out}}
<pre>Q: Is 1234687 prime? A: true
Line 5,360 ⟶ 6,453:
=={{header|Scala}}==
{{libheader|Scala}}<
object MillerRabinPrimalityTest extends App {
val (n, certainty )= (BigInt(args(0)), args(1).toInt)
println(s"$n is ${if (n.isProbablePrime(certainty)) "probably prime" else "composite"}")
}</
Direct implementation of algorithm:
<
import scala.annotation.tailrec
import scala.language.{implicitConversions, postfixOps}
Line 5,407 ⟶ 6,500:
}) != 1
}
}</
=={{header|Scheme}}==
<
(import (rnrs base (6))
(srfi :27 random-bits))
Line 5,452 ⟶ 6,545:
(and (> n 1)
(or (= n 2)
(pseudoprime? n 50))))</
=={{header|Seed7}}==
<
include "bigint.s7i";
Line 5,500 ⟶ 6,593:
end if;
end for;
end func;</
=={{header|Sidef}}==
<
return false if (n <= 1)
Line 5,530 ⟶ 6,623:
}
say miller_rabin.grep(^1000).join(', ')</
=={{header|Smalltalk}}==
{{works with|GNU Smalltalk}}
Smalltalk handles big numbers naturally and trasparently (the parent class <tt>Integer</tt> has many subclasses, and <cite>a subclass is picked according to the size</cite> of the integer that must be handled)
<
millerRabinTest: kl [ |k| k := kl.
self <= 3
Line 5,568 ⟶ 6,661:
]
]
].</
<
(n millerRabinTest: 10) ifTrue: [ n printNl ]
].</
=={{header|Standard ML}}==
<
val mr_iterations = Int.toLarge 20;
Line 5,622 ⟶ 6,715:
then (n,t)
else findPrime t end
in List.tabulate (10, fn e => findPrime 0) end;</
{{out|Sample run}}
<pre>
Line 5,649 ⟶ 6,742:
{{libheader|Attaswift BigInt}}
<
private let numTrails = 5
Line 5,685 ⟶ 6,778:
return true
}</
=={{header|Tcl}}==
Use Tcl 8.5 for large integer support
<
proc miller_rabin {n k} {
Line 5,722 ⟶ 6,815:
puts $i
}
}</
{{out}}
<pre>1
Line 5,741 ⟶ 6,834:
I've therefore used this method to check the same numbers as in my Kotlin entry.
<
var iters = 10
Line 5,757 ⟶ 6,850:
for (bi in bia) {
System.print("%(bi) is %(bi.isProbablePrime(iters) ? "probably prime" : "composite")")
}</
{{out}}
Line 5,770 ⟶ 6,863:
=={{header|zkl}}==
Using the Miller-Rabin primality test in GMP:
<
zkl: BN("4547337172376300111955330758342147474062293202868155909489").probablyPrime()
True
zkl: BN("4547337172376300111955330758342147474062293202868155909393").probablyPrime()
False</
| |||