FFT学习及简单应用(一点点详细)

什么是FFT#

既然打开了这篇博客，大家肯定都已经对FFT(Fast Fourier Transformation)有一点点了解了吧
FFT即为快速傅里叶变换，可以快速求卷积（当然不止这一些应用，但是我不会）

系数表示法与点值表示法#

我们通常表示一个(n-1)次多项式是利用系数表示法like this:(f(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1})
点值表示法即为将多项式用坐标系上的若干个点表示
我们对这个多项式代入不同的值{(x_1,x_2,...,x_n)}
我们就可以得到(n)个点((x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n)))
实际上只要保证代入的(n)个数互不相同，那么这(n)个点就对应了唯一的(f(x))（类似(n)元一次方程组？）
点值表示法它的优势在哪呢？
我们如果想要求出两个多项式相乘，系数表示法就很麻烦，而点值表示法却只需要将相同的(x)对应的点值相乘就行了

DFT与IDFT#

DFT##

我们若将一个多项式强行转为点值表示法则时间复杂度为(O(n^2))
自然有人表示：太慢啦，就不能快点吗！
于是我们开始思考，如果我们对于((x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n)))取一系列特殊的有关系的值，是不是能加速呢？
当然可以，要不然就不会有这篇博客了
在DFT中，我们可以将复平面中的单位园n等分，然后每一等分的顶点(omega _n^i)作为(x_i)，(omega _n^i)就是单位根
（如图为(8)等分）

其中(omega _n^i=cos (frac{{2pi i}}{n}) + isin (frac{{2pi i}}{n}))
当然DFT只是让((x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n)))变得特殊了一些，并没有优化复杂度的效果（要不然哪来的FFT）

IDFT##

IDFT(Inverse Discrete Fourier Transform)顾名思义就是DFT的逆变换（当然我不会，这个东西，可以会，但没必要）

关于单位根#

首先，显然(omega _{2n}^{2i}=omega _n^i)（结合图像可得）((n=2^k))
(omega _n^{i+frac{n}{2}}=-omega _n^i)为什么呢，因为这两个单位根在坐标中终点关于原点对称
最后一个没有多大用的东西，但是后面会用到(omega _n^0=omega _n^n=1)，这个显然吧，因为这个顶点是重合的，且虚部为0

FFT#

递归解法##

FFT(Fast Fourier Transformation)快速傅里叶变换，有了上面这些东西，我们就可以考虑对原多项式进行变形
(f(x)=a_0+a_1x+a_2x^2+...+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}))
是不是发现这两个东西很像啊
我们设(f_1(x)=(a_0+a_2x+...+a_{n-2}x^{frac{n}{2}-1}),f_2(x)=(a_1+a_3x+a_{n-1}x^{frac{n}{2}-1}))
所以(f(x)=f_1(x^2)+xf_2(x^2))
我们发现，(f_1(x),f_2(x))和(f(x))的形式一样诶！
也就是说我们现在这个FFT已经可以用递归解决了！
然而，递归自带sb大常数，时间上并不能通过
于是，厉害的人们又发现了新的操作，蝴蝶变换

蝴蝶变换##

首先我们令(i<frac{n}{2})
(f(omega _n^i)=f_1({(omega _n^i)}^2)+omega _n^if_2({(omega _n^i)}^2))
(=f_1(omega _n^{2i})+omega _n^if_2(omega _n^{2i}))
(=f_1(omega _{frac{n}{2}}^i)+omega _n^if_2(omega _{frac{n}{2}}^i))
(f(omega _n^{i+frac{n}{2}})=f_1({(omega _n^i)}^2omega _n^n)+omega _n^{i+frac{n}{2}}f_2({(omega _n^i)}^2omega _n^n))
(=f_1(omega _n^{2i})-omega _n^if_2(omega _n^{2i}))
(=f_1(omega _{frac{n}{2}}^i)-omega _n^if_2(omega _{frac{n}{2}}^i))
是不是很神奇，是不是！
就差了一个符号，是不是很蝴蝶啊(亲爱的，你慢慢飞，小心...)
有了蝴蝶变换，我们发现我们只要知道(f_1(omega _{frac{n}{2}}^i),f_2(omega _{frac{n}{2}}^i))就可以求出(n)等分下两个不同单位根带入后的点值了，是不是很厉害！

