Factorización de enteros

En este artículo listamos varios algoritmos para la factorización de enteros, cada uno de los cuales puede ser rápido o lento en distintos grados según su entrada.

Nótese que, si el número que se quiere factorizar es en realidad un número primo, la mayoría de los algoritmos se ejecutarán muy lentamente. Esto es especialmente cierto para los algoritmos de factorización de Fermat, p−1 de Pollard y rho de Pollard. Por lo tanto, lo más razonable es realizar un test de primalidad probabilístico (o uno determinista rápido) antes de intentar factorizar el número.

División por tentativa

Este es el algoritmo más básico para hallar una factorización prima.

Dividimos por cada divisor posible $d$. Se puede observar que es imposible que todos los factores primos de un número compuesto $n$ sean mayores que $\sqrt{n}$. Por lo tanto, solo hay que probar los divisores $2 \le d \le \sqrt{n}$, lo que nos da la factorización prima en $O(\sqrt{n})$. (Esto es tiempo seudopolinomial , es decir, polinomial en el valor de la entrada pero exponencial en la cantidad de bits de la entrada.)

El menor divisor debe ser un número primo. Quitamos el número factorizado y continuamos el proceso. Si no podemos hallar ningún divisor en el rango $[2; \sqrt{n}]$, entonces el número mismo tiene que ser primo.

vector<long long> trial_division1(long long n) {
    vector<long long> factorization;
    for (long long d = 2; d * d <= n; d++) {
        while (n % d == 0) {
            factorization.push_back(d);
            n /= d;
        }
    }
    if (n > 1)
        factorization.push_back(n);
    return factorization;
}

Factorización por rueda (wheel factorization)

Esta es una optimización de la división por tentativa. Una vez que sabemos que el número no es divisible por 2, no hace falta comprobar otros números pares. Esto nos deja solo el $50$ de los números para comprobar. Después de extraer el factor 2 y obtener un número impar, podemos empezar simplemente en 3 y contar solo los demás números impares.

vector<long long> trial_division2(long long n) {
    vector<long long> factorization;
    while (n % 2 == 0) {
        factorization.push_back(2);
        n /= 2;
    }
    for (long long d = 3; d * d <= n; d += 2) {
        while (n % d == 0) {
            factorization.push_back(d);
            n /= d;
        }
    }
    if (n > 1)
        factorization.push_back(n);
    return factorization;
}

Este método se puede extender más. Si el número no es divisible por 3, también podemos ignorar todos los demás múltiplos de 3 en los cálculos futuros. Así que solo hay que comprobar los números $5, 7, 11, 13, 17, 19, 23, \dots$. Podemos observar un patrón en estos números restantes. Hay que comprobar todos los números con $d \bmod 6 = 1$ y $d \bmod 6 = 5$. Así nos queda solo el $33.3$ de los números para comprobar. Podemos implementar esto extrayendo primero los primos 2 y 3, después de lo cual empezamos en 5 y solo contamos restos $1$ y $5$ módulo $6$.

Aquí hay una implementación para los números primos 2, 3 y 5. Es conveniente guardar los incrementos de salto en un arreglo.

vector<long long> trial_division3(long long n) {
    vector<long long> factorization;
    for (int d : {2, 3, 5}) {
        while (n % d == 0) {
            factorization.push_back(d);
            n /= d;
        }
    }
    static array<int, 8> increments = {4, 2, 4, 2, 4, 6, 2, 6};
    int i = 0;
    for (long long d = 7; d * d <= n; d += increments[i++]) {
        while (n % d == 0) {
            factorization.push_back(d);
            n /= d;
        }
        if (i == 8)
            i = 0;
    }
    if (n > 1)
        factorization.push_back(n);
    return factorization;
}

Si seguimos extendiendo este método para incluir aún más primos, se pueden alcanzar mejores porcentajes, pero las listas de saltos se volverán más grandes.

Primos precomputados

Extendiendo el método de factorización por rueda de forma indefinida, solo nos quedarán números primos para comprobar. Una buena forma de comprobar esto es precomputar todos los números primos con la Criba de Eratóstenes hasta $\sqrt{n}$, y probarlos individualmente.

vector<long long> primes;

vector<long long> trial_division4(long long n) {
    vector<long long> factorization;
    for (long long d : primes) {
        if (d * d > n)
            break;
        while (n % d == 0) {
            factorization.push_back(d);
            n /= d;
        }
    }
    if (n > 1)
        factorization.push_back(n);
    return factorization;
}

Método de factorización de Fermat

Podemos escribir un número compuesto impar $n = p \cdot q$ como la diferencia de dos cuadrados $n = a^2 - b^2$:

