

 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
   repeat s − 1 times:
      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)


 
\ modular multiplication and exponentiation
\
: 3rd  s" 2 pick" evaluate ; immediate
 
: mod*  ( a b m -- a*b {mod m} )
    >r um* r> ud/mod 2drop ;
 
: mod^  ( x n m -- x^n {mod m} )
    >r 1 swap
    begin ?dup while
        dup 1 and 1 =
        if
            swap 3rd r@ mod* swap 1-
        then dup 0>
        if
            rot dup r@ mod* -rot 2/
        then
    repeat nip rdrop ;
 
\ small divisor check: true => possibly prime; false => definitely not prime.
\
31 constant π-128
create maybe-prime?
      2 c,  3 c,  5 c,  7 c, 11 c,  13 c,  17 c,  19 c,  23 c,  29 c,
     31 c, 37 c, 41 c, 43 c, 47 c,  53 c,  59 c,  61 c,  67 c,  71 c,
     73 c, 79 c, 83 c, 89 c, 97 c, 101 c, 103 c, 107 c, 109 c, 113 c,
    127 c,
does>
    true -rot
    π-128 bounds do
        i c@ dup * over > if leave then
        dup i c@ mod 0= if 2drop false unloop exit then
    loop drop ;
 
\ actual Miller-Rabin test
\
: factor-2s  ( n -- s d )
    0 swap
    begin dup 1 and 0= while
        swap 1+  swap 2/
    repeat ;
 
: fermat-square-test ( n m s -- ? ) \ perform n = n^2 (mod m), s-1 times
    1- 0 ?do
        2dup - -1 =
        if   leave
        then >r dup r@ mod* r>
    loop
    - -1 = ;
 
: strong-fermat-pseudoprime? ( n a -- ? )
    over >r \ keep the modulus on the return stack
    >r 1- factor-2s r>  \ -- s d a
    swap r@ mod^        \ s d a -- s, a^d (mod n)
    dup 1 =             \ a^d == 1 (mod n) => Fermat pseudoprime
    if   2drop rdrop true
    else r> rot fermat-square-test
    then ;
 
4.759.123.141 drop constant mr-det-3   \ Deterministic threshold; 3 bases
 
create small-prime-bases 2 , 7 , 61 ,  \ deterministic up to mr-det-3
create large-prime-bases 2 , 325 , 9375 , 28178 , 450775 , 9780504 , 1795265022 , \ known to be deterministic for 64 bit integers.
 
: miler-rabin-bases ( n -- addr n )
    mr-det-3 <
    if    small-prime-bases 3
    else  large-prime-bases 7
    then ;
 
: miller-rabin-primality-test ( n -- f )
    dup miler-rabin-bases cells bounds do
        dup i @ strong-fermat-pseudoprime? invert
        if drop false unloop exit then
    cell +loop drop true ;
 
: prime? ( n -- f )
    dup 2 <
    if drop false
    else
        dup maybe-prime?
        if dup [ 127 dup * 1+ ] literal < 
            if   drop true
            else miller-rabin-primality-test
            then
        else drop false
        then
    then ;
 