【HDU5730】Shell Necklace（多项式运算，分治FFT）

【HDU5730】Shell Necklace（多项式运算，分治FFT）

题面

Vjudge
翻译：
有一个长度为(n)的序列
已知给连续的长度为(i)的序列装饰的方案数为(a[i])
求将(n)个位置全部装饰的总方案数。
答案(mod 313)

题解

很明显，是要求：
(f[n]=sum_{i=0}^na[i] imes f[n-i],f[0]=0)
卷积的形式啊。。

然后就可以开始搞了

忍不住的方法一

好明显啊，把生成函数(F,A)给搞出来
然后就有(F*A+1=F)
然后((1-A)F=1)
然后(F=frac{1}{1-A})
多项式求逆，没了。。。
但是这个(mod 313)很蛋疼啊。。
直接蒯模板了

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
using namespace std;
#define ll long long
#define ull unsigned long long
#define RG register
#define MAX 288888
#define MOD (313)
const double Pi=acos(-1);
const int m=sqrt(MOD);
inline int read()
{
    RG int x=0,t=1;RG char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
struct Complex{double a,b;}W[MAX],A[MAX],B[MAX],C[MAX],D[MAX];
Complex operator+(Complex a,Complex b){return (Complex){a.a+b.a,a.b+b.b};}
Complex operator-(Complex a,Complex b){return (Complex){a.a-b.a,a.b-b.b};}
Complex operator*(Complex a,Complex b){return (Complex){a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a};}
int r[MAX];
void FFT(Complex *P,int N,int opt)
{
	for(int i=0;i<N;++i)if(i<r[i])swap(P[i],P[r[i]]);
	for(int i=1;i<N;i<<=1)
		for(int p=i<<1,j=0;j<N;j+=p)
			for(int k=0;k<i;++k)
			{
				Complex w=(Complex){W[N/i*k].a,W[N/i*k].b*opt};
				Complex X=P[j+k],Y=P[i+j+k]*w;
				P[j+k]=X+Y;P[i+j+k]=X-Y;
			}
	if(opt==-1)for(int i=0;i<N;++i)P[i].a/=1.0*N;
}
void Multi(int *a,int *b,int len,int *ret)
{
	for(int i=0;i<(len<<1);++i)A[i]=B[i]=C[i]=D[i]=(Complex){0,0};
	for(int i=0;i<len;++i)
	{
		a[i]%=MOD;b[i]%=MOD;
		A[i]=(Complex){(a[i]/m)*1.0,0};
		B[i]=(Complex){(a[i]%m)*1.0,0};
		C[i]=(Complex){(b[i]/m)*1.0,0};
		D[i]=(Complex){(b[i]%m)*1.0,0};
	}
	int N,l=0;
	for(N=1;N<=len;N<<=1)++l;
	for(int i=0;i<N;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=1;i<N;i<<=1)
		for(int k=0;k<i;++k)W[N/i*k]=(Complex){cos(k*Pi/i),sin(k*Pi/i)};
	FFT(A,N,1);FFT(B,N,1);FFT(C,N,1);FFT(D,N,1);
	for(int i=0;i<N;++i)
	{
		Complex tmp=A[i]*C[i];
		C[i]=B[i]*C[i],B[i]=B[i]*D[i],D[i]=D[i]*A[i];
		A[i]=tmp;C[i]=C[i]+D[i];
	}
	FFT(A,N,-1);FFT(B,N,-1);FFT(C,N,-1);
	for(int i=0;i<len;++i)
	{
		ret[i]=0;
		ret[i]=(ret[i]+1ll*(ll)(A[i].a+0.5)%MOD*m%MOD*m%MOD)%MOD;
		ret[i]=(ret[i]+1ll*(ll)(C[i].a+0.5)%MOD*m%MOD)%MOD;
		ret[i]=(ret[i]+1ll*(ll)(B[i].a+0.5)%MOD)%MOD;
		ret[i]=(ret[i]+MOD)%MOD;
	}
}
int c[MAX],d[MAX];
void Inv(int *a,int *b,int len)
{
	if(len==1){b[0]=fpow(a[0],MOD-2);return;}
	Inv(a,b,len>>1);
	Multi(a,b,len,c);
	Multi(c,b,len,d);
	for(int i=0;i<len;++i)b[i]=(b[i]+b[i])%MOD;
	for(int i=0;i<len;++i)b[i]=(b[i]+MOD-d[i])%MOD;
}
int n,a[MAX],b[MAX];
int main()
{
	while(n=read())
	{
		for(int i=1;i<=n;++i)a[i]=read()%MOD;
		for(int i=1;i<=n;++i)a[i]=(MOD-a[i])%MOD;
		a[0]++;
		int N;for(N=1;N<=n;N<<=1);
		Inv(a,b,N);
		printf("%d
",b[n]);
		memset(a,0,sizeof(a));memset(b,0,sizeof(b));
	}
	return 0;
}

很生疏的方法二

这种东西显然也可以(CDQ)分治+(FFT)来做
简称分治(FFT)
怎么考虑？
式子长成这样：
(f[n]=sum_{i=0}^na[i] imes f[n-i],f[0]=0)
一定要求出前面的才能计算后面的。
每次考虑左半段对于右半段的贡献
直接拿左半段和右半段，（要求的东西）做卷积
这样，对于每一项的系数，都是右半段每个点的一部分贡献。累加即可。
复杂度(O(nlog^2n))
模数(313)依旧蛋疼。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
using namespace std;
#define ll long long
#define ull unsigned long long
#define RG register
#define MAX 288888
#define MOD (313)
const double Pi=acos(-1);
const int m=sqrt(MOD);
inline int read()
{
    RG int x=0,t=1;RG char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
struct Complex{double a,b;}W[MAX],A[MAX],B[MAX],C[MAX],D[MAX];
Complex operator+(Complex a,Complex b){return (Complex){a.a+b.a,a.b+b.b};}
Complex operator-(Complex a,Complex b){return (Complex){a.a-b.a,a.b-b.b};}
Complex operator*(Complex a,Complex b){return (Complex){a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a};}
int r[MAX];
void FFT(Complex *P,int N,int opt)
{
	for(int i=0;i<N;++i)if(i<r[i])swap(P[i],P[r[i]]);
	for(int i=1;i<N;i<<=1)
		for(int p=i<<1,j=0;j<N;j+=p)
			for(int k=0;k<i;++k)
			{
				Complex w=(Complex){W[N/i*k].a,W[N/i*k].b*opt};
				Complex X=P[j+k],Y=P[i+j+k]*w;
				P[j+k]=X+Y;P[i+j+k]=X-Y;
			}
	if(opt==-1)for(int i=0;i<N;++i)P[i].a/=1.0*N;
}
void Multi(int *a,int *b,int len,int *ret)
{
	for(int i=0;i<(len<<1);++i)A[i]=B[i]=C[i]=D[i]=(Complex){0,0};
	for(int i=0;i<len;++i)
	{
		a[i]%=MOD;b[i]%=MOD;
		A[i]=(Complex){(a[i]/m)*1.0,0};
		B[i]=(Complex){(a[i]%m)*1.0,0};
		C[i]=(Complex){(b[i]/m)*1.0,0};
		D[i]=(Complex){(b[i]%m)*1.0,0};
	}
	int N,l=0;
	for(N=1;N<=len;N<<=1)++l;
	for(int i=0;i<N;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=1;i<N;i<<=1)
		for(int k=0;k<i;++k)W[N/i*k]=(Complex){cos(k*Pi/i),sin(k*Pi/i)};
	FFT(A,N,1);FFT(B,N,1);FFT(C,N,1);FFT(D,N,1);
	for(int i=0;i<N;++i)
	{
		Complex tmp=A[i]*C[i];
		C[i]=B[i]*C[i],B[i]=B[i]*D[i],D[i]=D[i]*A[i];
		A[i]=tmp;C[i]=C[i]+D[i];
	}
	FFT(A,N,-1);FFT(B,N,-1);FFT(C,N,-1);
	for(int i=0;i<len;++i)
	{
		ret[i]=0;
		ret[i]=(ret[i]+1ll*(ll)(A[i].a+0.5)%MOD*m%MOD*m%MOD)%MOD;
		ret[i]=(ret[i]+1ll*(ll)(C[i].a+0.5)%MOD*m%MOD)%MOD;
		ret[i]=(ret[i]+1ll*(ll)(B[i].a+0.5)%MOD)%MOD;
		ret[i]=(ret[i]+MOD)%MOD;
	}
}
int n,a[MAX],x[MAX],y[MAX],f[MAX];
void CDQ(int l,int r)
{
	if(l==r){f[l]=(f[l]+a[l])%MOD;return;}
	int mid=(l+r)>>1,N;
	CDQ(l,mid);for(N=1;N<=(r-l+1);N<<=1);
	for(int i=0;i<=N;++i)x[i]=y[i]=0;
	for(int i=l;i<=mid;++i)x[i-l]=f[i];
	for(int i=0;i<=r-l;++i)y[i]=a[i];
	Multi(x,y,N,x);
	for(int i=mid+1;i<=r;++i)f[i]=(f[i]+x[i-l])%MOD;
	CDQ(mid+1,r);
}
int main()
{
	while(n=read())
	{
		for(int i=1;i<=n;++i)a[i]=read()%MOD;
		int N;for(N=1;N<=n;N<<=1);
		CDQ(1,n);
		printf("%d
",f[n]);
		memset(f,0,sizeof(f));
		memset(a,0,sizeof(a));
	}
	return 0;
}