Miller–Rabin primality test: Difference between revisions
Capra Hircus (talk | contribs) |
Capra Hircus (talk | contribs) |
||
| Line 1,215: | Line 1,215: | ||
miller_rabin(n) := block( |
miller_rabin(n) := block( |
||
[v: [2, 3, 5, 7, 11, 13, 17], s, d, q: true], |
[v: [2, 3, 5, 7, 11, 13, 17], s, d, q: true, a], |
||
if n < 19 then member(n, v) else ( |
if n < 19 then member(n, v) else ( |
||
[s, d]: facpow(n - 1, 2), |
[s, d]: facpow(n - 1, 2), |
||
Revision as of 19:00, 25 August 2012
You are encouraged to solve this task according to the task description, using any language you may know.
| This page uses content from Wikipedia. The original article was at Miller–Rabin primality test. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance) |
The Miller–Rabin primality test or Rabin–Miller primality test is a primality test: an algorithm which determines whether a given number is prime or not. The algorithm, as modified by Michael O. Rabin to avoid the generalized Riemann hypothesis, is a probabilistic algorithm.
The pseudocode, from Wikipedia is:
Input: n > 2, an odd integer to be tested for primality;
k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
write n − 1 as 2s·d with d odd by factoring powers of 2 from n − 1
LOOP: repeat k times:
pick a randomly in the range [2, n − 1]
x ← ad mod n
if x = 1 or x = n − 1 then do next LOOP
for r = 1 .. s − 1
x ← x2 mod n
if x = 1 then return composite
if x = n − 1 then do next LOOP
return composite
return probably prime
- The nature of the test involves big numbers, so the use of "big numbers" libraries (or similar features of the language of your choice) are suggested, but not mandatory.
- Deterministic variants of the test exist and can be implemented as extra (not mandatory to complete the task)
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. <lang Ada>with Ada.Text_IO; with Ada.Numerics.Discrete_Random;
procedure Miller_Rabin is
-- New Type for Number, for easy swapping with bigger type if necessary. -- Limit to 2**48 for now, since Discrete_Random in GNAT GPL 2010 doesn't -- guarantee statistical properties for size > 48. -- RM requires them to be guaranteed only up to 2**15. type Number is mod 2**48;
package Num_IO is new Ada.Text_IO.Modular_IO (Number); package Pos_IO is new Ada.Text_IO.Integer_IO (Positive);
type Result_Type is (Composite, Probably_Prime);
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;
-- start main procedure N : Number; K : Positive;
begin
for I in 2 .. 1000 loop
if Is_Prime (Number(I)) = Probably_Prime then
Ada.Text_IO.Put (Integer'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 Miller_Rabin;</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
<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.
(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
<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
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
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>
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
<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
<lang erlang>-module(miller_rabin).
-compile([export_all]).
basis(N) when N>2 ->
1 + random:uniform(N-2).
find_ds(D, S) ->
case D rem 2 == 0 of
true ->
find_ds(trunc(D/2), S+1);
false ->
{D, S}
end.
find_ds(N) ->
find_ds(N-1, 0).
pow_mod(B, E, M) ->
case E of
0 -> 1;
_ -> case trunc(E) rem 2 == 0 of
true -> trunc(math:pow(pow_mod(B, trunc(E/2), M), 2)) rem M;
false -> trunc(B*pow_mod(B, E-1, M)) rem M
end
end.
mr_series(N, A, D, S) when N rem 2 == 1 ->
Js = lists:seq(0, S), lists:map(fun(J) -> pow_mod(A, math:pow(2, J)*D, N) end, Js).
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).
proving_bases(N) when N < 1373653 ->
[2, 3];
proving_bases(N) when N < 25326001 ->
[2, 3, 5];
proving_bases(N) when N < 25000000000 ->
[2, 3, 5, 7];
proving_bases(N) when N < 2152302898747->
[2, 3, 5, 7, 11];
proving_bases(N) when N < 341550071728321 ->
[2, 3, 5, 7, 11, 13];
proving_bases(N) when N < 341550071728321 ->
[2, 3, 5, 7, 11, 13, 17].
random_bases(N, K) ->
[basis(N) || _ <- lists:seq(1, K)].
is_prime(1) -> false; is_prime(2) -> true; is_prime(N) when N rem 2 == 0 -> false; is_prime(N) when N < 341550071728321 ->
is_mr_prime(N, proving_bases(N)).
is_probable_prime(N) ->
is_mr_prime(N, random_bases(N, 20)).
first_1000() ->
L = lists:seq(1,1000),
lists:map(fun(X) ->
case is_prime(X) of
true ->
io:format("~w~n", [X]);
false ->
false
end
end,
L),
ok.</lang>
- Output:
42> miller_rabin:first_1000(). 2 5 7 11 13 17 ... 967 971 977 983 991 997 ok
Fortran
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
- 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
Another Miller Rabin test can be found in D. Amos's Haskell for Math module Primes.hs <lang Haskell>import System.Random import Data.List import Control.Monad import Control.Arrow
primesTo100 = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
powerMod :: (Integral a, Integral b) => a -> a -> b -> a powerMod m _ 0 = 1 powerMod m x n | n > 0 = join (flip f (n - 1)) x `rem` m where
f _ 0 y = y
f a d y = g a d where
g b i | even i = g (b*b `rem` m) (i `quot` 2)
| otherwise = f b (i-1) (b*y `rem` m)
witns :: (Num a, Ord a, Random a) => Int -> a -> IO [a] witns x y = do
g <- newStdGen
let r = [9080191, 4759123141, 2152302898747, 3474749600383, 341550071728321]
fs = [[31,73],[2,7,61],[2,3,5,7,11],[2,3,5,7,11,13],[2,3,5,7,11,13,17]]
if y >= 341550071728321
then return $ take x $ randomRs (2,y-1) g
else return $ snd.head.dropWhile ((<= y).fst) $ zip r fs
isMillerRabinPrime :: Integer -> IO Bool isMillerRabinPrime n | n `elem` primesTo100 = return True
| otherwise = do
let pn = pred n
e = uncurry (++) . second(take 1) . span even . iterate (`div` 2) $ pn
try = return . all (\a -> let c = map (powerMod n a) e in
pn `elem` c || last c == 1)
witns 100 n >>= try</lang>
Testing in GHCi:
*~> isMillerRabinPrime 4547337172376300111955330758342147474062293202868155909489 True *~> isMillerRabinPrime 4547337172376300111955330758342147474062293202868155909393 False *~> filterM isMillerRabinPrime [2..1000]>>= print. dropWhile(<=900) [907,911,919,929,937,941,946,947,953,967,971,977,983,991,997]
Icon and Unicon
The following code works in both languages: <lang unicon>procedure main()
every writes(primeTest(901 to 1000, 10)," ") write()
end
procedure primeTest(n, k)
if n = 2 then return n if n%2 = 0 then fail s := 0 d := n-1 while (d%2 ~= 0, s+:=1, d/:=2)
every (1 to k, x := ((2+?(n-2))^d)%n) do {
if x = (1 | (n-1)) then next
every (1 to s-1, x := (x^2)%n) do {
if x = 1 then fail
if x = n-1 then break next
}
fail
}
return n
end</lang>
- Sample run:
->mrpt 907 911 919 929 937 941 947 953 967 971 977 983 991 997 ->
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
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 264 (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
<lang Perl6>use v6;
- the modexp-function from: http://rosettacode.org/wiki/Modular_exponentiation
sub modexp(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 = modexp($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
<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>
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 postive. */ call suspenders /*suspenders now, belt later... */ primePi=# /*save the count of (real) primes*/ say "They're" 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 definately 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:
They're 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>def miller_rabin_prime?(n,k)
return true if n == 2
return false if n < 2 or n % 2 == 0
d = n - 1
s = 0
while d % 2 == 0
d /= 2
s += 1
end
k.times do
a = 2 + rand(n-4)
x = (a**d) % n
next if x == 1 or x == n-1
for r in (1 .. s-1)
x = (x**2) % n
return false if x == 1
break if x == n-1
end
return false if x != n-1
end
true # probably
end
p primes = (1..100).find_all {|i| miller_rabin_prime?(i,10)}</lang>
- Output:
[2, 3, 5, 7, 11, 13, 17, ..., 971, 977, 983, 991, 997]
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>
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:
1 2 3 5 7 11 ... 977 983 991 997