Fatoração de Inteiros

Neste artigo, listamos vários algoritmos para fatorar números inteiros, cada um deles pode ser rápido e lento (alguns mais lentos que outros), dependendo de suas entradas.

Se o número que você deseja fatorar for realmente um número primo, a maioria dos algoritmos, especialmente o algoritmo de fatoração de Fermat, o p-1 de Pollard, o algoritmo rho de Pollard, irão ser lentos. Portanto, faz sentido executar um teste de primalidade antes de tentar fatorar o número.

Trial division

Este é o algoritmo mais básico para encontrar a fatoração e os fatores primos do número.

Dividimos por cada divisor possível $d$. Podemos notar que é impossível que todos os fatores primos de um número composto $n$ sejam maiores que $\sqrt{n}$. Portanto, precisamos apenas testar os divisores $2 \le d \le \sqrt{n}$, o que nos fornece a fatoração em $O(\sqrt{n})$.

O menor divisor deve ser um número primo. Removemos o fator do número e repetimos o processo. Se não conseguirmos encontrar nenhum divisor no intervalo $[2; \sqrt{n}]$, então o próprio número deve 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;
}

Wheel factorization

Esta é uma otimização do trial division. A ideia é a seguinte. Uma vez que sabemos que o número não é divisível por 2, não precisamos checar os outros números pares. Isso nos deixa com apenas $50\%$ dos números para checar. Após marcar 2, podemos simplesmente começar com 3 e pular todos os outros números.

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 pode ser estendido. Se o número não for divisível por 3, também podemos ignorar todos os outros múltiplos de 3 nos cálculos futuros. Então precisamos checar apenas os números $5, 7, 11, 13, 17, 19, 23, \dots$. Podemos observar um padrão desses números restantes. Precisamos verificar todos os números com $d \bmod 6 = 1$ e $d \bmod 6 = 5$. Portanto, isso nos deixa com apenas $33.3\%$ dos números para checar. Podemos implementar isso verificando os números primos 2 e 3 primeiro e, em seguida, verificar com 5 e, alternativamente, pular 1 ou 3 números.

Podemos estender isso ainda mais. Aqui está uma implementação para os números primos 2, 3 e 5. É conveniente usar uma array para armazenar quanto precisamos pular.

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

Se estendermos ainda mais isso com mais números primos, podemos alcançar porcentagens melhores. No entanto, também as listas de pulos ficarão muito maiores.

Primos pré-calculados

Estendendo o "wheel factorization" com mais e mais primos deixará exatamente apenas primos para checar. Portanto, uma boa maneira de verificar é pré-calcular todos os números primos com o crivo de Eratóstenes até $\sqrt{n}$ e testá-los 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 fatoração de Fermat

Podemos escrever um número composto ímpar $n = p \cdot q$ como a diferença de dois quadrados $n = a^2 - b^2$: $$n = \left(\frac{p + q}{2}\right)^2 - \left(\frac{p - q}{2}\right)^2$$ O método de fatoração de Fermat tenta explorar o fato, adivinhando o primeiro quadrado $a^2$, e verificando se a parte restante $b^2 = a^2 - n$ também é um número quadrado. Se for, encontramos os fatores $a - b$ e $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;
}

Observe que esse método de fatoração pode ser muito rápido, se a diferença entre os dois fatores $p$ e $q$ for pequena. O algoritmo é executado em tempo $O(|p - q|)$. No entanto, como é muito lento, uma vez que os fatores estão distantes, raramente é usado na prática.

No entanto, ainda há um grande número de otimizações para essa abordagem. Por exemplo, olhando para os quadrados $a^2$ modulo um número fixo pequeno, você pode notar que não precisa olhar para determinados valores de $a$ já que eles não consegue formar um quadrado $a^2 - n$.

Método $p - 1$ de Pollard

É muito provável que pelo menos um fator de um número seja $B$. $B$-powersmooth significa que toda exponenciação $d^k$ de um primo $d$ que divide $p-1$ é no máximo $B$. Por exemplo, a fatoração de $4817191$ são $1303 \cdot 3697$. E os fatores máximo são $31$ e $16$ respectivamente, pois $1303 - 1 = 2 \cdot 3 \cdot 7 \cdot 31$ e $3697 - 1 = 2^4 \cdot 3 \cdot 7 \cdot 11$. Em 1974 John Pollard inventou um método para extrair esses fatores $B$ de um número composto.

