日常学习——FFT

FFT相关详解：

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

https://blog.csdn.net/ggn_2015/article/details/68922404

练习题：

bzoj2179 FFT快速傅立叶:

模板题

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define inf 0x7f7f7f7f
using namespace std;
const int maxn=2e5;
const double pi=acos(-1);
int n,m,l;
int sum[maxn*2+10],rev[maxn+10];

int read()
{
	char ch;int x=0,f=1;
	for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
	return x*f;
}

struct complex
{
	double x,y;
	complex(){};
	complex(double _x,double _y){x=_x,y=_y;}
	complex operator +(const complex &a){return complex(x+a.x,y+a.y);}
	complex operator -(const complex &a){return complex(x-a.x,y-a.y);}
	complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}A[maxn+10],B[maxn+10],Ans[maxn*2+10];

int get()
{
	int ch=getchar();
	for (;ch<'0'||ch>'9';ch=getchar());
	return ch-'0';
}

void fft(complex *a,int len,int flag)
{
	for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);
	for (int i=2;i<=len;i<<=1)
	{
		complex wn(cos(2*pi/i),sin(2*pi/i*flag));
		for (int j=0;j<len;j+=i)
		{
			complex nw(1,0);
			for (int k=0;k<(i>>1);k++,nw=nw*wn)
			{
				complex x=a[j+k],y=nw*a[j+k+(i>>1)];
				a[j+k]=x+y,a[j+k+(i>>1)]=x-y;
			}
		}
	}
	if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len;
}

int main()
{
	n=read();
	for (int i=n-1;~i;i--) A[i].x=get();
	for (int i=n-1;~i;i--) B[i].x=get();
	for (n*=2,m=1;m<=n;m<<=1)l++; 
	for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	fft(A,m,1),fft(B,m,1);
	for (int i=0;i<m;i++) Ans[i]=A[i]*B[i];
	fft(Ans,m,-1);
	for (int i=0;i<m;i++) sum[i]=(int)(Ans[i].x+0.5);
	for (int i=0;i<m;i++)
	{
		sum[i+1]+=sum[i]/10;
		sum[i]%=10;
	}
	while((!sum[n-1])&&n) n--;
	while(sum[n]) n++;
	for (int i=n-1;~i;i--) putchar(sum[i]+'0');
	printf("
");
	return 0;
}

bzoj2194 快速傅立叶之二:

模板题+1。计算(C_k=sum_{i=k}^{n-1}{a_i imes b_{i-k}}),将(b[])翻转，则(C_k=sum_{i=k}^{n-1}{a_i imes b_{n+1+k-i}})，即(C_k=sumlimits_{i+j=n+1+k,0<=i,j<=2 imes n}{a_i imes b_j})。直接FFT即可。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define inf 0x7f7f7f7f
using namespace std;
const int maxn=4e5;
const double pi=acos(-1);
int n,m,l;
int rev[(maxn<<2)+5];

int read()
{
	char ch;int x=0,f=1;
	for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
	return x*f;
}

struct complex
{
	double x,y;
	complex(){};
	complex(double _x,double _y){x=_x,y=_y;}
	complex operator +(const complex &a){return complex(x+a.x,y+a.y);};
	complex operator -(const complex &a){return complex(x-a.x,y-a.y);};
	complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);};
}a[maxn+10],b[maxn+10],ans[(maxn<<2)+5];

void fft(complex *a,int len,int flag)
{
	for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);
	for (int i=2;i<=len;i<<=1)
	{
		complex wn(cos(2*pi/i),sin(2*pi/i*flag));
		for (int j=0;j<len;j+=i)
		{
			complex nw(1,0);
			for (int k=0;k<(i>>1);k++,nw=nw*wn)
			{
				complex x=a[j+k],y=nw*a[j+k+(i>>1)];
				a[j+k]=x+y,a[j+k+(i>>1)]=x-y;
			}
		}
	}
	if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len;
}

int main()
{
	n=read();
	for (int i=0;i<n;i++) a[i].x=read(),b[n+1-i].x=read();
	for (m=1;m<=n*2;m<<=1)l++;
	for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	fft(a,m,1),fft(b,m,1);
	for (int i=0;i<m;i++) ans[i]=a[i]*b[i];
	fft(ans,m,-1);
	for (int i=n+1;i<=n*2;i++) printf("%d
",(int)(ans[i].x+0.5));
	return 0;
}

bzoj4827 [Hnoi2017]礼物:

假设移动a位，整体增加c，则答案为(sum_{i=1}^{n}{(x_{((i+a-1) \% n+1)}-y_{i}+c)^{2}})。把此式展开(sum _{i=1} ^{n} {x_{( (i+a-1) \% n+1 )} ^{2} +y_{i} ^{2} +c ^{2} +2 imes c imes x_{( (i+a-1) \% n+1 )} -2 imes y_{i} imes c -2 imes x_{( (i+a-1) \% n+1)} imes y_{i}})。

(sum _{i=1} ^{n} {x_{( (i+a-1) \% n+1 )} imes y_{i}})可以用FFT求，之后枚举a和c就行了。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
#define inf 0x7f7f7f7f
using namespace std;
const int maxn=5e4;
const double pi=acos(-1);
int n,m,k,l,s,ans=1e9;
int a[maxn+10],b[maxn+10],rev[maxn*4+10];

int read()
{
	char ch;int x=0,f=1;
	for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
	return x*f;
}

struct complex
{
	double x,y;
	complex(){};
	complex (double _x,double _y){x=_x,y=_y;};
	complex operator +(const complex &a){return complex(x+a.x,y+a.y);}
	complex operator -(const complex &a){return complex(x-a.x,y-a.y);}
	complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}x1[maxn*4+10],x2[maxn*4+10];

void fft(complex *a,int len,int flag)
{
	for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);
	for (int i=2;i<=len;i<<=1)
	{
		complex wn(cos(2*pi/i),sin(2*pi/i*flag));
		for (int j=0;j<len;j+=i)
		{
			complex nw(1,0);
			for (int k=0;k<(i>>1);k++,nw=nw*wn)
			{
				complex x=a[j+k],y=nw*a[j+k+(i>>1)];
				a[j+k]=x+y,a[j+k+(i>>1)]=x-y;
			}
		}
	}
	if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len;
}

