Multiplicação de polinômios usando FFT

Suponha que queiramos multiplicar dois polinômios $P, Q \in {\mathbb{C} }[x]$, de grau menor do que n, representados pelos seus coeficientes:
$$\begin{align}P(x) &= \sum_{0 \le j < n} a_j x^j &= a_0 + a_1 x + \cdots + a_{n-1... ...\le j < n} b_j x^j &= b_0 + b_1 x + \cdots + b_{n-1} x^{n-1}. \notag \end{align}$$
O método aprendido na escola exige n² multiplicações; o métode de Karatsuba pode ser adaptado para este problema e exige aproximadamente $n^\alpha$ multiplicações, com $\alpha = (\log 3)/(\log 2)$. Veremos agora como efetuar esta multiplicação com um número muito menor de operações.

Uma idéia é a de considerar os polinômios como representados não pelos seus coeficientes e sim pelos seus valores em npontos distintos $\xi_0, \ldots, \xi_{n-1}$. Temos evidentemente $(P \cdot Q)(\xi_j) = P(\xi_j) \cdot Q(\xi_j)$: se o produto PQ tem grau menor do que n então PQé o único polinômio que assume estes n valores. A dificuldade em usar este método está em calcular os valores de P e Qnos n pontos $\xi_0, \ldots, \xi_{n-1}$e em recuperar PQ a partir de seu valor nestes mesmos pontos. Se os valores $\xi_j$ forem escolhidos sem critério este método pode acabar sendo mais lento do que os outros que já apresentamos. Veremos que certas escolhas de n e $\xi_j$ tornam o algoritmo rápido: uma das mais simples é tomar n uma potência de 2 e $\xi_j = \omega^j$, onde $\omega = e^{2\pi i/n}$ é uma raiz da unidade de ordem n.

Suponha que $\xi_k = - \xi_j \ne 0$ então as potências pares de $\xi_j$ e $\xi_k$ coincidem, enquanto as potências ímpares diferem pelo sinal. Isto nos permite economizar multiplicações quando calculamos $P(\xi_j)$ e $P(\xi_k) = P(-\xi_j)$ simultaneamente. Se n é par, podemos escrever
$$\begin{align}P(\xi_j) &= P_+(\xi_j^2) + \xi_j P_-(\xi_j^2),\notag\\ P(-\xi_j) &= P_+(\xi_j^2) - \xi_j P_-(\xi_j^2),\notag \end{align}$$
onde
$$\begin{align}P_+(x) &= \sum_{0 \le j < n/2} a_{2j} x^j,\notag\\ P_-(x) &= \sum_{0 \le j < n/2} a_{2j+1} x^j,\notag. \end{align}$$
Ou seja, reduzimos o problema de calcular um polinômio de grau nem dois pontos ao problemas de calcular dois polinômios de grau n/2em um mesmo ponto, seguido de uma multiplicação, uma soma e uma subtração. Se os $\xi_j$ sempre ocorrerem aos pares, com por exemplo $\xi_{j + (n/2)} = -\xi_j$, o cálculo de $P(\xi_0), \ldots, P(\xi_n)$ reduz-se ao cálculo de $P_+(\xi_0^2), \ldots, P_+(\xi_{(n/2) - 1}^2), P_-(\xi_0^2), \ldots, P_-(\xi_{(n/2) - 1}^2)$seguido de 3n/2 operações.

O ideal é que pudéssemos repetir o processo acima, ou seja, que n seja múltiplo de 4 e que também no conjunto $\xi_0^2, \ldots, \xi_{(n/2) - 1}^2$os números ocorressem em pares diferindo apenas por sinal. Reordenando os termos, podemos reformular esta condição como $\xi_{j + (n/4)}^2 = - \xi_j^2$, ou, sem perda de generalidade, como $\xi_{j + (n/4)} = i \xi_j$. Para podermos repetir este processo um número máximo de vezes, devemos tomar n como uma potência de 2 e $\xi_{j+k} = \omega^k \xi_j$, onde $\omega = e^{2\pi i/n}$. Devemos assim tomar $\xi_j = \omega^j \xi_0$e a escolha $\xi_0 = 1$ parece particularmente simples.

Façamos agora uma estimativa de T(n), o número de operações usadas neste algoritmo para calcular $P(\xi_0), \ldots, P(\xi_n)$. Já vimos que T(n) = 2 T(n/2) + 3n/2; claramente T(1) = 0. Daí temos T(2) = 3, T(4) = 12 e, por uma indução simples, $T(2^k) = 3 k \cdot 2^{k-1}$. Assim, é possível calcular $P(1), \ldots, P(\omega^{n-1})$ muito rápidamente.