A ideia vem do pequeno Teorema de Fermat. Deixe que a fatoração de $n$ seja $n = p \cdot q$. Ele diz que se $a$ é coprimo para $p$(gcd(a,p)=1), a seguinte declaração é válida:

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

Isso também significa que

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

Portanto, para qualquer $M$ com $p - 1 ~|~ M$ sabemos que $a^M \equiv 1$. Isso significa que $a^M - 1 = p \cdot r$, e, por causa disso, também $p ~|~ \gcd(a^M - 1, n)$.

Portanto, se $p - 1$ para um fator $p$ de $n$ divide $M$, podemos extrair um fator usando algoritmo de Euclides.

O menor $M$ que é um múltiplo de todo número $B$ é $\text{lcm/mmc}(1,~2~,3~,4~,~\dots,~B)$. Ou alternativamente: $$M = \prod_{\text{primos } q \le B} q^{\lfloor \log_q B \rfloor}$$

Observe que, se $p-1$ divide $M$ para todos fatores primos $p$ de $n$, então $\gcd(a^M - 1, n)$ será apenas $n$. Nesse caso, não recebemos um fator. Portanto, tentaremos executar o $\gcd/mdc$ várias vezes, enquanto calculamos $M$.

Alguns números compostos não têm fatores $B$ para pequenos $B$. Por exemplo, os fatores do número composto $100.000.000.000.000.493 = 763.013 \cdot 131.059.365.961$ são $190.753$ e $1092161383$. Nós teríamos que escolher $B >= 190.753$ para decompor o número.

Na implementação seguinte. começamos com $B = 10$ e incrementamos $B$ depois de cada iteração.

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

Nota: Este é um algoritmo probabilístico, portanto, pode acontecer que o algoritmo não encontre um fator.

A complexidade é $O(B \log B \log^2 n)$ por iteração.

Algoritmo rho de Pollard

Outro algoritmo de fatoração de John Pollard.

Seja a fatoração primária de um número $n = p q$. O algoritmo analisa uma sequência pseudo-aleatória ${x_i} = {x_0,~f(x_0),~f(f(x_0)),~\dots}$ onde $f$ é uma função polinomial, usualmente $f(x) = x^2 + c \bmod n$ é escolhido com $c = 1$.

Na verdade, não estamos muito interessados ​​na sequência ${x_i}$, estamos mais interessados ​​na sequência ${x_i \bmod p}$. Como $f$ é uma função polinomial e todos os valores estão no range $[0;~p)$ essa sequência começará um ciclo mais cedo ou mais tarde. O birthday paradox sugere, na verdade, que o número esperado de elementos é $O(\sqrt{p})$ até que a repetição comece. Se $p$ for menor que $\sqrt{n}$, a repetição começará muito provavelmente em $O(\sqrt[4]{n})$.

Aqui está uma visualização dessa sequência ${x_i \bmod p}$ com $n = 2206637$, $p = 317$, $x_0 = 2$ e $f(x) = x^2 + 1$.

Pollard's rho visualization

Ainda existe uma grande questão em aberto. Ainda não sabemos $p$, então como podemos discutir sobre a sequência ${x_i \bmod p}$?

Na verdade, é bem fácil. Existe um ciclo na sequência $\{x_i \bmod p\}_{i \le j}$ se, e somente se, houver dois índices $s, t \le j$ e $t$ com $x_s \equiv x_t \bmod p$. Essa equação pode ser reescrita como $x_s - x_t \equiv 0 \bmod p$ que é igual a $p ~|~ \gcd(x_s - x_t, n)$.

Portanto, se encontrarmos dois índices $s$ e $t$ com $g = \gcd(x_s - x_t, n) > 1$, então encontramos um ciclo e também um fator $g$ de $n$. Note que é possível que $g = n$. Nesse caso, não encontramos um fator adequado e precisamos repetir o algoritmo com parâmetro diferente (valor inicial diferente $x_0$, constante diferente $c$ na função polinomial $f$).