rev数组##

但是我们又发现了一个新的问题，每次我们合并的时候，都是把要求的序列分为奇数位和偶数位做
可是我们原数组并不符合这个条件，直接做的话肯定死
我们可以研究一下每个数字和它的二进制表示
以8个数为例：

但是我们合并的顺序应该长成这样：
（读者可以自己模拟一下，段长每次为(2^k)，然后合并）
我们发现，其实如果把每个数的二进制位反过来一下，就变成顺序的了！

把二进制反过来其实也不复杂，就一行代码

rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1))

其实这段代码的意义大概是把除了最低位之外的位用i>>1翻转，然后再将最低位移到最高位
也就是说，做FFT的时候只要将系数重新排列，然后每段进行合并更新就好啦（这段看一下代码就会很清楚）

IFFT#

IFFT(Inverse Discrete Fourier Transform)快速傅立叶逆变换
实际上我们的傅立叶变换可以写成矩阵的形式：
(left[ {egin{array}{*{20}{c}} 1&1&1&{...}&1\ 1&{omega _n^1}&{omega _n^2}&{...}&{omega _n^{n - 1}}\ 1&{omega _n^2}&{omega _n^4}&{...}&{omega _n^{2n - 2}}\ {...}&{...}&{...}&{...}&{...}\ 1&{omega _n^{n - 1}}&{omega _n^{2n - 2}}&{...}&{omega _n^{(n - 1) imes (n - 1)}} end{array}} ight]left[ {egin{array}{*{20}{c}} {{a_0}}\ {{a_1}}\ {{a_2}}\ {...}\ {{a_{n - 1}}} end{array}} ight] = left[ {egin{array}{*{20}{c}} {f(omega _n^0)}\ {f(omega _n^1)}\ {f(omega _n^2)}\ {...}\ {f(omega _n^{n - 1})} end{array}} ight])
当然，看到这个形式，我们第一个想法就是对
(left[ {egin{array}{*{20}{c}} 1&1&1&{...}&1\ 1&{omega _n^1}&{omega _n^2}&{...}&{omega _n^{n - 1}}\ 1&{omega _n^2}&{omega _n^4}&{...}&{omega _n^{2n - 2}}\ {...}&{...}&{...}&{...}&{...}\ 1&{omega _n^{n - 1}}&{omega _n^{2n - 2}}&{...}&{omega _n^{(n - 1) imes (n - 1)}} end{array}} ight])
这个矩阵求逆，但是这样的话，复杂度就吃不消了，我们好不容易把(O(n^2))的DFT降为(O(nlogn))的FFT，结果逆不回去，这不是非常尴尬吗
但是我们发现，我们对这个矩阵求一下逆，可以得到这个形式是固定的
(left[ {egin{array}{*{20}{c}} {frac{1}{n}}&{frac{1}{n}}&{frac{1}{n}}&{...}&{frac{1}{n}}\ {frac{1}{n}}&{frac{1}{n}omega _n^{ - 1}}&{frac{1}{n}omega _n^{ - 2}}&{...}&{frac{1}{n}omega _n^{1 - n}}\ {frac{1}{n}}&{frac{1}{n}omega _n^{ - 2}}&{frac{1}{n}omega _n^{ - 4}}&{...}&{frac{1}{n}omega _n^{2 - 2n}}\ {...}&{...}&{...}&{...}&{...}\ {frac{1}{n}}&{frac{1}{n}omega _n^{1 - n}}&{frac{1}{n}omega _n^{2 - 2n}}&{...}&{frac{1}{n}omega _n^{ - (n - 1) imes (n - 1)}} end{array}} ight])
然后我们可以把(frac{1}{n})提出来，那么我们就可以先利用快速傅立叶变换，做出系数的(n)倍，然后再除掉
坑终于填完啦

FFT_Code#

