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 implementor can choose a suitable number of rounds to provide a trade off between the level of certainty and the running time of the algorithm.

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
      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($_[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

import java.math.BigInteger;
public class MillerRabin {
 
        public static final BigInteger ZERO = BigInteger.ZERO;  // declaring constants
        public static final BigInteger ONE  = BigInteger.ONE;
        public static final BigInteger TWO = BigInteger.valueOf(2);
        public static final int[] aValues = {2,3,5,7,11,13,17,19,23,29,31,37,41};  // values for bases
 
        public static boolean testPr (BigInteger n, BigInteger a, int s, BigInteger d){
                for (int i=0; i<s; i++){
                        BigInteger exp = TWO.pow(i);
                        exp = exp.multiply(d);
                        BigInteger res = a.modPow(exp, n);
                        if (res.equals(n.subtract(ONE)) || res.equals(ONE)){
                                return true;
                        }
                }
                return false;
        }
 
        public static boolean millerRabin (BigInteger n, int numValues){ // n = number to test
                BigInteger d = n.subtract(ONE);                          // numValues = # of bases to test
                int s = 0;
                while (d.mod(TWO).equals(ZERO)){
                        s++;
                        d = d.divide(TWO);
                }
                System.out.print("Base ");
                for (int i=0; i<numValues; i++){     // loops through the bases, terminating early if 
                        BigInteger a = BigInteger.valueOf(aValues[i]);  // composite
                        boolean r = testPr(n, a, s, d);
                        System.out.print(aValues[i] + " ");
                        if (!r){
                                return false;
                        }
                }
                return true;
        }
}