Transformada Rápida de Fourier - CP Algoritmos

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

CP Algoritmos
		
		
		
			
				
				
			
			Transformada rápida de Fourier
Neste artigo, discutiremos um algoritmo que nos permite multiplicar dois polinômios de tamanho nn em tempo O(nlogn)O(nlog⁡n), o que é melhor do que a multiplicação trivial que leva tempo O(n2)O(n2).
A multiplicação de dois números longos também pode ser reduzida a uma 
multiplicação de polinômios, portanto, dois números longos também podem 
ser multiplicados em tempo O(nlogn)O(nlog⁡n) (onde nn é a quantidade de dígitos nos números).
A descoberta da transformada rápida de Fourier (FFT)
 é atribuída a Cooley e Tukey, que publicaram um algoritmo em 1965.
Mas, na verdade, a FFT foi descoberta antes, mas sua importância não foi
 entendida antes das invenções dos computadores modernos.
Alguns pesquisadores atribuem a descoberta da FFT a Runge e König em 
1924.
Mas, na verdade, Gauss havia desenvolvido esse método em 1805, mas nunca
 o publicou.
Transformação discreta de Fourier
Seja um polinômio de grau n−1n−1:
A(x)=a0x0+a1x1+⋯+an−1xn−1A(x)=a0x0+a1x1+⋯+an−1xn−1
Assumimos que nn - o número de coeficientes - é uma potência de 22.
Se nn não for uma potência de 22, basta adicionar os termos ausentes aixiaixi e definir os coeficientes aiai para 00.
A teoria dos números complexos nos diz que a equação xn=1xn=1 tem nn soluções complexas, e as soluções são da forma wn,k=e2kπinwn,k=e2kπin com k=0…n−1k=0…n−1.
Além disso, esses números complexos têm propriedades interessantes: 
por exemplo, a principal nn-ésima raiz wn=wn,1=e2πinwn=wn,1=e2πin pode ser usado para descrever todas as outras nn-ésimas raízes da forma: wn,k=(wn)kwn,k=(wn)k.
A transformada discreta de Fourier (DFT) do polinômio A(x)A(x) (ou equivalentemente o vetor dos coeficientes (a0,a1,…,an−1)(a0,a1,…,an−1) é definida como os valores do polinômio nos pontos x=wn,kx=wn,k, ou seja, é o vetor:
DFT(a0,a1,…,an−1)=(y0,y1,…,yn−1)=(A(wn,0),A(wn,1),…,A(wn,n−1))=(A(w0n),A(w1n),…,A(wn−1n))DFT(a0,a1,…,an−1)=(y0,y1,…,yn−1)=(A(wn,0),A(wn,1),…,A(wn,n−1))=(A(wn0),A(wn1),…,A(wnn−1))
Similarmente a transformada discreta inversa de Fourier é definada como:
A DFT inversa dos valores do polinômio (y0,y1,…,yn−1)(y0,y1,…,yn−1) são os coeficientes do polinômio (a0,a1,…,an−1)(a0,a1,…,an−1).
InverseDFT(y0,y1,…,yn−1)=(a0,a1,…,an−1)InverseDFT(y0,y1,…,yn−1)=(a0,a1,…,an−1)
Assim, se uma DFT calcula os valores do polinômio nos pontos da nn-ésimas raízes, a DFT inversa pode restaurar os coeficientes do polinômio usando esses valores.
Aplicação da DFT: multiplicação rápida de polinômios
Sejam dois polinômios AA e BB.
Calculamos a DFT para cada um: DFT(A)DFT(A) e DFT(B)DFT(B).
O que acontece se multiplicarmos esses polinômios?
em cada ponto, os valores são simplesmente multiplicados, isto é
(A⋅B)(x)=A(x)⋅B(x).(A⋅B)(x)=A(x)⋅B(x).
Isso significa que, se multiplicarmos os vetores DFT(A)DFT(A) e DFT(B)DFT(B)
 - multiplicando cada elemento de um vetor pelo elemento correspondente 
do outro vetor - não obtemos nada diferente da DFT do polinômio DFT(A⋅B)DFT(A⋅B):
DFT(A⋅B)=DFT(A)⋅DFT(B)DFT(A⋅B)=DFT(A)⋅DFT(B)
Finalmente, aplicando a DFT inversa, obtemos:
A⋅B=InverseDFT(DFT(A)⋅DFT(B))A⋅B=InverseDFT(DFT(A)⋅DFT(B))
À direita, o produto das duas DFTs significa o produto em pares dos elementos dos vetores.
Isso pode ser calculado em tempo O(n)O(n).
Se pudermos calcular a DFT e a DFT inversa em O(nlogn)O(nlog⁡n),
 poderemos calcular o produto dos dois polinômios (e consequentemente 
também dois números longos) com a mesma complexidade de tempo.
Deve-se notar que os dois polinômios devem ter o mesmo grau.
Caso contrário, os doivetores resultantes da DFT irão ter tamanhos diferentes.
Podemos fazer isso adicionando coeficientes com o valor 00.
E também, como o resultado do produto de dois polinômios é um polinômio de grau 2(n−1)2(n−1), temos que dobrar os graus de cada polinômio (novamente preenchendo com 00s).
A partir de um vetor com nn valores não podemos reconstruir o polinômio desejado com 2n−12n−1 coeficientes.
Transformada rápida de Fourier
A transformada rápida de Fourier é um método que permite calcular a DFT em tempo O(nlogn)O(nlog⁡n).
A idéia básica da FFT é aplicar dividir e conquistar.
Dividimos o vetor coeficiente do polinômio em dois vetores, 
recursivamente calculamos a DFT para cada um deles, e combinamos os 
resultados para calcular a DFT do polinômio completo.
Assim, exista um polinômio A(x)A(x) com grau n−1n−1, onde nn é uma potência de 22, e n>1n>1:
A(x)=a0x0+a1x1+⋯+an−1xn−1A(x)=a0x0+a1x1+⋯+an−1xn−1
Dividimos em dois polinômios menores, o que contém apenas os 
coeficientes das posições pares e o que contém os coeficientes das 
posições ímpares:
A0(x)A1(x)=a0x0+a2x1+⋯+an−2xn2−1=a1x0+a3x1+⋯+an−1xn2−1A0(x)=a0x0+a2x1+⋯+an−2xn2−1A1(x)=a1x0+a3x1+⋯+an−1xn2−1
Vejamos que
A(x)=A0(x2)+xA1(x2).A(x)=A0(x2)+xA1(x2).
Os polinômios A0A0 e A1A1 são apenas metade dos coeficientes que o polinômio AA é.
Se pudermos calcular a DFT(A)DFT(A) em tempo linear usando a DFT(A0)DFT(A0) e DFT(A1)DFT(A1), obteremos a recorrência TDFT(n)=2TDFT(n2)+O(n)TDFT(n)=2TDFT(n2)+O(n) para a complexidade de tempo, o que resulta em TDFT(n)=O(nlogn)TDFT(n)=O(nlog⁡n).
Como podemos conseguir isso.
Suponha que tenhamos calculado os vetores (y0k)n/2−1k=0=DFT(A0)(yk0)k=0n/2−1=DFT(A0) and (y1k)n/2−1k=0=DFT(A1)(yk1)k=0n/2−1=DFT(A1).
Vamos encontrar uma expressão para (yk)n−1k=0=DFT(A)(yk)k=0n−1=DFT(A).
Para os primeiros valores de n2n2 podemos apenas usar a equação anterior A(x)=A0(x2)+xA1(x2)A(x)=A0(x2)+xA1(x2):
yk=y0k+wkny1k,k=0…n2−1.yk=yk0+wnkyk1,k=0…n2−1.
No entanto, para os segundos valores de n2n2 precisamos encontrar uma expressão diferente:
yk+n/2=A(wk+n/2n)=A0(w2k+nn)+wk+n/2nA1(w2k+nn)=A0(w2knwnn)+wknwn/2nA1(w2knwnn)=A0(w2kn)−wknA1(w2kn)=y0k−wkny1kyk+n/2=A(wnk+n/2)=A0(wn2k+n)+wnk+n/2A1(wn2k+n)=A0(wn2kwnn)+wnkwnn/2A1(wn2kwnn)=A0(wn2k)−wnkA1(wn2k)=yk0−wnkyk1
Aqui nós usado novamente A(x)=A0(x2)+xA1(x2)A(x)=A0(x2)+xA1(x2) e as duas identidades wnn=1wnn=1 e wn/2n=−1wnn/2=−1.
Portanto, obtemos as fórmulas desejadas para calcular todo o vetor (yk)(yk):
ykyk+n/2=y0k+wkny1k,=y0k−wkny1k,k=0…n2−1,k=0…n2−1.yk=yk0+wnkyk1,k=0…n2−1,yk+n/2=yk0−wnkyk1,k=0…n2−1.
(Esse padrão a+ba+b e a−ba−b é as vezes chamado de borboleta.)
Assim, aprendemos como calcular a DFT em tempo O(nlogn)O(nlog⁡n).
FFT Inversa
Seja o vetor (y0,y1,…yn−1)(y0,y1,…yn−1) - os valores do polinômio AA de grau n−1n−1 nos pontos x=wknx=wnk.
Queremos restaurar os coeficientes (a0,a1,…,an−1)(a0,a1,…,an−1) do polinômio.
Esse problema é chamado de interpolação, e existem 
algoritmos gerais para resolvê-lo.
Mas nesse caso (já que conhecemos os valores dos pontos nas raízes), 
podemos obter um algoritmo mais simples (que é praticamente o mesmo que o
 FFT).
Podemos escrever a DFT, de acordo com sua definição, na forma de matriz:
⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜w0nw0nw0nw0n⋮w0nw0nw1nw2nw3n⋮wn−1nw0nw2nw4nw6n⋮w2(n−1)nw0nw3nw6nw9n⋮w3(n−1)n⋯⋯⋯⋯⋱⋯w0nwn−1nw2(n−1)nw3(n−1)n⋮w(n−1)(n−1)n⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜a0a1a2a3⋮an−1⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜y0y1y2y3⋮yn−1⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟(wn0wn0wn0wn0⋯wn0wn0wn1wn2wn3⋯wnn−1wn0wn2wn4wn6⋯wn2(n−1)wn0wn3wn6wn9⋯wn3(n−1)⋮⋮⋮⋮⋱⋮wn0wnn−1wn2(n−1)wn3(n−1)⋯wn(n−1)(n−1))(a0a1a2a3⋮an−1)=(y0y1y2y3⋮yn−1)
Essa matriz é chamada de matriz de Vandermonde.
Assim, podemos calcular o vetor (a0,a1,…,an−1)(a0,a1,…,an−1) multiplicando o vetor (y0,y1,…yn−1)(y0,y1,…yn−1) da esquerda com o inverso da matriz:
⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜a0a1a2a3⋮an−1⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜w0nw0nw0nw0n⋮w0nw0nw1nw2nw3n⋮wn−1nw0nw2nw4nw6n⋮w2(n−1)nw0nw3nw6nw9n⋮w3(n−1)n⋯⋯⋯⋯⋱⋯w0nwn−1nw2(n−1)nw3(n−1)n⋮w(n−1)(n−1)n⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟−1⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜y0y1y2y3⋮yn−1⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟(a0a1a2a3⋮an−1)=(wn0wn0wn0wn0⋯wn0wn0wn1wn2wn3⋯wnn−1wn0wn2wn4wn6⋯wn2(n−1)wn0wn3wn6wn9⋯wn3(n−1)⋮⋮⋮⋮⋱⋮wn0wnn−1wn2(n−1)wn3(n−1)⋯wn(n−1)(n−1))−1(y0y1y2y3⋮yn−1)
A matriz inversa tem a seguinte forma:
1n⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜w0nw0nw0nw0n⋮w0nw0nw−1nw−2nw−3n⋮w−(n−1)nw0nw−2nw−4nw−6n⋮w−2(n−1)nw0nw−3nw−6nw−9n⋮w−3(n−1)n⋯⋯⋯⋯⋱⋯w0nw−(n−1)nw−2(n−1)nw−3(n−1)n⋮w−(n−1)(n−1)n⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟1n(wn0wn0wn0wn0⋯wn0wn0wn−1wn−2wn−3⋯wn−(n−1)wn0wn−2wn−4wn−6⋯wn−2(n−1)wn0wn−3wn−6wn−9⋯wn−3(n−1)⋮⋮⋮⋮⋱⋮wn0wn−(n−1)wn−2(n−1)wn−3(n−1)⋯wn−(n−1)(n−1))
Assim, obtemos a fórmula:
ak=1n∑j=0n−1yjw−kjnak=1n∑j=0n−1yjwn−kj
Comparando isso com a fórmula para ykyk
yk=∑j=0n−1ajwkjn,yk=∑j=0n−1ajwnkj,
observamos que esses problemas são quase os mesmos, portanto os coeficientes akak podem ser encontrados pelo mesmo algoritmo de divisão e conquista, bem como pela FFT direta, apenas em vez de wknwnk precisamos usar w−knwn−k, e no final precisamos dividir os coeficientes resultantes por nn.
Assim, o cálculo da DFT inversa é quase o mesmo que o cálculo da DFT direta, e também pode ser realizado em tempo O(nlogn)O(nlog⁡n).
Implementação
Aqui, apresentamos uma implementação recursiva da FFT
 direta e inversa, ambas em uma função, uma vez que a diferença entre a 
FFT direta e a inversa é tão mínima. Para armazenar os números 
complexos, usamos o tipo complexo da STL do C++.
using cd = complex<double>;
const double PI = acos(-1);
void fft(vector<cd> & a, bool invert) {
 int n = a.size();
 if (n == 1)
 return;
 vector<cd> a0(n / 2), a1(n / 2);
 for (int i = 0; 2 * i < n; i++) {
 a0[i] = a[2*i];
 a1[i] = a[2*i+1];
 }
 fft(a0, invert);
 fft(a1, invert);
 double ang = 2 * PI / n * (invert ? -1 : 1);
 cd w(1), wn(cos(ang), sin(ang));
 for (int i = 0; 2 * i < n; i++) {
 a[i] = a0[i] + w * a1[i];
 a[i + n/2] = a0[i] - w * a1[i];
 if (invert) {
 a[i] /= 2;
 a[i + n/2] /= 2;
 }
 w *= wn;
 }
}
A função recebe um vetor de coeficientes e calcula a DFT ou a DFT inversa e armazena o resultado novamente nesse vetor.
O argumento invertinvert mostra qual DFT(direta ou a inversa) deve ser calculada.
Dentro da função, primeiro verificamos se o comprimento do vetor é igual a um; se esse for o caso, não precisamos fazer nada.
Caso contrário, dividimos o vetor aa em dois vetores a0a0 e a1a1 e calculamos a DFT recursivamente para ambos.
Em seguida, inicializamos o valor wnwn e uma variável ww, que conterá a potência atual de wnwn.
Em seguida, os valores da DFT resultante são calculados usando as fórmulas.
Se o sinalizador invertinvert estiver definido, substituímos wnwn por wn−1wn−1, e cada um dos valores do resultado será dividido por 22 (como isso será feito em cada nível da recursão, então acabará dividindo os valores finais por nn).
Usando esta função, podemos criar uma função para multiplicar dois polinômios:
vector<int> multiply(vector<int> const& a, vector<int> const& b) {
 vector<cd> fa(a.begin(), a.end()), fb(b.begin(), b.end());
 int n = 1;
 while (n < a.size() + b.size()) 
 n <<= 1;
 fa.resize(n);
 fb.resize(n);
 fft(fa, false);
 fft(fb, false);
 for (int i = 0; i < n; i++)
 fa[i] *= fb[i];
 fft(fa, true);
 vector<int> result(n);
 for (int i = 0; i < n; i++)
 result[i] = round(fa[i].real());
 return result;
}
Essa função trabalha com polinômios com coeficientes inteiros, no 
entanto, você também pode ajustá-la para trabalhar com outros tipos.
Como existe algum erro ao trabalhar com números complexos, precisamos 
arredondar os coeficientes resultantes no final.
Finalmente, a função de multiplicar dois números longos praticamente não difere da função de multiplicar polinômios.
A única coisa que precisamos fazer depois é "normalizar" o número:
 int carry = 0;
 for (int i = 0; i < n; i++)
 result[i] += carry;
 carry = result[i] / 10;
 result[i] %= 10;
 }
