模板 lucas

void extend_gcd(ll a,ll &x,ll b,ll &y){
    if(b==0){
        x=1,y=0;
        return;
    }
    ll x1,y1;
    extend_gcd(b,x1,a%b,y1);
    x=y1;
    y=x1-(a/b)*y1;
}
ll inv(ll a,ll m){
    ll t1,t2;
    extend_gcd(a,t1,m,t2);
    return ( t1%m + m ) % m;
}
ll qpow(ll x,ll y,ll m){
    if(!y) return 1;
    ll ans=qpow(x,y>>1,m);
    ans= ans *ans %m;
    if(y&1) ans=ans *x %m;
    return ans;
}
ll nump(ll x,ll p){
    ll ans=0;
    while(x) ans += x/p, x /= p;
    return ans;
}
ll fac(ll n,ll p,ll pk){
    if(n==0) return 1;
    ll ans=1;
    for(ll i=1;i<=pk;i++){
        if(i%p==0) continue;
        ans= ans *i %pk;
    }
    ans = qpow(ans,n/pk,pk);
    ll to=n%pk;
    for(ll i=1;i<=to;i++){
        if(i%p==0) continue;
        ans = ans *i %pk;
    }
    return fac(n/p,p,pk) *ans %pk;
}
ll cal(ll n,ll m,ll p,ll pi,ll pk){
    ll a=fac(n,pi,pk),b=fac(m,pi,pk),c=fac(n-m,pi,pk);
    ll d=nump(n,pi)-nump(m,pi)-nump(n-m,pi);
    ll ans= a %pk * inv(b,pk) %pk * inv(c,pk) %pk * qpow(pi,d,pk) %pk;
    return ans * (p/pk) %p * inv(p/pk,pk) %p;
}
ll nCmmodp(ll n,ll m,ll p){
    if(n<m) return 0;
    ll ans=0;
    ll x=p;
    for(ll i=2;i*i<=x&&x>1;i++){
        ll k=0,pk=1;
        while(x%i==0){
            x /= i;
            k++;
            pk *= i;
        }
        if(k>0){
            ans =( ans + cal(n,m,p,i,pk) )%p;
        }
    }
    if(x>1) ans=( ans + cal(n,m,p,x,x) )%p;
    return ans;
}

 

posted @ 2015-09-17 15:52  qscqesze  阅读(264)  评论(0编辑  收藏  举报