组合数模任意数

本文内容是求组合数模合数的方法

本方法有局限性，假设求CKN mod P

其中P=Πmi=1Pqii,对于任意i,j∈[1,m],(Pi,Pj)=1

若要使用本方法，则Pqii不能太大

具体方法如下

我们知道，一般处理模合数的情况，我们要用到CRT（中国剩余定理），对于互质的模数来分别求答案，然后用CRT合并。

用M表示模的数，CRT的公式：ans=aitiMi

其中ai表示对于每一个模数mi得到的答案，Mi表示Mmi，ti表示Mi在模mi意义下的逆元

首先考虑组合数的通项公式

Ckn=n!k!(n−k)!

所以，我们就有一种思路，求n! mod mik!(n−k)! mod mi。在求的时候不乘上mi这个因子，然后再在外面得到mi的个数（详见代码）就可以算出来Ckn mod mqii了

#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;
}