卢卡斯定理

给定$n,m,p(1\le n,m,p\le 10^5)$

求 $C_{n+m}^{m}\ mod\ p$

保证P为prime

C表示组合数。

#include<iostream>
#include<algorithm>
#include<cstdio>
#define ed 100005
using namespace std;
int k,n,m,p;
long long a[ed],b[ed];
long long lucas(int x,int y)
{
    if(x<y) return 0;
    else if(x<p) return b[x]*a[y]*a[x-y]%p;
    else return lucas(x/p,y/p)*lucas(x%p,y%p)%p;
}
int main()
{
    scanf("%d",&k);
    while(k)
    {
        scanf("%d%d%d",&n,&m,&p);
        a[0]=a[1]=b[0]=b[1]=1;
        for(int i=2;i<=n+m;i++) b[i]=b[i-1]*i%p;
        for(int i=2;i<=n+m;i++) a[i]=(p-p/i)*a[p%i]%p;
        for(int i=2;i<=n+m;i++) a[i]=a[i-1]*a[i]%p;
        printf("%lld\n",lucas(n+m,m));
        k--;
    }
    return 0;
}

扩展卢卡斯

求$C_n^m \bmod{p}$

1≤m≤n≤10^{18}, 2≤p≤1000000，不保证 pp 是质数。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b){x=1;y=0;return a;}
    ll res=exgcd(b,a%b,x,y),t;
    t=x;x=y;y=t-a/b*y;
    return res;
}
ll p;
inline ll power(ll a,ll b,ll mod)
{
    ll sm;
    for(sm=1;b;b>>=1,a=a*a%mod)if(b&1)
        sm=sm*a%mod;
    return sm;
}
ll fac(ll n,ll pi,ll pk)
{
    if(!n)return 1;
    ll res=1;
    for(register ll i=2;i<=pk;++i)
        if(i%pi)(res*=i)%=pk;
    res=power(res,n/pk,pk);
    for(register ll i=2;i<=n%pk;++i)
        if(i%pi)(res*=i)%=pk;
    return res*fac(n/pi,pi,pk)%pk;
}
inline ll inv(ll n,ll mod)
{
    ll x,y;
    exgcd(n,mod,x,y);
    return (x+=mod)>mod?x-mod:x;
}
inline ll CRT(ll b,ll mod){return b*inv(p/mod,mod)%p*(p/mod)%p;}
const int MAXN=11;
static ll n,m;
static ll w[MAXN];
inline ll C(ll n,ll m,ll pi,ll pk)
{
    ll up=fac(n,pi,pk),d1=fac(m,pi,pk),d2=fac(n-m,pi,pk);
    ll k=0;
    for(register ll i=n;i;i/=pi)k+=i/pi;
    for(register ll i=m;i;i/=pi)k-=i/pi;
    for(register ll i=n-m;i;i/=pi)k-=i/pi;
    return up*inv(d1,pk)%pk*inv(d2,pk)%pk*power(pi,k,pk)%pk;
}
inline ll exlucus(ll n,ll m)
{
    ll res=0,tmp=p,pk;
    static int lim=sqrt(p)+5;
    for(register int i=2;i<=lim;++i)if(tmp%i==0)
    {
        pk=1;while(tmp%i==0)pk*=i,tmp/=i;
        (res+=CRT(C(n,m,i,pk),pk))%=p;
    }
    if(tmp>1)(res+=CRT(C(n,m,tmp,tmp),tmp))%=p;
    return res;
}
int main()
{
    scanf("%lld%lld%d",&n,&m,&p);
    printf("%d\n",exlucus(n,m));
    return 0;
}