FFT

最近刚学了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，快速傅里叶变换，是一种快速求卷积的方式

意义

两个数的竖式计算乘法，本质上就是在求他们的卷积，然后根据十进制进位，而FFT可以在\(O(nlogn)\)的时间内求卷积

不用卷积解释的话，那你可以理解为将两个函数转成点值表示法，\(O(n)\)相乘，然后转回乘积

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)\)
对于

和

乘积的点值表示法是

复数

复数肯定都会，就不讲了

蝴蝶变换

在复数平面上的单位圆上的点都是可以取的，因为易见它们的若干次方可以为1，读者自证不难，这里不证

当满足这个条件时，对于我们的计算，时空复杂度优化，码量以及IFFT大有裨益

假设\(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\)

所以易见

还有几条

DFT是代入这些值的变换，是\(O(n^2)\)的

但DFT可以分治来搞，这就是FFT

DFT很好写，就不贴代码了，把这些值直接带入多项式计算就行了

FFT

设

按\(A(x)\)下标的奇偶性把\(A(x)\)分成两半，右边再提一个\(x\)

两边非常相似，于是再设两个多项式\(A_1(x)\)和\(A_2(x)\)，令

比较明显得出

设\(k \leq \frac{n}{2}\)，把\(\omega_n^k\)带入，接下来几步变换要多想想单位根的性质

对于\(A(\omega_n^{k+\frac{n}{2}})\)

\(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})\)

当\(j=k\)时，易知

否则，通过等比数列求和可知

所以

证毕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 Problem

两个整数\(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})\)，但Karatsuba是可以压位的，实现的好甚至可以吊打不压位的fft

关于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;
}