12
3
2014
1

论逗逼的自我修养之FFT练习记

最近做大爷们出的训练题,基本每天一道FFT...但是似乎我对FFT的理解只停留在“可以求两个数组的卷积”+背代码阶段,感觉再不重新学一遍要跟不上时代的脚步了...

[12.7]差不多就到这儿吧。。感觉把教材pyx的博客上的每一块都看了一遍,但是感觉做题还是不大会用。。如果以后闲着无聊想自虐的话可能会去写一下策爷的[tex]e^{A_{(z)}}[/tex]。。(不过我感觉以我的智商估计是写不出来了)

1.首先是最基本的FFT:

重新看了一遍算法导论,感觉也不是特别难理解,果然当初学的时候还是太naive了

【BZOJ2194】FFT快速傅立叶

【BZOJ2179】快速傅立叶之二

void fft(atom *x,int n,int fl=1){
	for (int i=(n>>1),j=1;j<n;j++){
		if (i<j) swap(x[i],x[j]);
		int k;
		for (k=(n>>1);i&k;i^=k,k>>=1); i^=k;
	}
	for (int m=2;m<=n;m=(m<<1)){
		atom w=(atom){cos(pi*2*fl/m),sin(pi*2*fl/m)};
		for (int i=0;i<n;i+=m){
			atom cur=(atom){1,0};
			for (int j=i;j<i+(m>>1);j++){
				atom u=x[j],v=x[j+(m>>1)]*cur;
				x[j]=u+v; x[j+(m>>1)]=u-v; cur=cur*w;
			}
		}
	}
}

2.FNT

想起来感觉还是挺显然的,在模意义下用模意义下的单位根替换。