Reformulemos este problema na linguagem de álgebra linear. Temos

$$\begin{pmatrix} P(1) \\ P(\omega) \\ P(\omega^2) \\ \vdots \\... ...pmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{pmatrix};$$

a matriz de ordem n e coeficientes $\omega^{ij}$, $0 \le i, j < n$, é chamada de DFT_n, a transformada de Fourier discreta. O que aprendemos nos parágrafos acima foi a multiplicar DFT_npor um vetor rapidamente (pelo menos quando n é uma potência de 2). Em termos algébricos, aprendemos a escrever DFT_n como um produto de $\log_2 n$ matrizes esparsas cujos coeficientes não nulos são potências de $\omega$; cada matriz esparsa correspondendo a uma etapa do algoritmo FFT.

Falta aprender a recuperar os coeficientes de um polinômio P a partir de

$$P(1), \ldots, P(\omega^{n-1}),$$

ou seja, a multiplicar (DFT_n)^-1 por um vetor. Mas para isto basta observar que

$$(DFT_n)^2 = n \cdot \begin{pmatrix} 1 & 0 & 0 & \ldots & 0 & ... ...& \ldots & 0 & 0 \\ 0 & 1 & 0 & \ldots & 0 & 0 \end{pmatrix},$$

pois o coeficiente (i,k) de (DFT_n)² é

$$\sum_{j} \omega^{ij} \omega^{jk} = \sum_j \omega^{(i+k)j}$$

que é igual a n se $i+k \equiv 0 \pmod n$ e 0 caso contrário pois

$$\omega^{\ell n} - 1 = (\omega^\ell - 1) \sum_j \omega^{\ell j}.$$

Assim, o coeficiente (i,j) de (DFT_n)^-1 é $(1/n) \omega^{-ij}$e este processo de FFT inversa (ou interpolação) é tão fácil e rápido quanto FFT (ou evaluação). Temos portanto um algoritmo para multiplicar polinômios de grau nfazendo aproximadamente $C n \log n$ operações (onde C é uma constante positiva).

Reproduzimos abaixo o pseudo-código de [CB] para este algoritmo:

Input: O comprimento n (uma potência de 2); uma raiz primitiva da unidade $\omega$ de ordem n; um vetor $(a_0,\ldots,a_{n-1})$ de coeficientes complexos.

Output: O vetor $(A_0, \ldots, A_{n-1})^t = (DFT_n) (a_0, \ldots, a_{n-1})^t$.

procedure $FFT(n,\omega,a_0,a_1,\ldots,a_{n-1};A_0,A_1,\ldots,A_{n-1})$; begin         if n = 1 then             A₀ = a₀;         else              $FFT(n/2,\omega^2,a_0,a_2,\ldots,a_{n-2};E_0,\ldots,E_{n/2 - 1})$;              $FFT(n/2,\omega^2,a_1,a_3,\ldots,a_{n-1};O_0,\ldots,O_{n/2 - 1})$;             for k = 0 to n/2 - 1 do                  $A_k = E_k + \omega^k O_k$;                  $A_{k + n/2} = E_k - \omega^k O_k$; end

Até agora consideramos polinômios com coeficientes em ${\mathbb{C} }$mas o leitor atento já deve ter percebido que podemos usar o mesmo algoritmo para multiplicar polinômios sobre qualquer corpo Kdesde que exista em K um elemento $\omega$ que seja uma raiz da unidade de ordem n. Um exemplo de corpo onde existe um tal $\omega$ é ${\mathbb{Z} }/(p)$se $p \equiv 1 \pmod n$. Na verdade não é sequer necessário que os coeficientes estejam em um corpo: podemos trabalhar sobre qualquer anel A onde exista $\omega$com as seguintes propriedades:

1.

$\omega^n = 1$,

2.

n é inversível em A,

3.

se $0 < \ell < n$ então $\omega^\ell - 1$ é inversível em A.

Na próxima seção veremos uma situação onde será interessante trabalhar com $A = {\mathbb{Z} }/(2^K + 1)$.

Lembramos que este algoritmo calcula corretamente o produto dos polinômios P e Q desde que este produto tenha grau menor do que n. Mais geralmente, estaremos encontrando o único polinômio de grau menor que n que coincide com PQ em $\xi_0, \xi_1, \ldots, \xi_{n-1}$. Como estamos tomando sempre $\xi_j = \omega^j \xi_0$ temos

$$(x - \xi_0)(x - \xi_1)\cdots(x - \xi_{n-1}) = x^n - \xi_0^n$$

e nosso algoritmo calcula $PQ \bmod (x^n - \xi_0^n)$.