[51nod 1538] 一道难题 -- 特征多项式,多项式取模

问题相当于选择一些球排列起来使得总权值为n的方案数

显然有f[n]=∑f[n-ai]

这里ai会到23333,直接写O(n^2logn)的暴力解法会TLE 考虑用多项式取模+倍增求x^n mod h(x)

卡常 时限5s 一开始本机跑了20s 就是跑不过 注意快速幂的时候要从高位往低位算 这样只要算logn次取模而不是2logn次 还有取模的时候没有必要的数字就删掉 这样可以一直让长度在65536

预处理旋转因子会快很多

#include"bits/stdc++.h"

using namespace std;
typedef long long LL;
const int N=66666,M=1000005,P=104857601,g=3,L=23333;

int power(int a,int t){
	int r=1;
	while(t){
		if(t&1)r=(LL)r*a%P;
		a=(LL)a*a%P;t>>=1;
	}
	return r;
}

int _wn[25],z[25][N];

void NTT(int A[],int len,int dft){
	int i,j=len>>1,k,l,u,v,inv;
	for(i=1;i<len-1;i++){
		if(i<j)swap(A[i],A[j]);
		for(k=len>>1;(j^=k)<k;k>>=1);
	}
	if(!_wn[0]){
		int c=0,w,wn;
		for(_wn[i=22]=39193363;i;i--)_wn[i-1]=(LL)_wn[i]*_wn[i]%P;
		for(l=2;l<=(1<<16);l<<=1)for(i=l>>1,wn=_wn[++c],j=0;j<(1<<16);j+=l)for(w=1,k=j;k<j+i;k++,w=(LL)w*wn%P)z[c][k]=w;
	}
	int c=0,*w;
	for(l=2;l<=len;l<<=1)
		for(i=l>>1,w=z[++c],j=0;j<len;j+=l)
			for(k=j;k<j+i;k++)
				u=A[k],v=(LL)A[k+i]*w[k]%P,A[k]=(u+v)%P,A[k+i]=(u-v+P)%P;
	if(dft==-1){
		inv=power(len,P-2);
		for(i=0;i<len;i++)A[i]=(LL)A[i]*inv%P;
		for(i=1;i<len/2;i++)swap(A[i],A[len-i]); 
	}
}

int w1[N],w2[N],w3[N];

void convol(int A[],int B[],int R[],int len){
	memcpy(w1,A,len<<2);
	memcpy(w2,B,len<<2);
	NTT(w1,len,1);
	NTT(w2,len,1);
	for(int i=0;i<len;i++)w3[i]=(LL)w1[i]*w2[i]%P;
	NTT(w3,len,-1);
	memcpy(R,w3,len<<2);
}

int wi1[N],wi2[N],wi3[N],wi4[N];

void inverse(int A[],int T[],int len){
	if(len==1){T[0]=power(A[0],P-2);return;}
	inverse(A,wi1,len>>1);
	memset(wi2,0,len<<3);
	memset(wi3,0,len<<3);
	memcpy(wi2,wi1,len<<1);
	memcpy(wi3,A,len<<2);
	NTT(wi2,len<<1,1);
	NTT(wi3,len<<1,1);
	for(int i=0;i<(len<<1);i++)wi4[i]=((2*wi2[i]-(LL)wi2[i]*wi2[i]%P*wi3[i])%P+P)%P;
	NTT(wi4,len<<1,-1);
	memcpy(T,wi4,len<<2);
}

int wm1[N],wm2[N],wm3[N],wm4[N];

void modulo(int A[],int B[],int R[],int len){
	int degA,degB,degC;
	for(degA=len;degA&&A[degA]==0;--degA);
	for(degB=len;degB&&B[degB]==0;--degB);
	if(degA<degB){for(int i=0;i<len;i++)R[i]=A[i];return;}
	degC=degA-degB;
	memset(wm1,0,len<<2);
	memset(wm2,0,len<<2);
	memset(wm3,0,len<<2);
	memset(wm4,0,len<<2);
	for(int i=0;i<=degC;i++)wm1[i]=A[degA-i];
	for(int i=0;i<=degB;i++)wm2[i]=B[degB-i];
	inverse(wm2,wm3,len>>1);
	for(int i=degC+1;i<len;i++)wm3[i]=0;
	convol(wm1,wm3,wm4,len);
	memset(wm1,0,len<<2);
	for(int i=0;i<=degC;i++)wm1[degC-i]=wm4[i];
	memset(wm2,0,len<<2);
	memcpy(wm2,B,degB+1<<2);
	convol(wm1,wm2,wm3,len);
	for(int i=0;i<=degA;i++)R[i]=(A[i]-wm3[i]+P)%P;
}

int wl1[N],wl2[N];

void sqr(int A[],int D[],int len){
	memset(wl1,0,len<<2);
	memset(wl2,0,len<<2);
	memcpy(wl1,A,len<<1);
	NTT(wl1,len,1);
	for(int i=0;i<len;i++)wl1[i]=(LL)wl1[i]*wl1[i]%P;
	NTT(wl1,len,-1);
	modulo(wl1,D,wl2,len);
	memcpy(A,wl2,len<<1);
}

void work(int w[],int Mod[],int len,LL t){
	int l=0,i;
	while((1LL<<l)<=t)++l;
	while((t>>l)<=L)--l;
	w[t>>l+1]=1;
	while(l>=0){
		sqr(w,Mod,len);
		if((t>>l)&1){
			for(int i=L;i>=1;i--)w[i]=w[i-1];
			w[0]=0;
			for(int i=0;i<L;i++)w[i]=((w[i]-(LL)w[L]*Mod[i])%P+P)%P;
			w[L]=0;
		}
		--l;
	}
}

int n,a[M],A,B,w[N],mo[N],res[N];
LL m;

int main(){
	cin>>n>>m>>a[1]>>A>>B;
	for(int i=2;i<=n;i++)a[i]=((LL)a[i-1]*A+B)%L+1;
	for(int i=1;i<=n;i++)if(a[i]<=23333)w[a[i]]++;
	for(int i=mo[L]=1;i<=L;i++)mo[L-i]=(P-w[i])%P;
	work(res,mo,1<<16,m+L-1);
	printf("%d\n",res[L-1]);
	return 0;
}


评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值