多项式
对多项式进行乘法，就是对两个多项式进行加法卷积。其中卷积结果 \(C(k)=\sum_{i=0}^kA(i)B(k-i)\)。
分治乘法
将 \(A(x)\) 左右拆半，不足则末尾（最高位）补上 \(0\)，令 \(n=2^k\)。则
同理，\(B(x)\) 左右拆半，则卷积
这样可以转化为 \(3\) 个规模为 \(n/2\) 的子问题，总时间为 \(T(n)=O(n)+3T(n/2),T(n)=O(n^{\log_2 3})=O(n^{1.585})\)。
点值与插值
点值：对于一个 \(x\) 的 \(F(x)\) 的值。
插值：已知 \(F(x)\) 的若干点值，求其系数序列 \(G(x)\)。
根据定义，\(F(x)=\sum_{i=0}^{n-1}x^iG(i)\)。
（简记）拉格朗日插值
FFT 中的插值使用单位根反演实现。
复平面单位根
\(z=\cos\theta+i\sin\theta,\omega_n^k=\cos(\frac{2\pi k}{n})+i\sin(\frac{2\pi k}{n})\)，其点集为单位圆均分成 \(n\) 份的点集，且点集包含 \(1\)。
-
性质 1：\(\omega_n^k=\omega_n^{k+n}\)
-
性质 2：\(\omega_{nb}^{kb}=\omega_n^k\)。
-
性质 3：若 \(n\) 为偶数，\(\omega_n^{k+n/2}=-\omega_n^k\)
FFT
DFT
FFT 多项式乘法外层过程如下：
- 选定 \(n(n=2^{k'})\) 个数 \(\omega_n,\omega_n^{2},...\omega_n^{n-1}\)。
- 对于 \(i\in [0,n-1]\) 求出点值 \(A(\omega_n^i),B(\omega_n^i)\)（DFT）
- 根据 \(C(\omega_n^i)=A(\omega_n^i)B(\omega_n^i)\) 求出 \(C\) 的点值。
- 插值求出 \(C(x)\) 的系数。（IDFT，即 DFT 逆过程）
类似分治乘法，我们采用拆半，但是是奇偶拆半。
假设我们已经知道了一堆点值，代入 \(\omega_n^k,k\in[0,\frac{n}{2})\)。
考虑到 \(FL,FR\) 的计算，我们每次实际分治时只需要计算所有 \(k<\frac{n}{2}\) 的部分，\(\frac{n}{2}\le k<n\) 的部分根据性质 1 \(\omega_{n/2}^{k+n/2}=\omega_{n/2}^k\)，用到的 \(FL,FR\) 是和前面一样的，但是根据性质 3，我们在转移式子中用到的 \(\omega_n^{k+n/2}\) 需要变成 \(-\omega_n^k\)，即：
这是 DFT 的部分。
IDFT
根据上面的 DFT：
根据单位根反演：
所以求出点值以后，我们再跑一边 DFT 即可。
这样我们就可以用子问题求出原问题的答案了，时间复杂度 \(T(n)=T(n/2)+O(n),T(n)=O(n\log n)\)。
一些常数优化
迭代版神秘常数问题
- 如果去掉分治过程直接迭代，先枚举起始点再枚举增量，否则会变得异常缓慢。
蝴蝶变换
-
为了避免大规模数组拷贝，我们使用蝴蝶变换。注意到分治的过程中我们先不断地按照奇偶把数组切分然后分到两边，为了避免这个切分我们考虑一开始就把点放到它最下面那层分治所处的位置，然后合并和之前一样左边作为 \(FL\)，右边作为 \(FR\)。
我们发现一个点 \(i\) 的二进制形式决定了其最终所处位置，把 \(i\) 视为 \(n=2^{k'}\) 的 \(k'\) 位二进制数，那么 \(i\in[0,n-1]\)，从上至下每次奇偶分类相当于看二进制最低位，\(0\) 就分到左边，\(1\) 就分到右边，然后右移一位。不妨把这个二进制数
reverse（翻转）一下，那么每次的选择及其权重和就恰好对应了翻转后的二进制数。于是预处理 \(i\) 的reverse\(tr[i]\)，然后若 \(i<tr[i]\) 交换即可得到最下层位置（避免两次交换）。
三次变两次优化
构造出：
即把 \(F(x)\) 和 \(iG(x)\) 系数相加得到 \(P(x)\) 系数后做一次 DFT，得到的点值平方后做 IDFT，最后得到的系数虚部 \(/2\) 就是 \(F(x)G(x)\)。double 运算容易掉精度，承受能力为跨度上线 \(10^{12}\) 左右，如果相乘系数差太大（如一个 \(\in[10^{-6},10^{-5}]\)，一个 \([10^5,10^6]\)）在有平方项会导致 \(10^{24}\) 的精度跨度存在。解决办法就是乘上一个定值，做完 FFT 再除回去即可。
Code
分治版
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const double Pi=acos(-1);
const int N=1.35e6+5;
int n,m,tr[N<<1];
struct CP{CP(double xx=0,double yy=0){x=xx,y=yy;}double x,y;CP operator +(const CP &b)const{return (CP){x+b.x,y+b.y};}CP operator -(const CP &b)const{return (CP){x-b.x,y-b.y};}CP operator *(const CP &b)const{return (CP){x*b.x-y*b.y,x*b.y+y*b.x};}
}f[N<<1],g[N<<1];
void FFT(CP *f,int len,bool tf){if(len==1)return ;FFT(f,len/2,tf);FFT(f+len/2,len/2,tf);CP per(cos(2*Pi/len),sin(2*Pi/len)),now(1,0);if(!tf)per.y=-per.y;for(int k=0;k<len/2;k++){CP tt=now*f[k+len/2];f[k+len/2]=f[k]-tt;f[k]=f[k]+tt;now=now*per;}
}
int main(){ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n>>m;for(int i=0;i<=n;i++)cin>>f[i].x;for(int i=0;i<=m;i++)cin>>g[i].x;m=m+n;for(n=1;n<=m;n<<=1);for(int i=0;i<n;i++)tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);for(int i=0;i<n;i++)if(i<tr[i])swap(f[i],f[tr[i]]),swap(g[i],g[tr[i]]);FFT(f,n,1);FFT(g,n,1);for(int i=0;i<n;i++)f[i]=f[i]*g[i];for(int i=0;i<n;i++)if(i<tr[i])swap(f[i],f[tr[i]]);FFT(f,n,0);for(int i=0;i<=m;i++)cout<<(int)(f[i].x/n+0.49)<<' ';return 0;
}
迭代版
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const double Pi=acos(-1);
const int N=1.35e6+5;
int n,m,tr[N<<1];
struct CP{CP(double xx=0,double yy=0){x=xx,y=yy;}double x,y;CP operator +(const CP &b)const{return (CP){x+b.x,y+b.y};}CP operator -(const CP &b)const{return (CP){x-b.x,y-b.y};}CP operator *(const CP &b)const{return (CP){x*b.x-y*b.y,x*b.y+y*b.x};}
}f[N<<1],g[N<<1];
void FFT(CP *f,bool tf){if(n==1)return ;for(int i=0;i<n;i++)if(i<tr[i])swap(f[i],f[tr[i]]);for(int p=2;p<=n;p<<=1){int len=p>>1;CP per(cos(2*Pi/p),sin(2*Pi/p));if(!tf)per.y=-per.y;for(int l=0;l<n;l+=p){CP now(1,0);for(int k=l;k<l+len;k++){CP tt=now*f[k+len];f[k+len]=f[k]-tt;f[k]=f[k]+tt;now=now*per;}}}
}
int main(){ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n>>m;for(int i=0;i<=n;i++)cin>>f[i].x;for(int i=0;i<=m;i++)cin>>g[i].x;m=m+n;for(n=1;n<=m;n<<=1);for(int i=0;i<n;i++)tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);FFT(f,1);FFT(g,1);for(int i=0;i<n;i++)f[i]=f[i]*g[i];FFT(f,0);for(int i=0;i<=m;i++)cout<<(int)(f[i].x/n+0.49)<<' ';return 0;
}