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

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 !=- 1:
                  break  # could be strong liar, try another a
            else:
               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($_[0])->babs(),$_[1]);
  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);
	}
}