Miller-Rabin primality test

From Rosetta Code
Jump to: navigation, search
Task
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]
   xad mod n
   if x = 1 or x = n − 1 then do next LOOP
   for r = 1 .. s − 1
      xx2 mod n
      if x = 1 then return composite
      if x = n − 1 then do next LOOP
   return composite
return probably prime
  • The nature of the test involves big numbers, so the use of "big numbers" libraries (or similar features of the language of your choice) are suggested, but not mandatory.
  • Deterministic variants of the test exist and can be implemented as extra (not mandatory to complete the task)

Contents

[edit] Ada

[edit] 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 the Extensible prime generator, and the Emirp primes.

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;

The implementation of that package is as follows:

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;

Finally, the program itself:

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

[edit] using an external library to handle big integers

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

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;

Output:

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

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

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;

The output is the same.

[edit] ALGOL 68

Translation of: python
Works with: ALGOL 68 version Standard - with preludes manually inserted
Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386
MODE LINT=LONG INT;
MODE LOOPINT = INT;
 
MODE POWMODSTRUCT = LINT;
PR READ "prelude/pow_mod.a68" PR;
 
PROC miller rabin = (LINT n, LOOPINT k)BOOL: (
IF n<=3 THEN TRUE
ELIF NOT ODD n THEN FALSE
ELSE
LINT d := n - 1;
INT s := 0;
WHILE NOT ODD d DO
d := d OVER 2;
s +:= 1
OD;
TO k DO
LINT a := 2 + ENTIER (random*(n-3));
LINT x := pow mod(a, d, n);
IF x /= 1 THEN
TO s DO
IF x = n-1 THEN done FI;
x := x*x %* n
OD;
else: IF x /= n-1 THEN return false FI;
done: EMPTY
FI
OD;
TRUE EXIT
return false: FALSE
FI
);
 
FOR i FROM 937 TO 1000 DO
IF miller rabin(i, 10) THEN
print((" ",whole(i,0)))
FI
OD
Output:
 937 941 947 953 967 971 977 983 991 997

[edit] AutoHotkey

ahk forum: discussion

MsgBox % MillerRabin(999983,10) ; 1
MsgBox % MillerRabin(999809,10) ; 1
MsgBox % MillerRabin(999727,10) ; 1
MsgBox % MillerRabin(52633,10) ; 0
MsgBox % MillerRabin(60787,10) ; 0
MsgBox % MillerRabin(999999,10) ; 0
MsgBox % MillerRabin(999995,10) ; 0
MsgBox % MillerRabin(999991,10) ; 0
 
MillerRabin(n,k) { ; 0: composite, 1: probable prime (n < 2**31)
d := n-1, s := 0
While !(d&1)
d>>=1, s++
 
Loop %k% {
Random a, 2, n-2 ; if n < 4,759,123,141, it is enough to test a = 2, 7, and 61.
x := PowMod(a,d,n)
If (x=1 || x=n-1)
Continue
Cont := 0
Loop % s-1 {
x := PowMod(x,2,n)
If (x = 1)
Return 0
If (x = n-1) {
Cont = 1
Break
}
}
IfEqual Cont,1, Continue
Return 0
}
Return 1
}
 
PowMod(x,n,m) { ; x**n mod m
y := 1, i := n, z := x
While i>0
y := i&1 ? mod(y*z,m) : y, z := mod(z*z,m), i >>= 1
Return y
}

[edit] bc

Requires a bc with long names.

Works with: OpenBSD bc

