Miller–Rabin primality test

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.

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 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 2^s·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 ← a^d mod n
   if x = 1 or x = n − 1 then do next LOOP
   for r = 1 .. s − 1
      x ← x² 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)

Ada

ordinary integers

It's easy to get overflows doing exponential calculations. Therefore I implemented my own function for that.

For Number types >= 2**64 you may have to use an external library -- see below.

First, a package Miller_Rabin is specified. The same package is used else elsewhere in Rosetta Code, such as for the Carmichael 3 strong pseudoprimes.

<lang Ada>generic

  type Number is range <>;

package Miller_Rabin is

  type Result_Type is (Composite, Probably_Prime);

  function Is_Prime (N : Number; K : Positive := 10) return Result_Type;

end Miller_Rabin;</lang>

The implementation of that package is as follows:

<lang Ada>with Ada.Numerics.Discrete_Random;

package body Miller_Rabin is

  function Is_Prime (N : Number; K : Positive := 10)
                    return Result_Type
  is
     subtype Number_Range is Number range 2 .. N - 1;
     package Random is new Ada.Numerics.Discrete_Random (Number_Range);

     function Mod_Exp (Base, Exponent, Modulus : Number) return Number is
        Result : Number := 1;
     begin
        for E in 1 .. Exponent loop
           Result := Result * Base mod Modulus;
        end loop;
        return Result;
     end Mod_Exp;

     Generator : Random.Generator;
     D : Number := N - 1;
     S : Natural := 0;
     X : Number;
  begin
     -- exclude 2 and even numbers
     if N = 2 then
        return Probably_Prime;
     elsif N mod 2 = 0 then
        return Composite;
     end if;

     -- write N-1 as 2**S * D, with D mod 2 /= 0
     while D mod 2 = 0 loop
        D := D / 2;
        S := S + 1;
     end loop;

     -- initialize RNG
     Random.Reset (Generator);
     for Loops in 1 .. K loop
        X := Mod_Exp(Random.Random (Generator), D, N);
        if X /= 1 and X /= N - 1 then
       Inner : for R in 1 .. S - 1 loop
              X := Mod_Exp (X, 2, N);
              if X = 1 then return Composite; end if;
              exit Inner when X = N - 1;
           end loop Inner;
           if X /= N - 1 then return Composite; end if;
        end if;
     end loop;

     return Probably_Prime;
  end Is_Prime;

end Miller_Rabin;</lang>

Finally, the program itself:

<lang Ada>with Ada.Text_IO, Miller_Rabin;

procedure Mr_Tst is

  type Number is range 0 .. (2**48)-1;

  package Num_IO is new Ada.Text_IO.Integer_IO (Number);
  package Pos_IO is new Ada.Text_IO.Integer_IO (Positive);
  package MR     is new Miller_Rabin(Number); use MR;

  N : Number;
  K : Positive;