$$n = \left(\frac{p + q}{2}\right)^2 - \left(\frac{p - q}{2}\right)^2$$

El método de factorización de Fermat intenta aprovechar este hecho adivinando el primer cuadrado $a^2$, y comprobando si la parte restante, $b^2 = a^2 - n$, también es un cuadrado perfecto. Si lo es, entonces hemos hallado los factores $a - b$ y $a + b$ de $n$.

int fermat(int n) {
    int a = ceil(sqrt(n));
    int b2 = a*a - n;
    int b = round(sqrt(b2));
    while (b * b != b2) {
        a = a + 1;
        b2 = a*a - n;
        b = round(sqrt(b2));
    }
    return a - b;
}

Este método de factorización puede ser muy rápido si la diferencia entre los dos factores $p$ y $q$ es pequeña. El algoritmo se ejecuta en tiempo $O(|p - q|)$. En la práctica, sin embargo, este método se usa rara vez. Una vez que los factores se alejan más, es extremadamente lento.

Sin embargo, todavía hay una gran cantidad de opciones de optimización respecto a este enfoque. Mirando los cuadrados $a^2$ módulo un número pequeño fijo, se puede observar que ciertos valores $a$ no tienen que considerarse, ya que no pueden producir un cuadrado perfecto $a^2 - n$.

Método $p - 1$ de Pollard { data-toc-label=“Método $p - 1$ de Pollard” }

Es muy probable que un número $n$ tenga al menos un factor primo $p$ tal que $p - 1$ sea $\mathrm{B}$-powersmooth para un $\mathrm{B}$ pequeño. Se dice que un entero $m$ es $\mathrm{B}$-powersmooth (B-potencia-suave) si toda potencia de primo que divide a $m$ es a lo sumo $\mathrm{B}$. De forma formal, sea $\mathrm{B} \geqslant 1$ y sea $m$ cualquier entero positivo. Supongamos que la factorización prima de $m$ es $m = \prod {q_i}^{e_i}$, donde cada $q_i$ es un primo y $e_i \geqslant 1$. Entonces $m$ es $\mathrm{B}$-powersmooth si, para todo $i$, ${q_i}^{e_i} \leqslant \mathrm{B}$. P. ej. la factorización prima de $4817191$ es $1303 \cdot 3697$. Y los valores $1303 - 1$ y $3697 - 1$ son $31$-powersmooth y $16$-powersmooth respectivamente, porque $1303 - 1 = 2 \cdot 3 \cdot 7 \cdot 31$ y $3697 - 1 = 2^4 \cdot 3 \cdot 7 \cdot 11$. En 1974 John Pollard inventó un método para extraer factores $p$, t.q. $p-1$ es $\mathrm{B}$-powersmooth, de un número compuesto.

La idea viene del pequeño teorema de Fermat. Sea $n = p \cdot q$ una factorización de $n$. El teorema afirma que si $a$ es coprimo con $p$, se cumple lo siguiente:

$$a^{p - 1} \equiv 1 \pmod{p}$$

Esto también significa que

$${\left(a^{(p - 1)}\right)}^k \equiv a^{k \cdot (p - 1)} \equiv 1 \pmod{p}.$$

Así, para cualquier $M$ con $p - 1 | M$ sabemos que $a^M \equiv 1$. Esto significa que $a^M - 1 = p \cdot r$, y por eso también $p | \gcd(a^M - 1, n)$.

Por lo tanto, si $p - 1$ para un factor $p$ de $n$ divide a $M$, podemos extraer un factor usando el algoritmo de Euclides.

Está claro que el menor $M$ que es múltiplo de todo número $\mathrm{B}$-powersmooth es $\text{lcm}(1,2,3~,4~,~\dots,~B)$. O, de forma alternativa:

$$M = \prod_{\text{prime } q \le B} q^{\lfloor \log_q B \rfloor}$$

Nótese que, si $p-1$ divide a $M$ para todos los factores primos $p$ de $n$, entonces $\gcd(a^M - 1, n)$ será simplemente $n$. En este caso no obtenemos un factor. Por lo tanto, intentaremos realizar el $\gcd$ varias veces, mientras computamos $M$.

Algunos números compuestos no tienen factores $p$ t.q. $p-1$ sea $\mathrm{B}$-powersmooth para un $\mathrm{B}$ pequeño. Por ejemplo, para el número compuesto $100000000000000493 = 763013 \cdot 131059365961$, los valores $p-1$ son $190753$-powersmooth y $1092161383$-powersmooth respectivamente. Tendremos que elegir $B \geq 190753$ para factorizar el número.

En la siguiente implementación empezamos con $\mathrm{B} = 10$ y aumentamos $\mathrm{B}$ después de cada iteración.

