題目鏈接
->擴展Lucas

//求C_n^k%m 
#include <cstdio>
typedef long long LL;

LL FP(LL x,LL k,LL p)
{
    LL t=1ll;
    for(; k; k>>=1,x=x*x%p)
        if(k&1) t=t*x%p;
    return t;
}
void Exgcd(LL a,LL b,LL &x,LL &y)
{
    if(!b) x=1ll, y=0ll;
    else Exgcd(b,a%b,y,x),y-=a/b*x;
}
LL Inv(LL a,LL mod)
{
//  if(!a) return 0ll;//?
 

    LL x,y; Exgcd(a,mod,x,y);
    x=(x%mod+mod)%mod;//!
//  if(!x) x=mod;
    return x;
}
LL Fact(LL n,LL pi,LL pk)//factorial Calc n!%(pi^ki) (不計算pi因子 計算C()時提出)
{
    if(!n) return 1ll;
    LL ans=1ll;
    if(n/pk)
    {
        for(LL i=2; i<=pk; ++i)//每pi^ki一循環的部分 
            if(i%pi) (ans*=i)%=pk;
        ans=FP(ans,n/pk,pk);//一共n/pk個循環 
 

    }
    for(LL i=2; i<=n%pk; ++i)//pi^ki循環之外的部分 mod pk意義下所以i=2 to n%pk即可 
        if(i%pi) (ans*=i)%=pk;
    return ans*Fact(n/pi,pi,pk)%pk;//[n/pi]!部分 
}
LL C(LL n,LL m,LL mod,LL pi,LL pk)//Calc C_n^m%(pi^ki)
{
    if(n<m) return 0ll;
    LL a=Fact(n,pi,pk),b=Fact(m,pi,pk),c=Fact(n-m,pi,pk),k=0ll;//k:質因子pi的個數 
 

    for(LL i=n; i; i/=pi) k+=i/pi;//計算x!中pi因子個數:k=f(x)=f(x/pi)+x/pi 
    for(LL i=m; i; i/=pi) k-=i/pi;
    for(LL i=n-m; i; i/=pi) k-=i/pi;
    LL ans=a*Inv(b,pk)%pk*Inv(c,pk)%pk*FP(pi,k,pk)%pk;
    return ans*(mod/pk)%mod*Inv(mod/pk,pk)%mod;//這步不是很懂 大概是轉到模mod意義下？
}

int main()
{
    LL n,k,mod,ans=0ll;
    scanf("%I64d%I64d%I64d",&n,&k,&mod);
    for(LL now=mod,i=2; i<=mod; ++i)
        if(!(now%i))
        {
            LL pk=1ll;
            while(!(now%i)) pk*=i, now/=i;
            (ans+=C(n,k,mod,i,pk))%=mod;
        }
    printf("%I64d",ans);

    return 0;
}

CF.100633J.Ceizenpok's formula(擴展Lucas)