int main()
{
	n=read(),k=read();
	for (int i=0;i<n;i++) a[i]=read(),x1[i].x=a[i];
	for (int i=0;i<n;i++) b[i]=read(),x2[n-i-1].x=b[i],s+=b[i];
	for (m=1;m<=n*2;m<<=1)l++;
	for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	fft(x1,m,1),fft(x2,m,1);
	for (int i=0;i<m;i++) x1[i]=x1[i]*x2[i];
	fft(x1,m,-1);
	for (int c=0;c<=k;c++)
	{
		int sum=0;
		for (int i=0;i<n;i++) sum+=(a[i]+c)*(a[i]+c)+b[i]*b[i];
		ans=min(ans,sum-2*((int)(x1[n-1].x+0.5)+s*c));
		for (int i=1;i<n;i++)
		{
			int tmp=(int)(x1[n+i-1].x+0.5)+(int)(x1[i-1].x+0.5)+s*c;
			ans=min(ans,sum-2*tmp);
		}
	}
	printf("%d
",ans);
	return 0;
}

bzoj4555 [Tjoi2016&Heoi2016]求和

求(f_n=sum_{i=0}^{n}sum_{j=0}^{i}S_{i,j} imes 2^j imes j!),对答案(\%98244353)。其中(S_{i,j})为第二类斯特林数，意义为将n个不同的元素拆分成m个集合的方案数,所以当m>n时(S_{n,m})为0。而(S_{n,m} imes m!)就是将n个不同的元素拆分成m个不同的集合的方案数。用容斥可得到(S_{n,m} imes m!=sum_{i=0}^{m}{(-1)^i imes inom{m}{i} imes (m-i)^n})。

好，那我们开始推式子。
(f_n=sum_{i=0}^{n}sum_{j=0}^{i}S_{i,j} imes 2^j imes j!)
(f_n=sum_{i=0}^{n}sum_{j=0}^{n}{2^j imessum_{k=0}^{j}}(-1)^k imes inom{j}{k} imes (j-k)^i)
(f_n=sum_{j=0}^{n}2^j imes sum_{k=0}^{j}(-1)^k imes frac{j!}{k!(j-k)!} imes sum_{i=0}^{n} (j-k)^i​)
(f_n=sum_{j=0}^{n} 2^j imes j! imes sum_{k=0}^{j} frac{(-1)^k}{k!} imes frac{sum_{i=0}^{n}(j-k)^i}{(j-k)!})

后面的(sum_{k=0}^{j} frac{(-1)^k}{k!} imes frac{sum_{i=0}^{n}(j-k)^i}{(j-k)!})就是一个卷积的形式，因为要膜，所以直接上NTT即可。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
#define inf 0x7f7f7f7f
using namespace std;
const int mod=998244353;
const int maxn=1e5+5;
int ans,n,m,l;
int rev[maxn<<2],inv[maxn<<2],x[maxn<<2],y[maxn<<2];

int read()
{
	char ch;int x=0,f=1;
	for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
	return x*f;
}	

int power(int a,int k)
{
	int sum=1;
	for (;k;k>>=1,a=1ll*a*a%mod)
	if (k&1)
		sum=1ll*sum*a%mod;
	return sum;
}

void ntt(int *a,int len,int flag)
{
	for (int i=0;i<len;i++) if (rev[i]<i) swap(a[rev[i]],a[i]);
	for (int i=2;i<=len;i<<=1)
	{
		int gn=power(3,((mod-1)+flag*(mod-1)/i)%(mod-1));
		for (int j=0;j<len;j+=i)
		{
			int ng=1;
			for (int k=0;k<(i>>1);k++,ng=1ll*ng*gn%mod)
			{
				int x=a[j+k],y=1ll*ng*a[j+k+(i>>1)]%mod;
				a[j+k]=(x+y)%mod,a[j+k+(i>>1)]=1ll*((x-y)%mod+mod)%mod;
			}
		}
	}
	if (flag==-1)
	{
		int invlen=power(len,mod-2);
		for (int i=0;i<len;i++) x[i]=1ll*x[i]*invlen%mod;
	}
}

int main()
{
	n=read();
	for (m=1;m<n*2;m<<=1)l++;
	for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	inv[0]=1;
	for (int i=1;i<=n;i++) inv[i]=1ll*inv[i-1]*i%mod;
	inv[n]=power(inv[n],mod-2);
	for (int i=n-1;i;i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
	for (int i=0;i<=n;i++)  x[i]=(inv[i]*(i&1?-1:1)+mod)%mod;
	y[0]=1,y[1]=n+1;
	for (int i=2;i<=n;i++) y[i]=1ll*(power(i,n+1)-1)*power(i-1,mod-2)%mod*inv[i]%mod;
	ntt(x,m,1),ntt(y,m,1);
	for (int i=0;i<m;i++) x[i]=1ll*x[i]*y[i]%mod;
	ntt(x,m,-1);
	int two_power=1,fact=1;ans=x[0];
	for (int i=1;i<=n;i++)
	{
		two_power=1ll*two_power*2%mod,fact=1ll*fact*i%mod;
		ans=(ans+1ll*two_power*fact%mod*x[i]%mod)%mod;
	}
	printf("%d
",ans);
	return 0;
}