Miller–Rabin primality test
You are encouraged to solve this task according to the task description, using any language you may know.
This page uses content from Wikipedia. The original article was at Miller–Rabin primality test. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance) |
The Miller–Rabin primality test or Rabin–Miller primality test is a primality test: an algorithm which determines whether a given number is prime or not. The algorithm, as modified by Michael O. Rabin to avoid the generalized Riemann hypothesis, is a probabilistic algorithm.
The pseudocode, from Wikipedia is:
Input: n > 2, an odd integer to be tested for primality; k, a parameter that determines the accuracy of the test Output: composite if n is composite, otherwise probably prime write n − 1 as 2s·d with d odd by factoring powers of 2 from n − 1 LOOP: repeat k times: pick a randomly in the range [2, n − 2] x ← ad mod n if x = 1 or x = n − 1 then do next LOOP for r = 1 .. s − 1 x ← x2 mod n if x = 1 then return composite if x = n − 1 then do next LOOP return composite return probably prime
- The nature of the test involves big numbers, so the use of "big numbers" libraries (or similar features of the language of your choice) are suggested, but not mandatory.
- Deterministic variants of the test exist and can be implemented as extra (not mandatory to complete the task)
ALGOL 68
<lang algol>MODE LINT=LONG INT; MODE LOOPINT = INT;
MODE POWMODSTRUCT = LINT; PR READ "prelude/pow_mod.a68" PR;
PROC miller rabin = (LINT n, LOOPINT k)BOOL: (
IF n<=3 THEN TRUE ELIF NOT ODD n THEN FALSE ELSE LINT d := n - 1; INT s := 0; WHILE NOT ODD d DO d := d OVER 2; s +:= 1 OD; TO k DO LINT a := 2 + ENTIER (random*(n-3)); LINT x := pow mod(a, d, n); IF x /= 1 THEN TO s DO IF x = n-1 THEN done FI; x := x*x %* n OD; else: IF x /= n-1 THEN return false FI; done: EMPTY FI OD; TRUE EXIT return false: FALSE FI
);
FOR i FROM 937 TO 1000 DO
IF miller rabin(i, 10) THEN print((" ",whole(i,0))) FI
OD</lang> Sample output:
937 941 947 953 967 971 977 983 991 997
AutoHotkey
ahk forum: discussion <lang AutoHotkey> MsgBox % MillerRabin(999983,10) ; 1 MsgBox % MillerRabin(999809,10) ; 1 MsgBox % MillerRabin(999727,10) ; 1 MsgBox % MillerRabin(52633,10) ; 0 MsgBox % MillerRabin(60787,10) ; 0 MsgBox % MillerRabin(999999,10) ; 0 MsgBox % MillerRabin(999995,10) ; 0 MsgBox % MillerRabin(999991,10) ; 0
MillerRabin(n,k) { ; 0: composite, 1: probable prime (n < 2**31)
d := n-1, s := 0 While !(d&1) d>>=1, s++
Loop %k% { Random a, 2, n-2 ; if n < 4,759,123,141, it is enough to test a = 2, 7, and 61. x := PowMod(a,d,n) If (x=1 || x=n-1) Continue Cont := 0 Loop % s-1 { x := PowMod(x,2,n) If (x = 1) Return 0 If (x = n-1) { Cont = 1 Break } } IfEqual Cont,1, Continue Return 0 } Return 1
}
PowMod(x,n,m) { ; x**n mod m
y := 1, i := n, z := x While i>0 y := i&1 ? mod(y*z,m) : y, z := mod(z*z,m), i >>= 1 Return y
}</lang>
Bc
<lang bc>next = 1; # seed of the random generator scale = 0;
- random generator (according to POSIX...); since
- this gives number between 0 and 32767, we can
- test only numbers in this range.
define rand() {
next = next * 1103515245 + 12345; return ((next/65536) % 32768);
}
define rangerand(from, to) {
return (from + rand()%(to-from+1));
}
define miller_rabin_test(n, k)
{
auto d, r, a, x, s;
if ( n <= 3) { return(1); } if ( (n % 2) == 0) { return(0); }
# find s and d so that d*2^s = n-1 d = n-1; s = 0; while( (d%2) == 0 ) { d = d/2; s += 1; }
while(k-- > 0) { a = rangerand(2, n-2); x = (a^d) % n; if ( x != 1 ) { for(r=0; r < s; r++) { if ( x == (n-1) ) { break; } x = (x*x) % n; } if ( x != (n-1) ) { return (0); } } } return (1);
}
for(i=1; i < 1000; i++) {
if ( miller_rabin_test(i, 10) == 1 ) { i }
} quit;</lang>
C
miller-rabin.h
<lang c>#ifndef _MILLER_RABIN_H_
- define _MILLER_RABIN_H
- include <gmp.h>
bool miller_rabin_test(mpz_t n, int j);
- endif</lang>
miller-rabin.c
For decompose
(and header primedecompose.h), see Prime decomposition.
<lang c>#include <stdbool.h>
- include <gmp.h>
- include "primedecompose.h"
- define MAX_DECOMPOSE 100
bool miller_rabin_test(mpz_t n, int j) {
bool res; mpz_t f[MAX_DECOMPOSE]; mpz_t s, d, a, x, r; mpz_t n_1, n_3; gmp_randstate_t rs; int l=0, k;
res = false; gmp_randinit_default(rs);
mpz_init(s); mpz_init(d); mpz_init(a); mpz_init(x); mpz_init(r); mpz_init(n_1); mpz_init(n_3);
if ( mpz_cmp_si(n, 3) <= 0 ) { // let us consider 1, 2, 3 as prime gmp_randclear(rs); return true; } if ( mpz_odd_p(n) != 0 ) { mpz_sub_ui(n_1, n, 1); // n-1 mpz_sub_ui(n_3, n, 3); // n-3 l = decompose(n_1, f); mpz_set_ui(s, 0); mpz_set_ui(d, 1); for(k=0; k < l; k++) { if ( mpz_cmp_ui(f[k], 2) == 0 )
mpz_add_ui(s, s, 1);
else
mpz_mul(d, d, f[k]);
} // 2^s * d = n-1 while(j-- > 0) { mpz_urandomm(a, rs, n_3); // random from 0 to n-4 mpz_add_ui(a, a, 2); // random from 2 to n-2 mpz_powm(x, a, d, n); if ( mpz_cmp_ui(x, 1) == 0 ) continue; mpz_set_ui(r, 0); while( mpz_cmp(r, s) < 0 ) {
if ( mpz_cmp(x, n_1) == 0 ) break; mpz_powm_ui(x, x, 2, n); mpz_add_ui(r, r, 1);
} if ( mpz_cmp(x, n_1) == 0 ) continue; goto flush; // woops } res = true; }
flush:
for(k=0; k < l; k++) mpz_clear(f[k]); mpz_clear(s); mpz_clear(d); mpz_clear(a); mpz_clear(x); mpz_clear(r); mpz_clear(n_1); mpz_clear(n_3); gmp_randclear(rs); return res;
}</lang>
Testing
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <stdbool.h>
- include <gmp.h>
- include "miller-rabin.h"
- define PREC 10
- define TOP 4000
int main() {
mpz_t num;
mpz_init(num); mpz_set_ui(num, 1); while ( mpz_cmp_ui(num, TOP) < 0 ) { if ( miller_rabin_test(num, PREC) ) { gmp_printf("%Zd maybe prime\n", num); } /*else { gmp_printf("%Zd not prime\n", num); }*/ // remove the comment iff you're interested in // sure non-prime. mpz_add_ui(num, num, 1); }
mpz_clear(num); return EXIT_SUCCESS;
}</lang>
E
<lang e>def millerRabinPrimalityTest(n :(int > 0), k :int, random) :boolean {
if (n <=> 2 || n <=> 3) { return true } if (n <=> 1 || n %% 2 <=> 0) { return false } var d := n - 1 var s := 0 while (d %% 2 <=> 0) { d //= 2 s += 1 } for _ in 1..k { def nextTrial := __continue def a := random.nextInt(n - 3) + 2 # [2, n - 2] = [0, n - 4] + 2 = [0, n - 3) + 2 var x := a**d %% n # Note: Will do optimized modular exponentiation if (x <=> 1 || x <=> n - 1) { nextTrial() } for _ in 1 .. (s - 1) { x := x**2 %% n if (x <=> 1) { return false } if (x <=> n - 1) { nextTrial() } } return false } return true
}</lang>
<lang e>for i ? (millerRabinPrimalityTest(i, 1, entropy)) in 4..1000 {
print(i, " ")
} println()</lang>
Fortran
For the module PrimeDecompose, see Prime decomposition.
<lang fortran>module Miller_Rabin
use PrimeDecompose implicit none
integer, parameter :: max_decompose = 100
private :: int_rrand, max_decompose
contains
function int_rrand(from, to) integer(huge) :: int_rrand integer(huge), intent(in) :: from, to
real :: o call random_number(o) int_rrand = floor(from + o * real(max(from,to) - min(from, to))) end function int_rrand
function miller_rabin_test(n, k) result(res) logical :: res integer(huge), intent(in) :: n integer, intent(in) :: k integer(huge), dimension(max_decompose) :: f integer(huge) :: s, d, i, a, x, r
res = .true. f = 0
if ( (n <= 2) .and. (n > 0) ) return if ( mod(n, 2) == 0 ) then res = .false. return end if
call find_factors(n-1, f) s = count(f == 2) d = (n-1) / (2 ** s) loop: do i = 1, k a = int_rrand(2_huge, n-2) x = mod(a ** d, n) if ( x == 1 ) cycle do r = 0, s-1 if ( x == ( n - 1 ) ) cycle loop x = mod(x*x, n) end do if ( x == (n-1) ) cycle res = .false. return end do loop res = .true. end function miller_rabin_test
end module Miller_Rabin</lang>
Testing
<lang fortran>program TestMiller
use Miller_Rabin implicit none
integer, parameter :: prec = 30 integer(huge) :: i
! this is limited since we're not using a bignum lib call do_test( (/ (i, i=1, 29) /) )
contains
subroutine do_test(a) integer(huge), dimension(:), intent(in) :: a
integer :: i do i = 1, size(a,1) print *, a(i), miller_rabin_test(a(i), prec) end do
end subroutine do_test
end program TestMiller</lang>
Possible improvements: create bindings to the GMP library, change integer(huge)
into something like type(huge_integer)
, 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)
Python
<lang python>>>> import random >>> def miller_rabin (n,k):
if n<=3: return True if n % 2 == 0: return False d = n - 1 s = 0 while d%2 == 0: d = d / 2 s += 1 for i in xrange(k): a = random.randrange(2,n-1) x = pow(a, d, n) if x != 1: for r in range(s): if x == n-1: break x = x*x % n if x != n-1: return False return True
>>> [ i for i in range(1,1000) if miller_rabin(i, 10)][-10:] [937, 941, 947, 953, 967, 971, 977, 983, 991, 997] >>> </lang>
Ruby
<lang ruby>def miller_rabin_prime?(n,k)
return true if n == 2 return false if n < 2 or n % 2 == 0 d = n - 1 s = 0 while d % 2 == 0 d /= 2 s += 1 end k.times do a = 2 + rand(n-4) x = (a**d) % n next if x == 1 or x == n-1 for r in (1 .. s-1) x = (x**2) % n return false if x == 1 break if x == n-1 end return false if x != n-1 end true # probably
end
p primes = (1..100).find_all {|i| miller_rabin_prime?(i,10)}</lang>
[2, 3, 5, 7, 11, 13, 17, ..., 971, 977, 983, 991, 997]
Smalltalk
Smalltalk handles big numbers naturally and trasparently (the parent class Integer has many subclasses, and a subclass is picked according to the size of the integer that must be handled) <lang smalltalk>Integer extend [
millerRabinTest: kl [ |k| k := kl. self <= 3 ifTrue: [ ^true ] ifFalse: [ (self even) ifTrue: [ ^false ] ifFalse: [ |d s| d := self - 1. s := 0. [ (d rem: 2) == 0 ] whileTrue: [ d := d / 2. s := s + 1. ]. [ k:=k-1. k >= 0 ] whileTrue: [ |a x r| a := Random between: 2 and: (self - 2). x := (a raisedTo: d) rem: self. ( x = 1 ) ifFalse: [ |r|
r := -1.
[ r := r + 1. (r < s) & (x ~= (self - 1)) ] whileTrue: [ x := (x raisedTo: 2) rem: self ]. ( x ~= (self - 1) ) ifTrue: [ ^false ] ] ]. ^true ] ] ]
].</lang>
<lang smalltalk>1 to: 1000 do: [ :n |
(n millerRabinTest: 10) ifTrue: [ n printNl ]
].</lang>
Tcl
Use Tcl 8.5 for large integer support <lang tcl>package require Tcl 8.5
proc miller_rabin {n k} {
if {$n <= 3} {return true} if {$n % 2 == 0} {return false} # write n - 1 as 2^s·d with d odd by factoring powers of 2 from n − 1 set d [expr {$n - 1}] set s 0 while {$d % 2 == 0} { set d [expr {$d / 2}] incr s } while {$k > 0} { incr k -1 set a [expr {2 + int(rand()*($n - 4))}] set x [expr {($a ** $d) % $n}] if {$x == 1 || $x == $n - 1} continue for {set r 1} {$r < $s} {incr r} { set x [expr {($x ** 2) % $n}] if {$x == 1} {return false} if {$x == $n - 1} break }
if {$x != $n-1} {return false}
} return true
}
for {set i 1} {$i < 1000} {incr i} {
if {[miller_rabin $i 10]} { puts $i }
}</lang> produces
1 2 3 5 7 11 ... 977 983 991 997