(A previous version worked with GNU 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

[edit] Bracmat

Translation of: bc
( 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
);

output:

937 941 947 953 967 971 977 983 991 997

[edit] C

Library: GMP

miller-rabin.h

#ifndef _MILLER_RABIN_H_
#define _MILLER_RABIN_H
#include <gmp.h>
bool miller_rabin_test(mpz_t n, int j);
#endif

miller-rabin.c

Translation of: Fortran

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

#include <stdbool.h>
#include <gmp.h>
#include "primedecompose.h"
 
#define MAX_DECOMPOSE 100
 
bool miller_rabin_test(mpz_t n, int j)
{
bool res;
mpz_t f[MAX_DECOMPOSE];
mpz_t s, d, a, x, r;
mpz_t n_1, n_3;
gmp_randstate_t rs;
int l=0, k;
 
res = false;
gmp_randinit_default(rs);
 
mpz_init(s); mpz_init(d);
mpz_init(a); mpz_init(x); mpz_init(r);
mpz_init(n_1); mpz_init(n_3);
 
if ( mpz_cmp_si(n, 3) <= 0 ) { // let us consider 1, 2, 3 as prime
gmp_randclear(rs);
return true;
}
if ( mpz_odd_p(n) != 0 ) {
mpz_sub_ui(n_1, n, 1); // n-1
mpz_sub_ui(n_3, n, 3); // n-3
l = decompose(n_1, f);
mpz_set_ui(s, 0);
mpz_set_ui(d, 1);
for(k=0; k < l; k++) {
if ( mpz_cmp_ui(f[k], 2) == 0 )
mpz_add_ui(s, s, 1);
else
mpz_mul(d, d, f[k]);
} // 2^s * d = n-1
while(j-- > 0) {
mpz_urandomm(a, rs, n_3); // random from 0 to n-4
mpz_add_ui(a, a, 2); // random from 2 to n-2
mpz_powm(x, a, d, n);
if ( mpz_cmp_ui(x, 1) == 0 ) continue;
mpz_set_ui(r, 0);
while( mpz_cmp(r, s) < 0 ) {
if ( mpz_cmp(x, n_1) == 0 ) break;
mpz_powm_ui(x, x, 2, n);
mpz_add_ui(r, r, 1);
}
if ( mpz_cmp(x, n_1) == 0 ) continue;
goto flush; // woops
}
res = true;
}
 
flush:
for(k=0; k < l; k++) mpz_clear(f[k]);
mpz_clear(s); mpz_clear(d);
mpz_clear(a); mpz_clear(x); mpz_clear(r);
mpz_clear(n_1); mpz_clear(n_3);
gmp_randclear(rs);
return res;
}

Testing

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <gmp.h>
#include "miller-rabin.h"
 
#define PREC 10
#define TOP 4000
 
int main()
{
mpz_t num;
 
mpz_init(num);
mpz_set_ui(num, 1);
 
while ( mpz_cmp_ui(num, TOP) < 0 ) {
if ( miller_rabin_test(num, PREC) ) {
gmp_printf("%Zd maybe prime\n", num);
} /*else {
gmp_printf("%Zd not prime\n", num);
}*/
// remove the comment iff you're interested in
// sure non-prime.
mpz_add_ui(num, num, 1);
}
 
mpz_clear(num);
return EXIT_SUCCESS;
}


[edit] Deterministic up to 341,550,071,728,321

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

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


[edit] C#

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

[edit] Common Lisp

(defun factor-out (number divisor)
"Return two values R and E such that NUMBER = DIVISOR^E * R,
and R is not divisible by DIVISOR."

(do ((e 0 (1+ e))
(r number (/ r divisor)))
((/= (mod r divisor) 0) (values r e))))
 
(defun mult-mod (x y modulus) (mod (* x y) modulus))
 
(defun expt-mod (base exponent modulus)
"Fast modular exponentiation by repeated squaring."
(labels ((expt-mod-iter (b e p)
(cond ((= e 0) p)
((evenp e)
(expt-mod-iter (mult-mod b b modulus)
(/ e 2)
p))
(t
(expt-mod-iter b
(1- e)
(mult-mod b p modulus))))))
(expt-mod-iter base exponent 1)))
 
(defun random-in-range (lower upper)
"Return a random integer from the range [lower..upper]."
(+ lower (random (+ (- upper lower) 1))))
 
(defun miller-rabin-test (n k)
"Test N for primality by performing the Miller-Rabin test K times.
Return NIL if N is composite, and T if N is probably prime."

(cond ((= n 1) nil)
((< n 4) t)
((evenp n) nil)
(t
(multiple-value-bind (d s) (factor-out (- n 1) 2)
(labels ((strong-liar? (a)
(let ((x (expt-mod a d n)))
(or (= x 1)
(loop repeat s
for y = x then (mult-mod y y n)
thereis (= y (- n 1)))))))
(loop repeat k
always (strong-liar? (random-in-range 2 (- n 2)))))))))
CL-USER> (last (loop for i from 1 to 1000
                     when (miller-rabin-test i 10)
                     collect i)
               10)
(937 941 947 953 967 971 977 983 991 997)

[edit] D

Translation of: Ruby
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)));
}
Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