Como o comprimento do produto de dois números nunca excede o 
comprimento total de ambos os números, o tamanho do vetor é suficiente 
para executar todas as operações.
Implementação aprimorada
Para aumentar a eficiência, passaremos da implementação recursiva 
para a iterativa. Na implementação recursiva acima, separamos 
explicitamente o vetor aa
 em dois vetores - os elementos nas posi;'oes pares foram colocados em 
um vetor temporário assim como os elementos nas posições impares.
No entanto, se reordenarmos os elementos de uma certa maneira, não 
precisamos criar esses vetores temporários (ou seja, todos os cálculos 
podem ser feitos no próprio vetor AA).
Observe que, no primeiro nível de recursão, os elementos cujo bit mais baixo da posição era zero foram atribuídos ao vetor a0a0, e os que tinham um como o bit na posição mais baixa foram colocados no vetor a1a1.
No segundo nível de recursão, acontece o mesmo, mas com o segundo bit 
mais baixo, etc.
Portanto, se invertermos os bits da posição de cada coeficiente e os 
ordenarmos por esses valores invertidos, obteremos a ordem desejada ( 
chamada de permutação de reversão de bits).
Por exemplo, a ordem desejada para n=8n=8 tem a forma:
a={[(a0,a4),(a2,a6)],[(a1,a5),(a3,a7)]}a={[(a0,a4),(a2,a6)],[(a1,a5),(a3,a7)]}
No primeiro nível de recursão (cercado por chaves), o vetor é dividido em duas partes [a0,a2,a4,a6][a0,a2,a4,a6] e [a1,a3,a5,a7][a1,a3,a5,a7].
Na permutação de inversão de bits, isso corresponde simplesmente à divisão do vetor em duas metades: os primeiros n2n2 elementos e os últimos n2n2
 elementos.
