---
tags: fft, ntt, fwt, dft, polynomials, algebra, number-theory, combinatorics, convolution, interpolation, primitive-roots, generators
---
<style>
.markdown-body{
max-width: 1000px;
}
</style>
# FFT parte 1
## Polinomios
Un polinomio es una función de la forma $P(x) = a_0 + a_1x + a_2x^2 + \cdots + a_nx^n$. A los valores $a_0, a_1, \ldots$ los llamamos *coeficientes* y al entero no negativo finito $n$ el *grado*, es por eso que $a_n \neq 0$. Los polinomios son las funciones más básicas que puede haber, y nos van a servir para varios problemas de combinatoria debido a sus propiedades.
:::info
:information_source: **Ejemplo:**
- $P(x)=5+3x+2x^2$ es un polinomio de grado 2.
- $Q(x)=-2x+4-7x^3+\frac{5}{7}x^2$ es un polinomio de grado 3.
- $R(x) = 2 + 4x + 7x^{3/2}$ no es un polinomio.
- $e^x = 1 + x + \frac{1}{2}x^2 + \frac{1}{6}x^3 + \cdots$ no es un polinomio.
:::
Si $P(x)$ es de grado $n$ se considera que todos sus coeficientes a partir $a_{n+1}$ son cero, es decir, $a_{n+1}=a_{n+2}=\cdots=0$. Es por eso que para hallar el grado, solo hay que fijarnos en el coeficiente de la potencia más grande de $x$ que no sea cero.
### Representación en C++
La forma más cómoda de representar polinomios en algún lenguaje de programación es usando un arreglo: el coeficiente de $x^k$ se va a guardar en la posición $k$ de dicho arreglo. Por lo tanto, vemos que si $P(x)$ es de grado $n$, el tamaño del arreglo tiene que ser de tamaño $n+1$.
Es por eso que usualmente de forma matemática un polinomio también se representa en forma vectorial o de tupla como $P(x)=(a_0, a_1, \ldots, a_n)$.
:::info
:information_source: **Ejemplo:** el polinomio $P(x)=6-3x+12x^2-4x^3$ se puede representar en C++ como ``vector<int> P = {6, -3, 12, -4}``, y también como $P=(6,-3,12,-4)$.
:::
### Suma y resta de polinomios
Si tenemos dos polinomios $P(x)$ y $Q(x)$, hallar la suma $P(x)+Q(x)$ únicamente implica sumar los coeficientes correspondientes de cada potencia de $x$, y de forma análoga con la resta.
Más formalmente:
- Sea $\displaystyle P(x) = \sum_{k=0}^{m} p_k x^k$ de grado $m$.
- Sea $\displaystyle Q(x) = \sum_{k=0}^{n} q_k x^k$ de grado $n$.
Entonces la suma $R(x)=P(x)+Q(x)$ será otro polinomio de grado $\max(m,n)$, el cuál está dado por:
$$
R(x) = \sum_{k=0}^{\max(m,n)} (\underbrace{p_k + q_k}_{r_k}) x^k
$$
La expresión anterior igual funciona si los grados no son iguales, pues el polinomio de grado menor se "rellena" con ceros.
:::info
:information_source: **Ejemplo:** hallar la suma de $P(x)=7+2x+5x^2$ con $Q(x)=-4+3x-6x^2+12x^3$.
$R(x)=3+5x-x^2+12x^3$
:::
Una posible implementación quedaría como:
```cpp=
using poly = vector<T>; // T es un tipo de dato genérico, tal como int, long long o un número complejo
poly operator+(const poly& P, const poly& Q){
int m = P.size(), n = Q.size();
int g = max(m, n);
poly R(g);
for(int i = 0; i < g; i++){
if(i < m) R[i] += P[i];
if(i < n) R[i] += Q[i];
}
return R;
}
```
Y la resta queda de forma casi idéntica, solo cambiando `+` por `-`. La implementación tiene complejidad de $O(\max(m, n))$.
### Evaluación en un punto
Dado un polinomio $P(x)$ y un valor específico $x_0$ queremos hallar $P(x_0)$. Podemos hacerlo de tres maneras:
- Siguiendo la definición de $P(x_0) = p_0 + p_1x_0 + p_2x_0^2 + \cdots p_nx_0^n$, hallamos cada $x_0^k$ usando exponenciación binaria:
```cpp=
T eval(const poly& P, T x0){
T ans = 0;
for(int k = 0; k < P.size(); ++k){
ans += P[k] * power(x0, k);
}
return ans;
}
```
con una complejidad de $O(n \log n)$.
- Podemos hallar $x_0^k$ del paso actual a partir de $x_0^{k-1}$ del paso anterior al multiplicarlo por $x_0$, o sea, $x_0^k = x_0^{k-1} \cdot x_0$. Al principio comenzamos con $x_0=1$ y lo vamos actualizando:
```cpp=
T eval(const poly& P, T x0){
T ans = 0;
T x_k = 1;
for(int k = 0; k < P.size(); ++k){
ans += P[k] * x_k;
x_k = x_k * x0;
}
return ans;
}
```
Y así logramos una complejidad de $O(n)$.
- Usando el método de Horner podemos reducir el número de multiplicaciones al evaluar de reversa el polinomio. Está basado en la propiedad $P(x_0) = p_0 + p_1x_0 + p_2x_0^2 + \cdots + p_nx_0^n = p_0 + x_0 \left( p_1 + x_0 \left( p_2 + x_0 \left( p_3 + \cdots + x_0 (p_{n-1} + x_0 p_n) \right) \right) \right)$
```cpp=
T eval(const poly& P, T x0){
T ans = 0;
int n = P.size();
for(int i = n-1; i >= 0; i--){
ans *= x0;
ans += P[i];
}
return ans;
}
```
Con una complejidad de $O(n)$.
### Multiplicación de polinomios
Dados dos polinomios $A(x) = a_0 + a_1x + a_2x^2 + \cdots + a_mx^m$ y $B(x) = b_0 + b_1x + b_2x^2 + \cdots + b_nx^n$ de grados $m$ y $n$ respectivamente, queremos hallar su producto $A(x) * B(x)$ al que llamaremos $C(x)$.
La multiplicación clásica consiste en hacer uso repetido de la propiedad distributiva y combinar todos contra todos los coeficientes de cada polinomio, recordando que $x^i \cdot x^j = x^{i+j}$.
:::info
:information_source: **Ejemplo:** sean $A(x)=2+4x+3x^2+x^3$ y $B(x)=5+7x+2x^2+8x^3$. Entonces:
$C(x) = 10 + 34x + 47x^2 + 50x^3 + 45x^4 + 26x^5 + 8x^6$
:::
El grado de $C(x)$ es ahora $m+n$, ya que la máxima potencia de $x$ que podemos obtener del producto viene de $a_mx^m \cdot b_nx^n$.
Vemos que si al momento estamos multiplicando el monomio $a_i x^i$ con $b_j x^j$, obtenemos $a_i b_j x^{i+j}$, por lo que en el arreglo $C$ hay que sumarle $a_i b_j$ a la posición $i+j$. La implementación de esta idea queda como:
```cpp=
poly operator*(const poly& A, const poly& B){
int m = A.size(), n = B.size();
int sz = m+n-1; // grado = m+n-2, tamaño = m+n-1
poly C(sz);
for(int i = 0; i < m; ++i){
for(int j = 0; j < n; ++j){
C[i+j] += A[i]*B[j];
}
}
return C;
}
```
Cuya complejidad es $O(mn)$.
## Convolución lineal
Antes de intentar mejorar esta complejidad, interpretemos de otra forma la multiplicación de dos polinomios. Sea $C(x) = c_0 + c_1x + c_2x^2 + \cdots + c_{m+n}x^{m+n}$. El método anterior es muy similar a una cubeta, pues le va "aportando" términos a los coeficientes de $C$ pero no en orden, sino conforme va multiplicando los monomios.
Pero ahora hay que hacerlo en orden: ¿cómo hallamos un coeficiente específico, digamos $c_k$? Sabemos que el monomio es $c_k x^k$ y se obtuvo como resultado de multiplicar algunos monomios de $A$ y $B$ y sumarlos. De forma específica, todos los productos $a_i b_j$ aportarán a $c_k$ si $i+j=k$, o sea, todas las formas posibles de obtener $k$ como suma de un exponente de $x$ obtenido de $A$ y otro de $B$. Entonces, lo único que hay que hacer es iterar por todos los pares $(i,j)$ cuya suma sea $k$ y sumar $a_i b_j$, lo cual es equivalente a los pares $(i, k-i)$ y sumar $a_i b_{k-i}$, donde la $i$ comienza en $0$ y termina en $k$:
$$
c_k = \sum_{i=0}^{k} a_i b_{k-i}
$$
donde no olvidemos que asumimos que $a_i=0$ si $i>m$ y $b_i=0$ si $j>n$.
:::info
:information_source: **Ejemplo:** de nuevo sean $A(x)=2+4x+3x^2+x^3$ y $B(x)=5+7x+2x^2+8x^3$. Hallemos su producto coeficiente a coeficiente, sabemos que $C(x)$ tiene grado 6, por lo tanto:
- $c_0 = a_0 b_0 = (2)(5) = 10$
- $c_1 = a_0 b_1 + a_1 b_0 = (2)(7) + (4)(5) = 34$
- $c_2 = a_0 b_2 + a_1 b_1 + a_2 b_0 = (2)(2) + (4)(7) + (3)(5) = 47$
- $c_3 = a_0 b_3 + a_1 b_2 + a_2 b_1 + a_3 b_0 = (2)(8) + (4)(2) + (3)(7) + (1)(5) = 50$
- $c_4 = a_0 b_4 + a_1 b_3 + a_2 b_2 + a_3 b_1 + a_4 b_0 = (2)(0) + (4)(8) + (3)(2) + (1)(7) + (0)(5) = 45$
- ...
:::
Y la implementación queda como:
```cpp=
poly operator*(const poly& A, const poly& B){
int m = A.size(), n = B.size();
int sz = m+n-1;
poly C(sz);
for(int k = 0; k < sz; ++k){
// calcular C[k] de forma individual
for(int i = 0; i <= k; ++i){
if(i < m && k-i < n){
C[k] += A[i] * B[k-i];
}
}
}
return C;
}
```
con complejidad $O((n+m)^2)$. Todavía no hemos mejorado, pero interpretar así la multiplicación nos dará nuevas herramientas.
Si vemos los polinomios como vectores, decimos que $C$ es la *convolución lineal* o *acíclica* de $A$ con $B$, y escribimos $C = A * B$. Por lo tanto, multiplicar polinomios y realizar convoluciones son cosas equivalentes.
## Forma punto-valor
Existe una tercer forma de representar un polinomio: la *forma punto-valor*: podemos evaluar el polinomio $A(x)$ de grado $m$ en al menos $m+1$ valores distintos: $x_0, x_1, \ldots, x_m$ y quedarnos con los resultados guardados en un vector: $((x_0, A(x_0)), (x_1,A(x_1)), \ldots, (x_m, A(x_m)))$.
La elección de las $x$'s es arbitraria, la única condición es que sean todas distintas entre sí y al menos $m+1$. Esto es para que los coeficientes del polinomio se puedan recuperar sin ambigüedad: cada evaluación aporta una ecuación lineal involucrando a los coeficientes, y como hay $m+1$ coeficientes, necesitamos $m+1$ ecuaciones para que el sistema tenga solución única. Dicha recuperación se conoce como *interpolación*.
:::info
:information_source: **Ejemplo:** sea $A(x)=2+4x+3x^2+x^3$. Hallemos su forma punto-valor usando los valores $x=0,1,2,3$:
- $P(0)=2$
- $P(1)=10$
- $P(2)=30$
- $P(3)=68$
Entonces, la forma punto-valor es $((0,2), (1,10), (2,30), (3,68))$.
Y ahora recuperemos el polinomio original a partir de su forma punto-valor. Sea $A(x)=a_0 + a_1x + a_2x^2 + a_3x^3$, entonces tenemos el siguiente sistema de ecuaciones:
- $a_0 = 2$
- $a_0 + a_1 + a_2 + a_3 = 10$
- $a_0 + 2a_1 + 4a_2 + 8a_3 = 30$
- $a_0 + 3a_1 + 9a_2 + 27a_3 = 68$
Cuya solución es:
:::
Regresando al problema de la multiplicación, sea $C(x) = A(x) * B(x)$. Si ya conocemos la forma punto-valor de $A(x)$ y $B(x)$, ¿cuál será la forma punto-valor de $C(x)$ sin encontrar $C(x)$ directamente?
Primero vemos que $C(x)$ tiene grado $m+n$, por lo que necesitaremos evaluarlo en al menos $m+n+1$ valores distintos. Pero como $C(x)=A(x)B(x)$, entonces $C(x_k)=A(x_k)B(x_k)$, lo que nos dice que hay que hallar la forma punto-valor tanto de $A(x)$ como de $B(x)$ en $m+n+1$ valores (los mismos para ambos polinomios) y multiplicarlos uno a uno para hallar la forma punto-valor de $C(x)$.
:::info
:information_source: **Ejemplo:** sean $A(x)=2+4x+3x^2+x^3$ y $B(x)=5+7x+2x^2+8x^3$. $C(x)$ tendrá grado 6, entonces necesitaremos 7 valores para hallar la forma punto-valor de $C(x)$. Escojamos $x=0,1,2,3,4,5,6$:
- $C(0) = 2 \times 5 = 10$
- $C(1) = 10 \times 22 = 220$
- $C(2) = 30 \times 91 = 2730$
- $C(3) = 68 \times 260 = 17680$
- $C(4) = 130 \times 577 = 75010$
- $C(5) = 222 \times 1090 = 241980$
- $C(6) = 350 \times 1847 = 646450$
Y ahora regresemos $C(x)$ a su forma normal:
:::
La complejidad de evaluar $A(x)$ en los valores escogidos es $O(m(m+n))$ y evaluar $B(x)$ es $O(n(m+n))$ usando el método de Horner. Pero hallar las evaluaciones de $C(x)$ cuesta $O(m+n)$, por lo tanto en la forma punto-valor multiplicar es muy sencillo y en tiempo lineal. Sin embargo, para ir de la forma normal a la punto-valor necesitamos evaluar y de regreso necesitamos interpolar, lo que usualmente cuesta $O((m+n)^2)$ usando un método de Gauss-Jordan optimizado, el cual es la interpolación de Lagrange.
Un panorama general de lo que hemos visto hasta ahora:
Por lo que si encontramos una forma más eficiente de evaluar e interpolar, habremos logrado reducir la complejidad. La clave está en escoger con cuidado los valores en los cuáles evaluar.
## Convolución circular
Supongamos ahora que tenemos dos polinomios $A(x)$ y $B(x)$, pero ahora **ambos** de grado $n-1$, o sea, de **tamaño $n$**. Es decir, en su representación vectorial:
$$
\begin{align}
A &= (a_0, a_1, \ldots, a_{n-1}) \\
B &= (b_0, b_1, \ldots, b_{n-1})
\end{align}
$$
La convolución *circular* o *cíclica* de $A$ con $B$ se define como: $C = A * B = (c_0, c_1, c_2, \ldots, c_{n-1})$, donde: $\displaystyle c_k = \sum_{i=0}^{n-1} a_ib_{k-i}$, donde los índices de $a_i$ y $b_i$ siempre se toman módulo $n$. Es por esto que esta convolución tiene una naturaleza cíclica, pues a diferencia de la convolución lineal en donde hay ceros más allá del coeficiente más significativo, aquí se ciclan módulo $n$. Este comportamiento se da tanto en las entradas $A$ y $B$ como en la salida $C$.
:::info
:information_source: **Ejemplo:** sean de nuevo $A(x)=2+4x+3x^2+x^3$ y $B(x)=5+7x+2x^2+8x^3$. Hallemos su convolución circular:
- $c_0 = a_0 b_0 + a_1 b_{-1} + a_2 b_{-2} + a_3 b_{-3} = a_0 b_0 + a_1 b_3 + a_2 b_2 + a_3 b_1 = (2)(5) + (4)(8) + (3)(2) + (1)(7) = 55$
- $c_1 = a_0 b_1 + a_1 b_0 + a_2 b_3 + a_3 b_2 = (2)(7) + (4)(5) + (3)(8) + (1)(2) = 60$
- $c_2 = a_0 b_2 + a_1 b_1 + a_2 b_0 + a_3 b_3 = (2)(2) + (4)(7) + (3)(5) + (1)(8) = 55$
- $c_3 = a_0 b_3 + a_1 b_2 + a_2 b_1 + a_3 b_0 = (2)(8) + (4)(2) + (3)(7) + (1)(5) = 50$
Por lo tanto, $C=(55, 60, 55, 50)=55+60x+55x^2+50x^3$.
:::
La relación entre ambas convoluciones es que podemos obtener la circular a partir de la lineal, pero no siempre en sentido contrario. Como los índices de los coeficientes están en módulo $n$, también lo estarán los exponentes de $x$. Es decir, $x^n=1$, $x^{n+1}=x$, $x^{n+2}=x^2$ y así sucesivamente. Haciendo esos reemplazos en la convolución lineal se puede obtener la convolución circular.
:::info
:information_source: **Ejemplo:** Ya sabíamos que la convolución lineal de $A$ con $B$ es $A*B=10+34x+47x^2+50x^3+45x^4+26x^5+8x^6$.
Por lo tanto, la circular es igual a $A * B = (10+45) + (34+26)x + (47+8)x^2 + 50x^3 = 55+60x+55x^2+50x^3$.
:::
Hacerlo a la inversa no siempre es posible, pues la convolución circular tiende a traslapar las potencias de $x$ al hacer el módulo, sin embargo resulta que el método que estudiaremos en breve tiene una mejor complejidad pero solo funciona para convoluciones circulares.
¿Cómo podemos obtener la convolución lineal a partir de la circular? Vamos a rellenar con ceros $A$ y $B$ para que no haya traslape en las potencias de $x$. Si $A$ tiene tamaño $m$ y $B$ tiene tamaño $n$, entonces hay que rellenar con ceros hasta que ambos tengan el tamaño $m+n-1$.
## DFT (Discrete Fourier Transform)
### Raíces de la unidad
Resulta que los valores ideales para evaluar los polinomios y hacer más eficiente la multiplicación son las *raíces de la unidad*. Estas existen en varios campos, tanto en los números complejos como en los enteros módulo $p$.
Una raíz $n$-ésima de la unidad, donde $n$ es un entero positivo, es un número $z$ que satisface la ecuación $z^n=1$. Esta ecuación tendrá exactamente $n$ raíces distintas, por lo tanto hay $n$ raíces $n$-ésimas de la unidad.
Sea $\omega$ una raíz $n$-ésima de la unidad. Decimos que es *primitiva* si además cumple que $\omega^k \neq 1$ para toda $1 \leq k \leq n-1$. O sea, el mínimo exponente positivo al que hay que elevarla para que sea igual a $1$ es $n$. Además, todos los $n$ números $\omega^0$, $\omega^1$, $\omega^2$, $\ldots$, $\omega^{n-1}$ serán distintos y será un conjunto válido de raíces $n$-ésimas de la unidad.
Para cualquier raíz $n$-ésima de la unidad podemos tomar su exponente módulo $n$, es decir, $\omega^a = \omega^{a \mod n}$. Esto es debido a su naturaleza cíclica, es decir, si $a = nq+r$, entonces $\omega^{a} = \omega^{nq+r} = (\omega^n)^q \omega^r = \omega^r$.
Para la DFT usaremos números complejos. En los números complejos una raíz primitiva de la unidad puede ser $\omega = e^{2 \pi i / n} = \cos \frac{2\pi}{n} + i \sin \frac{2\pi}{n}$. Geométricamente, todas las raíces de la unidad graficadas en el plano complejo están en el perímetro del círculo unitario centrado en el origen, igualmente espaciadas respecto a su ángulo.
:::info
:information_source: **Ejemplo:** las raíces sextas de la unidad van a formar los vértices de un hexágono regular en el plano complejo:
Sea $\omega=e^{2 \pi i / 6}$ una raíz sexta primitiva de la unidad. En el dibujo de arriba se ve que al girarla de $60^\circ$ en $60^\circ$ en sentido antihorario eventualmente le pegamos a todas las raíces sextas, y como girar $60^\circ$ es equivalente a multiplicar por $\omega$, todas las raíces sextas son:
- $\omega^0=1$ (no giramos nada)
- $\omega^1$ (giramos $60^\circ$)
- $\omega^2$ (giramos $120^\circ$)
- $\omega^3$ (giramos $180^\circ$)
- $\omega^4$ (giramos $240^\circ$)
- $\omega^5$ (giramos $300^\circ$)
Al llegar a $\omega^6$ regresamos de nuevo al comienzo, esta es la interpretación geométrica de la periodicidad en los exponentes de $\omega$.
:::
### DFT
Entonces, definimos a la DFT de un polinomio/vector $A=(a_0, a_1, \ldots, a_{n-1})$ de tamaño $n$ como $\text{DFT}(A) = \hat{A} = (\hat{a_0}, \hat{a_1}, \ldots, \hat{a_{n-1}})$, donde:
$$
\hat{a_k} = \sum_{j=0}^{n-1} a_j (\omega^k)^j
$$
Si vemos al vector $A$ como un polinomio de grado $n-1$, o sea: $A(x) = a_0 + a_1x + \ldots + a_{n-1}x^{n-1}$, entonces $\text{DFT}(A) = (A(\omega^0), A(\omega^1), A(\omega^2), \ldots, A(\omega^{n-1}))$. Es decir, esta operación simplemente nos devuelve la forma punto-valor de $A(x)$ usando las raíces $n$-ésimas de la unidad.
:::info
:information_source: **Ejemplo:** hallemos la DFT de $A(x)=2+4x+3x^2+x^3$. Sea $\omega=e^{2\pi i/4}=i$, entonces:
- $\hat{a_0}=A(i^0)=A(1)=2+4+3+1=10$
- $\hat{a_1}=A(i^1)=A(i)=2+4i+3i^2+i^3=2+4i-3-i=-1+3i$
- $\hat{a_2}=A(i^2)=A(-1)=2-4+3-1=0$
- $\hat{a_3}=A(i^3)=A(-i)=2+4(-i)+3(-i)^2+(-i)^3=2-4i-3+i=-1-3i$
Por lo tanto, $\hat{A}=\text{DFT}(A)=(10,-1+3i,0,-1-3i)$.
:::
Una primera implementación queda como:
```cpp=
using T = complex<double>; // a partir de aquí T será un número complejo con componentes reales
using poly = vector<T>;
const double pi = acosl(-1);
poly DFT(const poly& A){
int n = A.size();
T w = polar(1.0, 2 * pi / n);
poly An(n);
for(int k = 0; k < n; k++){
for(int j = 0; j < n; j++){
An[k] += A[j] * pow(pow(w, k), j);
}
}
return An;
}
```
Al usar directamente la definición, calcular la DFT de $A$ nos cuesta $O(n^2)$, pues hallar un término nos cuesta $O(n)$ y hay $n$ términos.
----------
Ahora supongamos que tenemos el vector transformado $\hat{A}$ y queremos hallar el original $A$. Esta transformación siempre es invertible, y está dada por $\text{DFT}^{-1}\left(\hat{A}\right) = A$, donde:
$$
a_k = \dfrac{1}{n}\sum_{j=0}^{n-1} \hat{a_j} (\omega^{-k})^j
$$
Básicamente es casi lo mismo que la transformada directa, solo que ahora el exponente de $\omega$ es $-k$ y al final hay que dividir entre $n$.
:::spoiler **Demostración**
$$
\begin{align*}
a_k &= \dfrac{1}{n} \sum_{j=0}^{n-1} \hat{a_j} (\omega^{-k})^j \\
&= \dfrac{1}{n} \sum_{j=0}^{n-1} \left( \sum_{m=0}^{n-1} a_m (\omega^j)^m \right) (\omega^{-k})^j \\
&= \dfrac{1}{n} \sum_{m=0}^{n-1} \sum_{j=0}^{n-1} a_m (\omega^j)^{m-k} \\
&= \dfrac{1}{n} \sum_{m=0}^{n-1} a_m \sum_{j=0}^{n-1} (\omega^{m-k})^j \end{align*}
$$
Vemos que cuando $m=k$, $\omega^{m-k}=1$, entonces $\displaystyle \sum_{j=0}^{n-1} (\omega^{m-k})^j = \sum_{j=0}^{n-1} 1^j = \sum_{j=0}^{n-1} 1 = n$.
También vemos que si $m \neq k$, $\omega^{m-k} \neq 1$, entonces $\displaystyle \sum_{j=0}^{n-1} (\omega^{m-k})^j = \dfrac{(\omega^{m-k})^n - 1}{\omega^{m-k} - 1} = \dfrac{(\omega^n)^{m-k} - 1}{\omega^{m-k} - 1} = \dfrac{1^{m-k} - 1}{\omega^{m-k} - 1} = \dfrac{1 - 1}{\omega^{m-k} - 1} = 0$.
Finalmente concluimos que:
$$
\begin{align*} a_k &= \dfrac{1}{n} \sum_{m=0}^{n-1} a_m n \delta_{m,k} \\
&= \dfrac{1}{n} a_k n \\
&= a_k \quad \square
\end{align*}
$$
:::
Y una primera implementación queda como:
```cpp=
poly IDFT(const poly& A){
int n = A.size();
T w = polar(1.0, 2 * pi / n);
poly An(n);
for(int k = 0; k < n; k++){
for(int j = 0; j < n; j++){
An[k] += A[j] * pow(pow(w, -k), j);
}
An[k] /= n;
}
return An;
}
```
Si de nuevo vemos al vector $A$ como un polinomio $A(x)$ y al vector $\hat{A}$ como su forma punto-valor, esta operación de la transformada inversa nos devuelve los coeficientes a partir de esta forma punto-valor con las raíces $n$-ésimas de la unidad.
:::info
:information_source: **Ejemplo:** hallemos la DFT inversa de $\hat{A}=(10, -1+3i, 0, -1-3i)$.
- $a_0 = \frac{1}{4} [(10) + (-1+3i) + (0) + (-1-3i)] = \frac{1}{4}[8] = 2$
- $a_1 = \frac{1}{4} [(10) + (-1+3i)(-i) + (0)(-i)^2 + (-1-3i)(-i)^3] = \frac{1}{4} [10+i+3-i+3] = \frac{1}{4}[16] = 4$
- $a_2 = \frac{1}{4} [(10) + (-1+3i)(-1) + (0)(-1)^2 + (-1-3i)(-1)^3] = \frac{1}{4} [10+1-3i+1+3i] = 3$
- $a_3 = \frac{1}{4} [(10) + (-1+3i)(i) + (0)(i)^2 + (-1-3i)(i)^3] = \frac{1}{4} [10-i-3+i-3] = \frac{1}{4}[4] = 1$
Por lo tanto, $A = \text{DFT}^{-1}(\hat{A}) = (2,4,3,1)$, que es lo que ya esperábamos.
:::
De forma natural, la DFT asume que tanto la entrada como la salida son vectores periódicos, lo cual nos da una pista de que tiene que haber una relación con la convolución circular.
:::spoiler **Demostración de que la salida es periódica**
$$
\begin{align}
\hat{a}_{k+n} &= \sum_{j=0}^{n-1} a_j (\omega^{k+n})^j \\
&= \sum_{j=0}^{n-1} a_j (\omega^k \omega^n)^j \\
&= \sum_{j=0}^{n-1} a_j (\omega^k)^j \\
&= \hat{a_k}
\end{align}
$$
:::
## Teorema de la convolución
Si $\hat{A} = \text{DFT}(A)$ y $\hat{B} = \text{DFT}(B)$, entonces $\text{DFT}(A * B) = \hat{A} \cdot \hat{B}$, donde $\cdot$ es una multiplicación punto a punto (la misma que hacemos al multiplicar en forma punto-valor) y $*$ es una **convolución circular**, o sea, tanto $A$ como $B$ tienen el mismo tamaño $n$.
:::spoiler **Demostración**
Sea $C=(c_0, c_1, \ldots, c_{n-1})=A*B$, entonces:
$$
c_k = \sum_{i=0}^{n-1} a_i b_{k-i}
$$
Sea $\hat{C}=(\hat{c_0}, \hat{c_1}, \ldots, \hat{c_{n-1}})=\text{DFT}(C)$. Por un lado, por definición de DFT, tenemos que:
$$
\begin{align}
\hat{c_k} &= \sum_{m=0}^{n-1} c_m (\omega^k)^m \\
&= \sum_{m=0}^{n-1} \left( \sum_{i=0}^{n-1} a_i b_{m-i} \right) (\omega^k)^m
\end{align}
$$
Por otro lado, por la multiplicación punto a punto tenemos que:
$$
\begin{align}
\hat{c_k} &= \hat{a_k} \hat{b_k} \\
&= \left( \sum_{i=0}^{n-1} a_i (\omega^k)^i \right) \left( \sum_{j=0}^{n-1} b_j (\omega^k)^j \right) \\
&= \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} a_i b_j (\omega^k)^{i+j}
\end{align}
$$
Vemos que estamos haciendo una multiplicación de polinomios, pues la contribución $a_ib_j$ se la sumamos a la potencia $i+j$ de $w^k$. Hagamos algo muy similar a lo que hicimos en la convolución lineal para obtener en orden las potencias de $w^k$: sea $m=i+j$, entonces iteremos primero sobre $m$ y luego sobre $i$, y $j=m-i$:
$$
\begin{align}
\hat{c_k} &= \sum_{m=0}^{n-1} \left( \sum_{i=0}^{n-1} a_i b_{m-i} \right) (\omega^k)^m
\end{align}
$$
Obtuvimos el mismo valor de $\hat{c_k}$ de ambas formas, por lo tanto la demostración queda completa. $\square$
:::
Por lo tanto, despejamos la convolución deseada y obtenemos: $A * B = \text{DFT}^{-1}(\text{DFT}(A) \cdot \text{DFT}(B))$. Esto es la formalización de evaluar, multiplicar a través de la forma punto-valor e interpolar usando raíces de la unidad.
Como se vio anteriormente, es muy importante que $A$ y $B$ estén rellenos con los suficientes ceros si la convolución a realizar es lineal, de tal forma que se pueda recuperar de forma única el producto deseado $C = A * B$ y no ocurra el efecto indeseado de que varias potencias de $x$ se traslapen (efecto *aliasing*, *overlapping* o *wrap-around*), el cuál sí ocurre en la convolución circular.
:::info
:information_source: **Ejemplo completo con convolución circular:**
Hallemos la convolución circular de $A(x)=2+4x+3x^2+x^3$ y $B(x)=5+7x+2x^2+8x^3$.
Ya habíamos visto que $\hat{A} = (10, -1+3i, 0, -1-3i)$. Ahora hallemos $\hat{B}$:
- $\hat{b_0} = B(i^0) = B(1) = 5+7+2+8 = 22$
- $\hat{b_1} = B(i^1) = B(i) = 5+7i+2i^2+8i^3 = 5+7i-2-8i = 3-i$
- $\hat{b_2} = B(i^2) = B(-1) = 5-7+2-8 = -8$
- $\hat{b_3} = B(i^3) = B(-i) = 5+7(-i)+2(-i)^2+8(-i)^3 = 5-7i-2+8i=3+i$
Por lo tanto $\hat{B}=(22, 3-i, -8, 3+i)$.
Luego, para hallar $\hat{C}$, hay que hacer la multiplicación punto a punto:
- $\hat{c_0} = \hat{a_0} \hat{b_0} = (10)(22) = 220$
- $\hat{c_1} = \hat{a_1} \hat{b_1} = (-1+3i)(3-i) = 10i$
- $\hat{c_2} = \hat{a_2} \hat{b_2} = (0)(-8) = 0$
- $\hat{c_3} = \hat{a_3} \hat{b_3} = (-1-3i)(3+i) = -10i$
Finalmente hallamos $C$ tomando la DFT inversa de $\hat{C}$:
- $c_0 = \frac{1}{4} [(220) + (10i) + (0) + (-10i)] = \frac{1}{4} [220] = 55$
- $c_1 = \frac{1}{4} [(220) + (10i)(-i) + (0)(-i)^2 + (-10i)(-i)^3] = \frac{1}{4} [220+10+10] = 60$
- $c_2 = \frac{1}{4} [(220) + (10i)(-1) + (0)(-1)^2 + (-10i)(-1)^3] = \frac{1}{4} [220-10i+10i] = 55$
- $c_3 = \frac{1}{4} [(220) + (10i)(i) + (0)(i)^2 + (-10i)(i)^3] = \frac{1}{4} [220-10-10] = 50$
Y así obtenemos el resultado: $C=(55, 60, 55, 50)$, o sea, $C(x)=55+60x+55x^2+50x^3$.
:::
:::info
:information_source: **Ejemplo completo con convolución lineal:**
Hallemos la convolución lineal de $A(x)=2+4x+3x^2+x^3$ y $B(x)=5+7x+2x^2+8x^3$.
Como el grado de $C(x)$ es 6, necesitamos rellenar con ceros $A$ y $B$ hasta que tengan al menos 7 términos. Por lo tanto, $A=(2, 4, 3, 1, 0, 0, 0)$ y $B=(5, 7, 2, 8, 0, 0, 0)$. Comencemos hallando la DFT de ambos:
Luego hallemos la multiplicación punto a punto, que es $\hat{C}$:
Y finalmente hallemos $C$ usando la transformada inversa:
Por lo que $C(x)=10+34x+47x^2+50x^3+45x^4+26x^5+8x^6$, lo cual ya sabíamos.
:::
Implementación de la convolución lineal usando DFT y teorema de la convolución:
```cpp=
poly operator*(poly A, poly B){
int m = A.size(), n = B.size();
int sz = m+n-1; // tamaño de la convolución lineal
// rellenamos con ceros A y B
A.resize(sz);
B.resize(sz);
// transformamos A y B, o sea, los evaluamos en las raíces de la unidad
poly An = DFT(A);
poly Bn = DFT(B);
// hacemos la multiplicación punto a punto
poly Cn(sz);
for(int i = 0; i < sz; ++i){
Cn[i] = An[i] * Bn[i];
}
// recuperamos C a través de la transformada inversa, o sea, interpolamos usando raíces de la unidad
poly C = IDFT(Cn);
return C;
}
```
La complejidad sigue siendo $O((m+n)^2)$, lo único que resta por hacer es evaluar la DFT y la inversa de la DFT con una complejidad menor.
## Idea de la FFT recursiva
El objetivo será calcular $\text{DFT}(A)$ con complejidad $O(n \log_2 n)$, donde $A$ tiene tamaño $n$.
Para efectos de la convolución lineal podemos asumir que $n$ es una potencia de dos, pues si no lo es, rellenamos con ceros $A$ hasta que su tamaño sea potencia de dos, y esto no afectará el resultado. Esto nos permitirá proponer un algoritmo recursivo, usando divide y vencerás.
Sea $n=2^h$, entonces tendremos $h$ llamadas recursivas a la función. Por definición de DFT tenemos:
$$
\hat{a_k} = \sum_{j=0}^{n-1} a_j (\omega^k)^j
$$
para $0 \leq k \leq n-1$.
Dividamos esa suma en dos partes: una con los términos de índice par y otra con los impares. Cada suma tendrá la mitad de términos, es decir, $\frac{n}{2}=2^{h-1}$. Indexaremos los índices pares haciendo el cambio $j \to 2j$ y los impares con $j \to 2j+1$:
$$
\begin{align}
\hat{a_k} &= \sum_{j=0}^{n/2-1} a_{2j} (\omega^k)^{2j} + \sum_{j=0}^{n/2-1} a_{2j+1} (\omega^k)^{2j+1} \\
&= \sum_{j=0}^{n/2-1} a_{2j} (\omega^{2k})^{j} + \omega^k\sum_{j=0}^{n/2-1} a_{2j+1} (\omega^{2k})^{j}
\end{align}
$$
Ahora definamos $A_{even}=(a_0, a_2, \ldots, a_{n-2})$ y $A_{odd} = (a_1, a_3, \ldots, a_{n-1})$, ambos de tamaño $\frac{n}{2}$. Vemos que si $\omega$ es una raíz $n$-ésima primitiva de la unidad, entonces $\omega^2$ es una raíz $\frac{n}{2}$-ésima primitiva de la unidad. Entonces cada suma en la expresión anterior representa el $k$-ésimo elemento de la DFT de $A_{even}$ y $A_{odd}$ respectivamente, donde adicionalmente multiplicamos el segundo término por $\omega^k$.
Entonces:
$$
\hat{a_k} = \text{DFT}(A_{even})_k + \omega^k \text{DFT}(A_{odd})_k
$$
Pero hay un detalle, ambas DFT's solo contienen $\frac{n}{2}$ términos, pero nosotros necesitamos $n$ de ellos, pues $0 \leq k \leq n-1$. Sin embargo recordemos que la DFT es cíclica, entonces los índices de ambas DFT's se toman módulo $\frac{n}{2}$ permitiendo repetir la primera mitad de los términos y reusarlos para la segunda mitad. Es por esto que usamos las raíces de la unidad, ya que permiten realizar la recursión con solo la mitad de términos y combinarlos para obtener todos.
Por lo tanto, para la primera mitad de valores de $\hat{a_k}$ ($0 \leq k < \frac{n}{2}$) tenemos directamente que:
$$
\hat{a_k} = \text{DFT}(A_{even})_k + \omega^k \text{DFT}(A_{odd})_k
$$
Y para la segunda mitad de valores, los que indexaremos como $\hat{a}_{k+n/2}$, donde $0 \leq k < \frac{n}{2}$, tenemos que:
$$
\begin{align}
\hat{a}_{k+n/2} &= \text{DFT}(A_{even})_{k+n/2} + \omega^{k+n/2} \text{DFT}(A_{odd})_{k+n/2} \\
&= \text{DFT}(A_{even})_k + \omega^k \omega^{n/2} \text{DFT}(A_{odd})_k
\end{align}
$$
Sea $z=\omega^{n/2}$, entonces al elevar ambos lados al cuadrado tenemos que $z^2=\omega^n=1$, por lo que $z=\pm 1$, pero como $\omega$ es raíz $n$-ésima primitiva $\omega^{n/2}$ no puede ser $1$, por lo tanto $\omega^{n/2}=-1$, y:
$$
\hat{a}_{k+n/2} = \text{DFT}(A_{even})_k - \omega^k \text{DFT}(A_{odd})_k
$$
De aquí vemos que en un solo ciclo desde $k=0$ hasta $\frac{n}{2}-1$ podremos hallar todos los $n$ valores de $\hat{a_k}$, ya que en cada iteración hallamos tanto $\hat{a_k}$ como $\hat{a}_{k+n/2}$.
Únicamente nos queda ver el caso base, que se da cuando ya no podemos dividir a la mitad el número de elementos de $A$. Eso se da cuando $n=1$, y en ese caso, $\text{DFT}(A)=A$.
Al implementar la DFT así se le conoce como *Fast Fourier Transform* o *FFT*.
:::info
:information_source: **Ejemplo:**
:::
Sea $T(n)$ la complejidad de la función ``FFT``. Hace dos llamadas recursivas a sí misma con la mitad del tamaño, y combinar la respuesta toma tiempo lineal. Por lo tanto se cumple que $T(n) = 2T(n/2) + O(n)$, y por teorema maestro, $T(n) = O(n \log_2 n)$.
### FFT inversa
Recordemos que los únicos cambios para hallar la DFT inversa era cambiarle el exponente a $\omega$ y dividir todos los elementos entre $n$ al final. Por eso, en nuestra función ``FFT`` solo le cambiamos el exponente a $\omega$ y al final dividimos entre $2$ todos los elementos en cada llamada recursiva. ¿Por qué entre $2$ y no entre $n$? Recordemos que tenemos $h$ llamadas recursivas, entonces al final habremos dividido entre $2^h=n$.
### Implementación
La implementación, tanto de la FFT normal como la inversa en una sola función, queda como:
```cpp=
poly FFT(const poly& A, int inv = 1){ //inv=1 (normal) o inv=-1 (inversa)
int n = A.size(); //asumimos que n es potencia de 2
// caso base
if(n == 1) return A;
// dividir A en pares e impares
poly odd, even;
for(int i = 0; i < n; i++){
if(i&1) odd.push_back(A[i]);
else even.push_back(A[i]);
}
// transformar odd y even
odd = FFT(odd, inv);
even = FFT(even, inv);
// combinar las respuestas de mis hijos
T w = polar(1.0, 2*pi/n * inv);
poly ans(n);
for(int k = 0; k < n/2; k++){
ans[k] = even[k] + pow(w, k)*odd[k];
ans[k + n/2] = even[k] - pow(w, k)*odd[k];
}
if(inv == -1){
for(int i = 0; i < n; i++){
ans[i] /= 2;
}
}
return ans;
}
```
Con esto hemos logrado multiplicar dos polinomios de grado $m$ y $n$ en tiempo $O((m+n)\log_2(m+n))$, ahora solo hay que llamar a ``FFT`` en lugar de ``DFT`` y no olvidar que hay que rellenar con ceros ambos polinomios hasta que su tamaño sea la mínima potencia de 2 mayor o igual a $m+n-1$:
```cpp=
int nearestPowerOfTwo(int n){
int ans = 1;
while(ans < n) ans <<= 1;
return ans;
}
// convolución lineal en O(n log n)
poly operator*(poly A, poly B){
int m = A.size(), n = B.size();
int sz = m+n-1;
// rellenamos con ceros hasta la potencia de 2 más cercana
int SZ = nearestPowerOfTwo(sz);
A.resize(SZ);
B.resize(SZ);
// transformamos A y B
poly An = FFT(A, 1);
poly Bn = FFT(B, 1);
poly Cn(SZ);
// multiplicación punto a punto
for(int i = 0; i < SZ; ++i){
Cn[i] = An[i] * Bn[i];
}
// transformada inversa
poly C = FFT(Cn, -1);
C.resize(sz);
return C;
}
```
## FFT iterativa
Aunque la versión recursiva de la FFT funciona en $O(n \log_2 n)$, al hacerla iterativa tendremos dos ventajas: no existirá *overhead* por la recursión y la calcularemos *in-place*, es decir, sobre el mismo vector sin memoria extra.
Para ejemplificar, sea $A = (a_0, a_1, \ldots, a_{15})$ un vector de tamaño $16=2^4$. Sea $\omega$ una raíz primitiva dieciseisava de la unidad. Primero veamos cómo se van moviendo los elementos del arreglo conforme se hacen las llamadas recursivas:
```graphviz
digraph fft{
node [color=Red,shape=box]
edge [color=Blue, style=dashed]
layer00 [label="(a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9], a[10], a[11], a[12], a[13], a[14], a[15])"]
layer10 [label="(a[0], a[2], a[4], a[6], a[8], a[10], a[12], a[14])"]
layer11 [label="(a[1], a[3], a[5], a[7], a[9], a[11], a[13], a[15])"]
layer20 [label="(a[0], a[4], a[8], a[12])"]
layer21 [label="(a[2], a[6], a[10], a[14])"]
layer22 [label="(a[1], a[5], a[9], a[13])"]
layer23 [label="(a[3], a[7], a[11], a[15])"]
layer30 [label="(a[0], a[8])"]
layer31 [label="(a[4], a[12])"]
layer32 [label="(a[2], a[10])"]
layer33 [label="(a[6], a[14])"]
layer34 [label="(a[1], a[9])"]
layer35 [label="(a[5], a[13])"]
layer36 [label="(a[3], a[11])"]
layer37 [label="(a[7], a[15])"]
layer40 [label="(a[0])"]
layer41 [label="(a[8])"]
layer42 [label="(a[4])"]
layer43 [label="(a[12])"]
layer44 [label="(a[2])"]
layer45 [label="(a[10])"]
layer46 [label="(a[6])"]
layer47 [label="(a[14])"]
layer48 [label="(a[1])"]
layer49 [label="(a[9])"]
layer410 [label="(a[5])"]
layer411 [label="(a[13])"]
layer412 [label="(a[3])"]
layer413 [label="(a[11])"]
layer414 [label="(a[7])"]
layer415 [label="(a[15])"]
layer00 -> {layer10 layer11}
{layer10} -> {layer20 layer21}
{layer11} -> {layer22 layer23}
{layer20} -> {layer30 layer31}
{layer21} -> {layer32 layer33}
{layer22} -> {layer34 layer35}
{layer23} -> {layer36 layer37}
{layer30} -> {layer40 layer41}
{layer31} -> {layer42 layer43}
{layer32} -> {layer44 layer45}
{layer33} -> {layer46 layer47}
{layer34} -> {layer48 layer49}
{layer35} -> {layer410 layer411}
{layer36} -> {layer412 layer413}
{layer37} -> {layer414 layer415}
}
```
Si nos fijamos en la última capa recursiva, vemos que los elementos de $A$ se permutaron de una forma especial: en la $i$-ésima posición se encuentra $a[rev(i)]$, donde $rev(i)$ es el número $i$ con sus 4 bits de reversa.
Numeremos las capas de arriba hacia abajo, comenzando en la 0 y terminando en la 4. En cada capa tenemos un cierto número de bloques, cada bloque tiene un mismo tamaño y la raíz de la unidad usada es la misma:
| Capa | Número de bloques | Tamaño de bloque | Raíz de la unidad |
| - | - | - | - |
| 0 | 1 | 16 | Dieciseisava: $\omega$ |
| 1 | 2 | 8 | Octava: $\omega^2$ |
| 2 | 4 | 4 | Cuarta: $\omega^4$ |
| 3 | 8 | 2 | Cuadrada: $\omega^8$ |
| 4 | 16 | 1 | |
De forma general, si $n=2^h$, la $c$-ésima capa tendrá $2^c$ bloques, el tamaño de bloque será $2^{h-c}$ y la raíz de la unidad a usar será la $(2^{h-c})$-ésima.
El primer paso será permutar los elementos, asignemos $A \gets bitReversal(A)$, lo cual se puede hacer in-place. Después iremos de la penúltima capa hacia la primera:
```graphviz
digraph fft{
node [color=Red,shape=box]
edge [dir="back", color=Blue, style=dashed]
layer00 [label="(a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9], a[10], a[11], a[12], a[13], a[14], a[15])"]
layer10 [label="(a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7])"]
layer11 [label="(a[8], a[9], a[10], a[11], a[12], a[13], a[14], a[15])"]
layer20 [label="(a[0], a[1], a[2], a[3])"]
layer21 [label="(a[4], a[5], a[6], a[7])"]
layer22 [label="(a[8], a[9], a[10], a[11])"]
layer23 [label="(a[12], a[13], a[14], a[15])"]
layer30 [label="(a[0], a[1])"]
layer31 [label="(a[2], a[3])"]
layer32 [label="(a[4], a[5])"]
layer33 [label="(a[6], a[7])"]
layer34 [label="(a[8], a[9])"]
layer35 [label="(a[10], a[11])"]
layer36 [label="(a[12], a[13])"]
layer37 [label="(a[14], a[15])"]
layer40 [label="(a[0])"]
layer41 [label="(a[1])"]
layer42 [label="(a[2])"]
layer43 [label="(a[3])"]
layer44 [label="(a[4])"]
layer45 [label="(a[5])"]
layer46 [label="(a[6])"]
layer47 [label="(a[7])"]
layer48 [label="(a[8])"]
layer49 [label="(a[9])"]
layer410 [label="(a[10])"]
layer411 [label="(a[11])"]
layer412 [label="(a[12])"]
layer413 [label="(a[13])"]
layer414 [label="(a[14])"]
layer415 [label="(a[15])"]
layer00 -> {layer10 layer11}
{layer10} -> {layer20 layer21}
{layer11} -> {layer22 layer23}
{layer20} -> {layer30 layer31}
{layer21} -> {layer32 layer33}
{layer22} -> {layer34 layer35}
{layer23} -> {layer36 layer37}
{layer30} -> {layer40 layer41}
{layer31} -> {layer42 layer43}
{layer32} -> {layer44 layer45}
{layer33} -> {layer46 layer47}
{layer34} -> {layer48 layer49}
{layer35} -> {layer410 layer411}
{layer36} -> {layer412 layer413}
{layer37} -> {layer414 layer415}
}
```
- Vamos a iterar por todas las capas $c$ de abajo hacia arriba, donde $h-1 \geq c \geq 0$:
- Sea $sz \gets 2^{h-c}$.
- Sea $\omega$ una raíz $sz$-ésima de la unidad.
- Vamos a iterar por todos los bloques $i$ en orden, donde $0 \leq i \leq 2^c-1$:
- El $i$-ésimo bloque tiene tamaño $sz$ y su primer elemento comienza en la posición $sz \cdot i$. Por lo tanto, sea $start \gets sz \cdot i$. De esta forma sabemos que el bloque abarca las posiciones $[start, start+sz-1]$. Cada bloque tiene que transformarse con su DFT usando sus dos bloques hijos, que ya están transformados. Las posiciones del bloque izquierdo (even) son $[start, start+sz/2-1]$ y las del derecho (odd) son $[start+sz/2,start+sz-1]$. Usaremos de nuevo las ecuaciones de transformación, para $0 \leq j < sz/2$, donde ambas asignaciones se realizan en **paralelo**:
$$
\begin{align}
a_{start+j} &\gets a_{start+j} + \omega^j a_{start+sz/2+j} \\
a_{start+sz/2+j} &\gets a_{start+j} - \omega^j a_{start+sz/2+j}
\end{align}
$$
Notemos que en cada asignación solo se están involucrando dos posiciones diferentes, por lo que se puede hacer in-place con tan solo dos variables temporales.
Tenemos $h = \log_2(n)$ capas para procesar, y en cada capa procesamos cada elemento del arreglo, o sea, una capa tiene una complejidad de $O(n)$. Por lo tanto, la complejidad total es $O(n \log_2 n)$.
:::info
:information_source: **Ejemplo:** sea $A=(2,-8,-1,3,1,-7,9,6)$ de tamaño $8=2^3$. Hallemos $\text{DFT}(A)$ usando el algoritmo de la FFT iterativa.
Primero vamos a permutar las posiciones de $A$ reflejando los 3 bits de cada posición:
- $a_0$ se intercambia con $a_0$, pues $0=000_2$.
- $a_1$ se intercambia con $a_4$, pues $1=001_2$ y $4=100_2$.
- $a_2$ se intercambia con $a_2$, pues $2=010_2$.
- $a_3$ se intercambia con $a_6$, pues $3=011_2$ y $6=110_2$.
- $a_5$ se intercambia con $a_5$, pues $5=101_2$.
- $a_7$ se intercambia con $a_7$, pues $7=111_2$.
Una forma de visualizarlo es usando un diagrama en donde las capas ahora las veremos de izquierda a derecha, conocido como *butterfly diagram*:
:::
La implementación queda como:
```cpp=
int rev(int i, int h){
int r = 0;
while(h--){
r = (r<<1) | (i&1);
i >>= 1;
}
return r;
}
void FFT(poly& A, int inv = 1){
int n = A.size(); // n debe ser potencia de 2
int h = __builtin_ctz(n); //log2(n)
// permutar A de acuerdo al bit reversal de sus posiciones
for(int i = 0; i < n; i++){
int r = rev(i, h);
if(i < r) swap(A[i], A[r]);
}
// iteramos por todas las capas de abajo hacia arriba
for(int c = h-1; c >= 0; c--){
int sz = 1<<(h-c);
int blocks = 1<<c;
T w = polar(1.0, 2*pi*inv/sz);
// iteramos por todos los bloques de la capa actual
for(int i = 0; i < blocks; i++){
int start = sz*i;
// transformar el bloque actual
for(int j = 0; j < sz/2; j++){
T u = A[start + j];
T v = pow(w,j) * A[start + sz/2 + j];
A[start + j] = u + v;
A[start + sz/2 + j] = u - v;
}
}
}
if(inv == -1){
for(int i = 0; i < n; i++){
A[i] /= n;
}
}
}
```
Optimizaciones adicionales:
- ``pow`` es muy costoso, entonces una mejor opción es en cada capa precalcular todas potencias de $\omega$ a usar, entre la línea 22 y 23 del código anterior:
```cpp=
poly w(sz/2);
for(int j = 0; j < sz/2; j++){
w[j] = polar(1.0, 2*pi*j*inv/sz);
}
```
Y reemplazar `pow(w,j)` por `w[j]`. De hecho, calcular las raíces justo de esta forma minimiza el error numérico, por lo que la FFT es un algoritmo numéricamente estable al usar números reales.
- Cuando en un mismo programa se van a hacer varias llamadas a ``FFT``, todas con el mismo tamaño del vector, conviene precalcular el bit reversal y las raíces de la unidad desde antes.
- El orden de los dos ciclos for anidados (los que iteran en `i` y el `j`) así como está ya es óptimo, pues si se intercambian habrá muchos *cache-misses*, ya que estaremos barriendo el zig-zag los elementos del arreglo. Pero así como están anidados, los elementos del arreglo se van barriendo lo más ordenadamente posible.
- Implementación más corta (omitiendo comentarios, llaves y algunas variables de un solo uso): https://github.com/alaneos777/reference/blob/master/fft.cpp
## Referencias
- https://cp-algorithms.com/algebra/fft.html
- https://neerc.ifmo.ru/trains/toulouse/2017/fft.pdf
- http://neerc.ifmo.ru/trains/toulouse/2017/fft2.pdf
- https://drive.google.com/file/d/1B9BIfATnI_qL6rYiE5hY9bh20SMVmHZ7/view
- https://www.youtube.com/watch?v=h7apO7q16V0
- https://www.youtube.com/watch?v=iTMn0Kt18tg
Sign in with Wallet
Connect another wallet