卢卡斯Lucas & 扩展卢卡斯exLucas

卢卡斯定理

适用条件:

只算数量较少的C(n, m), 模数P较小(约1e5),且P为质数。

递推式:C(n, m) = C(n/P, m/P) * C(n%P, m%P) % P;

当 n < m 的时候, \(C_n^m = 0\)

当 n == 0 或者 m == 0 的时候就可以停止了。

Code:

//阶乘已经预处理好,逆元用快速幂
long long Lucas(long long n, long long m) {//n > m 
	if (m == 0 || n == 0)	return 1;
	long long res = Lucas(n / p, m / p) * C(n % p, m % p) % p;
	return res;
}

扩展卢卡斯

适用范围:

只求少量组合数,m, n巨大,但模数 \(M = \prod p^q\) ,且 \(p^q\) 不大。

复杂度:\(O(\log_Pn * (p^q))\)

exlucas

扩展卢卡斯实际上还可以搞任何阶乘除阶乘的类似问题。比如:P2183 [国家集训队]礼物。主要思路是对每个 \(p_q\) 单独考虑,对因子 \(p\) 和其它因子单独考虑。提出来 \(p\) 的倍数,剩下的利用暴力前缀积搞,\(p\) 的倍数部分递归搞。

Code:

//P2183礼物
int realP;
int n, m, w[10];
int p[44], c[44], ptot, rp[44];
inline void Div(int x) {
	for (int i = 2; i * i <= x; ++i) {
		if (x % i == 0) {
			p[++ptot] = i; rp[ptot] = 1;
			while (x % i == 0)	x /= i, ++c[ptot], rp[ptot] *= i;
		}
	}
	if (x > 1)	p[++ptot] = x, c[ptot] = 1, rp[ptot] = x;
}

inline int quickpow(int x, int k, int P)

int exgcd(int a, int b, int &x, int &y)

int calc(int a, int b) {//ax + by = 1 (gcd(a,b)=1), return x
	int x, y; exgcd(a, b, x, y);
	return (x % b + b) % b;
}

int nw;
int yu, ct;
int Prod[101000];
int fakeans[44];
void sol(int n) {
	if (!n)	return ;
	yu = yu * quickpow(Prod[rp[nw]], n / rp[nw], rp[nw]) % rp[nw] * Prod[n % rp[nw]] % rp[nw];
	ct += n / p[nw];
	sol(n / p[nw]);
}

inline void merge() {
	int ans = 0;
	for (int i = 1; i <= ptot; ++i) {
		int tmp = calc(realP / rp[i], rp[i]);
		ans = (ans + fakeans[i] * (realP / rp[i]) % realP * tmp) % realP;
	}
	printf("%lld\n", (ans % realP + realP) % realP);
}

signed main() {//calculate (n!)/(w_1! w_2!...)
	Div(realP);
	for (int i = 1; i <= ptot; ++i) {
		nw = i;
		yu = 1, ct = 0;
		Prod[0] = 1;
		for (int j = 1; j <= rp[i]; ++j)
			if (j % p[i])	Prod[j] = Prod[j - 1] * j % rp[i];
			else	Prod[j] = Prod[j - 1];
		sol(n);
		int memo = yu, memoct = ct; yu = 1; ct = 0;
		for (int j = 1; j <= m; ++j)	sol(w[j]);
		yu = calc(yu, rp[nw]) * memo % rp[i];
		fakeans[i] = yu * quickpow(p[i], memoct - ct, rp[i]) % rp[i];
		yu = 1; ct = 0;
	}
	merge();
}
posted @ 2020-10-22 17:49  JiaZP  阅读(194)  评论(0编辑  收藏  举报