bzoj 2445 最大团（阶乘取模+中国剩余定理CRT）

最新推荐文章于 2021-09-13 17:14:16 发布

lethalboy

最新推荐文章于 2021-09-13 17:14:16 发布

阅读量902

点赞数

分类专栏：数论文章标签：数论

本文链接：https://blog.csdn.net/qq_34564984/article/details/70224817

版权

数论专栏收录该内容

12 篇文章 0 订阅

订阅专栏

题意即求下式：

m∑x|nn!(x!)nx(nx!)m∑x|nn!(x!)nx(nx!)m∑x|nn!(x!)nx(nx!)

m∑x|nn!(x!)nx(nx!)

根据欧拉定理： 对于互质的正整数 a 和 n ，有 a^φ(n) ≡ 1 mod n 可知，求指数式子模1e9-402的值然后快速幂即可。

O(sqrt(n))分解质因数求值，注意到n非常大，无法预处理阶乘及逆元。

定理： (p-1)!%p=-1

对于 n!，将p提取出并记下p的指数k，fac(n)=(-1)^(n/p) * (n%p)! * fac(n/p)

做除法时求逆元并判断分子和分母中p的指数能否相消，不能为0，能则为fac(n)*inv( fac(x)^(n/x) )*fac(n/x)

(注：fac(n)表示n的阶乘，inv(n)表示n的逆元)

代码：

#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<queue>
#include<stack>
#include<map>
#include<set>
#define ll long long
#define mod 999999599
#define Mod 999999598
#define N 200200
using namespace std;
/*
2 1
13 1
5281 1
7283 1
*/
int n,m,cas;
int MOD[5];
ll ret[5],ans;
ll pre[N][5],inv[N][5];
int num[N][5];
ll qm(ll a,ll b,ll md)
{
	ll ret=1,t=a%md;
	while(b)
	{
		if(b&1)ret=ret*t%md;
		t=t*t%md;b>>=1;
	}return ret;
}
void ex_gcd(int a,int b,int &x,int &y)
{
	if(b==0){x=1,y=0;return ;}
	ex_gcd(b,a%b,x,y);
	int t=x;x=y,y=t-a/b*y;
}
void init()
{
	MOD[1]=2,MOD[2]=13;
	MOD[3]=5281,MOD[4]=7283;
	pre[0][1]=pre[0][2]=pre[0][3]=pre[0][4]=1;
	inv[0][1]=inv[0][2]=inv[0][3]=inv[0][4]=1;
	for(int i=1;i<=4;i++)
		for(int j=1;j<MOD[i];j++)
		{
			int t=j,cnt=0;
			while(t%MOD[i]==0)t/=MOD[i],cnt++;
			pre[j][i]=pre[j-1][i]*t%MOD[i];
			num[j][i]=num[j-1][i]+cnt;
			inv[j][i]=qm(pre[j][i],MOD[i]-2,MOD[i]);
		}
}
ll fac(int x,int p)
{
	if(x<MOD[p])return pre[x][p];
	ll t=x/MOD[p],ret=1;
	ret=((t&1)?-1:1);
	return ret*pre[x%MOD[p]][p]%MOD[p]*fac(x/MOD[p],p)%MOD[p];
}
int con(int x,int p)
{
	if(x<MOD[p])return 0;
	int ret=x/MOD[p];
	return con(x/MOD[p],p)+ret;
}
void cal(int x)
{
	//ret[1]=1;
	for(int i=1;i<=4;i++)
	{
		//if(n>=MOD[i])continue;
		ll tmp,txp,now;
		tmp=fac(x,i),now=con(x,i);
		tmp=qm(tmp,n/x,MOD[i]),now=now*(n/x);
		tmp=qm(tmp,MOD[i]-2,MOD[i]),now=-now;
		txp=fac(n/x,i),txp=qm(txp,MOD[i]-2,MOD[i]);
		tmp=tmp*txp%MOD[i],now-=con(n/x,i);
		tmp=tmp*fac(n,i)%MOD[i],now+=con(n,i);
		if(now!=0)tmp=0;
		ret[i]=(ret[i]+tmp)%MOD[i];
	}
}
void CRT()
{
	for(int i=1;i<=4;i++)
	{
		ll M=Mod/MOD[i];
		//int Mo,y;ex_gcd(M,MOD[i],Mo,y);
		ll Mo=qm(M,MOD[i]-2,MOD[i]);
		ans=(ans+M*Mo%Mod*ret[i]%Mod)%Mod;
	}ans=ans+Mod;
	ans=qm(m,ans,mod);
}
void solve()
{
	int x=sqrt(n);
	for(int i=1;i<=x;i++)
	{
		if(n%i==0)
		{
			int y=n/i;cal(i);
			if(y!=i)cal(y);
		}
	}CRT();
	printf("%lld\n",ans);
}
int main()
{
	scanf("%d",&cas);init();
	while(cas--)
	{
		memset(ret,0,sizeof(ret));
		scanf("%d%d",&n,&m);
		ans=0;solve();
	}
}