【某模拟题T1】循环计数 mx爷的题,求斯特林数。把斯特林数的那个式子用分治FNT搞就好了(我偷懒空间是[tex]O(n \log n)[/tex]的,稍微改一改就可以做到[tex]O(n)[/tex]

void fft(int *x,int n,int fl=1){
	for (int i=(n>>1),j=1;j<n;j++){
		if (i<j) swap(x[i],x[j]);
		int k;
		for (k=(n>>1);k&i;i^=k,k>>=1); i^=k;
	}
	int now=0;
	for (int m=2;m<=n;m<<=1){
		int w; now++; if (fl>0) w=G[now]; else w=nG[now]; 
		for (int i=0;i<n;i+=m){
			int cur=1;
			for (int j=i;j<i+(m>>1);j++){
				int u=x[j],v=1ll*x[j+(m>>1)]*cur%mo;
				x[j]=(u+v)%mo; x[j+(m>>1)]=(u-v+mo)%mo;
				cur=1ll*cur*w%mo;
			}
		}
	}
}
atom solve(int l,int r){
	if (l==r){
		C[++tot]=l; C[++tot]=1; return (atom){tot-1,tot};
	}
	int mid=(l+r)>>1; atom k1=solve(l,mid),k2=solve(mid+1,r);
	int n=max(mid-l+1,r-mid),len=1;
	while ((1<<len)<=(n<<1)) len++;
	for (int i=0;i<(1<<len);i++){A[i]=0; B[i]=0;}
	for (int i=k1.l;i<=k1.r;i++)A[i-k1.l]=C[i]; 
	for (int i=k2.l;i<=k2.r;i++)B[i-k2.l]=C[i];
	fft(A,1<<len); fft(B,1<<len);
	for (int i=0;i<(1<<len);i++) A[i]=1ll*A[i]*B[i]%mo;
	fft(A,1<<len,-1);
	for (int i=0;i<(1<<len);i++) A[i]=1ll*A[i]*two[len]%mo;
	atom ans; ans.l=tot+1;
	for (int i=0;i<=r-l+1;i++) C[++tot]=A[i]; ans.r=tot; 
	return ans;
}

3.分治FFT

利用CDQ分治的思想,FFT可以求解形如[tex]f_n=\sum_{i=1}^{n-1}f_i \times A_{n-i}[/tex]的递推式。

【BZOJ3456】城市规划 可以得到递推式[tex]f_n=2^{\binom{n}{2}}- \sum_{i=1}^{n-1}2^{\binom{n-i}{2}} \times \binom{n-1}{i-1} \times f_i[/tex],意义就是用所有图减去不连通的图的数量,强制1在一个联通块中枚举这个块的大小。之后转化成卷积的形式就可以分治FFT了。

void fft(int *x,int n,int fl=1){
	for (int i=(n>>1),j=1;j<n;j++){
		if (i<j) swap(x[i],x[j]);
		int k;
		for (k=(n>>1);k&i;i^=k,k>>=1); i^=k;
	}
	int now=0;
	for (int m=2;m<=n;m<<=1){
		int w; now++; if (fl==1) w=G[now]; else w=nG[now];
		for (int i=0;i<n;i+=m){
			int cur=1;
			for (int j=i;j<i+(m>>1);j++){
				int u=x[j],v=1ll*x[j+(m>>1)]*cur%mo;
				x[j]=(u+v)%mo; x[j+(m>>1)]=(u-v+mo)%mo;
				cur=1ll*cur*w%mo;
			}
		}
	}
}
void solve(int l,int r){
	if (l==r){
		f[l]=(T[l]-1ll*I[l-1]*f[l]%mo+mo)%mo; return;
	}
	int mid=l+r>>1; solve(l,mid);
	int n=max(mid-l+1,r-mid),len=1;
	while ((1<<len)<=(n<<1)) len++;
	for (int i=0;i<(1<<len);i++){A[i]=0; B[i]=0;}
	for (int i=l;i<=mid;i++)A[i-l]=1ll*nI[i-1]*f[i]%mo; 
	for (int i=1;i<=r-l;i++)B[i]=1ll*T[i]*nI[i]%mo;
	fft(A,1<<len); fft(B,1<<len);
	for (int i=0;i<(1<<len);i++) A[i]=1ll*A[i]*B[i]%mo;
	fft(A,1<<len,-1);
	for (int i=mid+1;i<=r;i++) f[i]=(f[i]+1ll*A[i-l]*two[len])%mo;
	solve(mid+1,r);
}

4.多项式求逆

看了pyx的博客(感觉pyx的博客已经可以成为中国多项式从入门到提高到进阶三合一通用教材了)后感觉其实如果想到了倍增的话剩下的还是不怎么难想的。。不过似乎这个模型比较奇怪?感觉转化起来有点费劲。似乎满足[tex]f_n=C_{n}-\sum_{i=1}^{n-1}f_i \times A_{n-i} \times B_{n}[/tex],如果[tex]B_n[/tex]可以把[tex]f_n[/tex]形式转化为[tex]A_{n-i}[/tex]相同的话就可以用这个了。

【BZOJ3456】城市规划

可以得到[tex] \sum_{i=1}^n 2^{\binom{n-i}{2}} \times (n-i)!^{-1} \times (i-1)!^{-1} \times f_i = 2^{\binom{n}{2}} \times (n-1)!^{-1}[/tex],然后左边可以看成卷积,就可以把一个多项式除过去(就是乘上逆元)。

一直听别人说这玩意常数巨大,跑得还没有分治FFT快。。实际写了一下速度还是快了那么一点点的,但是复杂度的优势完全没了。

void fft(int *x,int n,int fl=1){
	for (int i=(n>>1),j=1;j<n;j++){
		if (i<j) swap(x[i],x[j]);
		int k;
		for (k=(n>>1);i&k;i^=k,k>>=1); i^=k;
	}
	int now=0;
	for (int m=2;m<=n;m<<=1){
		int w; now++; if (fl==1) w=G[now]; else w=nG[now];
		for (int i=0;i<n;i+=m){
			int cur=1;
			for (int j=i;j<i+(m>>1);j++){
				int u=x[j],v=1ll*x[j+(m>>1)]*cur%mo;
				x[j]=(u+v)%mo; x[j+(m>>1)]=(u-v+mo)%mo;
				cur=1ll*cur*w%mo;
			}
		}
	}
}
void getrev(int n){
	if (n==1){
		A[0]=quick(C[0],mo-2); return;
	}
	getrev((n+1)/2); int len=1;
	while ((1<<len)<=((n+1)<<1)) len++;
	for (int i=0;i<n;i++) B[i]=C[i]; for (int i=n;i<(1<<len);i++) B[i]=0;
	fft(A,1<<len); fft(B,1<<len);
	for (int i=0;i<(1<<len);i++) B[i]=(2-1ll*A[i]*B[i]%mo+mo)%mo;
	for (int i=0;i<(1<<len);i++) A[i]=1ll*A[i]*B[i]%mo;
	fft(A,1<<len,-1);
	for (int i=0;i<n;i++) A[i]=1ll*A[i]*two[len]%mo;
	for (int i=n;i<(1<<len);i++) A[i]=0;
}

5.多项式除法

依然是pyx的博客...关键就是把多项式reverse,后面的就挺好推了。

练习嘛写了一个多项式多点求值愉♂悦一下,写出来之后妈呀怎么10W10W要跑10S..实在是被慢哭了。(策爷:这种常数巨大的东西没前途了


6.多项式开根

仍旧是pyx的博客。。依然是倍增的思想,考虑平方后把式子转化为[tex]A_{(x)}^2 \equiv B_{(x)}(\text{mod}\ x^{n})[/tex]的形式就好了,还是挺好推的。

【codeforces 250E】小朋友和二叉树 

接一个方程可以得到母函数的式子是[tex]F_{(x)}=\frac{2}{1+ \sqrt{1-4A_{(x)}}}[/tex],其中[tex]A_{(x)}[/tex]是那个集合中的数组成的函数。然后就用多项式开根加求逆算出来[tex]F_{(x)}[/tex]就好了。

void fft(int *x,int n,int fl=1){
	for (int i=(n>>1),j=1;j<n;j++){
		if (i<j) swap(x[i],x[j]);
		int k;
		for (k=(n>>1);i&k;i^=k,k>>=1); i^=k;
	}
	int now=0;
	for (int m=2;m<=n;m<<=1){
		int w; now++; if (fl==1) w=G[now]; else w=nG[now];
		for (int i=0;i<n;i+=m){
			int cur=1;
			for (int j=i;j<i+(m>>1);j++){
				int u=x[j],v=1ll*x[j+(m>>1)]*cur%mo;
				x[j]=(u+v)%mo; x[j+(m>>1)]=(u-v+mo)%mo;
				cur=1ll*cur*w%mo;
			}
		}
	}
}
void getrev(int n){
	if (n==1){
		rev[0]=quick(fir[0],mo-2); return;
	}
	getrev((n+1)>>1); int len=1;
	while ((1<<len)<=(n<<1)) len++;
	for (int i=((n+1)>>1);i<(1<<len);i++) rev[i]=0;
	for (int i=0;i<n;i++) A[i]=fir[i];
	for (int i=n;i<(1<<len);i++) A[i]=0;
	fft(rev,1<<len); fft(A,1<<len);
	for (int i=0;i<(1<<len);i++) A[i]=(2-1ll*A[i]*rev[i]%mo+mo)%mo;
	for (int i=0;i<(1<<len);i++) rev[i]=1ll*rev[i]*A[i]%mo;
	fft(rev,1<<len,-1);
	for (int i=0;i<n;i++) rev[i]=1ll*rev[i]*two[len]%mo;
	for (int i=n;i<(1<<len);i++) rev[i]=0;
}
void getroot(int n){
	if (n==1){
		root[0]=1; return;
	}
	getroot((n+1)>>1); int len=1,m=((n+1)>>1),n2=two[1];
	while ((1<<len)<=(n<<1)) len++;
	for (int i=0;i<(1<<len);i++){rev[i]=0; fir[i]=0;}
	for (int i=0;i<m;i++) fir[i]=root[i]; 
	getrev(n);
	for (int i=0;i<(1<<len);i++) A[i]=0;
	for (int i=0;i<n;i++) A[i]=C[i];
	fft(fir,1<<len); fft(A,1<<len); fft(rev,1<<len);
	for (int i=0;i<(1<<len);i++) A[i]=(1ll*fir[i]*n2%mo+1ll*A[i]*rev[i]%mo*n2%mo)%mo;
	fft(A,1<<len,-1);
	for (int i=0;i<n;i++) root[i]=1ll*A[i]*two[len]%mo;
	for (int i=n;i<(1<<len);i++) root[i]=0;
}

 

Category: 论逗逼的自我修养系列 | Tags: | Read Count: 10456
Avatar_small
ssali 说:
2022年5月22日 00:35

Great stuff. i am going to bookmark it, if you are looking for lyca beltegoed, lebara beltegoed, lyca bundel xs, samsung a52, samsung a51, samsung a72, iphone 13 pro max in Netherlands, get in touch with us.

 


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter

Host by is-Programmer.com | Power by Chito 1.3.3 beta | Theme: Aeros 2.0 by TheBuckmaker.com