Implementación de algoritmos de teoría de números/Test de primalidad de Miller-Rabin

El test de primalidad de Miller-Rabin es un test de primalidad, es decir, un algoritmo para determinar si un número dado es primo, similar al test de primalidad de Fermat. Su versión original fue propuesta por G. L. Miller, se trataba de un algoritmo determinista, pero basado en la no demostrada hipótesis generalizada de Riemann; Michael Oser Rabin modificó la propuesta de Miller para obtener un algoritmo probabilístico que no utiliza resultados no probados.

Pseudocódigo

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

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);
	}
}

Julia

En el lenguaje de programación Julia. Se encuentra en este repositorio.

function miller_rabin_primality_test(n, k = 40)

	if n % 2 == 0
		return false
	end
	if n == 2
		return true
	end

    m = n - 1
	r = 0
	while m % 2 == 0
		r += 1
		m = Int(m / 2)
	end

	for i = 1 : k - 1
		a = first(rand(2 : n - 1, 1))
		b = (a^m) % n
		if b == 1 || b == n - 1
			continue
		end
		check = false
		for i = 1 : r - 1
			b = (b^2) % n
			if b == n - 1
				check = true
				break
			end
		end
		if check
			return false
		end
	end
	return true
end