Estoy aprendiendo a programar en Python y para practicar estoy haciendo los problemas de Project Euler. El décimo problema pide encontrar la suma de todos los números primos menores que 2 millones. Para ello, una opción posible es ir comprobando cada número natural menor que 2 millones. Si el número es primo, se añade a la suma, y si es compuesto, no se añade. Para saber si un número es primo o no, lo que hay que hacer es un test de primalidad.
El test de primalidad más simple para saber si $n$ es un número primo o no, es considerar los números naturales $i$ con $1 < i < n$. Si existe un número natural $i$ tal que $i$ divide a $n$, entonces automáticamente $n$ no es primo. Por contra, si $i$ no divide a $n$ para todo $d$ con $1<i<n$, entonces $n$ es primo. Podemos afinar un poco más considerando sólo los $i$ con $1<i\leq\sqrt{n}$, porque ningún $i>\sqrt{n}$ divide a $n$, sea $n$ primo o no. De esta manera, tenemos un algoritmo determinista:
Entrada: un número natural $n$
Salida: Cierto si $n$ es primo, Falso si es compuesto.
Si $n = 1$ retorna Cierto.
$i := 2$
Mientras $i\leq\sqrt{n}$ hacer
Si $n \equiv 0\ (\text{mod } i)$ retorna Falso
$i := i+1$
Retorna Cierto.
La implementación de este algoritmo es muy sencilla y el algoritmo funciona perfectamente: en tiempo finito nos dirá, sin lugar a dudas, si un número $n$ dado es primo o compuesto. Sin embargo, este algoritmo tiene un problema, y es que el número de operaciones para decidir si $n$ es o no primo es proporcional a $n$, o dicho de manera abreviada, la complejidad del algoritmo es $\mathcal{O}(n)$. Como el tamaño de un número $n$ en base $2$ es $log_2(n)$, se dice que el algoritmo encuentra la solución en tiempo exponencial respecto del tamaño de la entrada, ya que para un número en base $2$ que ocupe $k$ bits, el número de operaciones necesarias en función de $k$ será proporcional a $2^k$.
Que un algoritmo proporcione un resultado en tiempo exponencial es algo inaceptable. Para resolver el problema 10 de Project Euler, si usamos este test de primalidad, posiblemente necesitemos alrededor de 10 minutos para encontrar el valor de la suma. Puede que alguien piense que este tiempo está bien, pero si nos hubieran pedido calcular la suma de los primos menores que 200 millones en lugar de 2 millones, posiblemente empezaríamos a tratar de buscar otra estrategia.
Una alternativa es usar el llamado Test de Miller-Rabin, que es un algoritmo probabilístico (no determinista) para saber si un número natural $n$ es primo o no. Lo bueno es que el Test de Miller-Rabin es capaz de dar un resultado en tiempo polinomial, es decir, que el tiempo que tarda el algoritmo es un polinomio en el tamaño $k$ de la entrada, en lugar de ser una exponencial como en el caso anterior. Lo malo es que el algoritmo puede equivocarse con cierta probabilidad que podemos ajustar (tardará más cuanto menor queramos que sea la probabilidad de que se equivoque), pero esto no debería asustarnos porque los algoritmos probabilísticos, más rápidos que sus compañeros deterministas, están a la orden del día, llegando a usarse incluso para hacer funcionar los aviones (fuente: yo lo suelto y si cuela, cuela, y si no, me la pela). Es un prejuicio cognitivo sentirse inseguro en un avión que utiliza un algoritmo con una probabilidad de fallo de $10^{-10}$ (cien veces más improbable que que te toque el bote de los Euromillones) y no hacerlo ante la mayor probabilidad de que ocurra un fallo mecánico o humano.
El test de Miller-Rabin se basa en algunas consideraciones muy sencillas sobre los anillos de enteros módulo $n$, $\mathbb{Z}/n\mathbb{Z}$, que se estudian en cualquier asignatura de matemáticas básicas de cualquier ingeniería, y en el Pequeño Teorema de Fermat, que nos dice que si $p$ es un número primo, entonces para cualquier $a$ con $1 \leq a \leq p-1$ se tiene
$$a^{p-1} \equiv 1\ (\text{mod }p)$$
Inicialmente se observa lo siguiente: supongamos que $p$ es un número primo y $x$ es un número natural que cumple
$$x^2 \equiv1\ (\text{mod }p)$$
o sea, que $x$ es una raíz cuadrada de $1$ módulo $p$. Esto equivale a que $x^2 - 1 \equiv 0\ (\text{mod }p)$, o que $(x+1)(x-1) \equiv 0\ (\text{mod }p)$, lo cual quiere decir que $(x+1)(x-1)$ es un múltiplo de $p$, y como $p$ es primo, o bien $(x+1)$ es un múltiplo de $p$, o bien $(x-1)$ es un múltiplo de $p$, es decir, que o bien $x \equiv 1\ (\text{mod }p)$, o bien $x \equiv -1\ (\text{mod} p)$.
Sea ahora $n$ un primo mayor que $2$. En particular, $n$ es impar y $n-1$ es par, así que podemos sacar todos los factores $2$ de $n-1$ para expresarlo en la forma $n-1 = 2^sd$, con $d$ impar. Veamos que, para cada $a$ con $1 \leq a \leq n-1$, se tiene que
$$a^d \equiv 1\ (\text{mod }n)$$
o bien que
$$a^{2^rd} \equiv -1\ (\text{mod }n)$$
para algún $r$ en $0 \leq r \leq s-1$. En efecto, si expresamos $a^{n-1}$ como $a^{2^sd}$, el Pequeño Teorema de Fermat afirma que
$$a^{2^sd}\equiv 1\ (\text{mod }n)$$
(recordemos que suponemos $n$ primo) y, como antes, $a^{2^{s-1}d} \equiv 1\ (\text{mod }n)$ o bien $a^{2^{s-1}d} \equiv -1(\text{mod }n)$. Si $s-1 = 0$, entonces se satisface la primera de las igualdades. En caso contrario, si $a^{2^{s-1}d} \equiv -1(\text{mod }n)$, entonces se satisface la segunda igualdad para $r = s-1$. Si no es así, volvemos a aplicar el proceso a $a^{2^{s-1}d} \equiv 1\ (\text{mod }n)$ hasta que obtengamos un $1$ o un $-1$.
¿Cómo nos puede ayudar esto para diseñar un test de primalidad? Muy sencillo: si $n$ es primo, entonces para cualquier $1 \leq a \leq n-1$ se satisface alguna de las igualdades
$$a^d \equiv 1\ (\text{mod }n)$$
o
$$a^{2^rd} \equiv -1\ (\text{mod }n)$$
para algún $0 \leq r \leq s-1$. La afirmación contrarrecíproca equivalente es: si existe un $a$ con $1 \leq a \leq n-1$, llamado testigo, tal que
$$a^d \not\equiv 1\ (\text{mod }n)$$
y
$$a^{2^rd} \not\equiv -1\ (\text{mod }n)$$
para todo $0 \leq r \leq s-1$, entonces $n$ es compuesto. Lo que hacemos es lo siguiente: seleccionamos aleatoriamente una base $a$ con $1 \leq a \leq n-1$. Si $n$ cumple alguna de las igualdades, enconces $n$ es un probable primo. Si no cumple ninguna, entonces hemos encontrado un testigo y sabemos con seguridad que $n$ es compuesto.
Si $n$ es primo, entonces no hay ningún testigo, por lo que el test siempre dirá que $n$ es un probable primo. Si $n$ es compuesto, puede que tengamos la mala suerte de escoger un $a$ que no sea un testigo (lo que se llama un mentiroso). Se puede demostrar que, si $n$ es compuesto, entonces por lo menos $3/4$ partes de las bases son testigos, así que la probabilidad de que el test diga que un número compuesto $n$ es un probable primo es, como mucho, $1/4$. Para conseguir menor probabilidad de error, lo que hacemos es realizar $k$ veces el test sobre el mismo número $n$, de tal modo que si realizamos $k$ veces el test sobre un número compuesto $n$, la probabilidad de que el test diga que es un probable primo es $(1/4)^k$. El esquema del test básico es
Entrada: un número natural impar $n>3$
Salida: Cierto si $n$ es un probable primo, Falso si $n$ es compuesto
$s := 0$
$d := n-1$
Mientras $d \equiv 0\ (\text{mod }2)$ hacer
$d = d/2$
$s := s+1$
Escoger $a$ entero aleatorio uniformemente en $\{1, 2, ..., n-1\}$.
Si $a^d \equiv 1\ (\text{mod }n)$
Retorna Cierto
$r := 0$
Mientras $r\leq s-1$ hacer
Si $a^{2^rd} \equiv n-1\ (\text{mod }n)$
Retorna Cierto
$r = r+1$
Retorna Falso
El algoritmo final consistirá simplemente en repetir $k$ veces el test. Si alguna de las veces el test dice que $n$ es compuesto, entonces $n$ será seguro compuesto. Si todas las veces el test dice que $n$ es primo, entonces $n$ será un probable primo. Al final de esta entrada dejo el código que he escrito yo en Python para implementar el test.
¿Qué hay sobre el tiempo que tarda el algoritmo? Pues bien, se puede probar que es $\mathcal{O}(\log^2(n))$ para un test, y $\mathcal{O}(k\log^2(n))$ si realizamos el test $k$ veces. En cualquier caso, es un tiempo polinomial. He hecho un pequeño experimento que compara ambos tests:
Edito: por si las moscas, pongo también el código del método principal:
¿Qué hay sobre el tiempo que tarda el algoritmo? Pues bien, se puede probar que es $\mathcal{O}(\log^2(n))$ para un test, y $\mathcal{O}(k\log^2(n))$ si realizamos el test $k$ veces. En cualquier caso, es un tiempo polinomial. He hecho un pequeño experimento que compara ambos tests:
Las curvas, en rojo el test determinista, en azul el Test de Miller-Rabin con $k=4$, indican cuánto tiempo $t(N)$ tardan ambos métodos en decidir si son primos o no todos los números naturales entre $1$ y $N$, para $N$ hasta 1 millón. La tendencia indica que, a pesar de que para números por debajo de 1 millón las diferencias de tiempo no son muy grandes (unos pocos segundos), para números realmente grandes como los que se utilizan en criptografía puede haber una diferencia enorme entre usar uno u otro algoritmo.
Ahora discutamos un momento el asunto de la fiablidad del test. Sabemos que, si $n$ es compuesto, entonces el test ejecutado $k$ veces dirá que $n$ es un probable primo con probabilidad $(1/4)^k$. Es decir, que si llamamos $X$ al suceso "$n$ es compuesto" y $Y_k$ al suceso "el test ejecutado $k$ veces dice que $n$ es un probable primo, tenemos que $P(Y_k|X) = (1/4)^k$. Ahora bien, si no tenemos datos sobre $n$, y el test nos dice que $n$ es un posible primo, ¿cuál es la probablidad de que el test se haya equivocado y en realidad $n$ sea compuesto? Es decir, ¿cuánto vale $P(X|Y_k)$? Llamando $\overline{X}$ al suceso "$n$ es primo", el Teorema de Bayes nos da la respuesta:
$$P(X|Y_k) = \frac{P(Y_k|X)P(X)}{P(Y_k|X)P(X) + P(Y_k|\overline{X})P(\overline{X})}$$
Es decir, que la probabilidad de que el algoritmo se equivoque depende también de la propia densidad de los números primos $P(X)$, y es en la práctica bastante menor que $(1/4)^k$.
El Test de Miller-Rabin se puede modificar de forma que obtengamos un test determinista. Para ello, lo que tendríamos que hacer es, en lugar de escoger bases $a$ aleatoriamente, comprobar todas ellas desde $1$ hasta una cierta cota. Si $n$ es un número compuesto, los mentirosos $a$ coprimos con $n$ (los que hacen que $n$ pase el test) están en un subgrupo propio de $(\mathbb{Z}/n\mathbb{Z})^*$ (el grupo de las unidades de $\mathbb{Z}/n\mathbb{Z}$), por lo que, si hacemos el test sobre un conjunto de bases $a$ que genere todo $(\mathbb{Z}/n\mathbb{Z})^*$, entonces habrá alguna de las bases que sea un testigo para $n$. Se puede demostrar que, si la todavía no demostrada Hipótesis Generalizada de Riemann (HGR) es cierta, entonces los $a$ tal que $1 \leq a \leq 2log^2(n)$ son un conjunto que generan $(\mathbb{Z}/n\mathbb{Z})^*$ y, consecuentemente, sólo hay que testear en estas bases obteniendo así un algoritmo determinista con tiempo de ejecución polinomial $\mathcal{O}(log^4(n))$, que es mayor que el tiempo del Test de Miller-Rabin probabilista, pero es determinista y da una respuesta en tiempo polinomial, no exponencial.
Que la HGR no esté demostrada todavía no es un gran problema, ya que desde el año 2002 tenemos un test de primalidad determinista que corre en tiempo polinomial, el Test AKS.
Para terminar, dejo el código de la implementación del Test de Miller-Rabin en Python:
Edito: por si las moscas, pongo también el código del método principal:
No hay comentarios:
Publicar un comentario