这是一个高精乘的模板

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
char s[60050],s1[60050];
struct r{
	double real,image;
}a[131090],b[131090];
double pi=acos(-1);
r operator +(r a,r b){return (r){a.real+b.real,b.image+a.image};}
r operator -(r a,r b){return (r){a.real-b.real,a.image-b.image};}
r operator *(r a,r b){return (r){a.real*b.real-b.image*a.image,a.real*b.image+a.image*b.real};}
int ans[131090],rev[131090],bit,n,l;
void fft(r *a,int n,int op){
	for (int i=0;i<n;i++)if (i<rev[i])std::swap(a[rev[i]],a[i]);
	for (int i=1;i<n;i<<=1){
		r wn=(r){cos(pi/i),op*sin(pi/i)};//将单位圆分成i*2个部分
		for (int j=0;j<n;j+=i<<1){
			r wnk=(r){1,0};
			for (int k=j;k<i+j;k++,wnk=wnk*wn){
				r x=a[k],y=wnk*a[i+k];
				a[k]=x+y;
				a[k+i]=x-y;//这里很蝴蝶
			}
		}
	}
	if (op==-1)for (int i=0;i<n;i++)a[i].real/=n;
}
int main(){
	scanf("%d",&l);
	scanf("%s",s);
	scanf("%s",s1);
	int n=2;
	for (bit=1;(1<<bit)<(l<<1)-1;bit++)n<<=1;
	for (int i=0;i<l;i++)a[i]=(r){(double)(s[l-i-1]-'0'),0};
	for (int i=0;i<l;i++)b[i]=(r){(double)(s1[l-i-1]-'0'),0};
	for (int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));//处理rev数组
	fft(a,n,1);
	fft(b,n,1);
	for (int i=0;i<n;i++)a[i]=a[i]*b[i];
	fft(a,n,-1); //IFFT
	for (int i=0;i<n;i++){
		ans[i]+=(a[i].real+0.5);
		ans[i+1]=ans[i]/10;
		ans[i]%=10;
	}
	int i;
	for (i=n;i>=0&&!ans[i];i--);
	if (i<0)putchar('0');
	else for (;i>=0;i--)printf("%d",ans[i]);
}

什么是NTT#

NTT(Number Theoretic Transform)快速数论变换（不要问我快速是哪来的，我不知道）
可以将类似FFT的做法用整数来做

与FFT区别##

1、可以取模（必须取模）
2、比FFT快这是一个遥远的传说，究竟哪个快还是和数据有关

原根#

对于一个正整数m，存在一个整数g，满足(g^{varphi (m)}equiv 1(mod) (m))，则g为m的一个原根
一个数存在原根的充要条件是该数可以表示为(2,4,p^k,2*p^k)((p)为奇质数,(k geq 1))

NTT#

有了原根我们就可以考虑用原根来代替单位根，因为原根也有类似的性质（一般此类题目模数为质数，所以以下用(p-1)代替(varphi(p))）
那么我们考虑用(G_n^i)=(g^{i(p-1)/n})来代替(omega _n^i)
为什么这样代替呢，因为这样代替就满足：
(G_n^0=G_n^n=1)，这个显然吧？
(G_n^i=-G_n^{i+n/2})，这个通过原根的定义稍微推一下就能得到
(G _{2n}^{2i}=G _n^i)，这个应该问题也不大吧？
然后我们就发现NTT它锅了，哪里锅了呢？
性质当然是没有问题的，锅就锅在(n mid (p-1))的时候，我们没办法做浮点数的快速幂
那我们该怎么办，难道抛弃原根去寻找新方法？
Of course not.我们就规定NTT的模数只能取(k*2^m+1)并且(n leq 2^m)就行啦（最常用的可以NTT的模数就是998244353，如果出现新的大家自己检验一下啦）
（顺便送给各位，998244353的一个原根是3是不是非常良心）

NTT_Code#

这还是一个高精乘的模板

