大组合数模100W以内任意整数

那么要求模数为质数的lucas定理就不能用了

我们只能直接算阶乘了

但是对组合数分母来说与模数不互质的话逆元是不存在的

记当前问题为C(n,k) % m

把m分解质因数后得到的m=p1^r1 * p2^r2 *......

C(n,k) % p1^r1

C(n,k) % p2^r2

......

之后把所有结果用中国剩余定理合并就得到原来的结果了

那么现在的问题就是怎么求 C(n,k) % p^r

首先还是要解决互质的问题我们才能算分母的逆元

我们计算出 n! k! (n-k)!分别含有多少个p因子，就可以算出该组合数总的p因子数，如果大于等于r说明是p^r的倍数，取模得0，如果不是那还要往下算

算一个阶乘含多少个p因子的方法是：n/p+n/p/p+n/p/p/p+......

然后就是要算剩下的与p互质的部分了，也就是要算一个不含p因子的阶乘

记mod=p^r

我们可以先预处理出到mod的阶乘，遇到p的倍数先暂时忽略，乘1就可以了，其他的照常乘进去

然后对于n!来说，每mod个数的乘积对mod取模是相等的，就是说 1*2*...*mod == (mod+1)*(mod+2)*(2mod)

当然我说的都是除去含p因子的数之后的，为了方便描述我就先直接这么说了

那么就可以看n能被mod分成多少份，就得到n!=(mod!)^(n/mod) * (n%mod)! * （含p因子的数除去p因子后剩下的东西）

现在问题又缩小到怎么求含p因子的数除去p因子后剩下的东西了

因为 p 除去p剩下1 ，2p除去p剩下2 ，。。。。。。

因为剩下的东西要乘起来，所以就相当于能被p分成多少份，就是它的阶乘，即(n/p)!

这个东西也是要除去含p的因子的，也就是那些p的高次方的数除去一个p后还会剩下一些p

那么就跟我们上面的过程类似，可以n/=p循环求解，把所有结果乘一起

现在得到了不含p因子的阶乘了，那么分母就可以通过欧拉定理求逆了

把分子和分母的逆和一开始求出来的该组合数总的p因子数全部乘一起就得到一条方程的结果了

好了之后剩下的就是中国剩余定理合并所有方程的结果了。

代码：

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<string>
#include<iomanip>
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<vector>
#include<set>
#include<map>
#include<queue>
#include<list>

using namespace std;
typedef long long LL;
typedef unsigned long long ULL;

#define rep(i,k,n) for(int i=(k);i<=(n);i++)
#define red(i,k,n) for(int i=(k);i>=(n);i--)
#define sqr(x) ((x)*(x))
#define clr(x,y) memset((x),(y),sizeof(x))
#define MAX(a,b) ((a)>(b)?(a):(b))
#define MIN(a,b) ((a)>(b)?(b):(a))

LL n,k,m;
LL p[25],r[25];
int pcnt;

LL a[25];

LL quipow(LL x,LL k,LL mod)
{
    LL ret=1;
    while(k)
    {
        if(k&1)ret=ret*x%mod;
        x=x*x%mod;
        k>>=1;
    }
    return ret;
}

int quipow(int x,int k)
{
    int ret=1;
    while(k)
    {
        if(k&1)ret=ret*x;
        x=x*x;
        k>>=1;
    }
    return ret;
}

void getprime(int m)
{
    clr(p,0);
    clr(r,0);
    pcnt=0;
    for(LL i=2;i*i<=m;i++)
    {
        if(m%i==0)
        {
            p[++pcnt]=i;
            while(m%i==0){m/=i;r[pcnt]++;}
        }
    }
    if(m>1)
    {
        p[++pcnt]=m;
        r[pcnt]=1;
    }
}

LL f[1000010];

void pref(int p,int mod)
{
    clr(f,0);
    f[0]=f[1]=1;
    rep(i,2,mod-1)f[i]=f[i-1]*(i%p?i:1)%mod;
}

LL getp(LL n,LL p)
{
    LL ret=0;
    while(n) ret+=(n/=p);
    return ret;
}

LL getf(LL n,LL mod,LL p)
{
    LL ret=1;
    while(n)
    {
        ret=ret*quipow(f[mod-1],n/mod,mod) %mod *f[n%mod] %mod ;
        n/=p;
    }
    return ret;
}

LL getc(LL n,LL m,LL p,LL r)
{
    LL mod=quipow(p,r);
    pref(p,mod);
    LL nump=getp(n,p)-getp(m,p)-getp(n-m,p);
    if(nump>=r)return 0;
    return quipow(p,nump,mod)*getf(n,mod,p)%mod *quipow(getf(m,mod,p)*getf(n-m,mod,p)%mod,(p-1)*quipow(p,r-1)-1,mod)%mod ;
}

LL exgcd(LL a,LL b,LL &x,LL &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;
    }
    else
    {
        LL g=exgcd(b,a%b,y,x);
        y-=x*(a/b);
        return g;
    }
}

int CR(int cnt)
{
    LL ret=0,x,y;
    rep(i,1,cnt)
    {
        exgcd(m/p[i],p[i],x,y);
        ret=(ret+m/p[i]*x*a[i])%m;
    }
    return (ret+m)%m;
}

int main()
{
//#define LOCAL
#ifdef LOCAL
	freopen("e:\\read.txt","r",stdin);
	//freopen("e:\\write.txt","w",stdout);
#endif
    while(cin>>n>>k>>m)
    {
        getprime(m);
        rep(i,1,pcnt)
        {
            a[i]=getc(n,k,p[i],r[i]);
            p[i]=quipow(p[i],r[i]);
        }
        cout<<CR(pcnt)<<endl;
    }


	return 0;
}