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

题面

题解

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

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

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

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

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

\[\left(\prod_{i=1}^p(x+i)\right)^a\times \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 1\pmod{p}\)恒成立，所以\(x(x^{p-1}-1)\)也有且仅有\(p\)个根\(0,1,2,...,p-1\)

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

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

于是式子就可以写成

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

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

\[\sum_{i=1}^ax^{(p-1)\times i}(-1)^{a-i}{a\choose i} \]

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

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

现在转过头回来考虑

\[(x^{p-1}-1)^a\times \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\n",g);
}
int main(){
//	freopen("testdata.in","r",stdin);
	freopen("math.in","r",stdin);
	freopen("math.out","w",stdout);
	solve();
	return 0;
}