FFT详解

用于计算一类朴素卷积问题。

笔者学习的是这份博客 内容中可能有很多相同之处，敬请谅解。

现在要计算两个一元\(n\)次多项式\(F(x)\)与\(G(x)\)的乘积，如何计算？前置知识：多项式的表示方法 一. 系数表示法 对于一个\(n\)次多项式\(F(x)\)，它可以被表示成\[F(x) = a_nx^n+a_{n-1}x^{n-1}+...+a_1x^1+a_0x^0.\] 更加形式化的来说，它可以表示成\[F(x) = \sum_{i=0}^{n} a_ix^i.\] 举个例子，2次多项式，其中\(a_0=1,a_1=2,a_2=1\) 那么\(F(x)=1x^2+2x+1.\)这样即可通俗的表示出一个\(n\)次多项式。 二. 点值表示法众所周知，两个点确定一个一次函数，三个点确定一个二次函数。 所以，\(n+1\)个点确定一个一元\(n\)次多项式。 所以我们可以通过\(n+1\)个点来表示它。 那么相乘之后的点值如何计算？ 比如说两个2次多项式\(F(x)=x^2+2x+1\)（红色）与\(G(x)=3x^2-4x-2\)（蓝色），它们的图像如图所示： 那么观察图中\(x=1\)时的情况。此时\(F(1)=4,G(1)=-3.\) 所以，显然，\(F(1)\times G(1)=-12\).也就是说在\(Z=F*G\)这一多项式内，带入\(1\)，得到的结果是\(-12\). 等等，好像有哪里不对。如果说\(Z=F*G\)的话，那么Z的次数应该是\(2n\). 那\(Z\)需要\(2n+1\)个点来确定。但是原来只需要\(n+1\)个点，咋办？ 很简单，在原来的多项式里每个都多加\(n\)个点即可。反正多项式已知。 这样就可以用点值来进行操作。也就是说先转成点值，再一乘，再转回来，就是计算流程。 但是好像还是很慢。那么如何优化呢?复数部分 复数，即形如\(a+bi\)的数，其中\(i^2=-1.\) \(a\)称为实部，\(bi\)称为虚部。 或者说：在一个数轴上（只有x轴），我们可以表示出任何实数。 那么，多加一维（y轴），也就是类似于平面直角坐标系一样，我们就可以表示出任意一个复数。 所以我们把这个坐标系叫做复平面，其中x轴称为实轴，y轴称为虚轴。复数运算 复数相加：实部相加，虚部相加，例如\[(a+bi)+(c+di)=(a+c)+(b+d)i.\] 复数相减：同理。\[(a+bi)-(c+di)=(a-c)+(b-d)i.\] 复数相乘：像一次多项式一样相乘。 注意\(i^2=-1\).\[(a+bi)(c+di)=ac+(ad+bc)i-bd=(ac-bd)+(ad+bc)i.\] 复数相除： 相信大家都学过共轭根式。同样的，复数也有共轭。 即：\(a+bi\)的共轭为\(a-bi\)。 这两个复数卡乘在一起一定是个实数。即\[(a+bi)(a-bi)=a^2-(bi)^2=a^2+b^2.\] 所以再除的时候，将分子分母同乘分母的共轭，就可以将分母有理化。 即\[\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{c^2+d^2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i.\] 复数逆元：\[\frac{1}{a+bi}=\frac{a}{a^2+b^2}-\frac{b}{a^2+b^2}i.\] 同样的，复数运算在复平面内也有一定规律可循。 考虑在复平面内的两个复数：（借张图）

表示的是复数\((1+4i)\)与复数\((3+2i)\)相乘所得的结果：\((-5+14i)\)。 设其中\((5,0)\)点为位置\(P\)，则\(\angle{POC}=\angle{BOA}.\) 还有：\(\overline{OB}\times\overline{OC}=\overline{OA}.\) 第二个证明：勾股定理。先把\(i\)消掉。\(\overline{OB}^2=a^2+b^2.\) \(\overline{OC}^2=c^2+d^2.\)\(\overline{OB}^2\times \overline{OC}^2=(a^2+b^2)(c^2+d^2)=a^2c^2+a^2d^2+b^2c^2+b^2d^2.\)\(\overline{OA}^2=(ac-bd)^2+(ad+bc)^2=a^2c^2+a^2d^2+b^2c^2+b^2d^2.\)得证。 我们将复数中，复数向量的长度称为模长，向量与x轴正方向的夹角称为幅角。 根据上面的东西：复数相乘时，模长相乘，幅角相加。单位根 一个n次的单位根即为方程\(x^n=1\)的复数解。 考虑这样一个图： 其中圆上的所有复数模长都是1，这个圆称为单位圆。 考虑\(|x|\)的取值范围： 如果\(|x|<1\)，那么\(|x^n|<1\)（模长*模长，越乘越短） 如果\(|x|>1\)，那么\(|x^n|>1\)（模长*模长，越乘越长） 所以一定有\(|x|=1\)，即点一定在这个单位圆上。 那么还有一条，就是一个n次的单位根共有n个，并且从(1,0)开始，每次逆时针走\(\frac{1}{n}\)个圆周。分别称为\(\omega_{n}^{0},\omega_{n}^{1},...\omega_{n}^{n-1}.\)（逆时针编号！！！） 其实还有\(k\geq n\)和\(k<0\)的单位根，本质上就模个n就可以了，这里不再多加考虑。 有关单位根也有很多性质：

\(\omega_{n}^{0}=1\)\(\omega_{n}^{k}=(\omega_{n}^{1})^k\)，很好理解，就是每次转\(\frac{1}{n}\)圆周。\(\omega_{n}^{j} \times \omega_{n}^{k}=\omega_{n}^{j+k}.\) 证明：很简单，代入性质2即可。 根据3，有：\[\omega_{n}^{j+k}=\omega_{n}^{j} \times \omega_{n}^{k}.\]\[(\omega_{n}^{k})^{-1}=\frac{1}{\omega_{n}^{k}}=\frac{\omega_{n}^{0}}{\omega_{n}^{k}}=\omega_{n}^{-k}=\omega_{n}^{n-k}.\] （根据上面模个n的说法）\[\omega_{pn}^{pk}=\omega_{n}^{k}.\]（将蛋糕分成pn份取pk份=将蛋糕分成n份取k份）\[\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}.\]（相当于绕了半圈\(\rightarrow\)取了个相反数）\[(\omega_{n}^{k})^j=(\omega_{n}^{j})^k.\]（大小为k，共有j个=大小为j，共有k个） ...差不多就这些了 ---\(FFT\)流程：\(DFT,IDFT\)\(DFT\)加速版本：分治 首先，你需要控制得当将整个n-1次多项式（保证n是2的正整数次幂）按系数奇偶拍成两部分：\[F(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+...+a_{n-1}x^{n-1}).\] 再设两个n/2项多项式：\[FL(x)=a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1}\]\[FR(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}\] 那么有\(F(x)=FL(x^2)+x*FR(x^2).\)（不懂可以自己设一下代入进去） 下面，我们来代入单位根，看看会有什么神奇的事情发生。 比如说我们代一个\(\omega_{n}^{k}(k<\frac{n}{2}):\)\[F(\omega_{n}^{k})=FL((\omega_{n}^{k})^2)+\omega_{n}^{k}\times FR((\omega_{n}^{k})^2)\] 其中，\((\omega_{n}^{k})^2=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k}.\)然后代入：\[F(\omega_{n}^{k})=FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}\times FR(\omega_{\frac{n}{2}}^{k}).\] 然后呢？有啥用？ 不难发现，如果我们知道了\(FL(x),FR(x)\)在x取\(\omega_{\frac{n}{2}}^{0},\omega_{\frac{n}{2}}^{1}...\omega_{\frac{n}{2}}^{\frac{n}{2}-1}\)下的所有值，我们就可以知道\(F(x)\)在x取\(\omega_{n}^{0},\omega_{n}^{1}...\omega_{n}^{\frac{n}{2}-1}\)下的所有值（代回原式）。 但是这还不是全部啊，剩下那个部分怎么办？ 没关系，考虑\(k<\frac{n}{2}\)，代入\(\omega_{n}^{k+\frac{n}{2}}\)时的情况。