begin

  for I in Number(2) .. 1000 loop
     if Is_Prime (I) = Probably_Prime then
        Ada.Text_IO.Put (Number'Image (I));
     end if;
  end loop;
  Ada.Text_IO.Put_Line (".");

  Ada.Text_IO.Put ("Enter a Number: ");           Num_IO.Get (N);
  Ada.Text_IO.Put ("Enter the count of loops: "); Pos_IO.Get (K);
  Ada.Text_IO.Put_Line ("What is it? " & Result_Type'Image (Is_Prime(N, K)));

end MR_Tst;</lang>

Output:

 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 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997.
Enter a Number: 1234567
Enter the count of loops: 20
What is it? COMPOSITE

using an external library to handle big integers

Using the big integer implementation from a cryptographic library [1].

<lang Ada>with Ada.Text_IO, Crypto.Types.Big_Numbers, Ada.Numerics.Discrete_Random;

procedure Miller_Rabin is

  Bound: constant Positive := 256; -- can be any multiple of 32

  package LN is new Crypto.Types.Big_Numbers (Bound);
  use type LN.Big_Unsigned; -- all computations "mod 2**Bound"

  function "+"(S: String) return LN.Big_Unsigned
    renames LN.Utils.To_Big_Unsigned;

  function Is_Prime (N : LN.Big_Unsigned; K : Positive := 10) return Boolean is

     subtype Mod_32 is Crypto.Types.Mod_Type;
     use type Mod_32;
     package R_32 is new Ada.Numerics.Discrete_Random (Mod_32);
     Generator : R_32.Generator;

     function Random return LN.Big_Unsigned is
        X: LN.Big_Unsigned := LN.Big_Unsigned_Zero;
     begin
        for I in 1 .. Bound/32 loop
           X := (X * 2**16) * 2**16;
           X := X + R_32.Random(Generator);
        end loop;
        return X;
     end Random;

     D:    LN.Big_Unsigned := N - LN.Big_Unsigned_One;
     S:    Natural := 0;
     A, X: LN.Big_Unsigned;
  begin
     -- exclude 2 and even numbers
     if N = 2 then
        return True;
     elsif N mod 2 = LN.Big_Unsigned_Zero then
        return False;
     else

        -- write N-1 as 2**S * D, with odd D
        while D mod 2 = LN.Big_Unsigned_Zero loop
           D := D / 2;
           S := S + 1;
        end loop;

        -- initialize RNG
        R_32.Reset (Generator);

        -- run the real test
        for Loops in 1 .. K loop
           loop
              A := Random;
              exit when (A > 1) and (A < (N - 1));
           end loop;
           X := LN.Mod_Utils.Pow(A, D, N); -- X := (Random**D) mod N
           if X /= 1 and X /= N - 1 then
              Inner:
              for R in 1 .. S - 1 loop
                 X := LN.Mod_Utils.Pow(X, LN.Big_Unsigned_Two, N);
                 if X = 1 then
                    return False;
                 end if;
                 exit Inner when X = N - 1;
              end loop Inner;
              if X /= N - 1 then
                 return False;
              end if;
           end if;
        end loop;
     end if;

     return True;
  end Is_Prime;

  S: constant String :=
    "4547337172376300111955330758342147474062293202868155909489";
  T: constant String :=
    "4547337172376300111955330758342147474062293202868155909393";

  K: constant Positive := 10;

begin

  Ada.Text_IO.Put_Line("Prime(" & S & ")=" & Boolean'Image(Is_Prime(+S, K)));
  Ada.Text_IO.Put_Line("Prime(" & T & ")=" & Boolean'Image(Is_Prime(+T, K)));

end Miller_Rabin;</lang>

Output:

Prime(4547337172376300111955330758342147474062293202868155909489)=TRUE
Prime(4547337172376300111955330758342147474062293202868155909393)=FALSE

Using the built-in Miller-Rabin test from the same library:

<lang Ada>with Ada.Text_IO, Crypto.Types.Big_Numbers, Ada.Numerics.Discrete_Random;

procedure Miller_Rabin is

  Bound: constant Positive := 256; -- can be any multiple of 32

  package LN is new Crypto.Types.Big_Numbers (Bound);
  use type LN.Big_Unsigned; -- all computations "mod 2**Bound"

  function "+"(S: String) return LN.Big_Unsigned
    renames LN.Utils.To_Big_Unsigned;

  S: constant String :=
    "4547337172376300111955330758342147474062293202868155909489";
  T: constant String :=
    "4547337172376300111955330758342147474062293202868155909393";

  K: constant Positive := 10;

begin

  Ada.Text_IO.Put_Line("Prime(" & S & ")="
      & Boolean'Image (LN.Mod_Utils.Passed_Miller_Rabin_Test(+S, K)));
  Ada.Text_IO.Put_Line("Prime(" & T & ")="
      & Boolean'Image (LN.Mod_Utils.Passed_Miller_Rabin_Test(+T, K)));

end Miller_Rabin;</lang>

The output is the same.

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

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

Requires a bc with long names.

Works with: OpenBSD bc

(A previous version worked with GNU bc.) <lang bc>seed = 1 /* seed of the random number generator */ scale = 0

/* Random number from 0 to 32767. */ define rand() {

 /* Cheap formula (from POSIX) for random numbers of low quality. */
 seed = (seed * 1103515245 + 12345) % 4294967296
 return ((seed / 65536) % 32768)

}

/* Random number in range [from, to]. */ define rangerand(from, to) {

 auto b, h, i, m, n, r

 m = to - from + 1
 h = length(m) / 2 + 1  /* want h iterations of rand() % 100 */
 b = 100 ^ h % m        /* want n >= b */
 while (1) {
   n = 0                /* pick n in range [b, 100 ^ h) */
   for (i = h; i > 0; i--) {
     r = rand()
     while (r < 68) { r = rand(); }  /* loop if the modulo bias */
     n = (n * 100) + (r % 100)       /* append 2 digits to n */
   }
   if (n >= b) { break; }  /* break unless the modulo bias */
 }
 return (from + (n % m))

}

/* n is probably prime? */ 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 /= 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>

Bracmat

Translation of: bc

<lang bracmat>( 1:?seed & ( rand

 =   
   .   mod$(!seed*1103515245+12345.4294967296):?seed
     & mod$(div$(!seed.65536).32768)
 )

& ( rangerand

 =   from to b h i m n r length
   .   !arg:(?from,?to)
     & !to+-1*!from+1:?m
     & @(!m:? [?length)
     & div$(!length+1.2)+1:?h
     & 100^mod$(!h.!m):?b
     &   whl
       ' ( 0:?n
         & !h+1:?i
         &   whl
           ' ( !i+-1:>0:?i
             & rand$:?r
             & whl'(!r:<68&rand$:?r)
             & !n*100+mod$(!r.100):?n
             )
         & !n:>!b
         )
     & !from+mod$(!n.!m)
 )

& ( miller-rabin-test

 =   n k d r a x s return
   .   !arg:(?n,?k)
     & ( !n:~>3&1
       | mod$(!n.2):0
       |   !n+-1:?d
         & 0:?s
         &   whl
           ' ( mod$(!d.2):0
             & !d*1/2:?d
             & 1+!s:?s
             )
         & 1:?return
         &   whl
           ' ( !k+-1:?k:~<0
             & rangerand$(2,!n+-2):?a
             & mod$(!a^!d.!n):?x
             & ( !x:1
               |   0:?r
                 &   whl
                   ' ( !r+1:~>!s:?r
                     & !n+-1:~!x
                     & mod$(!x*!x.!n):?x
                     )
                 & ( !n+-1:!x
                   | 0:?return&~
                   )
               )
             )
         & !return
       )
 )

& 0:?i & :?primes & whl

 ' ( 1+!i:<1000:?i
   & (   miller-rabin-test$(!i,10):1
       & !primes !i:?primes
     | 
     )
   )

& !primes:? [-11 ?last & out$!last );</lang> output:

937 941 947 953 967 971 977 983 991 997

C

Library: GMP

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

Translation of: Fortran

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>

Deterministic up to 341,550,071,728,321

<lang c>// calcul a^n%mod size_t power(size_t a, size_t n, size_t mod) {

   size_t power = a;
   size_t result = 1;

   while (n)
   {
       if (n & 1)
           result = (result * power) % mod;
       power = (power * power) % mod;
       n >>= 1;
   }
   return result;

}

// n−1 = 2^s * d with d odd by factoring powers of 2 from n−1 bool witness(size_t n, size_t s, size_t d, size_t a) {

   size_t x = power(a, d, n);
   size_t y;

   while (s) {
       y = (x * x) % n;
       if (y == 1 && x != 1 && x != n-1)
           return false;
       x = y;
       --s;
   }
   if (y != 1)
       return false;
   return true;

}

/*

* if n < 1,373,653, it is enough to test a = 2 and 3;
* if n < 9,080,191, it is enough to test a = 31 and 73;
* if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
* if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;
* if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;
* if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;
* if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.
*/

bool is_prime_mr(size_t n) {

   if (((!(n & 1)) && n != 2 ) || (n < 2) || (n % 3 == 0 && n != 3))
       return false;
   if (n <= 3)
       return true;

   size_t d = n / 2;
   size_t s = 1;
   while (!(d & 1)) {
       d /= 2;
       ++s;
   }

   if (n < 1373653)
       return witness(n, s, d, 2) && witness(n, s, d, 3);
   if (n < 9080191)
       return witness(n, s, d, 31) && witness(n, s, d, 73);
   if (n < 4759123141)
       return witness(n, s, d, 2) && witness(n, s, d, 7) && witness(n, s, d, 61);
   if (n < 1122004669633)
       return witness(n, s, d, 2) && witness(n, s, d, 13) && witness(n, s, d, 23) && witness(n, s, d, 1662803);
   if (n < 2152302898747)
       return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11);
   if (n < 3474749660383)
       return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13);
   return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13) && witness(n, s, d, 17);

}</lang> Inspiration from http://stackoverflow.com/questions/4424374/determining-if-a-number-is-prime

C#

<lang csharp>public static class RabinMiller {

   public static bool IsPrime(int n, int k)
   {

if(n < 2)

return false;

if(n != 2 && n % 2 == 0)

return false;

int s = n - 1; while(s % 2 == 0)

s >>= 1;

       }
       Random r = new Random();

for (int i = 0; i < k; i++)

       {
           double a = r.Next((int)n - 1) + 1;
           int temp = s;
           int mod = (int)Math.Pow(a, (double)temp) % n;
           while(temp != n - 1 && mod != 1 && mod != n - 1)
           {

mod = (mod * mod) % n; temp = temp * 2;

if(mod != n - 1 && temp % 2 == 0)

return false;

           }
       }

return true;

}</lang> <lang csharp>// Miller-Rabin primality test as an extension method on the BigInteger type. // Based on the Ruby implementation on this page. public static class BigIntegerExtensions {

 public static bool IsProbablePrime(this BigInteger source, int certainty)
 {
   if(source == 2 || source == 3)
     return true;
   if(source < 2 || source % 2 == 0)
     return false;

   BigInteger d = source - 1;
   int s = 0;

   while(d % 2 == 0)
   {
     d /= 2;
     s += 1;
   }

   // There is no built-in method for generating random BigInteger values.
   // Instead, random BigIntegers are constructed from randomly generated
   // byte arrays of the same length as the source.
   RandomNumberGenerator rng = RandomNumberGenerator.Create();
   byte[] bytes = new byte[source.ToByteArray().LongLength];
   BigInteger a;

   for(int i = 0; i < certainty; i++)
   {
     do
     {
       // This may raise an exception in Mono 2.10.8 and earlier.
       // http://bugzilla.xamarin.com/show_bug.cgi?id=2761
       rng.GetBytes(bytes);
       a = new BigInteger(bytes);
     }
     while(a < 2 || a >= source - 2);

     BigInteger x = BigInteger.ModPow(a, d, source);
     if(x == 1 || x == source - 1)
       continue;

     for(int r = 1; r < s; r++)
     {
       x = BigInteger.ModPow(x, 2, source);
       if(x == 1)
         return false;
       if(x == source - 1)
         break;
     }

     if(x != source - 1)
       return false;
   }

   return true;
 }

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

