Lucas定理及扩展

Lucas定理

不会证明。。。

ppp为质数
C(n,m)≡C(n/p,m/p)∗C(n%p,m%p)(mod p)C(n, m)\equiv C(n/p, m/p)*C(n\%p, m\%p)(mod\ p)C(n,m)C(n/p,m/p)C(n%p,m%p)(mod p)

扩展

C(n,m)C(n,m)C(n,m)MMM 意义下的值
M=∏piaiM=\prod p_i^{a_i}M=piai
那么就只要求出模 piaip_i^{a_i}piai 的值,然后 CRTCRTCRT 合并即可
考虑求 C(n,m)%piaiC(n, m) \% p_i^{a_i}C(n,m)%piai
C(n,m)=n!m!(n−m)!C(n,m)=\frac{n!}{m!(n-m)!}C(n,m)=m!(nm)!n!

  1. 首先可以把分子分母中 pip_ipi 的因子约分
    n!n!n!pip_ipi 的个数为
    ∑k=1ai⌊npik⌋\sum_{k=1}^{a_i}\lfloor\frac{n}{p_i^k}\rfloork=1aipikn
  2. 提出 pip_ipi 后,就只要求出 n!%piain!\% p_i^{a_i}n!%piai 就好了,逆元也可以直接 exgcdexgcdexgcd
    先把 n!n!n! 中含有 pip_ipi 这个因子的项单独拿出,那么
    n!=1×2×...×(pi−1)×(pi+1)×...×...(pi2+1)×...×pik(1×2×3×...)n!=1\times 2 \times ... \times (p_i-1) \times (p_i+1) \times ... \times ... (p_i^2+1) \times... \times p_i^k(1\times 2 \times 3 \times ...)n!=1×2×...×(pi1)×(pi+1)×...×...(pi2+1)×...×pik(1×2×3×...)
    对于 pikp_i^kpik 之前提出来算过了,所以递归处理后面的就好了
    考虑前面的求法,1×2×...×(pi−1)×(pi+1)×...×...(pi2+1)×...1\times 2 \times ... \times (p_i-1) \times (p_i+1) \times ... \times ... (p_i^2+1) \times...1×2×...×(pi1)×(pi+1)×...×...(pi2+1)×...
    由于是模 piaip_i^{a_i}piai 意义下的,所以这些东西被分成若干段相乘,每段值一样,直接预处理即可
# include <bits/stdc++.h>
using namespace std;
typedef long long ll;

namespace IO {
	const int maxn((1 << 21) + 1);

	char ibuf[maxn], *iS, *iT, c;
	int f;
	
	char Getc() {
		return (iS == iT ? (iT = (iS = ibuf) + fread(ibuf, 1, maxn, stdin), (iS == iT ? EOF : *iS++)) : *iS++);
	}
	
	template <class Int> void In(Int &x) {
		for (f = 1, c = Getc(); c < '0' || c > '9'; c = Getc()) f = c == '-' ? -1 : 1;
		for (x = 0; c <= '9' && c >= '0'; c = Getc()) x = (x << 3) + (x << 1) + (c ^ 48);
		x *= f;
	}
}

using IO :: In;

const int maxn(1e6 + 5);

int mod, p[20], a[20], x[20], b[20], num, fac[maxn];

inline int Pow(ll x, ll y, int m) {
	ll ret = 1;
	for (; y; y >>= 1, x = x * x % m)
		if (y & 1) ret = ret * x % m;
	return ret;
}

inline void ExGcd(int a, int b, int c, int &xx, int &yy, int m) {
	if (!b) {
		xx = (c / a + m) % m, yy = 0;
		return;
	}
	ExGcd(b, a % b, c, yy, xx, m);
	yy = (yy - 1LL * (a / b) * xx % m + m) % m;
}

inline ll Gcd(ll a, ll b) {
	return !b ? a : Gcd(b, a % b);
}

inline ll F(ll xx, int yy) {
	return xx < yy ? 0 : xx / yy + F(xx / yy, yy);
}

int ans, cur, xx, yy;

inline int Inv(int a, int m) {
	return ExGcd(a, m, Gcd(a, m), xx, yy, m), xx;
}

inline int Fac(ll n, int pi, int xi) {
	return n <= pi ? fac[n] : 1LL * Pow(fac[xi], n / xi, xi) * fac[n % xi] % xi * Fac(n / pi, pi, xi) % xi;
}

ll n, m;

inline ll Solve(int pi, int ai, int xi) {
	ll nn = F(n, pi) - F(m, pi) - F(n - m, pi);
	if (nn >= ai) return 0;
	nn = Pow(pi, nn, xi);
	int facn = Fac(n, pi, xi), im = Inv(Fac(m, pi, xi), xi), inm = Inv(Fac(n - m, pi, xi), xi);
	return 1LL * facn * im % xi * inm % xi * nn % xi;
}

int main() {
	In(n), In(m), In(mod), cur = mod, fac[0] = 1;
	for (int i = 2; i * i <= cur; ++i)
		if (cur % i == 0) {
			p[++num] = i, x[num] = 1;
			while (cur % i == 0) cur /= i, ++a[num], x[num] *= i;
		}
	if (cur > 1) p[++num] = cur, ++a[num], x[num] = cur;
	for (int i = 1; i <= num; ++i) {
		for (int j = 1; j <= x[i]; ++j)
			if (j % p[i]) fac[j] = 1LL * fac[j - 1] * j % x[i];
			else fac[j] = fac[j - 1];
		b[i] = Solve(p[i], a[i], x[i]);
	}
	for (int i = 2; i <= num; ++i) {
		int xx, yy, c = b[i] - b[1], lcm = x[1] * x[i];
		ExGcd(x[1], x[i], c, xx, yy, lcm);
		b[1] = (1LL * xx * x[1] % lcm + b[1]) % lcm, x[1] = lcm;
	}
	printf("%d\n", b[1]);
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值