组合数模任意数

本文内容是求组合数模合数的方法
本方法有局限性，假设求$C_N^K \ mod \ P$

$$其中P = \Pi_{i=1}^{m}P_i^{q_i},$$ $$对于任意i,j\in[1,m],(P_i, P_j) = 1$$
若要使用本方法，则 $P_i^{q_i}$不能太大
$$具体方法如下$$
我们知道，一般处理模合数的情况，我们要用到CRT（中国剩余定理），对于互质的模数来分别求答案，然后用CRT合并。
用M表示模的数，CRT的公式： $ans = a_it_iM_i$
其中 $a_i表示对于每一个模数m_i得到的答案$， $M_i表示\frac{M}{m_i}$， $t_i表示M_i在模m_i意义下的逆元$
首先考虑组合数的通项公式
$$C_n^k = \frac{n!}{k!(n-k)!}$$
所以，我们就有一种思路，求 $\frac{n! \ mod \ m_i}{k! (n-k)! \ mod \ m_i}$。在求的时候不乘上 $m_i$这个因子，然后再在外面得到 $m_i$的个数（详见代码）就可以算出来 $C_n^k \ mod \ m_i^{q_i}$了

#include <iostream>
#include <cstdio>
#define LL long long
using namespace std;

LL ksm(LL a, LL k, LL MOD)
{
    LL re = 1;
    for(; k; k >>= 1, a = a*a % MOD)
        if(k & 1) re = re*a % MOD;
    return re;
}

LL factor(LL N, LL P, LL K) /**< To Get N!%(P^K)*/
{
    LL Pk, t1, t2, t3, t4;
    Pk = ksm(P, K, 0x7fffffff);/**< Because P^K <= 100000, we can get the exact value of P^K */
    if(N < P)
    {
        for(t1 = t2 = 1; t1 <= N; ++ t1)
            t2 = t2 * t1 % P;
    }
    else
    {
        t2 = t3 = 1; t4 = N % Pk;
        for(t1 = 1; t1 <= t4; ++ t1)
            if(t1 % P != 0) t2 = t2*t1%Pk;
        for(t1 = 1; t1 < Pk;  ++ t1)
            if(t1 % P != 0) t3 = t3*t1%Pk;
        t3 = ksm(t3, N/Pk, Pk);
        t2 = t2*t3%Pk;
        t3 = factor(N/P, P, K);
        t2 = t2*t3%Pk;
    }
    return t2;
}

LL a[30], p[30];
LL C(LL N, LL K, LL P) /**< C(N, K) % P */
{
    int i, c = 0;
    LL t1, t2, t3, pa, ans;
    t1 = P;
    for(i = 2; i * i <= t1; ++ i) if(t1 % i == 0)
        for(p[++ c] = i; t1%i == 0; t1/=i, ++ a[c]);
    if(t1 > 1) p[++ c] = t1, a[c] = 1;
    ans = 0;
    for(i = 1; i <= c; ++ i)
    {
        pa = ksm(p[i], a[i], 0x7fffffff);
        t1 = factor(N, p[i], a[i]);
        t2 = factor(K, p[i], a[i]);
        t3 = factor(N-K,p[i],a[i]);
        t2 = t2*t3%pa;
        t2 = ksm(t2, pa/p[i]*(p[i]-1)-1, pa); /**< Eular's Theorem */
        t1 = t1*t2%pa;
        t3 = 0;
        for(t2 = p[i]; t2 <= N; t2*=p[i]) t3 += N/t2;
        for(t2 = p[i]; t2 <= K; t2*=p[i]) t3 -= K/t2;
        for(t2 = p[i]; t2 <=N-K;t2*=p[i]) t3 -= (N-K)/t2;
        t2 = ksm(p[i], t3, pa);
        t1 = t1 * t2 % pa;
        t2 = ksm(P/pa%pa, pa/p[i]*(p[i]-1)-1, pa);
        t2 = P/pa * t2 % P;
        t1 = t1 * t2 % P;
        ans = (ans + t1) % P;
    }
    return ans;
}

int main() {
    LL N, K, P;
    scanf("%I64d%I64d%I64d", &N, &K, &P);
    printf("%I64d\n", C(N, K, P));
    return 0;
}