[edit] E

def millerRabinPrimalityTest(n :(int > 0), k :int, random) :boolean {
if (n <=> 2 || n <=> 3) { return true }
if (n <=> 1 || n %% 2 <=> 0) { return false }
var d := n - 1
var s := 0
while (d %% 2 <=> 0) {
d //= 2
s += 1
}
for _ in 1..k {
def nextTrial := __continue
def a := random.nextInt(n - 3) + 2 # [2, n - 2] = [0, n - 4] + 2 = [0, n - 3) + 2
var x := a**d %% n # Note: Will do optimized modular exponentiation
if (x <=> 1 || x <=> n - 1) { nextTrial() }
for _ in 1 .. (s - 1) {
x := x**2 %% n
if (x <=> 1) { return false }
if (x <=> n - 1) { nextTrial() }
}
return false
}
return true
}
for i ? (millerRabinPrimalityTest(i, 1, entropy)) in 4..1000 {
print(i, " ")
}
println()

[edit] Erlang

-module(miller_rabin).
 
-export([is_prime/1, power/2]).
 
is_prime(1) -> false;
is_prime(2) -> true;
is_prime(3) -> true;
is_prime(N) when N > 3, ((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),
%% this is a test for compositeness; the two case patterns disprove
%% compositeness.
not lists:any(fun(A) ->
case mr_series(N, A, D, S) of
[1|_] -> false; % first elem of list = 1
L -> not lists:member(N-1, L) % some elem of list = N-1
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-3). % random:uniform returns a single random number in range 1 -> N-3, to which is added 1, shifting the range to 2 -> N-2
 
 
power(B, E) ->
power(B, E, 1).
 
power(_, 0, Acc) ->
Acc;
power(B, E, Acc) ->
power(B, E - 1, B * Acc).

[edit] Fortran

Works with: Fortran version 95

For the module PrimeDecompose, see Prime decomposition.

 
module Miller_Rabin
use PrimeDecompose
implicit none
 
integer, parameter :: max_decompose = 100
 
private :: int_rrand, max_decompose
 
contains
 
function int_rrand(from, to)
integer(huge) :: int_rrand
integer(huge), intent(in) :: from, to
 
real :: o
call random_number(o)
int_rrand = floor(from + o * real(max(from,to) - min(from, to)))
end function int_rrand
 
function miller_rabin_test(n, k) result(res)
logical :: res
integer(huge), intent(in) :: n
integer, intent(in) :: k
 
integer(huge), dimension(max_decompose) :: f
integer(huge) :: s, d, i, a, x, r
 
res = .true.
f = 0
 
if ( (n <= 2) .and. (n > 0) ) return
if ( mod(n, 2) == 0 ) then
res = .false.
return
end if
 
call find_factors(n-1, f)
s = count(f == 2)
d = (n-1) / (2 ** s)
loop: do i = 1, k
a = int_rrand(2_huge, n-2)
x = mod(a ** d, n)
 
if ( x == 1 ) cycle
do r = 0, s-1
if ( x == ( n - 1 ) ) cycle loop
x = mod(x*x, n)
end do
if ( x == (n-1) ) cycle
res = .false.
return
end do loop
res = .true.
end function miller_rabin_test
 
end module Miller_Rabin

Testing

program TestMiller
use Miller_Rabin
implicit none
 
integer, parameter :: prec = 30
integer(huge) :: i
 
! this is limited since we're not using a bignum lib
call do_test( (/ (i, i=1, 29) /) )
 
contains
 
subroutine do_test(a)
integer(huge), dimension(:), intent(in) :: a
 
integer :: i
 
do i = 1, size(a,1)
print *, a(i), miller_rabin_test(a(i), prec)
end do
 
end subroutine do_test
 
end program TestMiller

Possible improvements: create bindings to the GMP library, change integer(huge) into something like type(huge_integer), write a lots of interfaces to allow to use big nums naturally (so that the code will be unchanged, except for the changes said above)

[edit] FunL

Direct implementation of the task algorithm.

import util.rnd
 
def isProbablyPrimeMillerRabin( n, k ) =
d = n - 1
s = 0
 
while 2|d
s++
d /= 2
 
repeat k
a = rnd( 2, n )
x = a^d mod n
 
if x == 1 or x == n - 1 then continue
 
repeat s - 1
x = x^2 mod n
 
if x == 1 then return false
 
if x == n - 1 then break
else
return false
 
true
 
for i <- 3..100
if isProbablyPrimeMillerRabin( i, 5 )
println( i )
Output:
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97

[edit] 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.

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
}