Depois, há uma chamada recursiva para cada metade.
Deixe a DFT resultante para cada um deles retornar no lugar dos próprios
 elementos (ou seja, a primeira metade e a segunda metade do vetor aa respectivamente.
a={[y00,y01,y02,y03],[y10,y11,y12,y13]}a={[y00,y10,y20,y30],[y01,y11,y21,y31]}
Agora, queremos combinar as duas DFTs em uma para o vetor completo.
A ordem dos elementos é ideal, mas também podemos realizar a união diretamente nesse vetor.
Podemos pegar os elementos y00y00 e y10y01 executar a transformação de borboleta.
O local dos dois valores resultantes é igual ao local dos dois valores iniciais, portanto, obtemos:
a={[y00+w0ny10,y01,y02,y03],[y00−w0ny10,y11,y12,y13]}a={[y00+wn0y01,y10,y20,y30],[y00−wn0y01,y11,y21,y31]}
Da mesma forma, podemos calcular a transformação de borboleta de y01y10 e y11y11 e colocar os resultados no lugar deles, e assim por diante.
Como resultado, obtemos:
a={[y00+w0ny10,y01+w1ny11,y02+w2ny12,y03+w3ny13],[y00−w0ny10,y01−w1ny11,y02−w2ny12,y03−w3ny13]}a={[y00+wn0y01,y10+wn1y11,y20+wn2y21,y30+wn3y31],[y00−wn0y01,y10−wn1y11,y20−wn2y21,y30−wn3y31]}
Assim, calculamos a DFT necessária a partir do vetor aa.
Aqui, descrevemos o processo de computação da DFT apenas no primeiro 
nível de recursão, mas o mesmo funciona obviamente também para todos os 
outros níveis. Assim, após aplicar a permutação de reversão de bits, podemos calcular a DFT no mesmo vetor, sem nenhuma memória adicional.
Além disso, isso permite nos livrar da recursão. Apenas começamos no 
nível mais baixo, ou seja, dividimos o vetor em pares e aplicamos a 
transformação de borboleta neles.
Isso resulta com o vetor aa com o trabalho do último nível aplicado.
Na próxima etapa, dividimos o vetor em vetores de tamanho 44, e novamente aplicamos a transformação de borboleta, que nos fornece a DFT para cada bloco de tamanho 44.
E assim por diante.
Finalmente, na última etapa, obtemos o resultado das DFTs de ambas as metades de aa, e, aplicando a transformação de borboleta, obtemos o DFT para o vetor completo aa.
using cd = complex<double>;
const double PI = acos(-1);
int reverse(int num, int lg_n) {
 int res = 0;
 for (int i = 0; i < lg_n; i++) {
 if (num & (1 << i))
 res |= 1 << (lg_n - 1 - i);
 }
 return res;
}
void fft(vector<cd> & a, bool invert) {
 int n = a.size();
 int lg_n = 0;
 while ((1 << lg_n) < n)
 lg_n++;
 for (int i = 0; i < n; i++) {
 if (i < reverse(i, lg_n))
 swap(a[i], a[reverse(i, lg_n)]);
 }
 for (int len = 2; len <= n; len <<= 1) {
 double ang = 2 * PI / len * (invert ? -1 : 1);
 cd wlen(cos(ang), sin(ang));
 for (int i = 0; i < n; i += len) {
 cd w(1);
 for (int j = 0; j < len / 2; j++) {
 cd u = a[i+j], v = a[i+j+len/2] * w;
 a[i+j] = u + v;
 a[i+j+len/2] = u - v;
 w *= wlen;
 }
 }
 }
 if (invert) {
 for (cd &
x : a)
 x /= n;
 }
}
Inicialmente, aplicamos a permutação de reversão de bits(bit-reversal
 permutation) trocando cada elemento pelo elemento da posição reversa. 
Então na logn−1log⁡n−1 fase do algoritmo calculamos a DFT para cada bloco de tamanho correspondente lenlen.
Para todos esses blocos, temos a mesma raiz de unidade wlenwlen.
Nós iteramos todos os blocos e executamos a transformação de borboleta em cada um deles.
Podemos otimizar ainda mais a reversão dos bits. Na implementação 
anterior, iteramos todos os bits do índice e criamos o índice reverso 
bit a bit. No entanto, podemos inverter os bits de uma maneira 
diferente.
Suponha que jj já contenha o reverso de ii.
Então, para ir para i+1i+1, temos que incrementar ii, e também precisamos incrementar jj,
 mas em um sistema numérico "invertido".
Adicionar um no sistema binário convencional é equivalente a inverter 
todos os ums juntos em zeros e inverter o zero antes deles em um.
Equivalentemente, no sistema numérico "invertido", invertemos todos os 
principais ums, e também o próximo zero.
Implementação:
using cd = complex<double>;
const double PI = acos(-1);
void fft(vector<cd> & a, bool invert) {
 int n = a.size();
 for (int i = 1, j = 0; i < n; i++) {
 int bit = n >> 1;
 for (; j & bit; bit >>= 1)
 j ^= bit;
 j ^= bit;
 if (i < j)
 swap(a[i], a[j]);
 }
 for (int len = 2; len <= n; len <<= 1) {
 double ang = 2 * PI / len * (invert ? -1 : 1);
 cd wlen(cos(ang), sin(ang));
 for (int i = 0; i < n; i += len) {
 cd w(1);
 for (int j = 0; j < len / 2; j++) {
 cd u = a[i+j], v = a[i+j+len/2] * w;
 a[i+j] = u + v;
 a[i+j+len/2] = u - v;
 w *= wlen;
 }
 }
 }
 if (invert) {
 for (cd & x : a)
 x /= n;
 }
}
Além disso, podemos pré-calcular previamente a permutação de reversão de bits.
Isso é especialmente útil quando o tamanho nn é o mesmo para todas as chamadas.
Mas mesmo quando temos apenas três chamadas (necessárias para multiplicar dois polinômios), o efeito é perceptível.
Também podemos pré-calcular todas as raízes da unidade e suas potências.
Transformação teórica dos números(number theoretic transform NTT)
Agora vamos mudar um pouco o objetivo. Ainda queremos multiplicar dois polinômios em tempo O(nlogn)O(nlog⁡n), mas desta vez queremos calcular os coeficientes módulo algum número primo pp.
Obviamente, para esta tarefa, podemos usar a DFT normal e aplicar o operador de módulo ao resultado.
No entanto, isso pode levar a erros, especialmente quando se lida com números grandes.
A transformação teórica dos números (NTT) funciona apenas com números inteiros e, portanto, é garantido que o resultado seja correto.
A transformada discreta de Fourier é baseada em números complexos e nas raizes nn-ésima
 de unidade.
Para calculá-la eficientemente, usamos extensivamente as propriedades 
das raízes (por exemplo: que existe uma raiz que gera todas as outras 
raízes por exponenciação).
Mas as mesmas propriedades são válidas para as raízes nn-ésimas de unidade na arithmética modular.
Uma nn-ésima raiz de unidade é um número wnwn que satisfaz:
(wn)n(wn)k=1(modp),≠1(modp),1≤k<n.(wn)n=1(modp),(wn)k≠1(modp),1≤k<n.
As outras n−1n−1 raízes podem ser obtidas como potências da raiz wnwn.
Para aplicar no algoritmo de transformação rápida de Fourier, precisamos que uma raiz exista para algum nn, na qual seja uma potência de 22, e também para todas as potências menores.
Podemos notar uma seguinte propriedade interessante:
(w2n)m=wnn(w2n)k=w2kn=1(modp),with m=n2≠1(modp),1≤k<m.(wn2)m=wnn=1(modp),with m=n2(wn2)k=wn2k≠1(modp),1≤k<m.
Então, se wnwn for uma nn-ésima raiz da unidade, então w2nwn2 é uma n2n2-ésima
 raiz.
E, conseqüentemente, para todas as potências menores que dois, existem 
raízes do grau exigido e elas podem ser calculadas usando wnwn.
Para calcular a DFT inversa, precisamos do w−1nwn−1 inverso de wnwn.
Mas, para um módulo primo, o inverso sempre existe.
Assim, todas as propriedades que precisamos das raízes complexas 
também estão disponíveis na aritmética modular, desde que tenhamos um 
módulo grande o suficiente pp no qual uma nn-ésima raiz da unidade exista.
Por exemplo, podemos usar os seguintes valores: módulo p=7340033p=7340033, w220=5w220=5.
Se este módulo não for suficiente, precisamos encontrar um par diferente.
Podemos usar esse fato de que, para os módulos da forma p=c2k+1p=c2k+1 (pp é primo), sempre existe 2k2k-ésima raiz da unidade.
Pode-se mostrar que gcgc é uma 2k2k-ésima raiz, onde gg é uma raiz primitiva de pp.
const int mod = 7340033;
const int root = 5;
const int root_1 = 4404020;
const int root_pw = 1 << 20;
void fft(vector<int> & a, bool invert) {
 int n = a.size();
 for (int i = 1, j = 0; i < n; i++) {
 int bit = n >> 1;
 for (; j & bit; bit >>= 1)
 j ^= bit;
 j ^= bit;
 if (i < j)
 swap(a[i], a[j]);
 }
 for (int len = 2; len <= n; len <<= 1) {
 int wlen = invert ? root_1 : root;
 for (int i = len; i < root_pw; i <<= 1)
 wlen = (int)(1LL * wlen * wlen % mod);
 for (int i = 0; i < n; i += len) {
 int w = 1;
 for (int j = 0; j < len / 2; j++) {
 int u = a[i+j], v = (int)(1LL * a[i+j+len/2] * w % mod);
 a[i+j] = u + v < mod ? u + v : u + v - mod;
 a[i+j+len/2] = u - v >= 0 ? u - v : u - v + mod;
 w = (int)(1LL * w * wlen % mod);
 }
 }
 }
 if (invert) {
 int n_1 = inverse(n, mod);
 for (int & x : a)
 x = (int)(1LL * x * n_1 % mod);
 }
}
A função inverse calcula o (see Inverso Modular).
As constantes mod, root, root_pw determinam o módulo e a raiz, e root_1 é o inverso de root módulo mod.
Na prática, essa implementação é mais lenta que a implementação 
usando números complexos (devido ao grande número de operações com 
módulo), mas possui algumas vantagens, como menos uso de memória e 
nenhum erro de arredondamento.
Multiplicação com módulo arbitrário
Aqui queremos alcançar o mesmo objetivo da seção anterior. Multiplicando dois polinômios A(x)A(x) e B(x)B(x), e calcular os coeficientes módulo algum número MM.
A transformação teórica dos números funciona apenas para determinados 
números primos, ela não funciona quando o módulo não possui a forma 
desejada.
Uma opção seria realizar transformações teóricas de múltiplos números com diferentes números primos da forma c2k+1c2k+1, e, em seguida, aplicar o teorema chinês do resto para calcular os coeficientes finais.
Outra opção é distribuir os polinômios A(x)A(x) e B(x)B(x) em dois polinômios menores
A(x)B(x)=A1(x)+A2(x)⋅C=B1(x)+B2(x)⋅CA(x)=A1(x)+A2(x)⋅CB(x)=B1(x)+B2(x)⋅C
with C≈M−−√C≈M.
Em seguida, o produto de A(x)A(x) e B(x)B(x) pode ser representado como:
A(x)⋅B(x)=A1(x)⋅B1(x)+(A1(x)⋅B2(x)+A2(x)⋅B1(x))⋅C+(A2(x)⋅B2(x))⋅C2A(x)⋅B(x)=A1(x)⋅B1(x)+(A1(x)⋅B2(x)+A2(x)⋅B1(x))⋅C+(A2(x)⋅B2(x))⋅C2
Os polinômios A1(x)A1(x), A2(x)A2(x), B1(x)B1(x) e B2(x)B2(x) contêm apenas coeficientes menores que M−−√M, portanto, os coeficientes de todos os produtos exibidos são menores que M⋅nM⋅n, que geralmente é pequeno o suficiente para lidar com tipos float.
Essa abordagem, portanto, exige calcular os produtos de polinômios 
com coeficientes menores (usando a FFT normal e a FFT inversa) e, em 
seguida, o produto original pode ser restaurado usando adição e 
multiplicação modular em tempo O(n)O(n).
Aplicações
A DFT pode ser usada em uma enorme variedade de outros problemas, que
 de primeira não têm nada a ver com a multiplicação de polinômios.
Todas as somas possíveis
São dadas duas arrays a[]a[] e b[]b[].
Temos que encontrar todas as somas possíveis a[i]+b[j]a[i]+b[j], e, para cada soma, conte com que frequência ela aparece.
Por exemplo, para a=[1, 2, 3]a=[1, 2, 3] e b=[2, 4]b=[2, 4], teremos:
a
soma 33 pode ser obtida em 11 maneira, a soma 44 também em 11 maneira, 55 em 22, 66 em 11, 77 em 11.
Construímos para as arrays aa e bb dois polinômios AA e BB.
Os números da array atuarão como expoentes no polinômio (a[i]⇒xa[i]a[i]⇒xa[i]); e os coeficientes desse termo serão com que frequência o número aparece na array.
Então, multiplicando esses dois polinômios em tempo O(nlogn)O(nlog⁡n), conseguimos o polinômio CC, onde os expoentes nos dizem quais somas podem ser obtidas e os coeficientes nos dizem com que frequência elas acontecem.
Para demonstrar:
(1x1+1x2+1x3)(1x2+1x4)=1x3+1x4+2x5+1x6+1x7(1x1+1x2+1x3)(1x2+1x4)=1x3+1x4+2x5+1x6+1x7
Todos os possíveis produtos escalares
São dadas duas arrays a[]a[] e b[]b[] de tamanho nn.
Precisamos calcular os produtos de aa a cada mudança(shift) cíclico de bb.
Geramos duas novas arrays de tamanho 2n2n:
Revertemos aa e adicionamos nn zeros.
E apenas adicionamos bb a ele mesmo.
Quando multiplicamos essas duas arrays como polinômios, e observamos os coeficientes c[n−1], c[n], c[2n−2]c[n−1], c[n], c[2n−2] do produto cc, obtemos:
c[k]=∑i+j=ka[i]b[j]c[k]=∑i+j=ka[i]b[j]
Como a[i]=0a[i]=0 para i≥ni≥n:
c[k]=∑i=0n−1a[i]b[k−i]c[k]=∑i=0n−1a[i]b[k−i]
Essa soma é o produto escalar do vetor aa com o (k−(n−1))(k−(n−1))-ésimo ciclo de bb.
Portanto, esses coeficientes são a resposta para o problema.
Nota: o c[2n−1]c[2n−1] nos dá o nn-ésimo ciclo mas esse é o mesmo que o ciclo 00, portanto, não precisamos considerar isso separadamente em nossa resposta.
Duas faixas(two stripes)
São dados duas faixas de booleanos (arrays cíclicas de valores 00 e 11) aa e bb.
Queremos encontrar todas as maneiras de anexar a primeira faixa à segunda, de modo que em nenhuma posição tenha 11 da primeira faixa ao lado de 11 da segunda faixa.
O problema não difere muito do problema anterior. Anexar duas faixas 
significa apenas que realizamos uma mudança cíclica na segunda array e 
podemos anexar as duas faixas, se o produto escalar das duas arrays for 00.
Coincidência de strings
Nos são dadas duas strings, um texto TT e um padrão PP, consistindo em letras minúsculas. Temos que calcular todas as ocorrências do padrão p no texto t.
Criamos um polinômio para cada string (T[i]T[i] e P[I]P[I] são número entre 00 e 2525 correspondentes às 2626 letras do alfabeto):
A(x)=a0x0+a1x1+⋯+an−1xn−1,n=|T|A(x)=a0x0+a1x1+⋯+an−1xn−1,n=|T|
com
ai=cos(αi)+isin(αi),αi=2πT[i]26.ai=cos⁡(αi)+isin⁡(αi),αi=2πT[i]26.
e
B(x)=b0x0+b1x1+⋯+bm−1xm−1,m=|P|B(x)=b0x0+b1x1+⋯+bm−1xm−1,m=|P|
com
bi=cos(βi)−isin(βi),βi=2πP[m−i−1]26.bi=cos⁡(βi)−isin⁡(βi),βi=2πP[m−i−1]26.
Observe que com a expressão P[m−i−1]P[m−i−1] inverte explicitamente o padrão.
O (m−1+i)(m−1+i)ésimo coeficiente do produto de dois polinômios C(x)=A(x)⋅B(x)C(x)=A(x)⋅B(x) nos dirá se o padrão aparece no texto na posição ii.
cm−1+i=∑j=0m−1ai+j⋅bm−1−j=∑j=0m−1(cos(αi+j)+isin(αi+j))⋅(cos(βj)−isin(βj))cm−1+i=∑j=0m−1ai+j⋅bm−1−j=∑j=0m−1(cos⁡(αi+j)+isin⁡(αi+j))⋅(cos⁡(βj)−isin⁡(βj))
com αi+j=2πT[i+j]26αi+j=2πT[i+j]26 and βj=2πP[j]26βj=2πP[j]26
Se eles corresponderem, então T[i+j]=P[j]T[i+j]=P[j], e, portanto, αi+j=βjαi+j=βj.
Isso nos dá (usando a identidade trigonométrica pitagórica):
cm−1+i=∑j=0m−1(cos(αi+j)+isin(αi+j))⋅(cos(αi+j)−isin(αi+j))=∑j=0m−1cos(αi+j)2+sin(αi+j)2=∑j=0m−11=mcm−1+i=∑j=0m−1(cos⁡(αi+j)+isin⁡(αi+j))⋅(cos⁡(αi+j)−isin⁡(αi+j))=∑j=0m−1cos⁡(αi+j)2+sin⁡(αi+j)2=∑j=0m−11=m
Se não houver uma correspondência, pelo menos um caractere é diferente, o que significa que um dos produtos ai+1⋅bm−1−jai+1⋅bm−1−j não é igual a 11, nos levando ao coeficiente cm−1+i≠mcm−1+i≠m.
Coincidência de strings com curingas
Esta é uma extensão do problema anterior. Desta vez, permitimos que o padrão contenha o caractere curinga ∗∗, que pode corresponder a todas as letras possíveis. Por exemplo, o padrão a∗ca∗c aparece no texto como abccaaccabccaacc em exatamente três posições, no índice 00, índice 44 e índice 55.
Criamos exatamente os mesmos polinômios, exceto que definimos bi=0bi=0 se P[m−i−1]=∗P[m−i−1]=∗.
Se xx é o número de curingas em PP, então teremos uma correspondência de PP em TT no índice ii se cm−1+i=m−xcm−1+i=m−x.
Problemas
	SPOJ - POLYMUL
	SPOJ - MAXMATCH
	SPOJ - ADAMATCH
	Codeforces - Yet Another String Matching Problem
	Codeforces - Lightsabers (hard)
	Kattis - K-Inversions
	Codeforces - Dasha and cyclic table
Referências
	transformada teórica dos números
			
 © 2019 traduzido por http://github.com/samuraiwesleyjack