D

Translation of: Ruby

<lang d>import std.random;

bool isProbablePrime(in ulong n, in int k) {

   static long modPow(long b, long e, in long m)
   pure nothrow {
       long result = 1;
       while (e > 0) {
           if ((e & 1) == 1) {
               result = (result * b) % m;
           }
           b = (b * b) % m;
           e >>= 1;
       }
       return result;
   }

   if (n < 2 || n % 2 == 0) 
       return n == 2;

   ulong d = n - 1;
   ulong s = 0;
   while (d % 2 == 0) {
       d /= 2;
       s++;
   }
   assert(2 ^^ s * d == n - 1);

   outer:
   foreach (_; 0 .. k) {
       ulong a = uniform(2, n);
       ulong x = modPow(a, d, n);
       if (x == 1 || x == n - 1)
           continue;
       foreach (__; 1 .. s) {
           x = modPow(x, 2, n);
           if (x == 1) return false;
           if (x == n - 1) continue outer;
       }
       return false;
   }

   return true;

}

void main() { // demo code

   import std.stdio, std.range, std.algorithm;
   writeln(filter!(n => isProbablePrime(n, 10))(iota(2, 30)));

}</lang>

Output:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

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>

Erlang_Fortran

<lang erlang>-module(miller_rabin).

-export([is_prime/1, power/2]).

% Replaced floating point functions with integer functions and arithmetic, % and added an integer power function. % Increased number of probabilistic trials from 20 to 100, % and augmented number of deterministic proving bases. % Dogwood, January 18, 2014 @ 11:45am PST.

is_prime(1) -> false; is_prime(2) -> true; is_prime(N) when N > 2, ((N rem 2) == 0) -> false; is_prime(N) when ((N rem 2) ==1), N < 341550071728321 ->

			is_mr_prime(N, proving_bases(N));

is_prime(N) when ((N rem 2) == 1) -> is_mr_prime(N, random_bases(N, 100)).

proving_bases(N) when N < 1373653 -> [2, 3]; proving_bases(N) when N < 9080191 ->

   [31, 73];

proving_bases(N) when N < 25326001 -> [2, 3, 5]; proving_bases(N) when N < 3215031751 -> [2, 3, 5, 7]; proving_bases(N) when N < 4759123141 ->

   [2, 7, 61];

proving_bases(N) when N < 1122004669633 -> [2, 13, 23, 1662803]; proving_bases(N) when N < 2152302898747 -> [2, 3, 5, 7, 11]; proving_bases(N) when N < 3474749660383 ->

   [2, 3, 5, 7, 11, 13];

proving_bases(N) when N < 341550071728321 ->

   [2, 3, 5, 7, 11, 13, 17].

is_mr_prime(N, As) when N>2, N rem 2 == 1 ->

   {D, S} = find_ds(N),
   not lists:any(fun(A) ->
                         case mr_series(N, A, D, S) of
                             [1|_] -> false;
                             L     -> not lists:member(N-1, L)
                         end
                 end,
                 As).

find_ds(N) ->

   find_ds(N-1, 0).

find_ds(D, S) ->

   case D rem 2 == 0 of
       true ->
           find_ds(D div 2, S+1);
       false ->
           {D, S}
   end.

mr_series(N, A, D, S) when N rem 2 == 1 ->

   Js = lists:seq(0, S),
   lists:map(fun(J) -> pow_mod(A, power(2, J)*D, N) end, Js).

pow_mod(B, E, M) ->

   case E of
       0 -> 1;
       _ -> case ((E rem 2) == 0) of
                true  -> (power(pow_mod(B, (E div 2), M), 2)) rem M;
                false -> (B*pow_mod(B, E-1, M)) rem M
            end
   end.

random_bases(N, K) ->

   [basis(N) || _ <- lists:seq(1, K)].

basis(N) when N>2 ->

   1 + random:uniform(N-2).

power(B, E) ->

   power(B, E, 1).

power(_, 0, Acc) ->

   Acc;

power(B, E, Acc) ->

   power(B, E - 1, B * Acc).

%%%%%%% HAVE UNINTENTIONALLY MERGED FORTRAN PAGE AND UNABLE TO REPAIR %%%%%%%%

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)

Go

Library

Go has it in math/big in standard library as ProbablyPrime. The argument n to ProbablyPrime is the input k of the pseudocode in the task description.

Deterministic

Below is a deterministic test for 32 bit unsigned integers. Intermediate results in the code below include a 64 bit result from multiplying two 32 bit numbers. Since 64 bits is the largest fixed integer type in Go, a 32 bit number is the largest that is convenient to test.

The main difference between this algorithm and the pseudocode in the task description is that k numbers are not chosen randomly, but instead are the three numbers 2, 7, and 61. These numbers provide a deterministic primality test up to 2^32. <lang go>package main

import "log"

func main() {

   // max uint32 is not prime
   c := uint32(1<<32 - 1)
   // a few primes near the top of the range.  source: prime pages.
   for _, p := range []uint32{1<<32 - 5, 1<<32 - 17, 1<<32 - 65, 1<<32 - 99} {
       for ; c > p; c-- {
           if prime(c) {
               log.Fatalf("prime(%d) returned true", c)
           }
       }
       if !prime(p) {
           log.Fatalf("prime(%d) returned false", p)
       }
       c--
   }

}

func prime(n uint32) bool {

   // bases of 2, 7, 61 are sufficient to cover 2^32
   switch n {
   case 0, 1:
       return false
   case 2, 7, 61:
       return true
   }
   // compute s, d where 2^s * d = n-1
   nm1 := n - 1
   d := nm1
   s := 0
   for d&1 == 0 {
       d >>= 1
       s++
   }
   n64 := uint64(n)
   for _, a := range []uint32{2, 7, 61} {
       // compute x := a^d % n
       x := uint64(1)
       p := uint64(a)
       for dr := d; dr > 0; dr >>= 1 {
           if dr&1 != 0 {
               x = x * p % n64
           }
           p = p * p % n64
       }
       if x == 1 || uint32(x) == nm1 {
           continue
       }
       for r := 1; ; r++ {
           if r >= s {
               return false
           }
           x = x * x % n64
           if x == 1 {
               return false
           }
           if uint32(x) == nm1 {
               break
           }
       }
   }
   return true

}</lang>

