记录编号 337443 评测结果 AAAAAAAAAA
题目名称 [HZOI 2016]定约servant 最终得分 100
用户昵称 Gravatarkito 是否通过 通过
代码语言 C++ 运行时间 5.587 s
提交时间 2016-11-04 15:27:16 内存使用 0.31 MiB
显示代码纯文本
#include<cstdio>
#include<cmath>
#include<vector>
#include<ctime>
using namespace std;
#define fcl	fclose(stdin);	fclose(stdout);	return 0
	#define SUBMIT 2333
int m,p,phi,Siz,n;
vector<int> Mu[15],Inv[15];
int prime[15],pp[15],tot=0,mi[15];
int AA[15];

/***********************Matrix****************************/
struct Matrix{
	int a[4][4];
	void Init(int k){
		for(int i=1;i<3;++i){
			for(int j=1;j<3;++j){
				a[i][j]=0;
			}
		}
		for(int i=1;i<3;++i)	a[i][i]=k;
	}
}X,Q,F;

Matrix operator * (const Matrix& A,const Matrix& B){
	X.Init(0);
	for(int i=1;i<3;++i){
		for(int j=1;j<3;++j){
			for(int k=1;k<3;++k){
				X.a[i][j]=(X.a[i][j]+(A.a[i][k]*1ll*B.a[k][j])%phi)%phi;
			}
		}
	}
	return X;
}

Matrix QMatrixPow(Matrix a,int b){
	Matrix ans;
	ans.Init(1);
	while(b){
		if(b&1){
			ans=ans*a;
		}
		b>>=1;
		a=a*a;
	}
	return ans;
}

void GetMatrix(){
	Q.a[1][1]=0;	Q.a[1][2]=1;
	Q.a[2][1]=1;	Q.a[2][2]=1;
	Matrix P=QMatrixPow(Q,n-1);
	F.a[1][1]=(P.a[1][1]+P.a[2][1])%phi;
	F.a[1][2]=(P.a[1][2]+P.a[2][2])%phi;
	Siz=(F.a[1][1]*1ll*F.a[1][2])%phi;
}
/***********************Matrix****************************/

int Qpow(int a,int b,int mod){
	int c=1;
	while(b){
		if(b&1){
			c=(c*1ll*a)%mod;
		}
		b>>=1;
		a=(a*1ll*a)%mod;
	}
	return c;
}

void ExGcd(int a,int b,int& x,int& y){
	if(b==0){
		x=1;	y=0;
		return;
	}
	ExGcd(b,a%b,x,y);
	int t=x;
	x=y;	y=t-(a/b)*y;
}

int Gcd(int a,int b){
	return b==0?a:Gcd(b,a%b);
}

bool Divprime(){
	bool flag=false;
	phi=p;
	int x=(int)sqrt(double(p)),mod=p;
	for(int i=2;i<=x;++i){
		if(mod%i==0){
			phi/=i;
			phi*=(i-1);
			prime[++tot]=i;
			pp[tot]=1;
			while(mod%i==0){
				mi[tot]++;
				mod/=i;
				pp[tot]*=i;
			}
			if(mi[tot]>1)	flag=true;
		}
		if(mod==1)	break;
	}
	if(mod!=1){
		phi/=mod;
		phi*=(mod-1);
		prime[++tot]=mod;
		pp[tot]=mod;
		mi[tot]=1;
	}
	return flag;
}

void Init(){
	int x;
	for(int j=1;j<=tot;++j){
		if(n>p)	x=prime[j]-1;
		else x=n;
		Mu[j].push_back(1);	Inv[j].push_back(1);
		for(int i=1;i<=x;++i){
			Mu[j].push_back((Mu[j][i-1]*1ll*i)%p);
			Inv[j].push_back(0);
		}
		int y;
		ExGcd(Mu[j][x],prime[j],Inv[j][x],y);
		Inv[j][x]%=p;	Inv[j][x]=(Inv[j][x]+p)%p;
		for(int i=x-1;i>=0;--i){
			Inv[j][i]=(Inv[j][i+1]*1ll*(i+1))%p;
		}
	}
}

int C(int a,int b,int i){
	if(a<b)	return 0;
	int x=(Mu[i][a]*1ll*Inv[i][b])%prime[i];
	x=(x*1ll*Inv[i][a-b])%prime[i];
	return x;
}

int Lucas(int a,int b,int i){
	if(b==0)	return 1;
	return (Lucas(a/prime[i],b/prime[i],i)*1ll*C(a%prime[i],b%prime[i],i))%prime[i];
}

void Union(int j,int i){
	int a=pp[j],b=pp[i],c=AA[i]-AA[j],x,y;
	ExGcd(a,b,x,y);
	x=(x*1ll*c)%b;
	x=(x+b)%b;
	pp[i]*=pp[j];
	AA[i]=((x*1ll*pp[j])%pp[i]+AA[j])%pp[i];
}

int CRT(){
	for(int i=2;i<=tot;++i){
		Union(i-1,i);
	}
	return AA[tot];
}

void work1(){
	GetMatrix();
	//处理模数拆分后每个质因子的幂为1的情况
	Init();
	for(int i=1;i<=tot;++i){
		for(int j=1;j<=m+1;++j){
			if(Gcd(j,m)==1){
				AA[i]=(AA[i]+Lucas(n,j-1,i))%prime[i];
			}
		}
	}
	int x=CRT();
	printf("%d\n",Qpow(x,Siz,p));
}

int INV(int a,int mod){
	int x,y;
	ExGcd(a,mod,x,y);
	x%=mod;	x=(x+mod)%mod;
	return x;
}

int size;
vector<int> STD;
int Fact(int N,int i){
	if(N==0)	return 1;
	int x=Qpow(STD[pp[i]],N/pp[i],pp[i]);
	x=(x*1ll*STD[N%pp[i]])%pp[i];
	size+=N/prime[i];
	return (x*1ll*Fact(N/prime[i],i))%pp[i];
}

void work2(){
	GetMatrix();
	//处理p的素因子指数大于1的情况 
	int x,y;
	for(int i=1;i<=tot;++i){
		STD.clear();
		STD.push_back(1);
		for(int j=1;j<=pp[i];++j){
			if(j%prime[i]){
				STD.push_back((STD[j-1]*1ll*j)%pp[i]);
			}
			else	STD.push_back(STD[j-1]);
		}
		for(int j=1;j<=m+1;++j){
			if(Gcd(j,m)==1){
				x=1;	y=0;
				size=0;	x=(x*1ll*Fact(n,i))%pp[i];	y+=size;
				size=0;	x=(x*1ll*INV(Fact(j-1,i),pp[i]))%pp[i];	y-=size;
				size=0;	x=(x*1ll*INV(Fact(n-j+1,i),pp[i]))%pp[i];	y-=size;
				x=(x*1ll*Qpow(prime[i],y,pp[i]))%pp[i];
				AA[i]=(AA[i]+x)%pp[i];
			}
		}
	}
	x=CRT();
	printf("%d",Qpow(x,Siz,p));
}

int main(){
	#ifdef SUBMIT
	freopen("servant.in","r",stdin);
	freopen("servant.out","w",stdout);
	#endif
	
	scanf("%d%d%d",&n,&m,&p);
	if(!Divprime())	work1();
	else	work2();
	#ifndef SUBMIT
	getchar();	getchar();
	#endif
	fcl;
}