\[F(\omega_{n}^{k+\frac{n}{2}})=FL((\omega_{n}^{k+\frac{n}{2}})^2)+\omega_{n}^{k+\frac{n}{2}}\times FR((\omega_{n}^{k+\frac{n}{2}})^2).\]

其中，\((\omega_{n}^{k+\frac{n}{2}})^2=\omega_{n}^{2k+n}.\)

故\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{n}^{2k+n}).\]

\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{n}^{2k})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{n}^{2k}).\] （模性质）

\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{\frac{n}{2}}^{k}).\] （性质3 延展）

\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}\times FR(\omega_{\frac{n}{2}}^{k}).\]

（取反性质）

然后观察这个式子和前面式子的区别，发现就是把+号改成了-号。也就是说我们可以通过那组数据来算出当\(k>\frac{n}{2}\)时的情况。这就是加速DFT的原理。IDFT 至此，我们已经可以将其转化为点值，那么如何还原回来呢？使用IDFT插值。 考虑一个向量\(c_0,c_1...c_{n-1}\)，满足\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i\times y_i.\] 其中\(y_i\)表示乘出来的点的y坐标。 那么，这个式子也就代表一个以\(y\)为系数的方程在取\(x=\omega_{n}^{-k}\)时的复数解。 我们尝试着来化简一波。