Haskell

Works with: Haskell version 7.6.3

Ideas taken from Primality proving
Functions witns and isMillerRabinPrime follow closely the code outlined in J/Essays
A useful powerMod function is taken from Multiplicative order#Haskell
Original Rosetta code has been simplified to be easier to follow

Another Miller Rabin test can be found in D. Amos's Haskell for Math module Primes.hs <lang Haskell>module Primes where

import System.Random import System.IO.Unsafe

-- Miller-Rabin wrapped up as an (almost deterministic) pure function isPrime :: Integer -> Bool isPrime n = unsafePerformIO (isMillerRabinPrime 100 n)

isMillerRabinPrime :: Int -> Integer -> IO Bool isMillerRabinPrime k n

  | even n    = return (n==2)
  | n < 100   = return (n `elem` primesTo100)
  | otherwise = do ws <- witnesses k n
                   return $ and [test n (pred n) evens (head odds) a | a <- ws]
 where
   (evens,odds) = span even (iterate (`div` 2) (pred n))

test :: Integral nat => nat -> nat -> [nat] -> nat -> nat -> Bool test n n_1 evens d a = x `elem` [1,n_1] || n_1 `elem` powers

 where
   x = powerMod n a d
   powers = map (powerMod n a) evens

witnesses :: (Num a, Ord a, Random a) => Int -> a -> IO [a] witnesses k n

 | n < 9080191         = return [31,73]
 | n < 4759123141      = return [2,7,61]
 | n < 3474749660383   = return [2,3,5,7,11,13]
 | n < 341550071728321 = return [2,3,5,7,11,13,17]
 | otherwise           = do g <- newStdGen
                            return $ take k (randomRs (2,n-1) g)

primesTo100 :: [Integer] 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 m x n = x^n `mod` m powerMod :: Integral nat => nat -> nat -> nat -> nat powerMod m x n = f (n - 1) x x `rem` m

 where
 f d a y = if d==0 then y else g d a y 
 g i b y | even i    = g (i `quot` 2) (b*b `rem` m) y
         | otherwise = f (i-1) b (b*y `rem` m)

</lang>

Sample output:

Testing in GHCi:
~> isPrime 4547337172376300111955330758342147474062293202868155909489
True

*~> isPrime 4547337172376300111955330758342147474062293202868155909393
False

*~> dropWhile (<900) $ filter isPrime [2..1000]
[907,911,919,929,937,941,947,953,967,971,977,983,991,997]

Icon and Unicon

The following runs in both languages: <lang unicon>procedure main(A)

  every n := !A do write(n," is ",(mrp(n,5),"probably prime")|"composite")

end

procedure mrp(n, k)

   if n = 2 then return ""
   if n%2 = 0 then fail
   nm1 := decompose(n-1)
   s := nm1[1]
   d := nm1[2]
   every !k do {
       a := ?(n-2)+1
       x := (a^d)%n
       if x = (1|(n-1)) then next
       every !(s-1) do {
           x := (x*x)%n
           if x = 1 then fail
           if x = (n-1) then break next
           }
       fail
       }
   return ""

end

procedure decompose(nm1)

   s := 1
   d := nm1
   while d%2 = 0 do {
       d /:= 2
       s +:= 1
       }
   return [s,d]

end</lang>

Sample run:

->mrp 219 221 223 225 227 229
219 is composite
221 is composite
223 is probably prime
225 is composite
227 is probably prime
229 is probably prime
->

J

See Primality Tests essay on the J wiki.

Java

The Miller-Rabin primality test is part of the standard library (java.math.BigInteger) <lang java>import java.math.BigInteger;

public class MillerRabinPrimalityTest {

 public static void main(String[] args) {
   BigInteger n = new BigInteger(args[0]);
   int certainty = Integer.parseInt(args[1]);
   System.out.println(n.toString() + " is " + (n.isProbablePrime(certainty) ? "probably prime" : "composite"));
 }

}</lang>

Sample output:

java MillerRabinPrimalityTest 123456791234567891234567 1000000
123456791234567891234567 is probably prime

JavaScript

This covers (almost) all integers in JavaScript (~2^53).

<lang JavaScript>function modProd(a,b,n){

 if(b==0) return 0;
 if(b==1) return a%n;
 return (modProd(a,(b-b%10)/10,n)*10+(b%10)*a)%n;

} function modPow(a,b,n){

 if(b==0) return 1;
 if(b==1) return a%n;
 if(b%2==0){
   var c=modPow(a,b/2,n);
   return modProd(c,c,n);
 }
 return modProd(a,modPow(a,b-1,n),n);

} function isPrime(n){

 if(n==2||n==3||n==5) return true;
 if(n%2==0||n%3==0||n%5==0) return false;
 if(n<25) return true;
 for(var a=[2,3,5,7,11,13,17,19],b=n-1,d,t,i,x;b%2==0;b/=2);
 for(i=0;i<a.length;i++){
   x=modPow(a[i],b,n);
   if(x==1||x==n-1) continue;
   for(t=true,d=b;t&&d<n-1;d*=2){
     x=modProd(x,x,n); if(x==n-1) t=false;
   }
   if(t) return false;
 }
 return true;

}

for(var i=1;i<=1000;i++) if(isPrime(i)) console.log(i);</lang>

Julia

The built-in isprime function uses the Miller-Rabin primality test. Here is the implementation of isprime from the Julia standard library (Julia version 0.2): <lang julia> witnesses(n::Union(Uint8,Int8,Uint16,Int16)) = (2,3) witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? (2,3) : (2,7,61) witnesses(n::Union(Uint64,Int64)) =

       n < 1373653         ? (2,3) :
       n < 4759123141      ? (2,7,61) :
       n < 2152302898747   ? (2,3,5,7,11) :
       n < 3474749660383   ? (2,3,5,7,11,13) :
                             (2,325,9375,28178,450775,9780504,1795265022)

function isprime(n::Integer)

   n == 2 && return true
   (n < 2) | iseven(n) && return false
   s = trailing_zeros(n-1)
   d = (n-1) >>> s
   for a in witnesses(n)
       a < n || break
       x = powermod(a,d,n)
       x == 1 && continue
       t = s
       while x != n-1
           (t-=1) <= 0 && return false
           x = oftype(n, Base.widemul(x,x) % n)
           x == 1 && return false
       end
   end
   return true

end </lang>

Liberty BASIC

<lang lb> DIM mersenne(11) mersenne(1)=7 mersenne(2)=31 mersenne(3)=127 mersenne(4)=8191 mersenne(5)=131071 mersenne(6)=524287 mersenne(7)=2147483647 mersenne(8)=2305843009213693951 mersenne(9)=618970019642690137449562111 mersenne(10)=162259276829213363391578010288127 mersenne(11)=170141183460469231731687303715884105727