Para encontrar o ciclo, podemos usar qualquer algoritmo comum de detecção de ciclo.

Algoritmo para encontrar ciclos de Floyd

Esse algoritmo encontra um ciclo usando dois ponteiros. Esses ponteiros se movem sobre a sequência em velocidades diferentes. Em cada iteração, o primeiro ponteiro avança para o próximo elemento, mas o segundo ponteiro avança dois elementos. Não é difícil ver que, se existir um ciclo, o segundo ponteiro fará pelo menos um ciclo completo e, em seguida, encontrará o primeiro ponteiro durante os próximos loops do ciclo. Se o comprimento do ciclo é $\lambda$ e $\mu$ é o primeiro índice onde o ciclo começa, então o algoritmo será executado em tempo $O(\lambda + \mu)$.

Esse algoritmo também é conhecido como algoritmo da tartaruga e da lebre, com base na história em que uma tartaruga e uma lebre fazem uma corrida.

Na verdade, é possível determinar o parâmetro $\lambda$ e $\mu$ usando esse algoritmo (também em tempo $O(\lambda + \mu)$ e espaço $O(1)$), mas aqui está apenas a versão simplificada para encontrar o ciclo. O algoritmo e retorna verdadeiro assim que detecta um ciclo. Se a sequência não tiver um ciclo, a função nunca irá parar. No entanto, isso não pode acontecer durante o algoritmo de Pollard.

function floyd(f, x0):
    tartaruga = x0
    lebre = f(x0)
    while tartaruga != lebre:
        tartaruga = f(tartaruga)
        lebre = f(f(lebre))
    return true

Implementação

Primeiro, aqui está uma implementação usando o algoritmo de ciclos de Floyd. O algoritmo é executado (normalmente) em tempo $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;
}

A tabela a seguir mostra os valores de $x$ e $y$ durante a execução do algoritmo com $n = 2206637$, $x_0 = 2$ e $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}$$

A implementação usa uma função mult, na qual multiplica dois inteiros $\le 10^{18}$, sem overflow, usando o tipo do GCC: __int128 para inteiros 128-bits. Se o GCC não estiver disponível, você poderá usar uma ideia semelhante à exponenciação binária.

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, você também pode implementar a multiplicação modular rápida.

Como visto anteriormente: se $n$ é composto e o algoritmo retorna $n$ como um fator, você precisa repetir o procedimento com os parâmetros $x_0$ e $c$ diferentes. Por exemplo, a opção $x_0 = c = 1$ não vai fatorar $25 = 5 \cdot 5$. O algoritmo apenas retornará $25$. Entretanto, a escolha $x_0 = 1$, $c = 2$ irá fatorar.

Algoritmo de Brent

Brent usa um algoritmo similar ao de Floyd. Também utiliza dois ponteiros. Ao invés de avançar os ponteiros por um e dois respectivamente, avançamos eles em potências de 2. Assim que $2^i$ for maior que $\lambda$ e $\mu$, encontraremos o ciclo.

function floyd(f, x0):
    tartaruga = x0
    lebre = f(x0)
    l = 1
    while tartaruga != lebre:
        tartaruga = lebre
        repete 1 vez:
            lebre = f(lebre)
            if tartaruga == lebre:
                return true
        l *= 2
    return true

O algoritmo de Brent também é executado em tempo linear, mas geralmente é mais rápido que o algoritmo de Floyd, uma vez que utiliza menos avaliações da função $f$.

Implementação

A implementação direta usando os algoritmos de Brent pode ser acelerada com a observação de que podemos omitir os termos $x_l - x_k$ se $k < \frac{3 \cdot l}{2}$. Além disso, em vez de executar o cálculo do $\gcd/mdc$ em cada etapa, multiplicamos os termos e calculamos apenas a cada poucas etapas.

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

A combinação de um trial division para pequenos números primos, juntamente com a versão de Brent do algoritmo rho de Pollard, criará um algoritmo de fatoração muito eficiente.

Problemas