[edit] Haskell

Works with: Haskell version 7.6.3

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

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

[edit] Icon and Unicon

The following runs in both languages:

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

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

[edit] J

See Primality Tests essay on the J wiki.

[edit] Java

The Miller-Rabin primality test is part of the standard library (java.math.BigInteger)

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"));
}
}
Sample output:
java MillerRabinPrimalityTest 123456791234567891234567 1000000
123456791234567891234567 is probably prime

This is a translation of the Python solution for a deterministic test for n < 341,550,071,728,321:

import java.math.BigInteger;
 
public class Prime {
 
// this is the RabinMiller test, deterministically correct for n < 341,550,071,728,321
// http://rosettacode.org/wiki/Miller-Rabin_primality_test#Python:_Proved_correct_up_to_large_N
public static boolean isPrime(BigInteger n, int precision) {
 
if (n.compareTo(new BigInteger("341550071728321")) >= 0) {
return n.isProbablePrime(precision);
}
 
int intN = n.intValue();
if (intN == 1 || intN == 4 || intN == 6 || intN == 8) return false;
if (intN == 2 || intN == 3 || intN == 5 || intN == 7) return true;
 
int[] primesToTest = getPrimesToTest(n);
if (n.equals(new BigInteger("3215031751"))) {
return false;
}
BigInteger d = n.subtract(BigInteger.ONE);
BigInteger s = BigInteger.ZERO;
while (d.mod(BigInteger.valueOf(2)).equals(BigInteger.ZERO)) {
d = d.shiftRight(1);
s = s.add(BigInteger.ONE);
}
for (int a : primesToTest) {
if (try_composite(a, d, n, s)) {
return false;
}
}
return true;
}
 
public static boolean isPrime(BigInteger n) {
return isPrime(n, 100);
}
 
public static boolean isPrime(int n) {
return isPrime(BigInteger.valueOf(n), 100);
}
 
public static boolean isPrime(long n) {
return isPrime(BigInteger.valueOf(n), 100);
}
 
private static int[] getPrimesToTest(BigInteger n) {
if (n.compareTo(new BigInteger("3474749660383")) >= 0) {
return new int[]{2, 3, 5, 7, 11, 13, 17};
}
if (n.compareTo(new BigInteger("2152302898747")) >= 0) {
return new int[]{2, 3, 5, 7, 11, 13};
}
if (n.compareTo(new BigInteger("118670087467")) >= 0) {
return new int[]{2, 3, 5, 7, 11};
}
if (n.compareTo(new BigInteger("25326001")) >= 0) {
return new int[]{2, 3, 5, 7};
}
if (n.compareTo(new BigInteger("1373653")) >= 0) {
return new int[]{2, 3, 5};
}
return new int[]{2, 3};
}
 
private static boolean try_composite(int a, BigInteger d, BigInteger n, BigInteger s) {
BigInteger aB = BigInteger.valueOf(a);
if (aB.modPow(d, n).equals(BigInteger.ONE)) {
return false;
}
for (int i = 0; BigInteger.valueOf(i).compareTo(s) < 0; i++) {
// if pow(a, 2**i * d, n) == n-1
if (aB.modPow(BigInteger.valueOf(2).pow(i).multiply(d), n).equals(n.subtract(BigInteger.ONE))) {
return false;
}
}
return true;
}
}
 

