【模板】快速傅里叶变换
P3803
给定一个 \(n\) 次多项式 \(F(x)\)，和一个 \(m\) 次多项式 \(G(x)\)。 请求出 \(F(x)\) 和 \(G(x)\) 的卷积。
\(n, m \leq {10}^6\)
快速傅里叶变换
FFT变换
\(n=2^k\)，已知\(f(x)=\sum_{i\in[0,2^k)}a_ix^i\)，求所有\(y_j = f(\omega_{2^k}^{j})\)。
做法：由于
\[y_j=f(\omega_{2^k}^{j})=\sum_{i\in[0,2^k)}a_i(\omega_{2^k}^{j})^i\\
y_{j+2^{k-1}}=f(\omega_{2^k}^{j+2^{k-1}})=\sum_{i\in[0,2^k)}a_i(\omega_{2^k}^{j+2^{k-1}})^i=\sum_{i\in[0,2^k)}a_i(-\omega_{2^k}^j)^i
\]
令 \(f_i(x) = \sum_{j\in[0,2^{k-1})}a_{2j+i}x^j\quad(i\in\{0,1\})\)：
\[y_j = f_0((\omega_{2^k}^j)^2)+\omega_{2^k}^j f_1((\omega_{2^k}^j)^2) = f_0(\omega_{2^{k-1}}^j)+\omega_{2^k}^j f_1(\omega_{2^{k-1}}^j) \\
y_{j+2^{k-1}} = f_0((-\omega_{2^k}^j)^2)-\omega_{2^k}^{j}f_1((-\omega_{2^k}^j)^2) = f_0(\omega_{2^{k-1}}^j)-\omega_{2^k}^j f_1(\omega_{2^{k-1}}^j)
\]
故只需求出所有\(y'_j=f_1(\omega_{2^{k-1}}^j),y''_j=f_2(\omega_{2^{k-1}}^j)\)即可。（成功划分成了两个子问题）
此时\(y_j=y'_j+\omega_{2^k}^jy''_j, y_{j+2^{k-1}}=y'_j-\omega_{2^k}^jy''_j\)。
FFT逆变换
假设已有\(y_j=\sum_{i\in[0,n)}a_i\omega_n^{ij}\)，要求\(a_i\)。
考虑对\(y_j\)再做一个反的FFT，得到：
\[\begin{aligned}
x_i &= \sum_{j\in[0,n)}y_j\omega_n^{-ij} \\
&= \sum_{j\in[0,n)}\sum_{k\in[0,n)}a_k\omega_n^{kj}\omega_n^{-ij} \\
&= \sum_{j\in[0,n)}\sum_{k\in[0,n)}a_k\omega_n^{(k-i)j} \\
&= \sum_{k\in[0,n)}a_k\sum_{j\in[0,n)}(\omega_n^{k-i})^j \\
&= na_{i}
\end{aligned}
\]
即得\(a_i=x_i/n\)
递归实现[1]
void fast_fast_tle(int limit, complex *a, int type) {
if (limit == 1) return ;
complex a1[limit >> 1], a2[limit >> 1];
for (int i = 0; i <= limit; i += 2)
a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
fast_fast_tle(limit >> 1, a1, type);
fast_fast_tle(limit >> 1, a2, type);
complex Wn = complex(cos(2.0 * Pi / limit) , type * sin(2.0 * Pi / limit)), w = complex(1, 0);
for (int i = 0; i < (limit >> 1); i++, w = w * Wn)
a[i] = a1[i] + w * a2[i],
a[i + (limit >> 1)] = a1[i] - w * a2[i];
}
非递归实现[1:1]
预处理r[i]（i的二进制表示的反转）：
for (int i = 0; i < limit; i++)
r[i] = ( r[i >> 1] >> 1 ) | ( (i & 1) << (l - 1) ) ;
代码：
void fast_fast_tle(complex *A, int type) {
for (int i = 0; i < limit; i++)
if (i < r[i]) swap(A[i], A[r[i]]); //r[r[i]] == i
for (int mid = 1; mid < limit; mid <<= 1) { //自底向上，mid表示区间长度的一半
complex Wn( cos(Pi / mid) , type * sin(Pi / mid) ); //mid<<1次单位根。
for (int R = mid << 1, j = 0; j < limit; j += R) { //R为区间长度，j为当前的最右端点
complex w(1, 0);
for (int k = 0; k < mid; k++, w = w * Wn) {
complex x = A[j + k], y = w * A[j + mid + k];
A[j + k] = x + y;
A[j + mid + k] = x - y;
}
}
}
}
有关complex
C++的complex<T>类
手写struct complex[1:2]
const double Pi = acos(-1.0);
struct complex {
double x, y;
complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
} a[MAXN], b[MAXN];
complex operator + (complex a, complex b) { return complex(a.x + b.x , a.y + b.y);}
complex operator - (complex a, complex b) { return complex(a.x - b.x , a.y - b.y);}
complex operator * (complex a, complex b) { return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);}
