Multiplicação de inteiros usando FFT

Quando escrevemos um inteiro a na base d,

$$a = \sum_k a_k d_k,$$

podemos pensar que estamos escrevendo

$$a = P(d), \quad P(x) = \sum_k a_k x^k.$$

Se desejarmos calcular ab onde

$$b = Q(d), \quad Q(x) = \sum_k b_k x^k$$

podemos usar o algoritmo da seção anterior para calcular os coeficientes c_k do produto

$$(PQ)(x) = \sum_k c_k x^k, \quad c_k = \sum_j a_j b_{k-j}$$

e temos ab = (PQ)(d). Os c_k em geral não serão algarismos aceitáveis para uma expansão na base d do inteiro ab mas isto pode facilmente ser corrigido: escrevemos c'₀ = c₀ + d e₀, c'₁ = c₁ - e₀ + d e₁, ..., c'_k = c_k - e_k-1 + d e_k, ..., onde a cada passo tomamos c'_k como sendo um algarismo aceitável. Ao final, teremos

$$ab = \sum_k c'_k d_k,$$

a expansão de ab na base d.

A dificuldade maior reside no fato que as contas descritas na seção anterior envolvem números complexos, e as partes real e imaginária destes números complexos são irracionais. Uma alternativa é fazer as contas usando variáveis do tipo double; teremos inevitavelmente erros de truncamento mas o fato de sabermos que a resposta final é um inteiro nos dá uma oportunidade de corrigir estes erros. É claro que precisamos ter o cuidado de evitar que os erros de truncamento somem mais do que 0,5: neste caso acabaríamos arredondando a resposta final para o inteiro errado. Esta possibilidade desastrosa pode ser evitada tomando d pequeno (e portanto grau grande, o que implica em uma transformada de Fourier de comprimento maior); também ajuda muito tomar o conjunto dos algarismos aceitáveis simétrico em relação ao zero, pois assim os produtos a_j b_k-jserão menores e terão sinais diferentes, o que evita que os coeficientes c_k sejam grandes demais. Mesmo para inteiros bem maiores do que o maior primo conhecido existem valores de d que garantem o bom funcionamento deste método, um dos mais rápidos para multiplicar inteiros grandes (em parte porque a maioria dos computadores é capaz de multiplicar doubles com grande rapidez). Por isso, ele é usado pelo programa mprime-prime95, que encontrou os últimos 4 primos de Mersenne. Escrevemos um pequeno programa que usa o critério de Lucas-Lehmer para testar a primalidade de M_p e que multiplica usando FFT: as funções FFT estão em fft.c e o programa principal em fftest2.c.

Uma segunda alternativa é escolher um primo pe fazer a multiplicação de polinômios considerando os coeficientes como elementos de ${\mathbb{Z} }/(p)$. Para recuperarmos os verdadeiros coeficientes do produto (que são inteiros), precisamos ter o cuidado de garantir que |c_k| < p/2 onde $c_k = \sum a_j b_{k-j}$. Um primo usado em alguns programas ^4.2 é p = 2⁶⁴ - 2³² + 1, que tem aliás várias propriedades especiais que o tornam particularmente apropriado para nossa tarefa. Com este valor de p, como 2³² | p-1, podemos fazer FFTs de comprimento 2³² com d = 2¹⁶, o que permite (em princípio) multiplicar inteiros de módulo menor do que $2^{16 \cdot 2^{32} - 1}$, ou seja, inteiros com alguns bilhões de algarismos; o simples armazenamento de um tal inteiro exige memória maior do que a que tem a maioria dos computadores atuais.

Mas estas alternativas, apesar de computacionalmente atraentes, não satisfazem ao matemático puro pois funcionam para inteiros menores do que um certo tamanho fixo (apesar de muito grande). A segunda alternativa apresentada acima pode ser levada adiante tomando primos cada vez maiores, mas não será fácil provar que existem sempre primos com as propriedades desejadas. Veremos agora como multiplicar inteiros de tamanho arbitrário em tempo baixo fazendo as contas não em ${\mathbb{Z} }/(p)$, mas em ${\mathbb{Z} }/(2^K + 1)$ (mesmo 2^K + 1 não sendo primo) e assim evitaremos esta dificuldade. Uma outra vantagem deste método é que será muito fácil multiplicar por potências de $\omega$ (assim tornando rápidas as FFTs).

