最近刚学了FFT，给我这个蒟蒻的脑细胞干废了，写篇笔记浅浅记录一下
先讲一下卷积，卷积就是求\(\int_{-\infty}^{\infty} f(\tau)g(x-\tau)d\tau\)，这玩意的用处有一个是多项式乘法\(\sum_{i+j=x}^{n+m} a_ib_j(i\leq n,j\leq m)\)，别的就比如说高斯模糊什么的
FFT，快速傅里叶变换，是一种快速求卷积的方式
点值表示法
正常我们表达都是\(f(x)=\sum_{i=0}^{n}a_ix^i\)，就和两个点只能划一条线一样，一个\(n\)次多项式由\(n+1\)个点唯一确定，都说到这个了，那我挖个拉格朗日插值法的坑回头讲，就是\(f(x)=(x_0,f(x_0)),(x_1,f(x_1))...(x_n,f(x_n))\)
点值表示法意义下的乘法是\(O(n)\)
对于
和
\[g(x)=(x_0,g(x_0)),(x_1,g(x_1))...(x_n,g(x_n)) \]乘积的点值表示法是
\[f\cdot g=(x_0,f(x_0)\cdot g(x_0)),(x_1,f(x_1)\cdot g(x_1))...(x_n,f(x_n)\cdot g(x_n)) \]复数肯定都会，就不讲了
蝴蝶变换
在复数平面上的单位圆上的点都是可以取的，因为易见它们的若干次方可以为1，读者自证不难，这里不证
假设\(n=8\)，我们取上面等分的8个点，逆时针从1开始，标到7
记编号为\(k\)的点代表的复数值为\(\omega_n^k\)，易得\((\omega_n^1)^k=\omega_n^k\)
众所周知，\(e^{i\theta}=cos\theta+isin\theta\)并且\(e^{i\pi}+1=0\)
所以易见
\[\omega_n^k=cos\frac{k}{n}2\pi+isin\frac{k}{n}2\pi \]还有几条
\[\omega_n^k=\omega_{2n}^{2k} \]\[\omega_n^{k+\frac{n}{2}}=-\omega_n^k \]\[\omega_n^0=\omega_n^n=1 \]DFT是代入这些值的变换，是\(O(n^2)\)的
但DFT可以分治来搞，这就是FFT
FFT
设
\[A(x)=\sum_{i=0}^{n-1}a_ix^i=a_0+a_1x+...+a_{n-1}x^{n-1} \]按\(A(x)\)下标的奇偶性把\(A(x)\)分成两半，右边再提一个\(x\)
\[A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1}) \]\[=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2}) \]两边非常相似，于是再设两个多项式\(A_1(x)\)和\(A_2(x)\)，令
\[A_1(x)=a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1} \]\[A_2(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1} \]比较明显得出
\[A(x)=A_1(x^2)+xA_2(x^2) \]设\(k \leq \frac{n}{2}\)，把\(\omega_n^k\)带入，接下来几步变换要多想想单位根的性质
对于\(A(\omega_n^{k+\frac{n}{2}})\)
\[A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n}) \]\[=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k}) \]\[=A_1(\omega_\frac{n}{2}^k)-\omega_n^kA_2(\omega_\frac{n}{2}^k) \]\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)后面只有符号不同
就是只要已知\(A_1(\omega_\frac{n}{2}^k)\)和\(A_2(\omega_\frac{n}{2}^k)\)，就可以求\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)了
现在我们就可以递归分治来搞FFT了
每一次回溯时只扫当前前面一半的序列，即可得出后面一半序列的答案
\(n==1\)时就一个常数项，直接\(return\)
\(O(nlognloglogn)\)
优化-迭代版FFT
- 每个位置分治后的最终位置为其二进制翻转后得到的位置
读者自证不难，这里不证
这样的话我们可以先把原序列变换好，把每个数放在最终的位置上，再一步一步向上合并
一句话就是可以\(O(n)\)预处理第\(i\)位最终的位置\(rev[i]\)
\(O(nlogn)\)
IFFT
- 一个多项式在分治的过程中乘上单位根的共轭复数，分治完的每一项除以\(n\)即为原多项式的每一项系数
证明
设\((y_0,y_1,y_2,...,y_{n-1})\)为\(A(x)=a_0+a_1x+...+a_{n-1}x^{n-1}\)做完FFT的样子，设\(B(x)=y_0+y_1x+...+y_{n-1}x^{n-1}\)，把它们的倒数\(\omega_n^0,\omega_n^{-1},...,\omega_n^{1-n}\)带入做FFT，得到\((z_0,z_1,...,z_{n-1})\)
\[z_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i \]\[=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i \]\[=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^i)^{j-k}) \]当\(j=k\)时，易知
\[\sum_{i=0}^{n-1}(\omega_n^i)^{j-k}=n \]否则，通过等比数列求和可知
\[\sum_{i=0}^{n-1}(\omega_n^i)^{j-k}=\frac{(\omega_n^{j-k})^n-1}{\omega_n^{j-k}-1}=\frac{(\omega_n^n)^{j-k}-1}{\omega_n^{j-k}-1}=\frac{1-1}{\omega_n^{j-k}-1}=0 \]所以
\[z_k=na_k \]\[a_k=\frac{z_k}{n} \]证毕Q.E.D
代码
学会了FFT之后就可以自己手搓板子啦
注意：FFT的\(n\)只能是\(2^k\)
void fft(cp *a,int inv){//a是要改变的数组，inv是区分FFT和IFFT的
for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));//二进制翻转
for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);//防止交换两次，就是没交换
for(int mid=1;mid<n;mid<<=1){//二分长度
cp temp(cos(pi/mid),inv*sin(pi/mid));//这个就是单位根，遍历每一个omega的
for(int i=0;i<n;i+=(mid<<1)){
cp omega(1,0);//初始元
for(int j=0;j<mid;j++,omega*=temp){
cp x=a[i+j],y=omega*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;//这个就是蝴蝶变换什么的
}
}
}
if(inv==-1) for(int i=0;i<n;i++) a[i]/=n;//IFFT
}
来实践一下
两个整数\(a\)和\(b\)可以看做两个多项式\(x=10\)带入的结果，可以把每一位看成系数，\(a=a_na_{n-1}...a_1a_0\)可以看成多项式\(a=a_nx^n+a_{n-1}x^{n-1}...+a_1x+a_0\)当\(x=10\)时的数，看到乘法就是这两个多项式的乘积带回\(x=10\)，所以，开卷！！！
正常的运算卷积速度，或者说乘法，是\(O(n^2)\)的，即使DFT和IDFT也是\(O(n^2)\)的
Karatsuba乘法的时间复杂度是\(O(n^{log3})≈O(n^{1.585})\)
关于NTT，实际上就是利用一些特殊的质数（例如\(998244353\)，这些质数的特征都是可以可以表示 \(q\times 2^k + 1\)的形式）进行膜意义下的运算代替浮点复数运算来保证精度，原理是质数原根和复数单位根在DFT运算中具有相同的性质
这时候就可以用FFT啦！主要是因为我不会NTT，回头会了再发一篇，诶不对，我好像又挖了一个坑
在2019年之前，FFT的时间复杂度是\(O(nlognloglogn)\)，虽然\(loglogn\)这项系数极小，但是直到2019年才去除掉，变成\(O(n logn)\)
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define cp complex<double>
const double pi=acos(-1);
const int N=2097152+114514;
int len1,len2,n=1,k=0;
string a1,b1;
cp a[N],b[N];
int ans[N],rev[N];
void fft(cp *a,int inv){
for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int mid=1;mid<n;mid<<=1){
cp temp(cos(pi/mid),inv*sin(pi/mid));
for(int i=0;i<n;i+=(mid<<1)){
cp omega(1,0);
for(int j=0;j<mid;j++,omega*=temp){
cp x=a[i+j],y=omega*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;
}
}
}
if(inv==-1) for(int i=0;i<n;i++) a[i]/=n;
}
main(){
cin>>a1>>b1;
len1=a1.length(),len2=b1.length();
for(int i=0;i<len1;i++) a[i]=(double)a1[len1-i-1]-48;
for(int i=0;i<len2;i++) b[i]=(double)b1[len2-i-1]-48;
while(n<len1+len2-1) n<<=1,k++;
for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
fft(a,1),fft(b,1);
for(int i=0;i<n;i++) a[i]*=b[i];
fft(a,-1);
for(int i=0;i<n;i++) ans[i]+=(int)(a[i].real()+0.5),ans[i+1]+=ans[i]/10,ans[i]%=10;
while(!ans[n]&&n) n--;
for(int i=n;i>=0;i--) cout<<ans[i];
return 0;
}