[edit] JavaScript

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

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

[edit] 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):

 
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
 

[edit] Liberty BASIC

 
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
 
next
 
END
 
 
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
 

[edit] 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] ]
Example output (not using the PrimeQ builtin):
MillerRabin[17388,10] 
->False

[edit] 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
)
)
)$

[edit] PARI/GP

[edit] Built-in

MR(n,k)=ispseudoprime(n,k);

[edit] Custom

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
};

[edit] 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 264 (by a factor of about 250).

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

[edit] Perl

use bigint;
sub is_prime
{
my ($n,$k) = @_;
return 1 if $n == 2;
return 0 if $n < 2 or $n % 2 == 0;
 
$d = $n - 1;
$s = 0;
 
while(!($d % 2))
{
$d /= 2;
$s++;
}
 
LOOP: for(1..$k)
{
$a = 2 + int(rand($n-2));
 
$x = $a->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);

[edit] Perl 6

Works with: Rakudo version 2011.11
# 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(", ");

[edit] 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";
?>

[edit] PicoLisp

(de longRand (N)
(use (R D)
(while (=0 (setq R (abs (rand)))))
(until (> R N)
(unless (=0 (setq D (abs (rand))))
(setq R (* R D)) ) )
(% R N) ) )
 
(de **Mod (X Y N)
(let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )
 
(de _prim? (N D S)
(use (A X R)
(while (> 2 (setq A (longRand N))))
(setq R 0 X (**Mod A D N))
(loop
(T
(or
(and (=0 R) (= 1 X))
(= X (dec N)) )
T )
(T
(or
(and (> R 0) (= 1 X))
(>= (inc 'R) S) )
NIL )
(setq X (% (* X X) N)) ) ) )
 
(de prime? (N K)
(default K 50)
(and
(> N 1)
(bit? 1 N)
(let (D (dec N) S 0)
(until (bit? 1 D)
(setq
D (>> 1 D)
S (inc S) ) )
(do K
(NIL (_prim? N D S))
T ) ) ) )
Output:
: (filter '((I) (prime? I)) (range 937 1000))
-> (937 941 947 953 967 971 977 983 991 997)

: (prime? 4547337172376300111955330758342147474062293202868155909489)
-> T

: (prime? 4547337172376300111955330758342147474062293202868155909393)
-> NIL

[edit] Prolog

:- module(primality, [is_prime/2]).
 
% is_prime/2 returns false if N is composite, true if N probably prime
% implements a Miller-Rabin primality test and is deterministic for N < 3.415e+14,
% and is probabilistic for larger N. Adapted from the Erlang version.
is_prime(1, Ret) :- Ret = false, !. % 1 is non-prime
is_prime(2, Ret) :- Ret = true, !. % 2 is prime
is_prime(3, Ret) :- Ret = true, !. % 3 is prime
is_prime(N, Ret) :-
N > 3, (N mod 2 =:= 0), Ret = false, !. % even number > 3 is composite
is_prime(N, Ret) :-
N > 3, (N mod 2 =:= 1), % odd number > 3
N < 341550071728321,
deterministic_witnesses(N, L),
is_mr_prime(N, L, Ret), !. % deterministic test
is_prime(N, Ret) :-
random_witnesses(N, 100, [], Out),
is_mr_prime(N, Out, Ret), !. % probabilistic test
 
% returns list of deterministic witnesses
deterministic_witnesses(N, L) :- N < 1373653,
L = [2, 3].
deterministic_witnesses(N, L) :- N < 9080191,
L = [31, 73].
deterministic_witnesses(N, L) :- N < 25326001,
L = [2, 3, 5].
deterministic_witnesses(N, L) :- N < 3215031751,
L = [2, 3, 5, 7].
deterministic_witnesses(N, L) :- N < 4759123141,
L = [2, 7, 61].
deterministic_witnesses(N, L) :- N < 1122004669633,
L = [2, 13, 23, 1662803].
deterministic_witnesses(N, L) :- N < 2152302898747,
L = [2, 3, 5, 7, 11].
deterministic_witnesses(N, L) :- N < 3474749660383,
L = [2, 3, 5, 7, 11, 13].
deterministic_witnesses(N, L) :- N < 341550071728321,
L = [2, 3, 5, 7, 11, 13, 17].
 
% random_witnesses/4 returns a list of K witnesses selected at random with range 2 -> N-2
random_witnesses(_, 0, T, T).
random_witnesses(N, K, T, Out) :-
G is N - 2,
H is 1 + random(G),
I is K - 1,
random_witnesses(N, I, [H | T], Out), !.
 
% find_ds/2 receives odd integer N and returns [D, S] s.t. N-1 = 2^S * D
find_ds(N, L) :-
A is N - 1,
find_ds(A, 0, L), !.
 
find_ds(D, S, L) :-
D mod 2 =:= 0,
P is D // 2,
Q is S + 1,
find_ds(P, Q, L), !.
find_ds(D, S, L) :-
L = [D, S].
 
is_mr_prime(N, As, Ret) :-
find_ds(N, L),
L = [D | T],
T = [S | _],
outer_loop(N, As, D, S, Ret), !.
 
outer_loop(N, As, D, S, Ret) :-
As = [A | At],
Base is powm(A, D, N),
inner_loop(Base, N, 0, S, Result),
( Result == false -> Ret = false
; Result == true, At == [] -> Ret = true
; outer_loop(N, At, D, S, Ret)
).
 
inner_loop(Base, N, Loop, S, Result) :-
Next_Base is (Base * Base) mod N,
Next_Loop is Loop + 1,
( Loop =:= 0, Base =:= 1 -> Result = true
; Base =:= N-1 -> Result = true
; Next_Loop =:= S -> Result = false
; inner_loop(Next_Base, N, Next_Loop, S, Result)
).

[edit] PureBasic

Enumeration 
#Composite
#Probably_prime
EndEnumeration
 
Procedure Miller_Rabin(n, k)
Protected d=n-1, s, x, r
If n=2
ProcedureReturn #Probably_prime
ElseIf n%2=0 Or n<2
ProcedureReturn #Composite
EndIf
While d%2=0
d/2
s+1
Wend
While k>0
k-1
x=Int(Pow(2+Random(n-4),d))%n
If x=1 Or x=n-1: Continue: EndIf
For r=1 To s-1
x=(x*x)%n
If x=1: ProcedureReturn #Composite: EndIf
If x=n-1: Break: EndIf
Next
If x<>n-1: ProcedureReturn #Composite: EndIf
Wend
ProcedureReturn #Probably_prime
EndProcedure

[edit] Python

The correctness of this code is in question! See discussion

[edit] 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.

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

[edit] 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.

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

[edit] 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
 

[edit] 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.

/*REXX program puts the  Miller-Rabin primality test  through its paces.*/
parse arg limit accur . /*get some optional arguments*/
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, use absolute value of K.*/
call suspenders limit /*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 /*just a temp variable for speed.*/
nL=n-1 /*saves a bit of time, down below*/
s=0 /*teeter-toter variable (see-saw)*/
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 /*if so, try another power of A. */
do r=1 for s-1; x=(x*x)//n /*compute new X ≡ X² modulus 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: parse arg high; @.=0; !.=0 /*crank up the ole prime factory.*/
_='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'
do #=1 for words(_); p=word(_,#); @.#=p; s.#=p*p;  !.p=1; end /*#*/
#=#-1 /*adjust the # of primes so far. */
do j=@.#+2 by 2 to high /*just process the odd integers. */
if j//3==0 then iterate /*divisible by three? Yes, ¬prime*/
if right(j,1)==5 then iterate /* " " five? " " */
do k=4 while s.k<=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; @.#=j; s.#=j*j; !.j=1 /*bump prime ctr, prime, primeidx*/
end /*j*/ /*this comment not left blank. */
return /*whew! All done with the primes*/

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}

[edit] 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)}
 
Output:
[3, 5, 7, 11, 13, 17, ..., 971, 977, 983, 991, 997]

The following larger examples all produce true:

 
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)
 

[edit] Run BASIC

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

[edit] Scala

Library: Scala
import scala.math.BigInt
 
object MillerRabinPrimalityTest extends App {
val (n, certainty )= (BigInt(args(0)), args(1).toInt)
println(s"$n is ${if (n.isProbablePrime(certainty)) "probably prime" else "composite"}")
}

[edit] 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;
Original source: [2]

[edit] Smalltalk

Works with: GNU Smalltalk

Smalltalk handles big numbers naturally and trasparently (the parent class Integer has many subclasses, and a subclass is picked according to the size of the integer that must be handled)

Integer extend [
millerRabinTest: kl [ |k| k := kl.
self <= 3
ifTrue: [ ^true ]
ifFalse: [
(self even)
ifTrue: [ ^false ]
ifFalse: [ |d s|
d := self - 1.
s := 0.
[ (d rem: 2) == 0 ]
whileTrue: [
d := d / 2.
s := s + 1.
].
[ k:=k-1. k >= 0 ]
whileTrue: [ |a x r|
a := Random between: 2 and: (self - 2).
x := (a raisedTo: d) rem: self.
( x = 1 )
ifFalse: [ |r|
r := -1.
[ r := r + 1. (r < s) & (x ~= (self - 1)) ]
whileTrue: [
x := (x raisedTo: 2) rem: self
].
( x ~= (self - 1) ) ifTrue: [ ^false ]
]
].
^true
]
]
]
].
1 to: 1000 do: [ :n | 
(n millerRabinTest: 10) ifTrue: [ n printNl ]
].

[edit] Standard ML

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

[edit] Tcl

Use Tcl 8.5 for large integer support

package require Tcl 8.5
 
proc miller_rabin {n k} {
if {$n <= 3} {return true}
if {$n % 2 == 0} {return false}
 
# write n - 1 as 2^s·d with d odd by factoring powers of 2 from n − 1
set d [expr {$n - 1}]
set s 0
while {$d % 2 == 0} {
set d [expr {$d / 2}]
incr s
}
 
while {$k > 0} {
incr k -1
set a [expr {2 + int(rand()*($n - 4))}]
set x [expr {($a ** $d) % $n}]
if {$x == 1 || $x == $n - 1} continue
for {set r 1} {$r < $s} {incr r} {
set x [expr {($x ** 2) % $n}]
if {$x == 1} {return false}
if {$x == $n - 1} break
}
if {$x != $n-1} {return false}
}
return true
}
 
for {set i 1} {$i < 1000} {incr i} {
if {[miller_rabin $i 10]} {
puts $i
}
}
Output:
1
2
3
5
7
11
...
977
983
991
997

[edit] zkl

Using the Miller-Rabin primality test in GMP:

zkl: var BN=Import("zklBigNum");
zkl: BN("4547337172376300111955330758342147474062293202868155909489").probablyPrime()
True
zkl: BN("4547337172376300111955330758342147474062293202868155909393").probablyPrime()
False
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox