CRT, lucas及其扩展形式

exgcd

int exgcd(int a, int b, int &x, int &y) {
	if (b == 0) return a, x = 1, y = 0; 
	int y = exgcd(b, a % b, x, y), t;
	t = x, x = y, y = t - a / b * y;
}

证明:

gcd的过程中, 假设我们已经求出了(b * x + (a~\%~b) * y = gcd(a, b))推导到(a*x + b*y = gcd(a, b))的形式

[b * x + (a - frac ab*b) y = gcd(a, b)\ a * y + b * (x - frac ab * y) = gcd(a, b) ]

中国剩余定理(CRT)

问题: 求最小整数解x, 其中(P_1,P_2cdots P_n)互素

[left{egin{matrix} x equiv a_1 mod P_1\ x equiv a_2 mod P_2\ vdots \ x equiv a_n mod P_n\ end{matrix} ight. ]

设(M = Pi^n_{i=1} P_i~~~ M_i=frac M{P_i}~~~M_i^{i-1})为(modM_i)意义下的逆元那么

[X = (sum_{i=1}^na_iM_iM_i^{-1}) mod M ]

扩展中国剩余定理(exCRT)

问题同中国剩余定理但各模数P不互质

假设已经求出前(k-1)个方程组成的同余方程组的一个解为(x), 且有 (M=LCM_{i-1}^{k-1}P_i), 则满足前k个方程的解为(x + t * M)

加入第k个方程得(x + t * M = a_i mod P_k)

即 (t*M=a_i-x mod P_i)

(exgcd)求解即可, 同时可以判是否有解

Lucas定理

解决的问题, P是质数且不大, 而n, m可以很大

[{n choose m} mod P ]

证明咕咕咕, 其实没啥用

[{nchoose m} mod P = {nmod P choose mmod P} * { frac np choose frac mp} ]

其实可以理解为转化为p进制下

ll C(int n,int m) {
	if (m > n) return 0;
	return jie[n] * fpw(jie[m], P-2) % P * fpw(jie[n-m], P-2) % P;
}

ll lau(ll x, ll y) {
	if (!y) return 1;
	return C(x % P, y % P) * lau(x / P, y / P) % P;
}

exLucas定理

这个比较有用, 不要求P是质数了

设(P= prod p_i^k), 如果求出了(inom{n}{m}mod{p_i^{k_i}}) 就可以利用CRT合并答案了

[{n choose m}= frac{n!}{(n−m)!m!} ]

把p都提取出来, 可以先计算上面有多少p, 下面有多少p

然后考虑其他的数, 即求(K!mod P), 发现在模(p^k)的时候有循环节, 因为(p^k+i equiv p^k)

所以对于与p互质的数直接暴力扫描一个(p^k)然后快速幂即可

对与不与p互质的数, 去除所有p后, 仍会有其他因子, 不妨递归解决子问题, 即先提取出一个p来, p的倍数为(1p, 2p, 3p cdots frac npp) ,都去掉一个p然后相乘等于(frac np!) , 递归下去就好喽

#include<iostream>
#include<cstring>
#include<cmath>
#include<cstdio>
#define ll long long
using namespace std;
ll m, n, p;

ll exgcd(ll a, ll b, ll &x, ll &y) {
	if (!b) return x = 1, y = 0, a;
	ll gcd = exgcd(b, a % b, x, y), t;
	t = x, x = y, y = t - a / b * y;
	return gcd;
}

ll inline inv(ll n, ll mod) {
	ll x, y; exgcd(n, mod, x, y);
	return (x % mod + mod) % mod;
}

ll fpw(ll x, ll mi, ll mod) {
	ll res = 1;
	while (mi) {
		if (mi & 1) res = res * x % mod;
		mi >>= 1;
		x = x * x % mod;
	}
	return res;
}

ll fac(ll n, ll pi, ll pk) {
	if (!n) return 1;
	ll res = 1;
	for (ll i = 2;i <= pk; i++) 
		if (i % pi) res = res * i % pk;
	res = fpw(res, n / pk, pk);
	for (ll i = 2;i <= n % pk; i++) 
		if (i % pi) res = res * i % pk;
	return res * fac(n / pi, pi, pk) % pk;
}

inline ll ptime(ll x, ll pi) {
	ll res = 0;
	for (;x ; x /= pi) res += x / pi;
	return res;
}

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 = ptime(n, pi) - ptime(n-m, pi) - ptime(m, pi);
	return up * inv(d1, pk) % pk * inv(d2, pk) % pk * fpw(pi, k, pk) % pk;
}

ll CRT(ll b, ll mod) {
	return b * inv(p / mod, mod) % p * (p / mod) % p;
}

inline ll exlucus(ll n, ll m) {
	ll res = 0, tmp = p, pk;
	ll lim = sqrt(p);
	for (int i = 2;i <= lim; i++) {
		if (tmp % i == 0) {
			pk = 1;
			while(tmp % i == 0) pk *= i, tmp /= i;
			res = (res + CRT(C(n, m, i, pk), pk)) % p;
//			cout << tt << endl;
		}
	}
	if (tmp > 1) res = (res + CRT(C(n, m, tmp, tmp), tmp)) % p; 
	return res;
}

int main() {
	cin >> n >> m >> p;
	printf ("%lld
", exlucus(n, m));
	return 0;
}