Mais precisamente, mostraremos como multiplicar inteiros (dados por suas expansões binárias) módulo 2^N + 1; esta é a versão simplificada de Schönhage de um algoritmo devido a Schönhage e Strassen. Se N for tomado suficientemente grande este algoritmo multiplica inteiros. Consideraremos apenas valores de $N \ge 320$ da forma

$$N = \nu \cdot 2^n, \quad n - 1 \le \nu \le 2n, \quad n \ge 4;$$

estes valores de N serão chamados de aceitáveis. Supomos que já sabemos multiplicar módulo 2^K + 1, onde $K = \kappa \cdot 2^k < N$ também é um valor aceitável (a ser escolhido).

Para usar a multiplicação de polinômios, escrevemos os inteiros a e ba serem multiplicados na base d, i.e.,

$$a = \sum_{0 \le i < 2^m - 1} a_i d^i, \quad b = \sum_{0 \le j < 2^m - 1} b_j d^j, \quad 0 \le a_i, b_j < d,$$

onde $m = \lfloor n/2 \rfloor + 1$ e d = 2^N/2m. Temos $d^{2^m} = 2^N \equiv -1 \pmod{2^N + 1}$. Assim, podemos escrever $c \equiv ab \pmod{2^N + 1}$ com

$$c = \sum_{0 \le \mu < 2^m - 1} b_\mu d^\mu, \quad c_\mu = \sum_{i+j = \mu} a_i b_j - \sum_{i+j = \mu + 2^m} a_i b_j.$$

Pela seção anterior e por indução, sabemos efetuar estas contas módulo 2^K + 1mas novamente precisamos do valor de cada $c_\mu$ como inteiro, ou seja, precisamos escolher K de tal forma que possamos garantir que $\vert c_\mu\vert \le 2^{K-1}$. É fácil verificar que podemos escolher $\kappa = \lceil (\nu + 1)/2 \rceil$e $k = \lceil n/2 \rceil + 1$; observe que $K = \kappa \cdot 2^k$ é de fato aceitável.

Sejam $\tilde\omega = 2^{K/2^m}$ e $\omega = \tilde\omega^2$. Como $\omega^{2^{m-1}} \equiv -1 \pmod{2^K + 1}$temos que $\omega$ é uma raiz da unidade em ${\mathbb{Z} }/(2^K + 1)$ de ordem 2^m: este valor de $\omega$ pode ser usado para fazer FFT como na seção anterior; temos

$$c_\mu \equiv \tilde\omega^{-\mu} \sum_{i+j \equiv \mu \pmod {2^m}}% (\tilde\omega^i a_i)(\tilde\omega^j b_j) \pmod {2^K + 1}.$$

Note que podemos efetuar tanto FFT quanto FFT inversa pois 2^m e $\omega^i - 1$ são inversíveis módulo 2^K + 1(o que deixamos como exercício).

Falta apenas estimar o número de operações gasto por este algoritmo; note que por operação aqui queremos dizer uma operação sobre bits. Em todo o algoritmo, efetuamos duas FFTs de comprimento 2^m sobre ${\mathbb{Z} }/(2^K + 1)$, 2^m multiplicações ponto a ponto (também sobre ${\mathbb{Z} }/(2^K + 1)$) e uma FFT inversa de comprimento 2^m. Observe que como $\omega$ é uma potência de 2, as multiplicações por potências de $\omega$ que ocorrem nas FFTs são rápidas pois são apenas translações dos algarismos; mais precisamente, exigem no máximo CK operações cada uma (para alguma constante positiva C). Assim, cada FFT exige no máximo $C m \cdot 2^m K$ operações. O número total T(N) de operações satisfaz assim a recorrência

$$T(N) \le 2^m T(K) + C m \cdot 2^m K$$

donde podemos demonstrar que, para alguma constante positiva C,

$$T(n) \le C N \log N \log\log N.$$