\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i\times y_i\]\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i(\sum_{j=0}^{n-1}a_j(\omega_{n}^{i})^j)\]\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j})^i)\]\[c_k=\sum _{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j})^i(\omega_{n}^{-k})^i)\]\[c_k=\sum _{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j-k})^i)\]\[c_k=\sum _{j=0}^{n-1}a_j\color{red}{(\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i)}\] 发现后面标红的东西就是一个等比数列，那么我们可以直接计算。即

\[\sum_{i=0}^{n-1}X^i=\frac{X^n-1}{X-1}.\]

我们设这个东西为公式\(Z.\) 于是就有\[c_k=\sum _{j=0}^{n-1}a_jZ(\omega_{n}^{j-k}).\]

观察函数\(Z\)的相关性质。 赫然发现：我们带进去的是单位根！单位根啥性质来着？\(X^n=1!\) 也就是说当\(X≠1\)，也就是\(k≠0\)时，原式\(=0.\) （分子为0） 那么\(X=1\)时，也就是\(k=0\)时呢？这时候这个求和公式就不管用了，带回原式算，发现由于\(1^i=1\)（i是自然数），那么\(Z(x)=n.\) 所以说，什么情况下带进去的值\(k=0\)呢？显然是\(j-k\)的时候！此时单项\(=na_k\). 否则的话会乘一个0，答案固然也是0. 所以有：\[c_k=na_k\]（飞跃性的一步！！！） 回顾整个式子，发现我们算\(c\)的时候连\(a\)的影子都没有，现在我们竟然推出了一个和\(a\)有关的关系式！太神奇了！这就是\(IDFT\)的原理。 至此，\(\text{FFT}\)的重要部分已经全部介绍完毕，下面进入代码实现与细节处理部分。

代码实现 朴素DFT:

complex < double > omega(const int n , const int k) { // 单位根 complex < double > ret = {cos(2 * Pi / n * k) , sin(2 * Pi / n * k)}; if(!inv) return ret; else return conj(ret); } void DFT(complex < double > a , const int n) { if(n == 1) return ; static complex < double > buf[MAXN]; const int m = n / 2; for(int i = 0; i < m; i++) { buf[i] = a[(i << 1)]; buf[i + m] = a[(i << 1 | 1)]; } copy(buf , buf + n , a); complex < double > *a1 = a , *a2 = a + m; DFT(a1 , m); DFT(a2 , m); for(int i = 0; i < m; i++) { complex < double > w = omega(n , i); buf[i] = a1[i] + w * a2[i]; buf[i + m] = a1[i] - w * a2[i]; } copy(buf , buf + n , a); }

发现这样的效率显然不行，那么考虑迭代形式。 首先观察分治到最后的二进制位表示：

000 001 010 011 100 101 110 111 0 1 2 3 4 5 6 7 0 2 4 6 - 1 3 5 7 0 4 - 2 6 - 1 5 - 3 7 0 - 4 - 2 - 6 - 1 - 5 - 3 - 7 000 100 010 110 001 101 011 111

然后...我们惊奇的发现最后的二进制位就是之前的二进制位倒过来！TQL! 接着考虑空间优化。 不难发现，由于DFT的两个式子十分相近，所以我们可以只算一次，同时通过技巧性的操作直接覆盖原数组，而不再新开一个。这个操作被称为蝴蝶操作。 蝴蝶操作的代码：

for(int l = 2; l <= n; l <<= 1) { int m = l / 2; for(complex < double > *p = a; p != a + n; p += l) { for(int i = 0; i < m; i++) { complex < double > t = omega[n / l * i] * p[m + i]; p[m + i] = p[i] - t; p[i] += t; } } }

然后穿起来，与以前的FFT模板结合即可。

#include <bits/stdc++.h> #define dbg(x) cerr << #x " = " << x << "\n" #define INF 0x3f3f3f3f /* do not give up ! */ /* try your best! */ /* Read the meaning clearly! */ typedef long long LL; typedef long double ld; typedef unsigned long long ULL; using namespace std; template < typename T > inline void inp(T& t) { char c = getchar(); T x = 1; t = 0; while(!isdigit(c)) {if(c == '-') x = -1; c = getchar();} while(isdigit(c)) t = t * 10 + c - '0' , c = getchar();t *= x; } template < typename T , typename... Args > inline void inp(T& t , Args&... args) {inp(t); inp(args...);} template < typename T > inline void outp(T t) { if(t < 0) putchar('-') , t = -t; T y = 10 , len = 1; while(y <= t) y *= 10 , len++; while(len--) y /= 10 , putchar(t / y + 48) , t %= y; } template < typename T > inline T mul(T x , T y , T MOD) {x=x%MOD,y=y%MOD;return ((x*y-(T)(((ld)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;} const int MAXN = 2097152 + 10; const double Pi = acos(-1.0); int ans[MAXN]; struct FastFourierTransform { complex < double > omega[MAXN] , omegaInverse[MAXN]; void init(const int n) { for(int i = 0; i < n; i++) { omega[i] = complex < double > (cos(2 * Pi / n * i) , sin(2 * Pi / n * i)); omegaInverse[i] = conj(omega[i]); } } void Transform(complex < double > *a , const int n , const complex < double > *omega) { int k = 0; while((1 << k) < n) ++k; for(int i = 0; i < n; i++) { int t = 0; for(int j = 0; j < k; j++) if(i & (1 << j)) t |= (1 << (k - j - 1)); if(i < t) swap(a[i] , a[t]); } for(int l = 2; l <= n; l <<= 1) { int m = l / 2; for(complex < double > *p = a; p != a + n; p += l) { for(int i = 0; i < m; i++) { complex < double > t = omega[n / l * i] * p[m + i]; p[m + i] = p[i] - t; p[i] += t; } } } } void DFT(complex < double > *a , const int n) { Transform(a , n , omega); } void IDFT(complex < double > *a , const int n) { Transform(a , n , omegaInverse); for(int i = 0; i < n; i++) a[i] /= n; } }FFT; int main() { int nn , mm , x; inp(nn , mm); nn++ , mm++; int n = 1; while(n < nn + mm) n <<= 1; complex < double > a[MAXN] , b[MAXN]; for(int i = 0; i < nn; i++) { inp(x); a[i].real(x); } for(int i = 0; i < mm; i++) { inp(x); b[i].real(x); } FFT.init(n); FFT.DFT(a , n); FFT.DFT(b , n); for(int i = 0; i < n; i++) a[i] *= b[i]; FFT.IDFT(a , n); for(int i = 0; i < nn + mm - 1; i++) { ans[i] = static_cast < int > (floor(a[i].real() + 0.5)); cout << ans[i] << " "; } return 0; }

至此，该代码已经可以效率很高地通过P3803一题！完结撒花！ FFT例题将会另外放在一个帖子~

转载于:https://www.cnblogs.com/LiM-817/p/10887264.html

相关资源：DFT以及FFT详解