dim SmallPrimes(1000) data 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 data 31, 37, 41, 43, 47, 53, 59, 61, 67, 71 data 73, 79, 83, 89, 97, 101, 103, 107, 109, 113 data 127, 131, 137, 139, 149, 151, 157, 163, 167, 173 data 179, 181, 191, 193, 197, 199, 211, 223, 227, 229 data 233, 239, 241, 251, 257, 263, 269, 271, 277, 281 data 283, 293, 307, 311, 313, 317, 331, 337, 347, 349 data 353, 359, 367, 373, 379, 383, 389, 397, 401, 409 data 419, 421, 431, 433, 439, 443, 449, 457, 461, 463 data 467, 479, 487, 491, 499, 503, 509, 521, 523, 541 data 547, 557, 563, 569, 571, 577, 587, 593, 599, 601 data 607, 613, 617, 619, 631, 641, 643, 647, 653, 659 data 661, 673, 677, 683, 691, 701, 709, 719, 727, 733 data 739, 743, 751, 757, 761, 769, 773, 787, 797, 809 data 811, 821, 823, 827, 829, 839, 853, 857, 859, 863 data 877, 881, 883, 887, 907, 911, 919, 929, 937, 941 data 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013 data 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069 data 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151 data 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223 data 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291 data 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373 data 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451 data 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511 data 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583 data 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657 data 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733 data 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811 data 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889 data 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987 data 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053 data 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129 data 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213 data 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287 data 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357 data 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423 data 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531 data 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617 data 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687 data 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741 data 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819 data 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903 data 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999 data 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079 data 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181 data 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257 data 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331 data 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413 data 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511 data 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571 data 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643 data 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727 data 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821 data 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907 data 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989 data 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057 data 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139 data 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231 data 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297 data 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409 data 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493 data 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583 data 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657 data 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751 data 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831 data 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937 data 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003 data 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087 data 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179 data 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279 data 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387 data 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443 data 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521 data 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639 data 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693 data 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791 data 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857 data 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939 data 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053 data 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133 data 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221 data 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301 data 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367 data 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473 data 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571 data 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673 data 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761 data 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833 data 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917 data 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997 data 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103 data 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207 data 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297 data 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411 data 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499 data 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561 data 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643 data 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723 data 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829 data 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919

print "Liberty Miller Rabin Demonstration" print "Loading Small Primes" for i=1 to 1000: read x : SmallPrimes(i)=x :next :NoOfSmallPrimes=1000 print NoOfSmallPrimes;" Primes Loaded"

'Prompt "Enter number to test:";resp$ 'x=val(resp$) 'goto [Jump]

For i=1 to 11

x=mersenne(i)

t1=time$("ms")
[TryAnother]
print

iterations=0
[Loop]
   iterations=iterations+1

   if MillerRabin(x,7)=1 then
    t2=time$("ms")
    print "Composite, found in ";t2-t1;" milliseconds"
   else
    t2=time$("ms")
    print x;" Probably Prime. Tested in ";t2-t1;" milliseconds"
    playwave "tada.wav", async
end if
print

Function GCD( m,n ) ' Find greatest common divisor with Extend Euclidian Algorithm ' Knuth Vol 1 P.13 Algorithm E

ap =1 :b =1 :a =0 :bp =0: c =m :d =n

[StepE2] q = int(c/d) :r = c-q*d

if r<>0 then

   c=d :d=r :t=ap :ap=a :a=t-q*a :t=bp  :bp=b  :b=t-q*b
   'print ap;" ";b;" ";a;" ";bp;" ";c;" ";d;" ";t;" ";q
   goto [StepE2]

end if

GCD=a*m+b*n

'print ap;" ";b;" ";a;" ";bp;" ";c;" ";d;" ";t;" ";q

End Function 'Extended Euclidian GCD

function IsEven( x )
   if ( x MOD 2 )=0 then
       IsEven=1
   else
        IsEven=0
   end if

end function

function IsOdd( x )

   if ( x MOD 2 )=0 then
       IsOdd=0
   else
       IsOdd=1
   end if

end function

Function FastExp(x, y, N)

 if (y=1) then                  'MOD(x,N)
     FastExp=x-int(x/N)*N
     goto [ExitFunction]
 end if

 if ( y and 1) = 0  then

    dum1=y/2
    dum2=y-int(y/2)*2              'MOD(y,2)

    temp=FastExp(x,dum1,N)
    z=temp*temp
    FastExp=z-int(z/N)*N            'MOD(temp*temp,N)
    goto [ExitFunction]
 else

    dum1=y-1
    dum1=dum1/2
    temp=FastExp(x,dum1,N)
    dum2=temp*temp
    temp=dum2-int(dum2/N)*N            'MOD(dum2,N)

    z=temp*x
    FastExp=z-int(z/N)*N             'MOD(temp*x,N)
    goto [ExitFunction]
 end if
 [ExitFunction]

end function

Function MillerRabin(n,b)

'print "Miller Rabin" 't1=time$("ms")

 if IsEven(n) then
   MillerRabin=1
   goto [ExtFn]
 end if

 i=0
 [Loop]
   i=i+1
   if i>1000 then goto [Continue]
   if ( n MOD SmallPrimes(i) )=0 then
     MillerRabin=0
     goto [ExtFn]
   end if
 goto [Loop]
 [Continue]

 if GCD(n,b)>1 then
   MillerRabin=1
   goto [ExtFn]
 end if

 q=n-1

t=0

 while  (int(q) AND 1 )=0
  t=t+1
  q=int(q/2)
 wend

 r=FastExp(b, q, n)

 if ( r <> 1 ) then
   e=0
   while ( e < (t-1) )
     if ( r <> (n-1) ) then
       r=FastExp(r, r, n)
       else
       Exit While
     end if

     e=e+1
   wend
   [ExitLoop]
 end if

 if ( (r=1) OR (r=(n-1)) ) then
     MillerRabin=0
   else
     MillerRabin=1
 end if

[ExtFn]

End Function </lang>

Mathematica

<lang Mathematica>MillerRabin[n_,k_]:=Module[{d=n-1,s=0,test=True},While[Mod[d,2]==0 ,d/=2 ;s++] Do[

 a=RandomInteger[{2,n-1}]; x=PowerMod[a,d,n];
 If[x!=1,
  For[ r = 0, r < s, r++, If[x==n-1, Continue[]]; x = Mod[x*x, n]; ];
  If[ x != n-1, test=False ];
 ];

,{k}]; Print[test] ]</lang>

Example output (not using the PrimeQ builtin):

<lang mathematica>MillerRabin[17388,10] ->False</lang>

Maxima

<lang maxima>/* Miller-Rabin algorithm is builtin, see function primep. Here is another implementation */

/* find highest power of p, p^s, that divide n, and return s and n / p^s */

facpow(n, p) := block(

  [s: 0],
  while mod(n, p) = 0 do (s: s + 1, n: quotient(n, p)),
  [s, n]

)$

/* check whether n is a strong pseudoprime to base a; s and d are given by facpow(n - 1, 2) */

