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

Entrada: n > 3, un número entero impar que se probará en cuanto a su primalidad;
Entrada: k, un parámetro que determina la precisión de la prueba.
Salida: Compuesto si n es compuesto, de lo contrario Probablemente primo
escribir n − 1 a 2s·d con d impar al factorizar potencias de 2 a partir de n - 1
BUCLE: repitir k veces:
   elegir al azar a en el rango [2, n - 2]
   x ← ad modulo n
   si x = 1 o x = n - 1 entonces haz el siguiente BUCLE
   para r = 1 .. s - 1 (Bucle for con incremento = 1)
      x ← ad modulo n
      si x = 1 entonces devuelve compuesto
      si x = n - 1 entonces haz el siguiente BUCLE
   retornar Compuesto
retornar Probablemente primo

Implementación en distintos lenguajes de programación

Python

def es_probable_primo(n, k = 7):
   """usa el algoritmo Rabin-Miller para devolver Verdadero (n probablemente sea primo)
      o Falso (n es definitivamente compuesto)"""
   if n < 6:  # asumiendo n >= 0 en todos los casos ... pequeño atajo de casos aquí
      return [False, False, True, True, False, True][n]
   elif n & 1 == 0:  # Debería ser más rápido que n % 2
      return False
   else:
      s, d = 0, n - 1
      while d & 1 == 0:
         s, d = s + 1, d >> 1
      # Utilice random.randint (2, n-2) para números muy grandes
      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  # compuesto seguro
               elif x == n - 1:
                  a = 0  # entonces sabemos que el bucle no continuó hasta terminar
                  break  # Podría ser un Falso fuerte, prueba con otro a
            if a:
               return False  # Compuesto si llegamos al final de este bucle
      return True  # probablemente primo si se alcanza el final del bucle exterior

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