[JZOJ6054]【NOI2019模拟2019.3.12】Z的礼物【计数】【斯特林反演】【分治FFT】

Description

Solution

设$p_i=\sum _{j=1}^i{b}_j$
由题意可以列出式子

$a_i=\sum _{j=1}^{i}S(i,j)p_j$
其中S(i,j)为第二类斯特林数

根据斯特林反演的式子
$$F(n)=\sum _{i=0}^{n}S(i,j)G(n)$$
$$G(n)=\sum _{i=0}^n{s}_s(i,j)F(n)或者=\sum _{i=0}^n{s}_u(i,j)(-1{)}^{i-j}F(j)$$

$s_s(n,m),s_u(n,m)$分别表示带符号和无符号第一类斯特林数

可以得到$p_i=\sum _{j=1}^i{s}_s(i,j)a_j$

而根据第一类斯特林数的定义，又有$\prod _{i=0}^{n-1}(x-i)=\sum _{i=0}^n{s}_s(n,i)x^i$

可以直接分治FFT，也可以倍增，求出第l-1行的
然后暴力递推出后面r-l行的

注意由于此处模数不是NTT模数，因此需要使用拆系数的MTT

Code

#include <bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define N 200005
#define M 262144
#define L 18
#define LL long long
#define mo 1000000007
using namespace std;
const double pi=acos(-1);
const LL P=31622;
struct Z
{
	double x,y;
	Z(double _x=0,double _y=0):x(_x),y(_y){};
}wi[M+1],wg[M+1],u1[M+1],u2[M+1],u3[M+1],u4[M+1],ux[M+1],uy[M+1];
Z operator +(Z x,Z y) {return Z(x.x+y.x,x.y+y.y);}
Z operator -(Z x,Z y) {return Z(x.x-y.x,x.y-y.y);}
Z operator *(Z x,Z y) {return Z(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
LL a[M+1],f[N],pr[N],a1[N];
int n,l,r,bit[M+1],cf[21],l2[M+1],st[N],le[N],t[2][N];
void prp(int num)
{
	fo(i,0,num)
	{
		bit[i]=(bit[i>>1]>>1)|((i&1)<<(l2[num]-1));
		wi[i]=wg[i*(M/num)];
	} 
}
void DFT(Z *a,bool pd,int num)
{
	fo(i,0,num-1) if(i<bit[i]) swap(a[i],a[bit[i]]);
	for(int h=1,m=2,lim=num>>1;m<=num;h=m,m<<=1,lim>>=1)
	{
		for(int j=0;j<num;j+=m)
		{
			fo(i,0,h-1)
			{
				Z v=((!pd)?wi[i*lim]:wi[num-i*lim])*a[i+j+h];
				a[i+j+h]=a[i+j]-v;
				a[i+j]=a[i+j]+v;
			}
		}
	}
	if(pd) fo(i,0,num-1) a[i].x/=num,a[i].y/=num; 
}
void rec(Z *a,Z *x,Z *y,int num)
{
	fo(i,0,num-1)
	{
		x[i].x=(a[i].x+a[(num-i)%num].x)/2;
		x[i].y=(a[i].y-a[(num-i)%num].y)/2;
		y[i].x=(a[i].y+a[(num-i)%num].y)/2;
		y[i].y=(a[(num-i)%num].x-a[i].x)/2;
	}
}
void doit(int l,int r)
{
	if(l>=r) return;
	int mid=(l+r)>>1;
	doit(l,mid),doit(mid+1,r);
	int num=cf[l2[le[l]+le[mid+1]-1]];
	prp(num);
	fo(i,0,num-1) ux[i]=uy[i]=u1[i]=u2[i]=u3[i]=u4[i]=Z(0,0);
	fo(i,0,le[l]-1) ux[i].x=a[st[l]+i]/P,ux[i].y=a[st[l]+i]%P;
	fo(i,0,le[mid+1]-1) uy[i].x=a[st[mid+1]+i]/P,uy[i].y=a[st[mid+1]+i]%P;
	DFT(ux,0,num),DFT(uy,0,num);
	rec(ux,u1,u2,num),rec(uy,u3,u4,num);
	fo(i,0,num-1)
	{
		ux[i]=u1[i]*u3[i]+(u1[i]*u4[i])*Z(0,1);
		uy[i]=u2[i]*u4[i]+(u2[i]*u3[i])*Z(0,1);
	}
	DFT(ux,1,num),DFT(uy,1,num);
	le[l]=le[l]+le[mid+1]-1;
	fo(i,0,le[l]-1)
	{
		LL u=ux[i].x+mo+0.5,v=ux[i].y+mo+0.5,p=uy[i].y+mo+0.5,q=uy[i].x+mo+0.5;
		u%=mo,v%=mo,p%=mo,q%=mo;
		a[st[l]+i]=((u*P%mo*P%mo+(v+p)%mo*P%mo+q)%mo+mo)%mo;
	}
}
int main()
{
	cin>>n>>l>>r;
	fo(i,0,M) wg[i]=Z(cos(2*i*pi/M),sin(2*i*pi/M));
	fo(i,1,n) scanf("%lld",&a1[i]);
	cf[0]=1;
	fo(i,1,18) l2[cf[i]=cf[i-1]<<1]=i;
	fod(i,M-1,2) if(!l2[i]) l2[i]=l2[i+1];
	int len=0;
	fo(i,0,l-2)
	{
		st[i]=len;
		le[i]=2;
		a[len]=mo-i;
		a[len+1]=1;
		len+=2;
	}
	doit(0,l-2);
	if(l==1) a[0]=1;
	fo(i,0,l-1) pr[l-1]=(pr[l-1]+a[i]*a1[i]%mo+mo)%mo,t[1][i]=a[i];
	fo(i,l,r)
	{
		int i1=(i-l)&1;
		t[i1][0]=0;
		fo(j,1,i) t[i1][j]=(t[1-i1][j-1]-t[1-i1][j]*(LL)(i-1))%mo,pr[i]=(pr[i]+a1[j]*t[i1][j]%mo+mo)%mo;
		printf("%lld ",(pr[i]-pr[i-1]+mo+mo)%mo);
	}
}