sppp(n, a, s, d) := block(

  [x: power_mod(a, d, n), q: false],
  if x = 1 or x = n - 1 then true else (
     from 2 thru s do (
        x: mod(x * x, n),
        if x = 1 then return(q: false) elseif x = n - 1 then return(q: true)
     ),
     q
  )

)$

/* Miller-Rabin primality test. For n < 341550071728321, the test is deterministic;

  for larger n, the number of bases tested is given by the option variable
  primep_number_of_tests, which is used by Maxima in primep. The bound for deterministic
  test is also the same as in primep. */

miller_rabin(n) := block(

  [v: [2, 3, 5, 7, 11, 13, 17], s, d, q: true, a],
  if n < 19 then member(n, v) else (
     [s, d]: facpow(n - 1, 2),
     if n < 341550071728321 then (    /* see http://oeis.org/A014233 */
        for a in v do (
           if not sppp(n, a, s, d) then return(q: false)
        ),
        q
     ) else (
        thru primep_number_of_tests do (
           a: 2 + random(n - 3),
           if not sppp(n, a, s, d) then return(q: false)
        ),
        q
     )
  )

)$</lang>

PARI/GP

Built-in

<lang parigp>MR(n,k)=ispseudoprime(n,k);</lang>

Custom

<lang parigp>sprp(n,b)={ my(s = valuation(n-1, 2), d = Mod(b, n)^(n >> s)); if (d == 1, return(1)); for(i=1,s-1, if (d == -1, return(1)); d = d^2; ); d == -1 };

MR(n,k)={

 for(i=1,k,
   if(!sprp(n,random(n-2)+2), return(0))
 );
 1

};</lang>

Deterministic version

A basic deterministic test can be obtained by an appeal to the ERH (as proposed by Gary Miller) and a result of Eric Bach (improving on Joseph Oesterlé). Calculations of Jan Feitsma can be used to speed calculations below 2⁶⁴ (by a factor of about 250). <lang parigp>A006945=[9, 2047, 1373653, 25326001, 3215031751, 2152302898747, 3474749660383, 341550071728321, 341550071728321, 3825123056546413051]; Miller(n)={

 if (n%2 == 0, return(n == 2)); \\ Handle even numbers
 if (n < 3, return(0)); \\ Handle 0, 1, and negative numbers

 if (n < 1<<64,
   \\ Feitsma
   for(i=1,#A006945,
     if (n < A006945[i], return(1));
     if(!sprp(n, prime(i)), return(0));
   );
   sprp(n,31)&sprp(n,37)
 ,
   \\ Miller + Bach
   for(b=2,2*log(n)^2,
     if(!sprp(n, b), return(0))
   );
   1
 )

};</lang>

Perl

<lang 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->bmodpow($d, $n);
               next if $x == 1 or $x == $n-1;

               for(1..$s-1)
               {
                       $x = ($x*$x) % $n;
                       return 0 if $x == 1;
                       next LOOP if $x == $n-1;
               }
               return 0;
       }
       return 1;

}

print join ", ", grep { is_prime $_,10 }(1..1000);</lang>

Perl 6

Works with: Rakudo version 2011.11

<lang Perl6># the expmod-function from: http://rosettacode.org/wiki/Modular_exponentiation sub expmod(Int $a is copy, Int $b is copy, $n) { my $c = 1; repeat while $b div= 2 { ($c *= $a) %= $n if $b % 2; ($a *= $a) %= $n; } $c; }

subset PrimeCandidate of Int where { $_ > 2 and $_ % 2 };

my Bool multi sub is-prime(Int $n, Int $k) { return False; } my Bool multi sub is-prime(2, Int $k) { return True; } my Bool multi sub is-prime(PrimeCandidate $n, Int $k) { my Int $d = $n - 1; my Int $s = 0;

while $d %% 2 { $d div= 2; $s++; }

for (2 ..^ $n).pick($k) -> $a { my $x = expmod($a, $d, $n);

# one could just write "next if $x == 1 | $n - 1" # but this takes much more time in current rakudo/nom next if $x == 1 or $x == $n - 1;

for 1 ..^ $s { $x = $x ** 2 mod $n; return False if $x == 1; last if $x == $n - 1; } return False if $x !== $n - 1; }

return True; }

say (1..1000).grep({ is-prime($_, 10) }).join(", "); </lang>

PHP

<lang php><?php function is_prime($n, $k) {

   if ($n == 2)
       return true;
   if ($n < 2 || $n % 2 == 0)
       return false;

   $d = $n - 1;
   $s = 0;

   while ($d % 2 == 0) {
       $d /= 2;
       $s++;
   }

   for ($i = 0; $i < $k; $i++) {
       $a = rand(2, $n-1);

       $x = bcpowmod($a, $d, $n);
       if ($x == 1 || $x == $n-1)
           continue;

       for ($j = 1; $j < $s; $j++) {
           $x = bcmod(bcmul($x, $x), $n);
           if ($x == 1)
               return false;
           if ($x == $n-1)
               continue 2;
       }
       return false;
   }
   return true;

}

for ($i = 1; $i <= 1000; $i++)

   if (is_prime($i, 10))
       echo "$i, ";

echo "\n"; ?></lang>

PicoLisp

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

Output:

: (filter '((I) (prime? I)) (range 937 1000))
-> (937 941 947 953 967 971 977 983 991 997)

: (prime? 4547337172376300111955330758342147474062293202868155909489)
-> T

: (prime? 4547337172376300111955330758342147474062293202868155909393)
-> NIL

PureBasic

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

Python

Python: Probably correct answers

This versions will give answers with a very small probability of being false. That probability being dependent on _mrpt_num_trials and the random numbers used for name a passed to function try_composite.

<lang 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.randrange(2, n)
       if try_composite(a):
           return False

   return True # no base tested showed n as composite</lang>

Python: Proved correct up to large N

This versions will give correct answers for n less than 341550071728321 and then reverting to the probabilistic form of the first solution. By selecting a certain number of primes for name a instead of random values mathematicians have proved the general algorithm correct.
For 341550071728321 and beyond, I have followed the pattern in choosing a from the set of prime numbers.

<lang python>def _try_composite(a, d, n, s):

   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

def is_prime(n, _precision_for_huge_n=16):

   if n in _known_primes or n in (0, 1):
       return True
   if any((n % p) == 0 for p in _known_primes):
       return False
   d, s = n - 1, 0
   while not d % 2:
       d, s = d >> 1, s + 1
   # Returns exact according to http://primes.utm.edu/prove/prove2_3.html
   if n < 1373653: 
       return not any(_try_composite(a, d, n, s) for a in (2, 3))
   if n < 25326001: 
       return not any(_try_composite(a, d, n, s) for a in (2, 3, 5))
   if n < 118670087467: 
       if n == 3215031751: 
           return False
       return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7))
   if n < 2152302898747: 
       return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11))
   if n < 3474749660383: 
       return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13))
   if n < 341550071728321: 
       return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17))
   # otherwise
   return not any(_try_composite(a, d, n, s) 
                  for a in _known_primes[:_precision_for_huge_n])

