Some of Sunday's edits have been lost. The edits from Saturday that were reverted have been restored. Site is now hosted on prgmr.com. Thank you for your patience. This notice will be removed one week from posting. --Michael Mol 18:12, 7 March 2010 (UTC)
Miller-Rabin primality test
From Rosetta Code
(Redirected from Miller-Rabin test)
Miller-Rabin primality test is a programming task. Visitors like you are encouraged to solve it according to the task description, using any language they may happen to 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 Free Documentation License. |
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 − 1]
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)
Contents |
[edit] ALGOL 68
Translation of: python
Works with: ALGOL 68 version Standard - with preludes manually inserted
Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386
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
Sample output:
937 941 947 953 967 971 977 983 991 997
[edit] AutoHotkey
ahk forum: discussion
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
}
[edit] Bc
Works with: GNU 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;
[edit] C
Library: GMP miller-rabin.h
#ifndef _MILLER_RABIN_H_
#define _MILLER_RABIN_H
#include <gmp.h>
bool miller_rabin_test(mpz_t n, int j);
#endif
miller-rabin.c
Translation of: Fortran
For decompose (and header primedecompose.h), see Prime decomposition.
#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;
}
Testing
#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;
}
[edit] Common Lisp
(defun factor-out (number divisor)
"Return two values R and E such that NUMBER = DIVISOR^E * R,
and R is not divisible by DIVISOR."
(do ((e 0 (1+ e))
(r number (/ r divisor)))
((/= (mod r divisor) 0) (values r e))))
(defun mult-mod (x y modulus) (mod (* x y) modulus))
(defun expt-mod (base exponent modulus)
"Fast modular exponentiation by repeated squaring."
(labels ((expt-mod-iter (b e p)
(cond ((= e 0) p)
((evenp e)
(expt-mod-iter (mult-mod b b modulus)
(/ e 2)
p))
(t
(expt-mod-iter b
(1- e)
(mult-mod b p modulus))))))
(expt-mod-iter base exponent 1)))
(defun random-in-range (lower upper)
"Return a random integer from the range [lower..upper]."
(+ lower (random (+ (- upper lower) 1))))
(defun miller-rabin-test (n k)
"Test N for primality by performing the Miller-Rabin test K times.
Return NIL if N is composite, and T if N is probably prime."
(cond ((= n 1) nil)
((< n 4) t)
((evenp n) nil)
(t
(multiple-value-bind (d s) (factor-out (- n 1) 2)
(labels ((strong-liar? (a)
(let ((x (expt-mod a d n)))
(or (= x 1)
(loop repeat s
for y = x then (mult-mod y y n)
thereis (= y (- n 1)))))))
(loop repeat k
always (strong-liar? (random-in-range 2 (- n 2)))))))))
CL-USER> (last (loop for i from 1 to 1000
when (miller-rabin-test i 10)
collect i)
10)
(937 941 947 953 967 971 977 983 991 997)
[edit] 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
}
for i ? (millerRabinPrimalityTest(i, 1, entropy)) in 4..1000 {
print(i, " ")
}
println()
[edit] Fortran
Works with: Fortran version 95
For the module PrimeDecompose, see Prime decomposition.
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
Testing
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
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)
[edit] Haskell
Works with: Haskell version 6.10.4
- Ideas taken from Primality proving
- Functions witns and isMillerRabinPrime follow closely the code outlined in J/Essays
- A useful powerMod function is taken from [1]
Another Miller Rabin test can be found in D. Amos's Haskell for Math module Primes.hs
import System.Random
import Data.List
import Control.Monad
import Control.Arrow
primesTo100 = [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]
powerMod :: (Integral a, Integral b) => a -> a -> b -> a
powerMod m _ 0 = 1
powerMod m x n | n > 0 = join (flip f (n - 1)) x `rem` m where
f _ 0 y = y
f a d y = g a d where
g b i | even i = g (b*b `rem` m) (i `quot` 2)
| otherwise = f b (i-1) (b*y `rem` m)
witns :: (Num a, Ord a, Random a) => Int -> a -> IO [a]
witns x y = do
g <- newStdGen
let r = [9080191, 4759123141, 2152302898747, 3474749600383, 341550071728321]
fs = [[31,73],[2,7,61],[2,3,5,7,11],[2,3,5,7,11,13],[2,3,5,7,11,13,17]]
if y >= 341550071728321
then return $ take x $ randomRs (2,y-1) g
else return $ snd.head.dropWhile ((<= y).fst) $ zip r fs
isMillerRabinPrime :: Integer -> IO Bool
isMillerRabinPrime n | n `elem` primesTo100 = return True
| otherwise = do
let pn = pred n
e = uncurry (++) . second(take 1) . span even . iterate (`div` 2) $ pn
try = return . all (\a -> let c = map (powerMod n a) e in
pn `elem` c || last c == 1)
witns 100 n >>= try
Testing in GHCi:
*~> isMillerRabinPrime 4547337172376300111955330758342147474062293202868155909489 True *~> isMillerRabinPrime 4547337172376300111955330758342147474062293202868155909393 False *~> filterM isMillerRabinPrime [2..1000]>>= print. dropWhile(<=900) [907,911,919,929,937,941,946,947,953,967,971,977,983,991,997]
[edit] J
See Primality Tests essay on the J wiki.
[edit] Perl
use bigint;
sub is_prime
{
my ($n,$k) = @_;
return 1 if $n == 2;
return 0 if $n < 2 or $n % 2 == 0;
$d = $n - 1;
$s = 0;
while(!($d % 2))
{
$d /= 2;
$s++;
}
LOOP: for(1..$k)
{
$a = 2 + int(rand($n-2));
$x = ($a**$d) % $n;
next if $x == 1 or $x == $n-1;
for(1..$s-1)
{
$x = ($x**2) % $n;
return 0 if $x == 1;
next LOOP if $x == $n-1;
}
return 0;
}
return 1;
}
print join ", ", grep { is_prime $_,10 }(1..1000);
[edit] Python
import random
_mrpt_num_trials = 5 # number of bases to test
def is_probable_prime(n):
"""
Miller-Rabin primality test.
A return value of False means n is certainly not prime. A return value of
True means n is very likely a prime.
>>> is_probable_prime(1)
Traceback (most recent call last):
...
AssertionError
>>> is_probable_prime(2)
True
>>> is_probable_prime(3)
True
>>> is_probable_prime(4)
False
>>> is_probable_prime(5)
True
>>> is_probable_prime(123456789)
False
>>> primes_under_1000 = [i for i in range(2, 1000) if is_probable_prime(i)]
>>> len(primes_under_1000)
168
>>> primes_under_1000[-10:]
[937, 941, 947, 953, 967, 971, 977, 983, 991, 997]
>>> is_probable_prime(6438080068035544392301298549614926991513861075340134\
3291807343952413826484237063006136971539473913409092293733259038472039\
7133335969549256322620979036686633213903952966175107096769180017646161\
851573147596390153)
True
>>> is_probable_prime(7438080068035544392301298549614926991513861075340134\
3291807343952413826484237063006136971539473913409092293733259038472039\
7133335969549256322620979036686633213903952966175107096769180017646161\
851573147596390153)
False
"""
assert n >= 2
# special case 2
if n == 2:
return True
# ensure n is odd
if n % 2 == 0:
return False
# write n-1 as 2**s * d
# repeatedly try to divide n-1 by 2
s = 0
d = n-1
while True:
quotient, remainder = divmod(d, 2)
if remainder == 1:
break
s += 1
d = quotient
assert(2**s * d == n-1)
# test the base a to see whether it is a witness for the compositeness of n
def try_composite(a):
if pow(a, d, n) == 1:
return False
for i in range(s):
if pow(a, 2**i * d, n) == n-1:
return False
return True # n is definitely composite
for i in range(_mrpt_num_trials):
a = random.randint(2, n-1)
if try_composite(a):
return False
return True # no base tested showed n as composite
[edit] 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)}
[2, 3, 5, 7, 11, 13, 17, ..., 971, 977, 983, 991, 997]
[edit] Smalltalk
Works with: GNU 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)
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
]
]
]
].
1 to: 1000 do: [ :n |
(n millerRabinTest: 10) ifTrue: [ n printNl ]
].
[edit] Tcl
Use Tcl 8.5 for large integer support
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
}
}
produces
1 2 3 5 7 11 ... 977 983 991 997
Categories: Programming Tasks | Prime Numbers | ALGOL 68 | AutoHotkey | Bc | C | GMP | Common Lisp | E | Fortran | Haskell | J | Perl | Python | Ruby | Smalltalk | Tcl
