[模板]扩展卢卡斯定理


#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#define WR WinterRain
using namespace std;
const long long WR=10010;
long long n0,m0,w[WR],mod;
long long a[WR],b[WR];
long long res=1;
long long read(){
    long long s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch<='9'&&ch>='0'){
        s=(s<<3)+(s<<1)+ch-48;
        ch=getchar();
    }
    return s*w;
}
long long quick_pow(long long a,long long b,long long md){
    long long base=a,ans=1;
    while(b>0){
        if(b&1) ans=ans*base%md;
        base=base*base%md;
        b>>=1;
    }
    return ans%md;
}
long long ex_gcd(long long &x,long long &y,long long a,long long b){
    if(b==0){
        x=1,y=0;
        return a;
    }
    long long d=ex_gcd(x,y,b,a%b);
    long long tmp=x;
    x=y,y=tmp-a/b*y;
    return d;
}
long long CRT(long long n,long long *a,long long *b){
    long long md=1,ans=0;
    for(long long i=1;i<=n;i++) md=md*b[i];
    for(long long i=1;i<=n;i++){
        long long x,y;
        ex_gcd(x,y,md/b[i],b[i]);
        ans=(ans+md/b[i]*a[i]*x)%md;
    }
    return (ans+md)%md;
}
//=============EX_LUCAS板子===============
long long calc(long long n,long long p,long long md){
    //printf("calculating %lld %lld %lld/n",n,p,md);
    if(!n) return 1;
    long long ans=1;
    for(long long i=1;i<=md;i++)
        if(i%p!=0) ans=ans*i%md;
    ans=quick_pow(ans,n/md,md);
    for(long long i=1;i<=n%md;i++)
        if(i%p!=0) ans=ans*i%md;
    ans=ans*calc(n/p,p,md)%md;
    return ans;
}
long long ex_Lucas(long long n,long long m,long long p,long long md){
    //printf("EX_Lucasing %lld %lld %lld %lld/n",n,m,p,md);
    long long cnt=0;
    for(long long i=n;i;i/=p) cnt+=i/p;
    for(long long i=m;i;i/=p) cnt-=i/p;
    for(long long i=n-m;i;i/=p) cnt-=i/p;
    return quick_pow(p,cnt,md)*calc(n,p,md)%md
           *quick_pow(calc(m,p,md),md/p*(p-1)-1,md)%md
           *quick_pow(calc(n-m,p,md),md/p*(p-1)-1,md)%md;
}
long long comb(long long n,long long m,long long md){
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    //printf("combing %lld %lld %lld/n",n,m,md);
    long long tot=0;
    for(long long i=2;i*i<=md;i++){
        if(md%i==0){
            b[++tot]=1;
            while(md%i==0) md/=i,b[tot]*=i;
            a[tot]=ex_Lucas(n,m,i,b[tot]);
            //printf("%lld calculate completed/n",i);
        }
    }
    if(md>1) b[++tot]=md,a[tot]=ex_Lucas(n,m,md,b[tot]);
    //for(long long i=1;i<=tot;i++){
    //    printf("%lld %lld/n",a[i],b[i]);
    //}
    return CRT(tot,a,b);
}
int main(){
    n0=read(),m0=read(),mod=read();
    printf("%lld",comb(n0,m0,mod));
    return 0;
}

 

原创文章,作者:dweifng,如若转载,请注明出处:https://blog.ytso.com/276418.html

(0)
上一篇 2022年7月23日
下一篇 2022年7月23日

相关推荐

发表回复

登录后才能评论