Miller-Rabin primality test
From Rosetta Code
(Redirected from Miller-Rabin 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 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] PicoLisp
(de longRand (N)
(use (R D)
(while (=0 (setq R (abs (rand)))))
(until (> R N)
(unless (=0 (setq D (abs (rand))))
(setq R (* R D)) ) )
(% R N) ) )
(de **Mod (X Y N)
(let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )
(de _prim? (N D S)
(use (A X R)
(while (> 2 (setq A (longRand N))))
(setq R 0 X (**Mod A D N))
(loop
(T
(or
(and (=0 R) (= 1 X))
(= X (dec N)) )
T )
(T
(or
(and (> R 0) (= 1 X))
(>= (inc 'R) S) )
NIL )
(setq X (% (* X X) N)) ) ) )
(de prime? (N K)
(default K 50)
(and
(> N 1)
(bit? 1 N)
(let (D (dec N) S 0)
(until (bit? 1 D)
(setq
D (>> 1 D)
S (inc S) ) )
(do K
(NIL (_prim? N D S))
T ) ) ) )
Output:
: (filter '((I) (prime? I)) (range 937 1000)) -> (937 941 947 953 967 971 977 983 991 997) : (prime? 4547337172376300111955330758342147474062293202868155909489) -> T : (prime? 4547337172376300111955330758342147474062293202868155909393) -> NIL
[edit] PureBasic
Enumeration
#Composite
#Probably_prime
EndEnumeration
Procedure Miller_Rabin(n, k)
Protected d=n-1, s, x, r
If n=2
ProcedureReturn #Probably_prime
ElseIf n%2=0 Or n<2
ProcedureReturn #Composite
EndIf
While d%2=0
d/2
s+1
Wend
While k>0
k-1
x=Int(Pow(2+Random(n-4),d))%n
If x=1 Or x=n-1: Continue: EndIf
For r=1 To s-1
x=(x*x)%n
If x=1: ProcedureReturn #Composite: EndIf
If x=n-1: Break: EndIf
Next
If x<>n-1: ProcedureReturn #Composite: EndIf
Wend
ProcedureReturn #Probably_prime
EndProcedure
[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 | WikipediaSourced | ALGOL 68 | AutoHotkey | Bc | C | GMP | Common Lisp | E | Fortran | Haskell | J | Perl | PicoLisp | PureBasic | Python | Ruby | Smalltalk | Tcl
