exLucas

参考博客

exLucas:

Cnmmodd ( d 不一定为质数)

1.

d 质因数分解为 d=p1c1×p2c2××pkck

i,j[1,k] , picipjcj 互质,所以可以构造出如下同余方程:

{a1Cnm (mod p1c1)a2Cnm (mod p2c2)akCnm (mod pkck)

在求出所有的 ai 后,就可以利用中国剩余定理求解出 Cnm

所以问题转化为求出 Cnmmodpα 的值。

2.

但是同余方程 ax1(mod p) (乘法逆元) 有解的充要条件为 gcd(a,p)=1

然而题目显然无法保证有解,所以无法直接求 invm!inv(nm)!

所以可将原式转化为:

n!pxm!py×(nm)!pz×pxyzmodpα

其中 xn! 中包含 p 的个数,y,z 同理。

3.

若对于一个 n 可以求出 n!px ,那么就可以求逆元了,所以现在要求 n!pxmodpk

n! 进行变形可得:

n!=1×2××n=(p×2p×)×i=1,i0(modp)ni=pnp×(np)!×i=1,i0(modp)ni=pnp×(np)!×(i=1,i0(modp)pki)npk×(i=pknpk,i0(modp)ni)

f(n)=n!px (f(0)=1)

f(n)=f(np)×(i=1,i0(modp)pki)npk×(i=pknpk,i0(modp)ni)

=f(n)f(m)f(nm)pxyzmodpk

这样就可以用 exgcd 求逆元了。

最后利用中国剩余定理得到答案。单次时间复杂度 O(plogp)

代码
#include<bits/stdc++.h>
#define ll long long
using namespace std;

inline ll read(){
	ll s=0,k=1;
	char c=getchar();
	while(c>'9'||c<'0'){
		if(c=='-')k=-1;
		c=getchar();
	}
	while(c>='0'&&c<='9'){
		s=(s<<3)+(s<<1)+(c^48);
		c=getchar();
	}
	return s*k;
}

ll n,m,p; 

ll ksm(ll a,ll b,ll mod){
	ll t=1;
	for(;b;b>>=1,a=a*a%mod)
		if(b&1) t=t*a%mod;
	return t;
}

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

ll inv(ll n,ll mod){
	ll x,y;
	exgcd(n,mod,x,y);
	return (x%mod+mod)%mod; 
}

ll CRT(int n,ll a[],ll P[]){
	ll ans=0;
	for(int i=1;i<=n;i++) ans=(ans+a[i]*(p/P[i])%p*inv(p/P[i],P[i])%p)%p;
	return ans;
}

ll f(ll n,ll pi,ll pk){
	if(!n) return 1;
	ll t=1;
	for(ll i=2;i<=pk;i++) if(i%pi) (t*=i)%=pk;
	t=ksm(t,n/pk,pk);
	for(ll i=n/pk*pk+1;i<=n;i++) if(i%pi) (t*=(i%pk))%=pk;
	return t*f(n/pi,pi,pk)%pk;
}

ll C(ll n,ll m,ll pi,ll pk){
	int k=0;
	for(ll i=n;i;i/=pi) k+=i/pi;
	for(ll i=m;i;i/=pi) k-=i/pi;
	for(ll i=n-m;i;i/=pi) k-=i/pi;
	return f(n,pi,pk)*inv(f(m,pi,pk),pk)%pk*inv(f(n-m,pi,pk),pk)%pk*ksm(pi,k,pk)%pk;
}

ll exLucas(ll n,ll m,ll p){
	int cnt=0;
	ll a[20],P[20],tmp=p;
	for(int i=2;i*i<=p;i++)
		if(tmp%i==0){
			P[++cnt]=1;
			while(tmp%i==0) P[cnt]*=i,tmp/=i;
			a[cnt]=C(n,m,i,P[cnt]);
		} 
	if(tmp>1) P[++cnt]=tmp,a[cnt]=C(n,m,tmp,tmp);
	return CRT(cnt,a,P);
}

int main(){
	n=read();m=read();p=read();
	printf("%lld\n",exLucas(n,m,p));
	return 0;
}
posted @   programmingysx  阅读(7)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· NetPad:一个.NET开源、跨平台的C#编辑器
· PowerShell开发游戏 · 打蜜蜂
· 凌晨三点救火实录:Java内存泄漏的七个神坑,你至少踩过三个!
点击右上角即可分享
微信分享提示