快速傅里叶变换FFT
快速傅里叶变换FFT
用途
\(\operatorname{FFT}\)算法支持在\(O(n log n)\)时间内计算两个\(n\)度的多项式的乘法。也可以用来加速大整数乘法运算。
前置知识
系数表示法
用一个多项式的各项系数来表达这个多项式（升幂顺序）：
则它的计算公式为：
时间复杂度：\(O(n^2)\)
点值表示法
把一个多项式看成一个坐标系中的函数图像，从图像上选取\(n +1\)个点，利用这\(n + 1\)个点来唯一地表示这个函数。
为什么\(n + 1\)个点可以唯一地表示\(f(n)\)
知道了正确性，那么我们可以设：
那么用点值表示法表示\(f(x)\)如下：
它的计算公式为：
时间复杂度：\(O(n^2)\)
单位负根
定义
\(x^n = 1\)在复数意义下的解释\(n\)次负根。这样的解有\(n\)个。
设：\(\large \omega_n = e^{\frac{2 \pi i}{n}}\)
则：\(x^n\)的解集表示为：
称\(\omega_n\)是\(n\)次单位负根，其他负根均可以用单位负根的幂表示。
由于我数学太菜了所以不会证明。
根据欧拉公式可得：
柿子很好记，可是我不知道怎么证，可是记住就行了。
性质
单位负根。对于任意正整数\(n\)和整数\(k\)：
快速傅里叶变换FFT
\(\operatorname{FFT}\)的基本思想是分治。
先看\(\operatorname{DFT}\)，它分治地来求当\(x = \omega^k_n\)的时候\(f(x)\)的值，分治思想体现在将多项式分为奇次项和偶次项处理。
令\(f(x)\)为一个\(n\)次的多项式。
令\(G(x)\)表示偶次项系数建立的新函数：
令\(H(x)\)表示奇次项系数建立的新函数：
那么原来的\(f(x)\)用新的两个函数表示为：
利用单位负根的性质推柿子：
同理：
因此，我们求出\(\operatorname{DFT}(G(\omega^k_{\frac{n}{2}}))\)和\(DFT(H(\omega^k_{\frac{n}{2}}))\)后，就可以同时求出\(\operatorname{DFT}(f(\omega^k_n)\)和\(\operatorname{DFT}(f(\omega^{k + \frac{n}{2}}_n))\)，于是对\(G\)和\(H\)分别递归\(\operatorname {DFT}\)即可。
注意：考虑到分治\(\operatorname {DFT}\)能处理的多项式长度只能是\(2^m(m \in N^*)\)，否则在分治的时候左式和右式长度不等，右式就取不到系数了。所以要在第一次\(\operatorname{DFT}\)之前就把序列向上补成长度为\(2^m(m \in N^*)\)，最高项次数为\(2^m - 1\)的多项式（高次系数补0）。
在代入值的时候，因为要代入\(n\)个不同值，所以我们代入\(\omega^0_n,\omega^1_n,\omega^2_n\ldots,\omega^{n - 1}_n(n = 2^m(m\in N^*))\)一共\(2^m\)个不同值。
快速傅里叶逆变换
作用
把点值表示法转化为系数表示法。
方法与推导
考虑原本的多项式：
考虑构造法。
我们已知\(y_i = f(\omega^i_n),i\in\{0,1,\ldots,n - 1\}\)，要求\(\{a_0,a_1,\ldots,a_{n - 1}\}\)
构造多项式：
也就是把\(\{y_0,y_1,\ldots,y_{n - 1}\}\)当做多项式\(A\)的系数表示法。
设\(b_i = \omega^{-i}_n\)，则多项式\(A\)在\(x = b_0,b_1,\ldots,b_{n-1}\)处的点值表示法为：\(\{A(b_0),A(b_1),\ldots,A(b_{n - 1})\}\)
对\(A(x)\)的定义式做一下变换，可以将\(A(b_k)\)表示为：
记\(S( \omega^a_n)=\sum^{n-1}{i=0}(\omega^a_n)^i\)
当\(a=0(\bmod n)\)时，\(S(\omega^a_n) = n\)
当\(a\ne0(\bmod n)\)时，错位相减：
也就是：
代回原式得：
那么多项式\(A\)的点值表示法为：
总结：我们取单位根为其倒数，对\(\{y_0,y_1,\ldots,y_{n-1}\}\)跑一遍\(\operatorname{FFT}\)，然后除以\(n\)即可得到\(f(x)\)的系数表示。
代码实现：luogu P3803[模板]多项式乘法(FFT)
#include <bits/stdc++.h>
using namespace std;
//#define ll long long
#define ri register int
const int maxn = 5e6 + 5;
const double pi = acos(-1.0);
int n,m,l,r[maxn],lim = 1;
inline int read() {
int s = 0,w = 1; char ch = getchar();
while (ch < '0' || ch > '9') {if (ch == '-') w = -1; ch = getchar();}
while (ch >= '0' && ch <= '9') s = s * 10 + ch - '0',ch = getchar();
return s * w;
}
struct cpx {
double x,y;
cpx(double a = 0,double b = 0) {x = a; y = b;}
cpx operator + (const cpx &a) {return cpx(x + a.x,y + a.y);}
cpx operator - (const cpx &a) {return cpx(x - a.x,y - a.y);}
cpx operator * (const cpx &a) {return cpx(x * a.x - y * a.y,x * a.y + y * a.x);}
}a[maxn],b[maxn];
inline void fft(cpx *f,int type) {
for (ri i = 0;i < lim;i++) if (i < r[i]) swap(f[i],f[r[i]]);
for (ri mid = 1;mid < lim;mid <<= 1) {
cpx W(cos(pi / mid),type * sin(pi / mid));
int len = (mid << 1);
for (ri j = 0;j < lim;j += len) {
cpx pw(1,0);
for (ri k = 0;k < mid;k++,pw = pw * W) {
cpx x = f[j + k],y = pw * f[j + k + mid];
f[j + k] = x + y;
f[j + k + mid] = x - y;
}
}
}
}
int main() {
n = read(); m = read();
for (ri i = 0;i <= n;i++) a[i].x = read();
for (ri i = 0;i <= m;i++) b[i].x = read();
while (lim <= n + m) {lim <<= 1; l++;}
for (ri i = 0;i < lim;i++) r[i] = (r[i >> 1] >> 1) | (i & 1) << (l - 1);
fft(a,1); fft(b,1);
// for (int i = 0;i <= lim;i++) cout << a[i].x << " " << a[i].y << endl;
for (ri i = 0;i <= lim;i++) a[i] = a[i] * b[i];
fft(a,-1);
for (ri i = 0;i <= n + m;i++) printf("%d ",(int)(a[i].x / lim + 0.5));
return 0;
}
本文大部分内容来自：oi-wiki
浙公网安备 33010602011771号