jzoj6005. 【PKUWC2019模拟2019.1.17】数学 （生成函数+FFT+抽代+高精）

题面

题解

幸好咱不是在晚上做的否则咱就不用睡觉了……都什么年代了居然还会出高精的题……

先考虑如果暴力怎么做，令(G(x))为(F(n,k))的生成函数，那么不难发现

[G^R(x)=prod_{i=1}^n(x+i) ]

也就是说如果把(G(x))的系数反过来就是后面那个东西，所以对于(nleq 100000)的数据直接分治(FFT)就行了。不过因为这里的模数不一定满足原根性质，所以要用三模数(NTT)或拆系数(FFT)（所以咱为了这题还特地去学了一下拆系数……）

我们要求的就是这个生成函数有多少项的系数在模(p)意义下不为(0)，设(n=a imes p+b)，那么生成函数可以写成

[left(prod_{i=1}^p(x+i) ight)^a imes prod_{i=1}^b(x+i) ]

然后是一个比较神仙的结论：

[prod_{i=1}^p(x+i)equiv x(x^{p-1}-1)pmod{p} ]

证：对于(prod_{i=1}^p(x+i))，在模(p)意义下有且仅有(p)个根(0,1,2,...,p-1)

根据费马小定理，对于(0<x<p)，(x^{p-1}equiv 1pmod{p})恒成立，所以(x(x^{p-1}-1))也有且仅有(p)个根(0,1,2,...,p-1)

因为(Z_p[x])是唯一分解整环（咱也不知道这是个啥），所以这两个多项式相等

详细的别问咱咱连抽代是啥都不知道

于是式子就可以写成

[x^a(x^{p-1}-1)^a imes prod_{i=1}^b(x+i) ]

把((x^{p-1}-1)^a)用二项式定理展开

[sum_{i=1}^ax^{(p-1) imes i}(-1)^{a-i}{achoose i} ]

考虑({achoose i})，(Lucas)的过程中，本质上就是对(a)和(i)进行(p)进制分解，如果在某一步中(a_k<i_k)，那么({a_kchoose i_k}=0)，所以({achoose i}equiv 0)就代表着(a,i)的(p)进制分解过程中某一位上(i)大于(a)

设(a=(a_0,a_1,...,a_k)_p)，那么满足({achoose i})在模(p)意义下不为(0)的(i)的个数就是(prod_{j=0}^k(a_j+1))

现在转过头回来考虑

[(x^{p-1}-1)^a imes prod_{i=1}^b(x+i) ]

如果(b<p-1)，那么左边式子里的非(0)项的次数都是(p-1)的倍数，那么任意两个这样的项乘上(prod_{i=1}^b(x+i))都不会重复，非(0)项的个数就是左右两边非(0)项个数的积，左边按上面的来，右边用分治(FFT)就行了

如果(b=p-1)，那么(prod_{i=1}^b(x+i)=x^{p-1}-1)，只要计算((x^{p-1}-1)^{a+1})就行了

//minamoto
#include<bits/stdc++.h>
#define R register
#define ll long long
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
const int N=5e5+5,P=998244353;
const double Pi=acos(-1.0);
struct cp{
	double x,y;
	cp(double xx=0,double yy=0){x=xx,y=yy;}
	inline cp operator +(const cp &b)const{return cp(x+b.x,y+b.y);}
	inline cp operator -(const cp &b)const{return cp(x-b.x,y-b.y);}
	inline cp operator *(const cp &b)const{return cp(x*b.x-y*b.y,x*b.y+y*b.x);}
	inline cp operator *(const double &b)const{return cp(x*b,y*b);}
}F[19][N],A[N],B[N],C[N],D[N],X[N],Y[N],Z[N],w[N];
int r[N],n,m,p,res;
void FFT(cp *A,int ty,int lim){
	fp(i,0,lim-1)if(i<r[i])swap(A[i],A[r[i]]);
	for(R int mid=1;mid<lim;mid<<=1)
		for(R int j=0;j<lim;j+=(mid<<1))
			for(R int k=0;k<mid;++k){
				cp x=A[j+k],y=w[mid+k]*A[j+k+mid];
				A[j+k]=x+y,A[j+k+mid]=x-y;
			}
	if(ty==-1){
		reverse(A+1,A+lim);
		double k=1.0/lim;fp(i,0,lim-1)A[i]=A[i]*k;
	}
}
void calc(int ql,int qr,int d){
	if(ql==qr)return F[d][0]=ql,F[d][1]=1,void();
	int mid=(ql+qr)>>1;
	calc(ql,mid,d),calc(mid+1,qr,d+1);
	int lim=1,l=0;while(lim<=qr-ql+1)lim<<=1,++l;
	fp(i,0,lim-1)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(R int i=1;i<lim;i<<=1)fp(k,0,i-1)w[i+k]=cp(cos(Pi*k/i),sin(Pi*k/i));
	fp(i,mid-ql+2,lim-1)F[d][i]=0;
	fp(i,qr-mid+1,lim-1)F[d+1][i]=0;
	fp(i,0,lim-1){
		A[i].x=(ll)(F[d][i].x+0.5)>>15,B[i].x=(ll)(F[d][i].x+0.5)&32767;
		C[i].x=(ll)(F[d+1][i].x+0.5)>>15,D[i].x=(ll)(F[d+1][i].x+0.5)&32767;
		A[i].y=B[i].y=C[i].y=D[i].y=0;
	}
	FFT(A,1,lim),FFT(B,1,lim),FFT(C,1,lim),FFT(D,1,lim);
	fp(i,0,lim-1)
		X[i]=A[i]*C[i],Y[i]=A[i]*D[i]+B[i]*C[i],Z[i]=B[i]*D[i];
	
	FFT(X,-1,lim),FFT(Y,-1,lim),FFT(Z,-1,lim);
	fp(i,0,lim-1){
		F[d][i].x=((((ll)(X[i].x+0.5))%p<<30)+((ll)(Y[i].x+0.5)<<15)+((ll)(Z[i].x+0.5)))%p,F[d][i].y=0;
	}
}
char s[N];int a[N],b[N],ans[N],top,st,ss,len,g;
void change(){
	//高精转进制 
	while(top){
		int res=0;st=top;
		fp(i,1,st){
			res=res*10+a[i];
			b[i]=res/p,res%=p;
		}
		ans[++ss]=res;
		//本轮的余数 
		int i=1;
		while(i<=st&&!b[i])++i;
		top=0;
		fp(j,i,st)a[++top]=b[j];
	}
}
void solve(){
	scanf("%s%d",s+1,&p),len=strlen(s+1);
	int k=1,res=0;
	while(k<=len){
		a[++top]=res/p,res%=p,res=res*10+s[k]-'0',++k;
	}a[++top]=res/p,res%=p;
	if(res!=p-1){
		calc(1,res,0),g=0;
		fp(i,0,res)if((ll)(F[0][i].x+0.5)%p!=0)++g;
	}else ++a[top],g=1;
	change();
	fp(i,1,ss)g=1ll*g*(ans[i]+1)%P;
	printf("%d
",g);
}
int main(){
//	freopen("testdata.in","r",stdin);
	freopen("math.in","r",stdin);
	freopen("math.out","w",stdout);
	solve();
	return 0;
}