long long pollards_p_minus_1(long long n) {
    int B = 10;
    long long g = 1;
    while (B <= 1000000 && g < n) {
        long long a = 2 + rand() %  (n - 3);
        g = gcd(a, n);
        if (g > 1)
            return g;

        // calcula a^M
        for (int p : primes) {
            if (p >= B)
                continue;
            long long p_power = 1;
            while (p_power * p <= B)
                p_power *= p;
            a = power(a, p_power, n);

            g = gcd(a - 1, n);
            if (g > 1 && g < n)
                return g;
        }
        B *= 2;
    }
    return 1;
}

Obsérvese que este es un algoritmo probabilístico. Una consecuencia de esto es que existe la posibilidad de que el algoritmo no sea capaz de hallar un factor en absoluto.

La complejidad es $O(B \log B \log^2 n)$ por iteración.

Algoritmo rho de Pollard

El algoritmo rho de Pollard es otro algoritmo de factorización de John Pollard.

Sea $n = p q$ la factorización prima de un número. El algoritmo observa una secuencia seudoaleatoria ${x_i} = {x_0,~f(x_0),f(f(x_0)),\dots}$ donde $f$ es una función polinómica; usualmente se elige $f(x) = (x^2 + c) \bmod n$ con $c = 1$.

En este caso, no nos interesa la secuencia ${x_i}$. Nos interesa más la secuencia ${x_i \bmod p}$. Como $f$ es una función polinómica, y todos los valores están en el rango $[0;~p)$, esta secuencia eventualmente convergerá a un bucle. La paradoja del cumpleaños de hecho sugiere que la cantidad esperada de elementos es $O(\sqrt{p})$ hasta que empieza la repetición. Si $p$ es menor que $\sqrt{n}$, la repetición probablemente empezará en $O(\sqrt[4]{n})$.

Aquí hay una visualización de una secuencia de ese tipo ${x_i \bmod p}$ con $n = 2206637$, $p = 317$, $x_0 = 2$ y $f(x) = x^2 + 1$. Por la forma de la secuencia se ve con mucha claridad por qué el algoritmo se llama algoritmo $\rho$ de Pollard.

Aun así, queda una pregunta abierta. ¿Cómo podemos aprovechar las propiedades de la secuencia ${x_i \bmod p}$ a nuestro favor sin siquiera conocer el número $p$ mismo?

En realidad es bastante fácil. Hay un ciclo en la secuencia ${x_i \bmod p}_{i \le j}$ si y solo si hay dos índices $s, t \le j$ tales que $x_s \equiv x_t \bmod p$. Esta ecuación se puede reescribir como $x_s - x_t \equiv 0 \bmod p$, que es lo mismo que $p | \gcd(x_s - x_t, n)$.

Por lo tanto, si hallamos dos índices $s$ y $t$ con $g = \gcd(x_s - x_t, n) > 1$, hemos hallado un ciclo y también un factor $g$ de $n$. Es posible que $g = n$. En este caso no hemos hallado un factor propio, así que hay que repetir el algoritmo con un parámetro distinto (distinto valor inicial $x_0$, distinta constante $c$ en la función polinómica $f$).

Para hallar el ciclo, podemos usar cualquier algoritmo común de detección de ciclos.

Algoritmo de detección de ciclos de Floyd

Este algoritmo halla un ciclo usando dos punteros que recorren la secuencia a distintas velocidades. En cada iteración, el primer puntero avanzará un elemento, mientras que el segundo puntero avanza de a dos elementos. Usando esta idea es fácil observar que, si hay un ciclo, en algún momento el segundo puntero dará la vuelta y se encontrará con el primero durante los bucles. Si la longitud del ciclo es $\lambda$ y $\mu$ es el primer índice en el que empieza el ciclo, entonces el algoritmo se ejecutará en tiempo $O(\lambda + \mu)$.

Este algoritmo también se conoce como el algoritmo de la tortuga y la liebre, basado en el cuento en el que una tortuga (el puntero lento) y una liebre (el puntero más rápido) tienen una carrera.

De hecho es posible determinar los parámetros $\lambda$ y $\mu$ usando este algoritmo (también en tiempo $O(\lambda + \mu)$ y espacio $O(1)$). Cuando se detecta un ciclo, el algoritmo devolverá ‘True’. Si la secuencia no tiene un ciclo, entonces la función se quedará en un bucle infinito. Sin embargo, usando el algoritmo rho de Pollard, esto se puede evitar.

function floyd(f, x0):
    tortoise = x0
    hare = f(x0)
    while tortoise != hare:
        tortoise = f(tortoise)
        hare = f(f(hare))
    return true

Implementación

Primero, aquí hay una implementación que usa el algoritmo de detección de ciclos de Floyd. El algoritmo en general se ejecuta en tiempo $O(\sqrt[4]{n} \log(n))$.

