Miller–Rabin primality test

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 − 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)

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

<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

Works with: GNU bc

<lang bc>next = 1; # seed of the random generator scale = 0;

1. random generator (according to POSIX...); since
2. this gives number between 0 and 32767, we can
3. 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

Library: GMP

miller-rabin.h

<lang c>#ifndef _MILLER_RABIN_H_

1. define _MILLER_RABIN_H
2. include <gmp.h>

bool miller_rabin_test(mpz_t n, int j);

1. endif</lang>

miller-rabin.c

Translation of: Fortran

For decompose (and header primedecompose.h), see Prime decomposition.

<lang c>#include <stdbool.h>

1. include <gmp.h>
2. include "primedecompose.h"

1. 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>

1. include <stdlib.h>
2. include <stdbool.h>
3. include <gmp.h>
4. include "miller-rabin.h"

1. define PREC 10
2. 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>

Common Lisp

<lang 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)))))))))

</lang>

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)

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

Works with: Fortran version 95

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)

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

<lang Haskell> 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 $ (fs!!).head.findIndices (> y) $ r

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</lang>

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]

J

See Primality Tests essay on the J wiki.

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

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) <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