Miller–Rabin primality test
You are encouraged to solve this task according to the task description, using any language you may know.
| This page uses content from Wikipedia. The original article was at Miller–Rabin primality test. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance) |
The Miller–Rabin primality test or Rabin–Miller primality test is a primality test: an algorithm which determines whether a given number is prime or not. The algorithm, as modified by Michael O. Rabin to avoid the generalized Riemann hypothesis, is a probabilistic algorithm.
The pseudocode, from Wikipedia is:
Input: n > 2, an odd integer to be tested for primality;
k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
write n − 1 as 2s·d with d odd by factoring powers of 2 from n − 1
LOOP: repeat k times:
pick a randomly in the range [2, n − 1]
x ← ad mod n
if x = 1 or x = n − 1 then do next LOOP
for r = 1 .. s − 1
x ← x2 mod n
if x = 1 then return composite
if x = n − 1 then do next LOOP
return composite
return probably prime
- The nature of the test involves big numbers, so the use of "big numbers" libraries (or similar features of the language of your choice) are suggested, but not mandatory.
- Deterministic variants of the test exist and can be implemented as extra (not mandatory to complete the task)
Ada
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.
<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
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> Sample output:
937 941 947 953 967 971 977 983 991 997
AutoHotkey
ahk forum: discussion <lang AutoHotkey>MsgBox % MillerRabin(999983,10) ; 1 MsgBox % MillerRabin(999809,10) ; 1 MsgBox % MillerRabin(999727,10) ; 1 MsgBox % MillerRabin(52633,10) ; 0 MsgBox % MillerRabin(60787,10) ; 0 MsgBox % MillerRabin(999999,10) ; 0 MsgBox % MillerRabin(999995,10) ; 0 MsgBox % MillerRabin(999991,10) ; 0
MillerRabin(n,k) { ; 0: composite, 1: probable prime (n < 2**31)
d := n-1, s := 0
While !(d&1)
d>>=1, s++
Loop %k% {
Random a, 2, n-2 ; if n < 4,759,123,141, it is enough to test a = 2, 7, and 61.
x := PowMod(a,d,n)
If (x=1 || x=n-1)
Continue
Cont := 0
Loop % s-1 {
x := PowMod(x,2,n)
If (x = 1)
Return 0
If (x = n-1) {
Cont = 1
Break
}
}
IfEqual Cont,1, Continue
Return 0
}
Return 1
}
PowMod(x,n,m) { ; x**n mod m
y := 1, i := n, z := x
While i>0
y := i&1 ? mod(y*z,m) : y, z := mod(z*z,m), i >>= 1
Return y
}</lang>
bc
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>
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 C sharp|C#> public static class RabinMiller {
public static bool isPrime(int p)
{
if(p<2)
{
return false;
}
if(p!=2 && p%2==0)
{
return false;
}
int s=p-1; while(s%2==0)
{
s>>=1;
}
for (int i = 1; i < 11;i++)
{
Random r = new Random();
double a = r.Next((int)p - 1) + 1;
int temp = s;
int mod = (int)Math.Pow(a, (double)temp) % p;
while(temp!=p-1 && mod!=1 && mod!=p-1)
{
mod=(mod*mod)%p; temp=temp*2;
}
if(mod!=p-1 && temp%2==0)
{
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(long n, int k) {
if (n < 2 || n % 2 == 0)
return n == 2;
long d = n - 1;
long s = 0;
while (d % 2 == 0) {
d /= 2;
s++;
}
bool isComposite(long a) {
if (modpow(a, d, n) == 1)
return false;
foreach (i; 0 .. s + 1) {
if (modpow(a, 2^^i * d, n) == n - 1)
return false;
}
return true;
}
foreach (_; 0 .. k + 1) {
if (isComposite(uniform(2, n)))
return false;
}
return true;
}
long modpow(long b, long e, long m) {
long result = 1;
while (e > 0) {
if ((e & 1) == 1) {
result = (result * b) % m;
}
b = (b * b) % m;
e >>= 1;
}
return result;
}</lang>
E
<lang e>def millerRabinPrimalityTest(n :(int > 0), k :int, random) :boolean {
if (n <=> 2 || n <=> 3) { return true }
if (n <=> 1 || n %% 2 <=> 0) { return false }
var d := n - 1
var s := 0
while (d %% 2 <=> 0) {
d //= 2
s += 1
}
for _ in 1..k {
def nextTrial := __continue
def a := random.nextInt(n - 3) + 2 # [2, n - 2] = [0, n - 4] + 2 = [0, n - 3) + 2
var x := a**d %% n # Note: Will do optimized modular exponentiation
if (x <=> 1 || x <=> n - 1) { nextTrial() }
for _ in 1 .. (s - 1) {
x := x**2 %% n
if (x <=> 1) { return false }
if (x <=> n - 1) { nextTrial() }
}
return false
}
return true
}</lang>
<lang e>for i ? (millerRabinPrimalityTest(i, 1, entropy)) in 4..1000 {
print(i, " ")
} println()</lang>
Fortran
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
Go has it in the standard library. ProbablyPrime
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 [1]
Another Miller Rabin test can be found in D. Amos's Haskell for Math module Primes.hs
<lang Haskell>import System.Random import Data.List import Control.Monad import Control.Arrow
primesTo100 = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
powerMod :: (Integral a, Integral b) => a -> a -> b -> a powerMod m _ 0 = 1 powerMod m x n | n > 0 = join (flip f (n - 1)) x `rem` m where
f _ 0 y = y
f a d y = g a d where
g b i | even i = g (b*b `rem` m) (i `quot` 2)
| otherwise = f b (i-1) (b*y `rem` m)
witns :: (Num a, Ord a, Random a) => Int -> a -> IO [a] witns x y = do
g <- newStdGen
let r = [9080191, 4759123141, 2152302898747, 3474749600383, 341550071728321]
fs = [[31,73],[2,7,61],[2,3,5,7,11],[2,3,5,7,11,13],[2,3,5,7,11,13,17]]
if y >= 341550071728321
then return $ take x $ randomRs (2,y-1) g
else return $ 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
PARI/GP
Built-in test: <lang parigp>MR(n,k)=ispseudoprime(n,k);</lang>
Custom function: <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>
A basic deterministic test can be obtained by an appeal to the ERH (as proposed by Gary Miller) and a result of Eric Bach. 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**$d) % $n;
next if $x == 1 or $x == $n-1;
for(1..$s-1)
{
$x = ($x**2) % $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>
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='k"] " z "isn't prime !" /*oopsy-doopsy and whoopsy-daisy!*/ end
say 'for 1-->'limit", K="k', Miller-Rabin primality test found' mrp,
'primes {out of' primePi"}"
end
exit
/*─────────────────────────────────────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
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 /*definately it's not prime. */
if x==nL then leave
end
if x\==nL then return 0
end
/*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>
[2, 3, 5, 7, 11, 13, 17, ..., 971, 977, 983, 991, 997]
Smalltalk
Smalltalk handles big numbers naturally and trasparently (the parent class Integer has many subclasses, and a subclass is picked according to the size of the integer that must be handled) <lang smalltalk>Integer extend [
millerRabinTest: kl [ |k| k := kl.
self <= 3
ifTrue: [ ^true ]
ifFalse: [
(self even)
ifTrue: [ ^false ]
ifFalse: [ |d s|
d := self - 1.
s := 0.
[ (d rem: 2) == 0 ]
whileTrue: [
d := d / 2.
s := s + 1.
].
[ k:=k-1. k >= 0 ]
whileTrue: [ |a x r|
a := Random between: 2 and: (self - 2).
x := (a raisedTo: d) rem: self.
( x = 1 )
ifFalse: [ |r|
r := -1.
[ r := r + 1. (r < s) & (x ~= (self - 1)) ]
whileTrue: [
x := (x raisedTo: 2) rem: self
].
( x ~= (self - 1) ) ifTrue: [ ^false ]
]
].
^true
]
]
]
].</lang>
<lang smalltalk>1 to: 1000 do: [ :n |
(n millerRabinTest: 10) ifTrue: [ n printNl ]
].</lang>
Tcl
Use Tcl 8.5 for large integer support <lang tcl>package require Tcl 8.5
proc miller_rabin {n k} {
if {$n <= 3} {return true}
if {$n % 2 == 0} {return false}
# write n - 1 as 2^s·d with d odd by factoring powers of 2 from n − 1
set d [expr {$n - 1}]
set s 0
while {$d % 2 == 0} {
set d [expr {$d / 2}]
incr s
}
while {$k > 0} {
incr k -1
set a [expr {2 + int(rand()*($n - 4))}]
set x [expr {($a ** $d) % $n}]
if {$x == 1 || $x == $n - 1} continue
for {set r 1} {$r < $s} {incr r} {
set x [expr {($x ** 2) % $n}]
if {$x == 1} {return false}
if {$x == $n - 1} break
}
if {$x != $n-1} {return false}
} return true
}
for {set i 1} {$i < 1000} {incr i} {
if {[miller_rabin $i 10]} {
puts $i
}
}</lang> produces
1 2 3 5 7 11 ... 977 983 991 997