long long mult(long long a, long long b, long long mod) {
    return (__int128)a * b % mod;
}

long long f(long long x, long long c, long long mod) {
    return (mult(x, x, mod) + c) % mod;
}

long long rho(long long n, long long x0=2, long long c=1) {
    long long x = x0;
    long long y = x0;
    long long g = 1;
    while (g == 1) {
        x = f(x, c, n);
        y = f(y, c, n);
        y = f(y, c, n);
        g = gcd(abs(x - y), n);
    }
    return g;
}

La siguiente tabla muestra los valores de $x$ e $y$ durante el algoritmo para $n = 2206637$, $x_0 = 2$ y $c = 1$.

$$\newcommand\T{\Rule{0pt}{1em}{.3em}}$$$$\begin{array}{|l|l|l|l|l|l|} \hline i & x_i \bmod n & x_{2i} \bmod n & x_i \bmod 317 & x_{2i} \bmod 317 & \gcd(x_i - x_{2i}, n) \ \hline 0 & 2 & 2 & 2 & 2 & - \ 1 & 5 & 26 & 5 & 26 & 1 \ 2 & 26 & 458330 & 26 & 265 & 1 \ 3 & 677 & 1671573 & 43 & 32 & 1 \ 4 & 458330 & 641379 & 265 & 88 & 1 \ 5 & 1166412 & 351937 & 169 & 67 & 1 \ 6 & 1671573 & 1264682 & 32 & 169 & 1 \ 7 & 2193080 & 2088470 & 74 & 74 & 317 \ \hline \end{array}$$

La implementación usa una función mult, que multiplica dos enteros $\le 10^{18}$ sin desbordamiento usando el tipo __int128 de GCC para enteros de 128 bits. Si GCC no está disponible, se puede usar una idea similar a la exponenciación binaria.

long long mult(long long a, long long b, long long mod) {
    long long result = 0;
    while (b) {
        if (b & 1)
            result = (result + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return result;
}

Como alternativa, también se puede implementar la multiplicación de Montgomery.

Como se indicó antes, si $n$ es compuesto y el algoritmo devuelve $n$ como factor, hay que repetir el procedimiento con distintos parámetros $x_0$ y $c$. P. ej. la elección $x_0 = c = 1$ no factorizará $25 = 5 \cdot 5$. El algoritmo devolverá $25$. Sin embargo, la elección $x_0 = 1$, $c = 2$ sí lo factorizará.

Algoritmo de Brent

Brent implementa un método similar al de Floyd, usando dos punteros. La diferencia es que, en lugar de avanzar los punteros de a uno y de a dos lugares respectivamente, se avanzan por potencias de dos. En cuanto $2^i$ es mayor que $\lambda$ y $\mu$, hallaremos el ciclo.

function floyd(f, x0):
    tortoise = x0
    hare = f(x0)
    l = 1
    while tortoise != hare:
        tortoise = hare
        repeat l times:
            hare = f(hare)
            if tortoise == hare:
                return true
        l *= 2
    return true

El algoritmo de Brent también se ejecuta en tiempo lineal, pero en general es más rápido que el de Floyd, ya que usa menos evaluaciones de la función $f$.

Implementación

La implementación directa del algoritmo de Brent se puede acelerar omitiendo los términos $x_l - x_k$ si $k < \frac{3 \cdot l}{2}$. Además, en lugar de realizar el cálculo del $\gcd$ en cada paso, multiplicamos los términos y solo comprobamos realmente el $\gcd$ cada ciertos pasos, y retrocedemos si nos pasamos.

long long brent(long long n, long long x0=2, long long c=1) {
    long long x = x0;
    long long g = 1;
    long long q = 1;
    long long xs, y;

    int m = 128;
    int l = 1;
    while (g == 1) {
        y = x;
        for (int i = 1; i < l; i++)
            x = f(x, c, n);
        int k = 0;
        while (k < l && g == 1) {
            xs = x;
            for (int i = 0; i < m && i < l - k; i++) {
                x = f(x, c, n);
                q = mult(q, abs(y - x), n);
            }
            g = gcd(q, n);
            k += m;
        }
        l *= 2;
    }
    if (g == n) {
        do {
            xs = f(xs, c, n);
            g = gcd(abs(xs - y), n);
        } while (g == 1);
    }
    return g;
}

La combinación de una división por tentativa para números primos pequeños junto con la versión de Brent del algoritmo rho de Pollard da un algoritmo de factorización muy potente.

Problemas de práctica

• SPOJ - FACT0 
• SPOJ - FACT1 
• SPOJ - FACT2 
• GCPC 15 - Divisions