# Algorithm Implementation/Mathematics/Primality Testing

In practice, primality testing for numbers of a size suitable for cryptographic applications has to be done in a probabilistic way. Such an algorithm can tell whether a given number is prime with extremely high probability, but cannot provide a certain proof. Many such algorithms operate with a variable number of rounds of testing; each additional round of testing makes it less likely to falsely classify the number, so the implementer can choose a suitable number of rounds to provide a trade off between the level of certainty and the running time of the algorithm.

## Miller-Rabin (or Rabin-Miller) primality test

### Pseudocode

```Input: n > 3, an odd integer to be tested for primality;
Input: 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 − 2]
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
```

### Python

```def is_probable_prime(n, k = 7):
"""use Rabin-Miller algorithm to return True (n is probably prime)
or False (n is definitely composite)"""
if n < 6:  # assuming n >= 0 in all cases... shortcut small cases here
return [False, False, True, True, False, True][n]
elif n & 1 == 0:  # should be faster than n % 2
return False
else:
s, d = 0, n - 1
while d & 1 == 0:
s, d = s + 1, d >> 1
# Use random.randint(2, n-2) for very large numbers
for a in random.sample(xrange(2, min(n - 2, sys.maxint)), min(n - 4, k)):
x = pow(a, d, n)
if x != 1 and x + 1 != n:
for r in xrange(1, s):
x = pow(x, 2, n)
if x == 1:
return False  # composite for sure
elif x == n - 1:
a = 0  # so we know loop didn't continue to end
break  # could be strong liar, try another a
if a:
return False  # composite if we reached end of this loop
return True  # probably prime if reached end of outer loop
```

### Perl

```use Math::BigInt;
use POSIX qw/INT_MAX/;

sub is_probable_prime(\$\$){
my (\$n,\$k) = (Math::BigInt->new(\$_)->babs(),\$_);
return 1 if \$n == 2 || \$n == 3;
return 0 if \$n < 5 || (\$n & 1) == 0;
my (\$s,\$d) = (0,\$n-1);
(\$s,\$d) = (\$s+1,\$d>>1) while (\$d & 1) == 0;
my \$maxa = (\$n-1 < int(INT_MAX) ? \$n-1 : int(INT_MAX));
for my \$i (0..\$k){
my \$a = Math::BigInt->new(int(rand(\$maxa-2)+2));
my \$x = \$a->bmodpow(\$d,\$n);
if(\$x != 1 && \$x+1 != \$n){
for my \$j (1..\$s){
\$x = \$x->bmodpow(2,\$n);
return 0 if \$x == 1;
if (\$x == \$n - 1){
\$a = 0;
last;
}
}
return 0 if \$a;
}
}
return 1;
}
```

### Java

For educational purpose only. See the official source code of BigInteger.isProbablePrime for a proper implementation on e.g. OpenJDK 7 or GNU classpath.

```import java.math.BigInteger;
import java.util.Random;

public class MillerRabin {

private static final BigInteger ZERO = BigInteger.ZERO;
private static final BigInteger ONE = BigInteger.ONE;
private static final BigInteger TWO = new BigInteger("2");
private static final BigInteger THREE = new BigInteger("3");

public static boolean isProbablePrime(BigInteger n, int k) {
if (n.compareTo(ONE) == 0)
return false;
if (n.compareTo(THREE) < 0)
return true;
int s = 0;
BigInteger d = n.subtract(ONE);
while (d.mod(TWO).equals(ZERO)) {
s++;
d = d.divide(TWO);
}
for (int i = 0; i < k; i++) {
BigInteger a = uniformRandom(TWO, n.subtract(ONE));
BigInteger x = a.modPow(d, n);
if (x.equals(ONE) || x.equals(n.subtract(ONE)))
continue;
int r = 0;
for (; r < s; r++) {
x = x.modPow(TWO, n);
if (x.equals(ONE))
return false;
if (x.equals(n.subtract(ONE)))
break;
}
if (r == s) // None of the steps made x equal n-1.
return false;
}
return true;
}

private static BigInteger uniformRandom(BigInteger bottom, BigInteger top) {
Random rnd = new Random();
BigInteger res;
do {
res = new BigInteger(top.bitLength(), rnd);
} while (res.compareTo(bottom) < 0 || res.compareTo(top) > 0);
return res;
}

public static void main(String[] args) {
// run with -ea to enable assertions
String[] primes = { "1", "3", "3613", "7297",
"226673591177742970257407", "2932031007403" };
String[] nonPrimes = { "3341", "2932021007403",
"226673591177742970257405" };
int k = 40;
for (String p : primes)
assert isProbablePrime(new BigInteger(p), k);
for (String n : nonPrimes)
assert !isProbablePrime(new BigInteger(n), k);
}
}
```