_known_primes = [2, 3] _known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]</lang>

Testing

Includes test values from other examples:

>>> is_prime(4547337172376300111955330758342147474062293202868155909489)
True
>>> is_prime(4547337172376300111955330758342147474062293202868155909393)
False
>>> [x for x in range(901, 1000) if is_prime(x)]
[907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]
>>> is_prime(643808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153)
True
>>> is_prime(743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153)
False
>>>

Racket

<lang Racket>#lang racket (define (miller-rabin-expmod base exp m)

 (define (squaremod-with-check x)
   (define (check-nontrivial-sqrt1 x square)
     (if (and (= square 1)
              (not (= x 1))
              (not (= x (- m 1))))
         0
         square))
   (check-nontrivial-sqrt1 x (remainder (expt x 2) m)))
 (cond ((= exp 0) 1)
       ((even? exp) (squaremod-with-check
                     (miller-rabin-expmod base (/ exp 2) m)))
       (else
        (remainder (* base (miller-rabin-expmod base (- exp 1) m))
                   m))))

(define (miller-rabin-test n)

 (define (try-it a)
   (define (check-it x)
     (and (not (= x 0)) (= x 1)))
   (check-it (miller-rabin-expmod a (- n 1) n)))
 (try-it (+ 1 (random (remainder (- n 1) 4294967087)))))

(define (fast-prime? n times)

 (for/and ((i (in-range times)))
   (miller-rabin-test n)))

(define (prime? n(times 100))

 (fast-prime? n times))

(prime? 4547337172376300111955330758342147474062293202868155909489) ;-> outputs true </lang>

REXX

With a K of 1, there seems to be a not uncommon number of failures, but

with a K ≥ 2, the failures are rare,

with a K ≥ 3, rare as hen's teeth.

This would be in the same vein as: 3 is prime, 5 is prime, 7 is prime, all odd numbers are prime. <lang rexx>/*REXX program puts Miller-Rabin primality test through its paces. */

arg limit accur . /*get some arguments (if any). */ if limit== | limit==',' then limit=1000 /*maybe assume LIMIT default*/ if accur== | accur==',' then accur=10 /* " " ACCUR " */ numeric digits max(200,2*limit) /*we're dealing with some biggies*/ tell=accur<0 /*show primes if K is negative.*/ accur=abs(accur) /*now, make K positive. */ call suspenders /*suspenders now, belt later... */ primePi=# /*save the count of (real) primes*/ say "There are" primePi 'primes ≤' limit /*might as well crow a wee bit.*/ say /*nothing wrong with whitespace. */

     do a=2 to accur                  /*(skipping 1)  do range of  K's.*/
     say copies('─',79)               /*show separator for the eyeballs*/
     mrp=0                            /*prime counter for this pass.   */

       do z=1 for limit               /*now, let's get busy and crank. */
       p=Miller_Rabin(z,a)            /*invoke and pray...             */
       if p==0 then iterate           /*Not prime?   Then try another. */
       mrp=mrp+1                      /*well, found another one, by gum*/

       if tell then say z,            /*maybe should do a show & tell ?*/
             'is prime according to Miller-Rabin primality test with K='a

       if !.z\==0 then iterate
       say '[K='a"] " z "isn't prime !" /*oopsy-doopsy & whoopsy-daisy!*/
       end   /*z*/

     say 'for 1──►'limit", K="a', Miller-Rabin primality test found' mrp,
         'primes {out of' primePi"}"
     end     /*a*/

exit /*stick a fork in it, we're done.*/ /*─────────────────────────────────────Miller─Rabin primality test.─────*/ /*─────────────────────────────────────Rabin─Miller (also known as)─────*/ Miller_Rabin: procedure; parse arg n,k if n==2 then return 1 /*special case of an even prime. */ if n<2 | n//2==0 then return 0 /*check for low, or even number.*/ d=n-1 nL=n-1 /*saves a bit of time, down below*/ s=0

  do while d//2==0;   d=d%2;   s=s+1;   end    /*while d//2==0 */

         do k
         a=random(2,nL)
         x=(a**d) // n                /*this number can get big fast.  */
         if x==1 | x==nL then iterate

             do r=1  for s-1
             x=(x*x) // n
             if x==1  then return 0   /*it's definitely not prime.     */
             if x==nL then leave
             end   /*r*/

         if x\==nL then return 0      /*nope, it ain't prime nohows.   */
         end       /*k*/
                                      /*maybe it is, maybe it ain't ...*/

return 1 /*coulda/woulda/shoulda be prime.*/ /*──────────────────────────────────SUSPENDERS subroutine───────────────*/ suspenders: @.=0; !.=0 /*crank up the ole prime factory.*/ @.1=2; @.2=3; @.3=5; #=3 /*prime the pump with low primes.*/ !.2=1; !.3=1; !.5=1 /*and don't forget the water jar.*/

      do j  =@.#+2 by 2 to limit      /*just process the odd integers. */
          do k=2 while @.k**2<=j      /*let's do the ole primality test*/
          if j//@.k==0 then iterate j /*the Greek way, in days of yore.*/
          end   /*k*/                 /*a useless comment, but hey!!   */
      #=#+1                           /*bump the prime counter.        */
      @.#=j                           /*keep priming the prime pump.   */
      !.j=1                           /*and keep filling the water jar.*/
      end     /*j*/                   /*this comment not left blank.   */

return /*whew! All done with the primes*/</lang>

Output when using the input of: 10000 10:

There are 1229 primes ≤ 10000

───────────────────────────────────────────────────────────────────────────────
[K=2]  2701 isn't prime !
for 1──►10000, K=2, Miller─Rabin primality test found 1230 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=2, Miller─Rabin primality test found 1229 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=3, Miller─Rabin primality test found 1229 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=4, Miller─Rabin primality test found 1229 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=5, Miller─Rabin primality test found 1229 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=6, Miller─Rabin primality test found 1229 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=7, Miller─Rabin primality test found 1229 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=8, Miller─Rabin primality test found 1229 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=9, Miller─Rabin primality test found 1229 primes {out of 1229}
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=10, Miller─Rabin primality test found 1229 primes {out of 1229}

Ruby

<lang ruby> require 'openssl' def miller_rabin_prime?(n,g)

 d = n - 1
 s = 0
 while d % 2 == 0
   d /= 2
   s += 1
 end
 g.times do
   a = 2 + rand(n-4)
   x = OpenSSL::BN::new(a.to_s).mod_exp(d,n) #x = (a**d) % n
   next if x == 1 or x == n-1
   for r in (1 .. s-1)
     x = x.mod_exp(2,n)  #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 = (3..1000).step(2).find_all {|i| miller_rabin_prime?(i,10)} </lang>

Output:

[3, 5, 7, 11, 13, 17, ..., 971, 977, 983, 991, 997]