#include<cstdio>
#include<algorithm>
int g=3,mo=998244353,a[131090],b[131090],rev[131090],l,bit;
char s[60050],s1[60050];
int ksm(int p,int k){
	int ret=1;
	while (k){
		if ((k&1))ret=(ret*1ll*p)%mo;
		p=(p*1ll*p)%mo;k>>=1;
	}
	return ret;
}
void ntt(int *a,int n,int f){
	for (int i=0;i<n;i++)if (i<rev[i])std::swap(a[i],a[rev[i]]);
	for (int i=1;i<n;i<<=1){
		int gn=ksm(g,(mo-1)/(i<<1));if (f==-1)gn=ksm(gn,mo-2);
		for(int j=0;j<n;j+=i<<1){
			int gqn=1;
			for (int k=j;k<j+i;k++){
				int x=a[k],y=gqn*1ll*a[k+i]%mo;
				a[k]=(x+y)%mo;
				a[k+i]=(x-y+mo)%mo;
				gqn=(gqn*1ll*gn)%mo;
			}
		}
	}
}
int main(){
	scanf("%d",&l);
	scanf("%s",s);
	scanf("%s",s1);
	int n=2;
	for (bit=1;(1<<bit)<(l<<1)-1;bit++)n<<=1;
	for (int i=0;i<l;i++)a[i]=s[l-i-1]-'0';
	for (int i=0;i<l;i++)b[i]=s1[l-i-1]-'0';
	for (int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
	ntt(a,n,1);
	ntt(b,n,1);
	for (int i=0;i<n;i++)a[i]=(a[i]*1ll*b[i])%mo;
	ntt(a,n,-1);
	int ny=ksm(n,mo-2);
	for (int i=0;i<n;i++)a[i]=(a[i]*1ll*ny)%mo;
	for (int i=0;i<n;i++)a[i+1]+=a[i]/10,a[i]%=10;
	int i;
	for (i=n;i>=0&&!a[i];i--);
	if (i==-1)printf("0");
	else for (;i>=0;i--)printf("%d",a[i]);
}

一点小应用#

多项式求逆##

多项式求逆是个啥呢，就是我们有多项式(f(x))我们希望求出多项式(g(x))满足
(f(x)g(x)equiv 1(mod) (x^n))
假设我们现在已经得到了(p(x))满足(f(x)p(x)equiv 1(mod) (x^{left lceil frac{n}{2} ight ceil}))
因为(g(x))满足(f(x)g(x)equiv 1(mod) (x^n))，则必然满足(f(x)g(x)equiv 1(mod) (x^{left lceil frac{n}{2} ight ceil}))
所以(f(x)(g(x)-p(x))equiv 0(mod) (x^{left lceil frac{n}{2} ight ceil}))
因为(f(x) otequiv 0(mod) (x^{left lceil frac{n}{2} ight ceil}))（否则不存在(g(x))，可以证明这样的多项式不存在逆多项式）
所以(g(x)-p(x)equiv 0(mod) (x^{left lceil frac{n}{2} ight ceil}))
将其平方一下得到：
(g^2(x)-2g(x)p(x)+p^2(x)equiv 0(mod) (x^n))（为什么后面mod的东西也平方了？因为(x)次多项式和(y)次多项式相乘后得到的是(x+y)次的多项式）
(g^2(x)equiv 2g(x)p(x)-p^2(x)(mod) (x^n))
(g(x)equiv 2p(x)-frac{p^2(x)}{g(x)}(mod) (x^n))
当然(f(x))与(g(x))是互逆的，所以(g(x)equiv 2p(x)-f(x)p^2(x)(mod) (x^n))
所以(g(x)equiv p(x)(2-f(x)p(x))(mod) (x^n))
到这里，我们发现，我们只要能够求出(p(x))就可以求出(g(x))，那递归求解就好啦

void inv(int *b,int n){
	if (n==1){a[0]=ksm(b[0],mo-2);return;}
	inv(b,(n+1)>>1);
	int bit,len=2;
	for (bit=1;(1<<bit)<(n<<1);bit++)len<<=1; 
	for (int i=0;i<n;i++)d[i]=b[i];
	for (int i=n;i<len;i++)d[i]=0;
	for (int i=0;i<len;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
	ntt(a,len,1);ntt(d,len,1);
	for (int i=0;i<len;i++)a[i]=(2ll-d[i]*1ll*a[i]%mo+mo)%mo*1ll*a[i]%mo;
	ntt(a,len,-1);
	for (int i=n;i<len;i++)a[i]=0;
}

多项式取ln##

相信大家已经熟练掌握了FFT（可能我还没？）我们就来学习一下多项式取ln的艺术吧
对于一个多项式(f(x))，我们想求出(ln(f(x)))
令(g(x)=ln(f(x)))，两边同时求导得：
(g'(x)=frac{f'(x)}{f(x)})
那么(g(x)=int frac{f'(x)}{f(x)})
那我们就只要做一下多项式求导，多项式求逆，多项式积分就好啦