FFT理论与习题
参考了 这篇博客，并且自己重新证了一下这篇博客中，笔者认为有误的 IDFT 中 \(j \neq k\) 的部分。
第零部分·原理
FFT 是一种 DFT 的高效算法，称为快速傅里叶变换（fast Fourier transform），当然也可以用来算 IDFT。
这俩玩意前者可以把多项式从系数表达式转化成点值表达式，后者可以把点值表达式转化成系数表达式。
点值表达式有很多非常好性质，比如两个多项式相乘，可以直接用点值表达式相乘。
这俩玩意的复杂度都是 \(O(n \log n)\)，非常厉害。
第一部分·\(n\) 次单位根
这里的 \(n\) 次单位根指的是满足 \(x^n=1\) 中 \(x\) 的 \(n\) 个解，我们按照他们在复平面上逆时针的顺序，依次记为 \(\omega_n^0,\omega_n^1……\omega_n^{n-1}\)。
这个东西有很好多条非常好的性质，在这里进行介绍：
-
\(\omega_n^0=1\)。
-
$\omega_n^1 = \cos{\frac{2\pi}{n} + i \sin{\frac{2\pi}{n}}} $，几何意义是把单位圆平局分成 \(n\) 份，从 \(x\) 轴正方向逆时针转遇到的第一个分割点就是 \(\omega_n^1\)，用三角函数算一下就理解了。
-
\(\omega_n^i\times \omega_n^j=\omega_n^{i+j}\)，即满足一个类似幂的相乘的东西。
-
\((\omega_n^i)^j=\omega_n^{ij}\)，即满足一个类似幂的乘方的东西。
-
\(\omega_n^{k}=\omega_n^{k\bmod n}\)，即单位圆上多转几圈。
-
\(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)，这条性质使用时需满足 \(n\bmod 2=0\)，即 \(n\) 为偶数。它的集合含义就相当于单位圆上转半圈。
-
\(\omega_{2n}^{2k}=\omega_{n}^{k}\)，即单位圆的分割分得稀疏一半，结果不会变。
差不多一共会用到这些性质。
第二部分·DFT
这里我们讲如何把一个 \(n-1\) 次多项式从系数表达式转化为点值表达式。
我们默认这里 \(n\) 是 \(2\) 的整次幂，如果原多项式不满足，可以补几个系数为 \(0\) 的项。
我们记
我们就有 \(F(x)=Fl(x^2) + x\cdot Fr(x^2)\)。
接下来，就可以考虑将 \(\omega_n^0,\omega_n^1……\omega_n^{n-1}\) 带入 \(F(x)\) 中了！我们记 \(k\in[0,\frac{n}{2})\)。
1.带入 \(\omega_n^k\)。
2.带入 \(\omega_n^{k+\frac{n}{2}}\)。
观察这两个等式，我们惊讶地发现：只要知道了 \(Fl(x)\) 和 \(Fr(x)\) 在 \(x=\omega_{\frac{n}{2}}^k)\) 处的取值，就可以 \(O(n)\) 算出 \(F(\omega_n^0,\omega_n^1……\omega_n^{n-1})\) ！
而由于 \(k \in [0,\frac{n}{2})\)，且 \(Fl(x)\) 和 \(Fr(x)\) 都是项数为 \(\frac{n}{2}\) 的多项式，我们就会发现计算 \(Fl(x)\) 和 \(Fr(x)\) 在 \(\omega_{\frac{n}{2}}^0,\omega_{\frac{n}{2}}^1……\omega_{\frac{n}{2}}^{\frac{n}{2}-1}\) 处的值，本质上是一个递归、分治的过程！
而分治的层数显然是 \(O(\log n)\) 级别的，每一层都要 \(O(n)\) 进行遍历，总复杂度就是 \(O(n\log n)\)！
这样子，我们就成功的在 \(O(n\log n)\) 的复杂度内完成了多项式系数表达式转点值表达式！
第三部分·IDFT
这里我们讲如何把点值表达式转化回系数表达式。
记 \(g_k=F(\omega_n^k)=\sum\limits_{i=0}^{n-1} f_i(\omega_n^k)^i\)。
则有 \(nf_k=\sum\limits_{i=0}^{n-1} g_i(\omega_n^{-k})^i\)。
下面给出证明：
这里 \(t_1= \sum\limits_{i=0}^{n-1}f_k\omega_n^{i(k-k)}\)，$t2=\sum\limits_{j=0}^{n-1} [j \ne k]\sum\limits_{i=0}^{n-1} f_j\omega_n^{i(j-k)} $，即将最外层关于 \(j\) 的求和拆成 \(j=k\) 和 \(j\neq k\) 的两部分，前者为 \(t_1\)，后者为 \(t_2\)。
所以 \(RHS=t1+t2=nf_k+0=nf_k=LHS\)。证毕。
在这里重新摆一下结论：\(nf_k=\sum\limits_{i=0}^{n-1} g_i(\omega_n^{-k})^i\)。
我们记 \(G(x)=\sum\limits_{i=0}^{n-1} g_ix^i\)，则 \(\sum\limits_{i=0}^{n-1} g_i(\omega_n^{-k})^i\) 可以看作将 \(\omega_n^{-k}\) 带入 \(G(x)\) 所得的值。
于是我们如果能求的 \(G(x)\) 在 \(\omega_n^0,\omega_n^{-1}……\omega_n^{-(n-1)}\) 处的点值，就可以 \(O(n)\) 求出 \(f_0,f_1……f_{n-1}\)了。
而这个点值表达式的求职过程和 DFT 的分治一模一样！就是把带进去的单位根的“指数”变了个相反数！
于是我们就可以把 DFT 和 IDFT 封装进一个函数里头，只需要加一个 \(flag\)，来判定要不要带一个负号就行了！
第四部分·优化
写一个分治的做法常数很大，我们这里有更好的写法：
第一次优化：减少数组拷贝。
我们发现在分治的过程中，每一层都会把奇数处的系数全放到左边，偶数处的系数全放到右边，而这样的操作放在每一层都要进行的递归了就会增大很大的常数。
我们能否一次性知道每一个数在最底层都应该放在哪里呢？
观察下列分治过程：
0 1 2 3 4 5 6 7 第1层
0 2 4 6|1 3 5 7 第2层
0 4|2 6|1 5|3 7 第3层
0|4|2|6|1|5|3|7 第4层
会发现位于 \(i\) 位置上的数，其实是 \(f_{rev_i}\)，其中 \(rev_i\) 指将 \(i\) 的二进制翻转后得到的数，比如 \(3=011\)，反转后 \(rev_3=110=6\)。
关于如何 \(O(n)\) 预处理出 \(rev\) 数组，可以使用一个类似 dp 的思想：
rev[i]=(rev[i>>1]>>1) | (i&1 ? (lim>>1) : 0)
假设此时二进制串长度为 \(m\)，则其中 \(lim=2^m\)。
借助一个例子 i=110100 ，这句话的含义如下：
i>>1=011011，则 rev[i>>1]=rev[001101]=101100，然后 rev[i>>1]>>1=(101100)>>1=010110，具体含义就是把 \(i\) 的最末位去掉，剩余部分翻转后的东西。
| (i&1 ? (lim>>1) : 0) 指的是如果原本最末尾是 \(1\)，就需要把 lim>>1 拼到刚才的出来的 010110最前面，得到 110110，就是所求的答案。而最末尾是 \(0\) 时，什么也不用干。
这样就得到了 \(rev\) 数组了。
第二次优化：递归改为枚举。
我们可以不递归，而直接通过枚举长度进行分治，这样的好处在于可以少算很多次三角函数，减少很多常数。
这样子，我们最常用的 fft 模板就成型了！
下面是非常优雅的模板：
//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize("Ofast","inline")
#include<bits/stdc++.h>
#define frr(a) freopen(a,"r",stdin)
#define fww(a) freopen(a,"w",stdout)
#define MP(a,b) make_pair(a,b)
#define DEBUG
using namespace std;
typedef long long ll;
const int N=2e6+10;
const double Pi=acos(-1);
ll n,m,tr[2*N];
struct cplx{
cplx (double xin=0,double yin=0){x=xin,y=yin;}
double x,y;
friend cplx operator+(cplx a,cplx b){return cplx(a.x+b.x,a.y+b.y);}
friend cplx operator-(cplx a,cplx b){return cplx(a.x-b.x,a.y-b.y);}
friend cplx operator*(cplx a,cplx b){return cplx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}f[2*N],g[2*N],res[2*N];
void fft(cplx *f,ll n,bool flag){//第三版·迭代实现
for(int i=0;i<n;i++) if(i<tr[i]) swap(f[i],f[tr[i]]);
for(int len=2;len<=n;len<<=1){
cplx tmp(cos(2*Pi/len),(flag==0?1.0:-1.0)*sin(2*Pi/len));
for(int k=0;k<n;k+=len){
cplx now(1,0);
for(int i=0;i<(len>>1);i++){
cplx tt=now*f[k+i+(len>>1)];
f[k+i+(len>>1)]=f[k+i] - tt;
f[k+i]=f[k+i] + tt;
now = now*tmp;
}
}
}return;
}
int main(){
scanf("%lld%lld",&n,&m);
for(int i=0;i<=n;i++) scanf("%lf",&f[i].x);
for(int i=0;i<=m;i++) scanf("%lf",&g[i].x);
ll tmp=1;while((1ll<<tmp) < n+m) tmp++; tmp=1ll<<tmp;
for(int i=0;i<tmp;i++) tr[i]=(tr[i>>1]>>1) | (i&1?(tmp>>1):0);
fft(f,tmp,1),fft(g,tmp,1);
for(int i=0;i<tmp;i++) res[i] = f[i] * g[i];
fft(res,tmp,0);
for(int i=0;i<=n+m;i++) printf("%lld ",(ll)(res[i].x/tmp + 0.5));
return 0;
}
第五部分·习题
一.P3338 [ZJOI2014] 力
此题咱都从 \(0\) 开始编号。
给定一个长为 \(n\) 数组 \(q\)，求一个数组 \(E\)，满足：
可以分成两部分考虑，先考虑前半部分。
我们令 \(g_i=\dfrac{1}{i^2}\)，特殊地，令 \(g_0=0\)，则：
这是一个卷积形式，直接 fft 求就行了。
而后半部分也是一个道理，只不过需要把 \(q\) 数组翻转一下，变成卷积形式就可以了。
注意写代码时为了保证精度，\(g\) 数组应写为 g[i]=(1.0/i)/i 。