The following larger examples all produce true: <lang ruby> puts miller_rabin_prime?(94366396730334173383107353049414959521528815310548187030165936229578960209523421808912459795329035203510284576187160076386643700441216547732914250578934261891510827140267043592007225160798348913639472564715055445201512461359359488795427875530231001298552452230535485049737222714000227878890892901228389026881,1000) puts miller_rabin_prime?(138028649176899647846076023812164793645371887571371559091892986639999096471811910222267538577825033963552683101137782650479906670021895135954212738694784814783986671046107023185842481502719762055887490765764329237651328922972514308635045190654896041748716218441926626988737664133219271115413563418353821396401,1000) puts miller_rabin_prime?(123301261697053560451930527879636974557474268923771832437126939266601921428796348203611050423256894847735769138870460373141723679005090549101566289920247264982095246187318303659027201708559916949810035265951104246512008259674244307851578647894027803356820480862664695522389066327012330793517771435385653616841,1000) puts miller_rabin_prime?(119432521682023078841121052226157857003721669633106050345198988740042219728400958282159638484144822421840470442893056822510584029066504295892189315912923804894933736660559950053226576719285711831138657839435060908151231090715952576998400120335346005544083959311246562842277496260598128781581003807229557518839,1000) puts miller_rabin_prime?(132082885240291678440073580124226578272473600569147812319294626601995619845059779715619475871419551319029519794232989255381829366374647864619189704922722431776563860747714706040922215308646535910589305924065089149684429555813953571007126408164577035854428632242206880193165045777949624510896312005014225526731,1000) puts miller_rabin_prime?(153410708946188157980279532372610756837706984448408515364579602515073276538040155990230789600191915021209039203172105094957316552912585741177975853552299222501069267567888742458519569317286299134843250075228359900070009684517875782331709619287588451883575354340318132216817231993558066067063143257425853927599,1000) puts miller_rabin_prime?(103130593592068072608023213244858971741946977638988649427937324034014356815504971087381663169829571046157738503075005527471064224791270584831779395959349442093395294980019731027051356344056416276026592333932610954020105156667883269888206386119513058400355612571198438511950152690467372712488391425876725831041,1000) </lang>

Run BASIC

<lang runbasic>input "Input a number:";n input "Input test:";k

test = millerRabin(n,k) if test = 0 then

 print "Probably Prime"
else
 print "Composite"

end if wait

' ---------------------------------------- ' Returns ' Composite = 1 ' Probably Prime = 0 ' ----------------------------------------

FUNCTION millerRabin(n, k)
 if n = 2 then
   millerRabin = 0 'probablyPrime
   goto [funEnd]
 end if

 if n mod 2 = 0 or n < 2 then
   millerRabin = 1 'composite
   goto [funEnd]
 end if

d = n - 1 while d mod 2 = 0

 d = d / 2
 s = s + 1

wend

while k > 0

 k = k - 1
 x = (int(rnd(1) * (n-4))^d) mod n
 if x = 1 or x = n-1 then
   for r=1 To s-1
     x =(x * x) mod n
     if x=1 then
      millerRabin = 1 ' composite 
      goto [funEnd]
     end if
     if x = n-1 then exit for
   next r
   
   if x<>n-1 then
     millerRabin =  1 ' composite 
     goto [funEnd]
   end if
 end if

wend [funEnd] END FUNCTION</lang>

Seed7

<lang seed7>$ include "seed7_05.s7i";

 include "bigint.s7i";

const func boolean: millerRabin (in bigInteger: n, in integer: k) is func

 result
   var boolean: probablyPrime is TRUE;
 local
   var bigInteger: d is 0_;
   var integer: r is 0;
   var integer: s is 0;
   var bigInteger: a is 0_;
   var bigInteger: x is 0_;
   var integer: tests is 0;
 begin
   if n < 2_ or (n > 2_ and not odd(n)) then
     probablyPrime := FALSE;
   elsif n > 3_ then
     d := pred(n);
     s := lowestSetBit(d);
     d >>:= s;
     while tests < k and probablyPrime do
       a := rand(2_, pred(n));
       x := modPow(a, d, n);
       if x <> 1_ and x <> pred(n) then
         r := 1;
         while r < s and x <> 1_ and x <> pred(n) do
           x := modPow(x, 2_, n);
           incr(r);
         end while;
         probablyPrime := x = pred(n);
       end if;
       incr(tests);
     end while;
   end if;
 end func;

const proc: main is func

 local
   var bigInteger: number is 0_;
 begin
   for number range 2_ to 1000_ do
     if millerRabin(number, 10) then
       writeln(number);
     end if;
   end for;
 end func;</lang>

Original source: [2]

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>

Standard ML

<lang sml>open LargeInt;

val mr_iterations = Int.toLarge 20; val rng = Random.rand (557216670, 13504100); (* arbitrary pair to seed RNG *)

fun expmod base 0 m = 1

 | expmod base exp m =
     if exp mod 2 = 0
     then let val rt = expmod base (exp div 2) m;
              val sq = (rt*rt) mod m
          in if sq = 1
                andalso rt <> 1     (* ignore the two *)
                andalso rt <> (m-1) (* 'trivial' roots *)
             then 0
             else sq
          end
     else (base*(expmod base (exp-1) m)) mod m;

(* arbitrary precision random number [0,n) *) fun rand n =

 let val base = Int.toLarge(valOf Int.maxInt)+1;
     fun step r lim =
       if lim < n then step (Int.toLarge(Random.randNat rng) + r*base) (lim*base)
                  else r mod n
 in step 0 1 end;

fun miller_rabin n =

 let fun trial n 0 = true
       | trial n t = let val a = 1+rand(n-1)
                     in (expmod a (n-1) n) = 1
                        andalso trial n (t-1)
                     end
 in trial n mr_iterations end;

fun trylist label lst = (label, ListPair.zip (lst, map miller_rabin lst));

trylist "test the first six Carmichael numbers"

       [561, 1105, 1729, 2465, 2821, 6601];

trylist "test some known primes"

       [7369, 7393, 7411, 27367, 27397, 27407];

(* find ten random 30 digit primes (according to Miller-Rabin) *) let fun findPrime trials = let val t = trials+1;

                              val n = 2*rand(500000000000000000000000000000)+1
                          in if miller_rabin n
                             then (n,t)
                             else findPrime t end

in List.tabulate (10, fn e => findPrime 0) end;</lang>

Sample run:

...
val it =
  ("test the first six Carmichael numbers",
   [(561,false),(1105,false),(1729,false),(2465,false),(2821,false),
    (6601,false)]) : string * (int * bool) list
val it =
  ("test some known primes",
   [(7369,true),(7393,true),(7411,true),(27367,true),(27397,true),
    (27407,true)]) : string * (int * bool) list
[autoloading]
[autoloading done]
val it =
  [(505776511533674858497882481471,8),(668742242620107711631417930007,111),
   (831749124005136073184150011961,24),(159858916052323079037919394483,14),
   (810857757001516064878680795563,43),(903375242242638088171051457359,6),
   (506008872035764637556989600477,91),(105574439115200786396150347661,29),
   (349239056313926786302179212509,7),(565349019043144709861293116613,126)]
  : (int * int) list

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>

Output: