Miller–Rabin primality test: Difference between revisions
m →{{header|REXX}}: added/changed comments, optimized the SUSPENDERS subroutine. |
m →{{header|REXX}}: reinstated original output introduction text. |
||
Line 2,308: | Line 2,308: | ||
end /*j*/ /*this comment not left blank. */ |
end /*j*/ /*this comment not left blank. */ |
||
return /*whew! All done with the primes*/</lang> |
return /*whew! All done with the primes*/</lang> |
||
'''output''' when using the input of: <tt> 10000 10 </tt> |
|||
<pre style="height:30ex"> |
<pre style="height:30ex"> |
||
There are 1229 primes ≤ 10000 |
There are 1229 primes ≤ 10000 |
Revision as of 05:16, 11 May 2014
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.
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.
<lang Ada>generic
type Number is range <>;
package Miller_Rabin is
type Result_Type is (Composite, Probably_Prime);
function Is_Prime (N : Number; K : Positive := 10) return Result_Type;
end Miller_Rabin;</lang>
The implementation of that package is as follows:
<lang Ada>with Ada.Numerics.Discrete_Random;
package body Miller_Rabin is
function Is_Prime (N : Number; K : Positive := 10) return Result_Type is subtype Number_Range is Number range 2 .. N - 1; package Random is new Ada.Numerics.Discrete_Random (Number_Range);
function Mod_Exp (Base, Exponent, Modulus : Number) return Number is Result : Number := 1; begin for E in 1 .. Exponent loop Result := Result * Base mod Modulus; end loop; return Result; end Mod_Exp;
Generator : Random.Generator; D : Number := N - 1; S : Natural := 0; X : Number; begin -- exclude 2 and even numbers if N = 2 then return Probably_Prime; elsif N mod 2 = 0 then return Composite; end if;
-- write N-1 as 2**S * D, with D mod 2 /= 0 while D mod 2 = 0 loop D := D / 2; S := S + 1; end loop;
-- initialize RNG Random.Reset (Generator); for Loops in 1 .. K loop X := Mod_Exp(Random.Random (Generator), D, N); if X /= 1 and X /= N - 1 then Inner : for R in 1 .. S - 1 loop X := Mod_Exp (X, 2, N); if X = 1 then return Composite; end if; exit Inner when X = N - 1; end loop Inner; if X /= N - 1 then return Composite; end if; end if; end loop;
return Probably_Prime; end Is_Prime;
end Miller_Rabin;</lang>
Finally, the program itself:
<lang Ada>with Ada.Text_IO, Miller_Rabin;
procedure Mr_Tst is
type Number is range 0 .. (2**48)-1;
package Num_IO is new Ada.Text_IO.Integer_IO (Number); package Pos_IO is new Ada.Text_IO.Integer_IO (Positive); package MR is new Miller_Rabin(Number); use MR;
N : Number; K : Positive;
begin
for I in Number(2) .. 1000 loop if Is_Prime (I) = Probably_Prime then Ada.Text_IO.Put (Number'Image (I)); end if; end loop; Ada.Text_IO.Put_Line (".");
Ada.Text_IO.Put ("Enter a Number: "); Num_IO.Get (N); Ada.Text_IO.Put ("Enter the count of loops: "); Pos_IO.Get (K); Ada.Text_IO.Put_Line ("What is it? " & Result_Type'Image (Is_Prime(N, K)));
end MR_Tst;</lang>
- Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997. Enter a Number: 1234567 Enter the count of loops: 20 What is it? COMPOSITE
using an external library to handle big integers
Using the big integer implementation from a cryptographic library [1].
<lang Ada>with Ada.Text_IO, Crypto.Types.Big_Numbers, Ada.Numerics.Discrete_Random;
procedure Miller_Rabin is
Bound: constant Positive := 256; -- can be any multiple of 32
package LN is new Crypto.Types.Big_Numbers (Bound); use type LN.Big_Unsigned; -- all computations "mod 2**Bound"
function "+"(S: String) return LN.Big_Unsigned renames LN.Utils.To_Big_Unsigned;
function Is_Prime (N : LN.Big_Unsigned; K : Positive := 10) return Boolean is
subtype Mod_32 is Crypto.Types.Mod_Type; use type Mod_32; package R_32 is new Ada.Numerics.Discrete_Random (Mod_32); Generator : R_32.Generator;
function Random return LN.Big_Unsigned is X: LN.Big_Unsigned := LN.Big_Unsigned_Zero; begin for I in 1 .. Bound/32 loop X := (X * 2**16) * 2**16; X := X + R_32.Random(Generator); end loop; return X; end Random;
D: LN.Big_Unsigned := N - LN.Big_Unsigned_One; S: Natural := 0; A, X: LN.Big_Unsigned; begin -- exclude 2 and even numbers if N = 2 then return True; elsif N mod 2 = LN.Big_Unsigned_Zero then return False; else
-- write N-1 as 2**S * D, with odd D while D mod 2 = LN.Big_Unsigned_Zero loop D := D / 2; S := S + 1; end loop;
-- initialize RNG R_32.Reset (Generator);
-- run the real test for Loops in 1 .. K loop loop A := Random; exit when (A > 1) and (A < (N - 1)); end loop; X := LN.Mod_Utils.Pow(A, D, N); -- X := (Random**D) mod N if X /= 1 and X /= N - 1 then Inner: for R in 1 .. S - 1 loop X := LN.Mod_Utils.Pow(X, LN.Big_Unsigned_Two, N); if X = 1 then return False; end if; exit Inner when X = N - 1; end loop Inner; if X /= N - 1 then return False; end if; end if; end loop; end if;
return True; end Is_Prime;
S: constant String := "4547337172376300111955330758342147474062293202868155909489"; T: constant String := "4547337172376300111955330758342147474062293202868155909393";
K: constant Positive := 10;
begin
Ada.Text_IO.Put_Line("Prime(" & S & ")=" & Boolean'Image(Is_Prime(+S, K))); Ada.Text_IO.Put_Line("Prime(" & T & ")=" & Boolean'Image(Is_Prime(+T, K)));
end Miller_Rabin;</lang>
Output:
Prime(4547337172376300111955330758342147474062293202868155909489)=TRUE Prime(4547337172376300111955330758342147474062293202868155909393)=FALSE
Using the built-in Miller-Rabin test from the same library:
<lang Ada>with Ada.Text_IO, Crypto.Types.Big_Numbers, Ada.Numerics.Discrete_Random;
procedure Miller_Rabin is
Bound: constant Positive := 256; -- can be any multiple of 32
package LN is new Crypto.Types.Big_Numbers (Bound); use type LN.Big_Unsigned; -- all computations "mod 2**Bound"
function "+"(S: String) return LN.Big_Unsigned renames LN.Utils.To_Big_Unsigned;
S: constant String := "4547337172376300111955330758342147474062293202868155909489"; T: constant String := "4547337172376300111955330758342147474062293202868155909393";
K: constant Positive := 10;
begin
Ada.Text_IO.Put_Line("Prime(" & S & ")=" & Boolean'Image (LN.Mod_Utils.Passed_Miller_Rabin_Test(+S, K))); Ada.Text_IO.Put_Line("Prime(" & T & ")=" & Boolean'Image (LN.Mod_Utils.Passed_Miller_Rabin_Test(+T, K)));
end Miller_Rabin;</lang>
The output is the same.
ALGOL 68
<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>
Deterministic up to 341,550,071,728,321
<lang c>// calcul a^n%mod size_t power(size_t a, size_t n, size_t mod) {
size_t power = a; size_t result = 1;
while (n) { if (n & 1) result = (result * power) % mod; power = (power * power) % mod; n >>= 1; } return result;
}
// n−1 = 2^s * d with d odd by factoring powers of 2 from n−1 bool witness(size_t n, size_t s, size_t d, size_t a) {
size_t x = power(a, d, n); size_t y;
while (s) { y = (x * x) % n; if (y == 1 && x != 1 && x != n-1) return false; x = y; --s; } if (y != 1) return false; return true;
}
/*
* if n < 1,373,653, it is enough to test a = 2 and 3; * if n < 9,080,191, it is enough to test a = 31 and 73; * if n < 4,759,123,141, it is enough to test a = 2, 7, and 61; * if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803; * if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11; * if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13; * if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17. */
bool is_prime_mr(size_t n) {
if (((!(n & 1)) && n != 2 ) || (n < 2) || (n % 3 == 0 && n != 3)) return false; if (n <= 3) return true;
size_t d = n / 2; size_t s = 1; while (!(d & 1)) { d /= 2; ++s; }
if (n < 1373653) return witness(n, s, d, 2) && witness(n, s, d, 3); if (n < 9080191) return witness(n, s, d, 31) && witness(n, s, d, 73); if (n < 4759123141) return witness(n, s, d, 2) && witness(n, s, d, 7) && witness(n, s, d, 61); if (n < 1122004669633) return witness(n, s, d, 2) && witness(n, s, d, 13) && witness(n, s, d, 23) && witness(n, s, d, 1662803); if (n < 2152302898747) return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11); if (n < 3474749660383) return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13); return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13) && witness(n, s, d, 17);
}</lang> Inspiration from http://stackoverflow.com/questions/4424374/determining-if-a-number-is-prime
C#
<lang csharp>public static class RabinMiller {
public static bool IsPrime(int n, int k) {
if(n < 2)
{
return false;
}
if(n != 2 && n % 2 == 0)
{
return false;
}
int s = n - 1; while(s % 2 == 0)
{
s >>= 1;
} Random r = new Random();
for (int i = 0; i < k; i++)
{ double a = r.Next((int)n - 1) + 1; int temp = s; int mod = (int)Math.Pow(a, (double)temp) % n; while(temp != n - 1 && mod != 1 && mod != n - 1) {
mod = (mod * mod) % n; temp = temp * 2;
}
if(mod != n - 1 && temp % 2 == 0)
{
return false;
} }
return true;
}
}</lang> <lang csharp>// Miller-Rabin primality test as an extension method on the BigInteger type. // Based on the Ruby implementation on this page. public static class BigIntegerExtensions {
public static bool IsProbablePrime(this BigInteger source, int certainty) { if(source == 2 || source == 3) return true; if(source < 2 || source % 2 == 0) return false;
BigInteger d = source - 1; int s = 0;
while(d % 2 == 0) { d /= 2; s += 1; }
// There is no built-in method for generating random BigInteger values. // Instead, random BigIntegers are constructed from randomly generated // byte arrays of the same length as the source. RandomNumberGenerator rng = RandomNumberGenerator.Create(); byte[] bytes = new byte[source.ToByteArray().LongLength]; BigInteger a;
for(int i = 0; i < certainty; i++) { do { // This may raise an exception in Mono 2.10.8 and earlier. // http://bugzilla.xamarin.com/show_bug.cgi?id=2761 rng.GetBytes(bytes); a = new BigInteger(bytes); } while(a < 2 || a >= source - 2);
BigInteger x = BigInteger.ModPow(a, d, source); if(x == 1 || x == source - 1) continue;
for(int r = 1; r < s; r++) { x = BigInteger.ModPow(x, 2, source); if(x == 1) return false; if(x == source - 1) break; }
if(x != source - 1) return false; }
return true; }
}</lang>
Common Lisp
<lang lisp>(defun factor-out (number divisor)
"Return two values R and E such that NUMBER = DIVISOR^E * R, and R is not divisible by DIVISOR." (do ((e 0 (1+ e)) (r number (/ r divisor))) ((/= (mod r divisor) 0) (values r e))))
(defun mult-mod (x y modulus) (mod (* x y) modulus))
(defun expt-mod (base exponent modulus)
"Fast modular exponentiation by repeated squaring." (labels ((expt-mod-iter (b e p) (cond ((= e 0) p) ((evenp e) (expt-mod-iter (mult-mod b b modulus) (/ e 2) p)) (t (expt-mod-iter b (1- e) (mult-mod b p modulus)))))) (expt-mod-iter base exponent 1)))
(defun random-in-range (lower upper)
"Return a random integer from the range [lower..upper]." (+ lower (random (+ (- upper lower) 1))))
(defun miller-rabin-test (n k)
"Test N for primality by performing the Miller-Rabin test K times. Return NIL if N is composite, and T if N is probably prime." (cond ((= n 1) nil) ((< n 4) t) ((evenp n) nil) (t (multiple-value-bind (d s) (factor-out (- n 1) 2) (labels ((strong-liar? (a) (let ((x (expt-mod a d n))) (or (= x 1) (loop repeat s for y = x then (mult-mod y y n) thereis (= y (- n 1))))))) (loop repeat k always (strong-liar? (random-in-range 2 (- n 2)))))))))</lang>
CL-USER> (last (loop for i from 1 to 1000 when (miller-rabin-test i 10) collect i) 10) (937 941 947 953 967 971 977 983 991 997)
D
<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).
-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).</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
- 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
- Original Rosetta code has been simplified to be easier to follow
Another Miller Rabin test can be found in D. Amos's Haskell for Math module Primes.hs <lang Haskell>module Primes where
import System.Random import System.IO.Unsafe
-- Miller-Rabin wrapped up as an (almost deterministic) pure function isPrime :: Integer -> Bool isPrime n = unsafePerformIO (isMillerRabinPrime 100 n)
isMillerRabinPrime :: Int -> Integer -> IO Bool
isMillerRabinPrime k n
| even n = return (n==2) | n < 100 = return (n `elem` primesTo100) | otherwise = do ws <- witnesses k n return $ and [test n (pred n) evens (head odds) a | a <- ws] where (evens,odds) = span even (iterate (`div` 2) (pred n))
test :: Integral nat => nat -> nat -> [nat] -> nat -> nat -> Bool test n n_1 evens d a = x `elem` [1,n_1] || n_1 `elem` powers
where x = powerMod n a d powers = map (powerMod n a) evens
witnesses :: (Num a, Ord a, Random a) => Int -> a -> IO [a] witnesses k n
| n < 9080191 = return [31,73] | n < 4759123141 = return [2,7,61] | n < 3474749660383 = return [2,3,5,7,11,13] | n < 341550071728321 = return [2,3,5,7,11,13,17] | otherwise = do g <- newStdGen return $ take k (randomRs (2,n-1) g)
primesTo100 :: [Integer] primesTo100 = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
-- powerMod m x n = x^n `mod` m powerMod :: Integral nat => nat -> nat -> nat -> nat powerMod m x n = f (n - 1) x x `rem` m
where f d a y = if d==0 then y else g d a y g i b y | even i = g (i `quot` 2) (b*b `rem` m) y | otherwise = f (i-1) b (b*y `rem` m)
</lang>
- Sample output:
Testing in GHCi: ~> isPrime 4547337172376300111955330758342147474062293202868155909489 True *~> isPrime 4547337172376300111955330758342147474062293202868155909393 False *~> dropWhile (<900) $ filter isPrime [2..1000] [907,911,919,929,937,941,947,953,967,971,977,983,991,997]
Icon and Unicon
The following runs in both languages: <lang unicon>procedure main(A)
every n := !A do write(n," is ",(mrp(n,5),"probably prime")|"composite")
end
procedure mrp(n, k)
if n = 2 then return "" if n%2 = 0 then fail nm1 := decompose(n-1) s := nm1[1] d := nm1[2] every !k do { a := ?(n-2)+1 x := (a^d)%n if x = (1|(n-1)) then next every !(s-1) do { x := (x*x)%n if x = 1 then fail if x = (n-1) then break next } fail } return ""
end
procedure decompose(nm1)
s := 1 d := nm1 while d%2 = 0 do { d /:= 2 s +:= 1 } return [s,d]
end</lang>
Sample run:
->mrp 219 221 223 225 227 229 219 is composite 221 is composite 223 is probably prime 225 is composite 227 is probably prime 229 is probably prime ->
J
See Primality Tests essay on the J wiki.
Java
The Miller-Rabin primality test is part of the standard library (java.math.BigInteger) <lang java>import java.math.BigInteger;
public class MillerRabinPrimalityTest {
public static void main(String[] args) { BigInteger n = new BigInteger(args[0]); int certainty = Integer.parseInt(args[1]); System.out.println(n.toString() + " is " + (n.isProbablePrime(certainty) ? "probably prime" : "composite")); }
}</lang>
- Sample output:
java MillerRabinPrimalityTest 123456791234567891234567 1000000 123456791234567891234567 is probably prime
JavaScript
This covers (almost) all integers in JavaScript (~2^53).
<lang JavaScript>function modProd(a,b,n){
if(b==0) return 0; if(b==1) return a%n; return (modProd(a,(b-b%10)/10,n)*10+(b%10)*a)%n;
} function modPow(a,b,n){
if(b==0) return 1; if(b==1) return a%n; if(b%2==0){ var c=modPow(a,b/2,n); return modProd(c,c,n); } return modProd(a,modPow(a,b-1,n),n);
} function isPrime(n){
if(n==2||n==3||n==5) return true; if(n%2==0||n%3==0||n%5==0) return false; if(n<25) return true; for(var a=[2,3,5,7,11,13,17,19],b=n-1,d,t,i,x;b%2==0;b/=2); for(i=0;i<a.length;i++){ x=modPow(a[i],b,n); if(x==1||x==n-1) continue; for(t=true,d=b;t&&d<n-1;d*=2){ x=modProd(x,x,n); if(x==n-1) t=false; } if(t) return false; } return true;
}
for(var i=1;i<=1000;i++) if(isPrime(i)) console.log(i);</lang>
Julia
The built-in isprime
function uses the Miller-Rabin primality test. Here is the implementation of isprime
from the Julia standard library (Julia version 0.2):
<lang julia>
witnesses(n::Union(Uint8,Int8,Uint16,Int16)) = (2,3)
witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? (2,3) : (2,7,61)
witnesses(n::Union(Uint64,Int64)) =
n < 1373653 ? (2,3) : n < 4759123141 ? (2,7,61) : n < 2152302898747 ? (2,3,5,7,11) : n < 3474749660383 ? (2,3,5,7,11,13) : (2,325,9375,28178,450775,9780504,1795265022)
function isprime(n::Integer)
n == 2 && return true (n < 2) | iseven(n) && return false s = trailing_zeros(n-1) d = (n-1) >>> s for a in witnesses(n) a < n || break x = powermod(a,d,n) x == 1 && continue t = s while x != n-1 (t-=1) <= 0 && return false x = oftype(n, Base.widemul(x,x) % n) x == 1 && return false end end return true
end </lang>
Liberty BASIC
<lang lb> DIM mersenne(11) mersenne(1)=7 mersenne(2)=31 mersenne(3)=127 mersenne(4)=8191 mersenne(5)=131071 mersenne(6)=524287 mersenne(7)=2147483647 mersenne(8)=2305843009213693951 mersenne(9)=618970019642690137449562111 mersenne(10)=162259276829213363391578010288127 mersenne(11)=170141183460469231731687303715884105727
dim SmallPrimes(1000)
data 2, 3, 5, 7, 11, 13, 17, 19, 23, 29
data 31, 37, 41, 43, 47, 53, 59, 61, 67, 71
data 73, 79, 83, 89, 97, 101, 103, 107, 109, 113
data 127, 131, 137, 139, 149, 151, 157, 163, 167, 173
data 179, 181, 191, 193, 197, 199, 211, 223, 227, 229
data 233, 239, 241, 251, 257, 263, 269, 271, 277, 281
data 283, 293, 307, 311, 313, 317, 331, 337, 347, 349
data 353, 359, 367, 373, 379, 383, 389, 397, 401, 409
data 419, 421, 431, 433, 439, 443, 449, 457, 461, 463
data 467, 479, 487, 491, 499, 503, 509, 521, 523, 541
data 547, 557, 563, 569, 571, 577, 587, 593, 599, 601
data 607, 613, 617, 619, 631, 641, 643, 647, 653, 659
data 661, 673, 677, 683, 691, 701, 709, 719, 727, 733
data 739, 743, 751, 757, 761, 769, 773, 787, 797, 809
data 811, 821, 823, 827, 829, 839, 853, 857, 859, 863
data 877, 881, 883, 887, 907, 911, 919, 929, 937, 941
data 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013
data 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069
data 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151
data 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223
data 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291
data 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373
data 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451
data 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511
data 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583
data 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657
data 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733
data 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811
data 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889
data 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987
data 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053
data 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129
data 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213
data 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287
data 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357
data 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423
data 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531
data 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617
data 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687
data 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741
data 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819
data 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903
data 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999
data 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079
data 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181
data 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257
data 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331
data 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413
data 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511
data 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571
data 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643
data 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727
data 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821
data 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907
data 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989
data 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057
data 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139
data 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231
data 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297
data 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409
data 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493
data 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583
data 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657
data 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751
data 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831
data 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937
data 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003
data 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087
data 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179
data 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279
data 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387
data 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443
data 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521
data 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639
data 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693
data 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791
data 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857
data 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939
data 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053
data 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133
data 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221
data 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301
data 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367
data 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473
data 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571
data 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673
data 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761
data 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833
data 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917
data 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997
data 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103
data 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207
data 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297
data 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411
data 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499
data 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561
data 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643
data 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723
data 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829
data 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919
print "Liberty Miller Rabin Demonstration"
print "Loading Small Primes"
for i=1 to 1000: read x : SmallPrimes(i)=x :next :NoOfSmallPrimes=1000
print NoOfSmallPrimes;" Primes Loaded"
'Prompt "Enter number to test:";resp$ 'x=val(resp$) 'goto [Jump]
For i=1 to 11
x=mersenne(i)
t1=time$("ms") [TryAnother] print
iterations=0 [Loop] iterations=iterations+1
if MillerRabin(x,7)=1 then t2=time$("ms") print "Composite, found in ";t2-t1;" milliseconds" else t2=time$("ms") print x;" Probably Prime. Tested in ";t2-t1;" milliseconds" playwave "tada.wav", async end if print
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 </lang>
Mathematica
<lang Mathematica>MillerRabin[n_,k_]:=Module[{d=n-1,s=0,test=True},While[Mod[d,2]==0 ,d/=2 ;s++] Do[
a=RandomInteger[{2,n-1}]; x=PowerMod[a,d,n]; If[x!=1, For[ r = 0, r < s, r++, If[x==n-1, Continue[]]; x = Mod[x*x, n]; ]; If[ x != n-1, test=False ]; ];
,{k}]; Print[test] ]</lang>
- Example output (not using the PrimeQ builtin):
<lang mathematica>MillerRabin[17388,10] ->False</lang>
Maxima
<lang maxima>/* Miller-Rabin algorithm is builtin, see function primep. Here is another implementation */
/* find highest power of p, p^s, that divide n, and return s and n / p^s */
facpow(n, p) := block(
[s: 0], while mod(n, p) = 0 do (s: s + 1, n: quotient(n, p)), [s, n]
)$
/* check whether n is a strong pseudoprime to base a; s and d are given by facpow(n - 1, 2) */
sppp(n, a, s, d) := block(
[x: power_mod(a, d, n), q: false], if x = 1 or x = n - 1 then true else ( from 2 thru s do ( x: mod(x * x, n), if x = 1 then return(q: false) elseif x = n - 1 then return(q: true) ), q )
)$
/* Miller-Rabin primality test. For n < 341550071728321, the test is deterministic;
for larger n, the number of bases tested is given by the option variable primep_number_of_tests, which is used by Maxima in primep. The bound for deterministic test is also the same as in primep. */
miller_rabin(n) := block(
[v: [2, 3, 5, 7, 11, 13, 17], s, d, q: true, a], if n < 19 then member(n, v) else ( [s, d]: facpow(n - 1, 2), if n < 341550071728321 then ( /* see http://oeis.org/A014233 */ for a in v do ( if not sppp(n, a, s, d) then return(q: false) ), q ) else ( thru primep_number_of_tests do ( a: 2 + random(n - 3), if not sppp(n, a, s, d) then return(q: false) ), q ) )
)$</lang>
PARI/GP
Built-in
<lang parigp>MR(n,k)=ispseudoprime(n,k);</lang>
Custom
<lang parigp>sprp(n,b)={ my(s = valuation(n-1, 2), d = Mod(b, n)^(n >> s)); if (d == 1, return(1)); for(i=1,s-1, if (d == -1, return(1)); d = d^2; ); d == -1 };
MR(n,k)={
for(i=1,k, if(!sprp(n,random(n-2)+2), return(0)) ); 1
};</lang>
Deterministic version
A basic deterministic test can be obtained by an appeal to the ERH (as proposed by Gary Miller) and a result of Eric Bach (improving on Joseph Oesterlé). Calculations of Jan Feitsma can be used to speed calculations below 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># the expmod-function from: http://rosettacode.org/wiki/Modular_exponentiation sub expmod(Int $a is copy, Int $b is copy, $n) { my $c = 1; repeat while $b div= 2 { ($c *= $a) %= $n if $b % 2; ($a *= $a) %= $n; } $c; }
subset PrimeCandidate of Int where { $_ > 2 and $_ % 2 };
my Bool multi sub is-prime(Int $n, Int $k) { return False; } my Bool multi sub is-prime(2, Int $k) { return True; } my Bool multi sub is-prime(PrimeCandidate $n, Int $k) { my Int $d = $n - 1; my Int $s = 0;
while $d %% 2 { $d div= 2; $s++; }
for (2 ..^ $n).pick($k) -> $a { my $x = expmod($a, $d, $n);
# one could just write "next if $x == 1 | $n - 1" # but this takes much more time in current rakudo/nom next if $x == 1 or $x == $n - 1;
for 1 ..^ $s { $x = $x ** 2 mod $n; return False if $x == 1; last if $x == $n - 1; } return False if $x !== $n - 1; }
return True; }
say (1..1000).grep({ is-prime($_, 10) }).join(", "); </lang>
PHP
<lang php><?php function is_prime($n, $k) {
if ($n == 2) return true; if ($n < 2 || $n % 2 == 0) return false;
$d = $n - 1; $s = 0;
while ($d % 2 == 0) { $d /= 2; $s++; }
for ($i = 0; $i < $k; $i++) { $a = rand(2, $n-1);
$x = bcpowmod($a, $d, $n); if ($x == 1 || $x == $n-1) continue;
for ($j = 1; $j < $s; $j++) { $x = bcmod(bcmul($x, $x), $n); if ($x == 1) return false; if ($x == $n-1) continue 2; } return false; } return true;
}
for ($i = 1; $i <= 1000; $i++)
if (is_prime($i, 10)) echo "$i, ";
echo "\n"; ?></lang>
PicoLisp
<lang PicoLisp>(de longRand (N)
(use (R D) (while (=0 (setq R (abs (rand))))) (until (> R N) (unless (=0 (setq D (abs (rand)))) (setq R (* R D)) ) ) (% R N) ) )
(de **Mod (X Y N)
(let M 1 (loop (when (bit? 1 Y) (setq M (% (* M X) N)) ) (T (=0 (setq Y (>> 1 Y))) M ) (setq X (% (* X X) N)) ) ) )
(de _prim? (N D S)
(use (A X R) (while (> 2 (setq A (longRand N)))) (setq R 0 X (**Mod A D N)) (loop (T (or (and (=0 R) (= 1 X)) (= X (dec N)) ) T ) (T (or (and (> R 0) (= 1 X)) (>= (inc 'R) S) ) NIL ) (setq X (% (* X X) N)) ) ) )
(de prime? (N K)
(default K 50) (and (> N 1) (bit? 1 N) (let (D (dec N) S 0) (until (bit? 1 D) (setq D (>> 1 D) S (inc S) ) ) (do K (NIL (_prim? N D S)) T ) ) ) )</lang>
- Output:
: (filter '((I) (prime? I)) (range 937 1000)) -> (937 941 947 953 967 971 977 983 991 997) : (prime? 4547337172376300111955330758342147474062293202868155909489) -> T : (prime? 4547337172376300111955330758342147474062293202868155909393) -> NIL
Prolog
<lang 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) ).</lang>
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
The correctness of this code is in question! See discussion
Python: Probably correct answers
This versions will give answers with a very small probability of being false. That probability being dependent on _mrpt_num_trials and the random numbers used for name a
passed to function try_composite.
<lang python>import random
_mrpt_num_trials = 5 # number of bases to test
def is_probable_prime(n):
""" Miller-Rabin primality test.
A return value of False means n is certainly not prime. A return value of True means n is very likely a prime.
>>> is_probable_prime(1) Traceback (most recent call last): ... AssertionError >>> is_probable_prime(2) True >>> is_probable_prime(3) True >>> is_probable_prime(4) False >>> is_probable_prime(5) True >>> is_probable_prime(123456789) False
>>> primes_under_1000 = [i for i in range(2, 1000) if is_probable_prime(i)] >>> len(primes_under_1000) 168 >>> primes_under_1000[-10:] [937, 941, 947, 953, 967, 971, 977, 983, 991, 997]
>>> is_probable_prime(6438080068035544392301298549614926991513861075340134\
3291807343952413826484237063006136971539473913409092293733259038472039\ 7133335969549256322620979036686633213903952966175107096769180017646161\ 851573147596390153)
True
>>> is_probable_prime(7438080068035544392301298549614926991513861075340134\
3291807343952413826484237063006136971539473913409092293733259038472039\ 7133335969549256322620979036686633213903952966175107096769180017646161\ 851573147596390153)
False """ assert n >= 2 # special case 2 if n == 2: return True # ensure n is odd if n % 2 == 0: return False # write n-1 as 2**s * d # repeatedly try to divide n-1 by 2 s = 0 d = n-1 while True: quotient, remainder = divmod(d, 2) if remainder == 1: break s += 1 d = quotient assert(2**s * d == n-1)
# test the base a to see whether it is a witness for the compositeness of n def try_composite(a): if pow(a, d, n) == 1: return False for i in range(s): if pow(a, 2**i * d, n) == n-1: return False return True # n is definitely composite
for i in range(_mrpt_num_trials): a = random.randrange(2, n) if try_composite(a): return False
return True # no base tested showed n as composite</lang>
Python: Proved correct up to large N
This versions will give correct answers for n
less than 341550071728321 and then reverting to the probabilistic form of the first solution. By selecting a certain number of primes for name a
instead of random values mathematicians have proved the general algorithm correct.
For 341550071728321 and beyond, I have followed the pattern in choosing a
from the set of prime numbers.
<lang python>def _try_composite(a, d, n, s):
if pow(a, d, n) == 1: return False for i in range(s): if pow(a, 2**i * d, n) == n-1: return False return True # n is definitely composite
def is_prime(n, _precision_for_huge_n=16):
if n in _known_primes or n in (0, 1): return True if any((n % p) == 0 for p in _known_primes): return False d, s = n - 1, 0 while not d % 2: d, s = d >> 1, s + 1 # Returns exact according to http://primes.utm.edu/prove/prove2_3.html if n < 1373653: return not any(_try_composite(a, d, n, s) for a in (2, 3)) if n < 25326001: return not any(_try_composite(a, d, n, s) for a in (2, 3, 5)) if n < 118670087467: if n == 3215031751: return False return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7)) if n < 2152302898747: return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11)) if n < 3474749660383: return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13)) if n < 341550071728321: return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17)) # otherwise return not any(_try_composite(a, d, n, s) for a in _known_primes[:_precision_for_huge_n])
_known_primes = [2, 3] _known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]</lang>
- Testing
Includes test values from other examples:
>>> is_prime(4547337172376300111955330758342147474062293202868155909489) True >>> is_prime(4547337172376300111955330758342147474062293202868155909393) False >>> [x for x in range(901, 1000) if is_prime(x)] [907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997] >>> is_prime(643808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153) True >>> is_prime(743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153) False >>>
Racket
<lang Racket>#lang racket (define (miller-rabin-expmod base exp m)
(define (squaremod-with-check x) (define (check-nontrivial-sqrt1 x square) (if (and (= square 1) (not (= x 1)) (not (= x (- m 1)))) 0 square)) (check-nontrivial-sqrt1 x (remainder (expt x 2) m))) (cond ((= exp 0) 1) ((even? exp) (squaremod-with-check (miller-rabin-expmod base (/ exp 2) m))) (else (remainder (* base (miller-rabin-expmod base (- exp 1) m)) m))))
(define (miller-rabin-test n)
(define (try-it a) (define (check-it x) (and (not (= x 0)) (= x 1))) (check-it (miller-rabin-expmod a (- n 1) n))) (try-it (+ 1 (random (remainder (- n 1) 4294967087)))))
(define (fast-prime? n times)
(for/and ((i (in-range times))) (miller-rabin-test n)))
(define (prime? n(times 100))
(fast-prime? n times))
(prime? 4547337172376300111955330758342147474062293202868155909489) ;-> outputs true </lang>
REXX
With a K of 1, there seems to be a not uncommon number of failures, but
- with a K ≥ 2, the failures are rare,
- with a K ≥ 3, rare as hen's teeth.
This would be in the same vein as: 3 is prime, 5 is prime, 7 is prime, all odd numbers are prime. <lang rexx>/*REXX program puts 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*/</lang> output when using the input of: 10000 10
There are 1229 primes ≤ 10000 ─────────────────────────────────────────────────────────────────────────────── [K=2] 2701 isn't prime ! for 1──►10000, K=2, Miller─Rabin primality test found 1230 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=2, Miller─Rabin primality test found 1229 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=3, Miller─Rabin primality test found 1229 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=4, Miller─Rabin primality test found 1229 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=5, Miller─Rabin primality test found 1229 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=6, Miller─Rabin primality test found 1229 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=7, Miller─Rabin primality test found 1229 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=8, Miller─Rabin primality test found 1229 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=9, Miller─Rabin primality test found 1229 primes {out of 1229} ─────────────────────────────────────────────────────────────────────────────── for 1──►10000, K=10, Miller─Rabin primality test found 1229 primes {out of 1229}
Ruby
<lang ruby> require 'openssl' def miller_rabin_prime?(n,g)
d = n - 1 s = 0 while d % 2 == 0 d /= 2 s += 1 end g.times do a = 2 + rand(n-4) x = OpenSSL::BN::new(a.to_s).mod_exp(d,n) #x = (a**d) % n next if x == 1 or x == n-1 for r in (1 .. s-1) x = x.mod_exp(2,n) #x = (x**2) % n return false if x == 1 break if x == n-1 end return false if x != n-1 end true # probably
end
p primes = (3..1000).step(2).find_all {|i| miller_rabin_prime?(i,10)} </lang>
- Output:
[3, 5, 7, 11, 13, 17, ..., 971, 977, 983, 991, 997]
The following larger examples all produce true: <lang ruby> puts miller_rabin_prime?(94366396730334173383107353049414959521528815310548187030165936229578960209523421808912459795329035203510284576187160076386643700441216547732914250578934261891510827140267043592007225160798348913639472564715055445201512461359359488795427875530231001298552452230535485049737222714000227878890892901228389026881,1000) puts miller_rabin_prime?(138028649176899647846076023812164793645371887571371559091892986639999096471811910222267538577825033963552683101137782650479906670021895135954212738694784814783986671046107023185842481502719762055887490765764329237651328922972514308635045190654896041748716218441926626988737664133219271115413563418353821396401,1000) puts miller_rabin_prime?(123301261697053560451930527879636974557474268923771832437126939266601921428796348203611050423256894847735769138870460373141723679005090549101566289920247264982095246187318303659027201708559916949810035265951104246512008259674244307851578647894027803356820480862664695522389066327012330793517771435385653616841,1000) puts miller_rabin_prime?(119432521682023078841121052226157857003721669633106050345198988740042219728400958282159638484144822421840470442893056822510584029066504295892189315912923804894933736660559950053226576719285711831138657839435060908151231090715952576998400120335346005544083959311246562842277496260598128781581003807229557518839,1000) puts miller_rabin_prime?(132082885240291678440073580124226578272473600569147812319294626601995619845059779715619475871419551319029519794232989255381829366374647864619189704922722431776563860747714706040922215308646535910589305924065089149684429555813953571007126408164577035854428632242206880193165045777949624510896312005014225526731,1000) puts miller_rabin_prime?(153410708946188157980279532372610756837706984448408515364579602515073276538040155990230789600191915021209039203172105094957316552912585741177975853552299222501069267567888742458519569317286299134843250075228359900070009684517875782331709619287588451883575354340318132216817231993558066067063143257425853927599,1000) puts miller_rabin_prime?(103130593592068072608023213244858971741946977638988649427937324034014356815504971087381663169829571046157738503075005527471064224791270584831779395959349442093395294980019731027051356344056416276026592333932610954020105156667883269888206386119513058400355612571198438511950152690467372712488391425876725831041,1000) </lang>
Run BASIC
<lang runbasic>input "Input a number:";n input "Input test:";k
test = millerRabin(n,k) if test = 0 then
print "Probably Prime" else print "Composite"
end if wait
' ---------------------------------------- ' Returns ' Composite = 1 ' Probably Prime = 0 ' ----------------------------------------
FUNCTION millerRabin(n, k) if n = 2 then millerRabin = 0 'probablyPrime goto [funEnd] end if
if n mod 2 = 0 or n < 2 then millerRabin = 1 'composite goto [funEnd] end if
d = n - 1 while d mod 2 = 0
d = d / 2 s = s + 1
wend
while k > 0
k = k - 1 x = (int(rnd(1) * (n-4))^d) mod n if x = 1 or x = n-1 then for r=1 To s-1 x =(x * x) mod n if x=1 then millerRabin = 1 ' composite goto [funEnd] end if if x = n-1 then exit for next r if x<>n-1 then millerRabin = 1 ' composite goto [funEnd] end if end if
wend [funEnd] END FUNCTION</lang>
Seed7
<lang seed7>$ include "seed7_05.s7i";
include "bigint.s7i";
const func boolean: millerRabin (in bigInteger: n, in integer: k) is func
result var boolean: probablyPrime is TRUE; local var bigInteger: d is 0_; var integer: r is 0; var integer: s is 0; var bigInteger: a is 0_; var bigInteger: x is 0_; var integer: tests is 0; begin if n < 2_ or (n > 2_ and not odd(n)) then probablyPrime := FALSE; elsif n > 3_ then d := pred(n); s := lowestSetBit(d); d >>:= s; while tests < k and probablyPrime do a := rand(2_, pred(n)); x := modPow(a, d, n); if x <> 1_ and x <> pred(n) then r := 1; while r < s and x <> 1_ and x <> pred(n) do x := modPow(x, 2_, n); incr(r); end while; probablyPrime := x = pred(n); end if; incr(tests); end while; end if; end func;
const proc: main is func
local var bigInteger: number is 0_; begin for number range 2_ to 1000_ do if millerRabin(number, 10) then writeln(number); end if; end for; end func;</lang>
Original source: [2]
Smalltalk
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