lucas及其拓展

模板题洛谷 P4720

本文侧向结论和代码实现，

推导请转至lucas定理及其拓展的推导 https://blog.csdn.net/yuyilahanbao/article/details/100568285

lucas定理

请阅读lucas定理 https://blog.csdn.net/yuyilahanbao/article/details/100550317

拓展lucas的结论

(nm)≡pk1−k2−k3×b1×b2−1×b3−1×cu1−u2−u3×k1!k2!×k3!(modpw)\tbinom{n}{m} \equiv
{
p^{k_1-k_2-k_3} \times b_1 \times {b_2}^{-1} \times {b_3}^{-1}
\times c^{u_1-u_2-u_3}
\times \frac{{k_1}!}{{k_2}! \times {k_3}!}
} \pmod{p^w}(mn)≡pk1−k2−k3×b1×b2−1×b3−1×cu1−u2−u3×k2!×k3!k1!(modpw).

当r1≥r2r_1 \geq r_2r1≥r2时，k1=k2+k3k_1=k_2+k_3k1=k2+k3,故上面这个式子最后分式的部分k1!k2!×k3!=(k1k2)\frac{{k_1}!}{{k_2}! \times {k_3}!}=\tbinom{k_1}{k_2}k2!×k3!k1!=(k2k1).
当r1<r2r_1 < r_2r1<r2时，k1=k2+k3+1k_1=k_2+k_3+1k1=k2+k3+1,最后的那个分式无法直接变成组合数，但是我们只需要分子分母同时乘以k1−k2k_1-k_2k1−k2，即可变成组合数。k1!k2!×k3!=(k1−k2)×(k1k2)\frac{{k_1}!}{{k_2}! \times {k_3}!}=(k_1-k_2) \times \tbinom{k_1}{k_2}k2!×k3!k1!=(k1−k2)×(k2k1)

各个字母代表的含义。

与n,m,n-m有关的量k,r,u,v分别用下标1,2,3区分.

k,r是除以ppp的商与余数，u,v是除以模数pwp^wpw的商与余数。

b是n!n!n!(m!m!m!或(n−m)!(n-m)!(n−m)!)最后剩下的v个数中不是p的倍数的数的乘积。

c是[1,pw]\left[1,p^w\right][1,pw]中不是p的倍数的数的乘积。

从结论中的式子可以看到b,c我们只关注模pkp^kpk意义下的值，因此可以预先求出[1..i] (i<pw)[1..i] \; (i < p^w)[1..i](i<pw)中不是p的倍数的数的乘积f(i)f(i)f(i)(模pkp^kpk意义下的)。

(nm)&VeryThinSpace;mod&VeryThinSpace;N\tbinom{n}{m} \bmod N(mn)modN的求取

N是任意正整数。对NNN进行素数分解。N=∏i=1qpikiN=\prod\limits_{i=1}^{q}p_i^{k_i}N=i=1∏qpiki.

对(nm)&VeryThinSpace;mod&VeryThinSpace;piki\tbinom{n}{m} \bmod p_i^{k_i}(mn)modpiki问题，可以通过上一小节的拓展lucas求得，记答案是cic_ici.

于是得到了qqq个线性同余方程,即线性同余方程组(nm)≡ci(modpiki)(1≤i≤q)\tbinom{n}{m} \equiv c_i \pmod{p_i^{k_i}} \quad (1 \leq i \leq q)(mn)≡ci(modpiki)(1≤i≤q).

对于线性同余方程组，并且注意到模数pikip_i^{k_i}piki两两互质，可以用中国剩余定理（也可以用拓欧）解出其通解x=x0+ktx=x_0+ktx=x0+kt。并且由于模数互质，k=lcm(piki)=N(1≤i≤q)k=lcm(p_i^{k_i})=N \quad (1 \leq i \leq q)k=lcm(piki)=N(1≤i≤q).所以在[0,N)[0,N)[0,N)内只有一个特解x0x_0x0,而这个特解就是(nm)&VeryThinSpace;mod&VeryThinSpace;N\tbinom{n}{m} \bmod N(mn)modN.

核心代码

ll pow(ll a, ll n)

{

	if (n == 0) return 1;

	// 始终维持要求的数可以表示为(a)^n*t

	ll t = 1;

	while (n > 1)

	{

		if (n&1) t = t*a;

		n >>= 1; a = a*a;

	}

	return a*t; // now n = 1

}

// 质因数分解，p_i^{k_i} 共q项 返回q

int factor(ll n, vector<ll>&p, vector<int>&k) {

	p.clear(); k.clear();

	if (n <= 1) return 0;

	int q = 0;

	for (ll i = 2; i*i <= n; ++i) { // 不必担心溢出，因为会溢出说明肯定会tle

		if (!(n%i)) {

			p.push_back(i);

			k.push_back(0);

			do {n /= i; ++k[q];} while (!(n%i));

			++q;

		}

	}

	if (n > 1) {

		p.push_back(n);

		k.push_back(1);

		++q;

	}

	return q;

}

// 求C(n,m)%(p^k)

const ll kMaxPk = 1000000;

// f[i]表示1..i中不是p的倍数的数的乘积(%pk) inv_f则是相应的逆元

ll f[kMaxPk],inv_f[kMaxPk];

ll ex_lucas(ll n, ll m, ll p, ll k, ll pk) {

	ll k1,k2,k3,r1,r2,r3,u1,u2,u3,v1,v2,v3;

	ll ans = 1;

	f[0] = 1; inv_f[0] = 1;

	for (ll j = 1; j < pk; ++j) {

		if (j%p) {

			f[j] = (f[j-1]*j)%pk;

			multiplicative_inverse(f[j],pk,inv_f[j]); // 肯定存在逆元

		} else {

			f[j] = f[j-1];

			inv_f[j] = inv_f[j-1];

		}

	}

	while(1) {

		if (m == 0) return ans;

		k1 = n/p, r1 = n%p;

		k2 = m/p, r2 = m%p;

		k3 = (n-m)/p, r3 = (n-m)%p;

		u1 = n/pk, v1 = n%pk;

		u2 = m/pk, v2 = m%pk;

		u3 = (n-m)/pk, v3 = (n-m)%pk;

		if (k1-k2-k3) { // == 1

			ans = (ans*p)%pk;

		} // else == 0

		ans = (ans*f[v1])%pk;

		ans = (ans*inv_f[v2])%pk;

		ans = (ans*inv_f[v3])%pk;

		if (u1-u2-u3) { // == 1

			ans = (ans*f[pk-1])%pk;

		} // else == 0

		if (r1 < r2) ans = ans*((k1-k2)%pk)%pk;

		n = k1; m = k2;

	}

}

const int kMaxQ = 30;

ll a[kMaxQ];

ll c[kMaxQ];

ll mm[kMaxQ];

// 返回C(n,m)%N N的分解因式后最大的x=p^k

// 可以开x大小的数组

ll ex_lucas_N(ll n, ll m, ll N) {

	vector<ll>p;

	vector<int>k;

	int q = factor(N,p,k);

	for (int i = 0; i < q; ++i) {

		a[i] = 1;

		mm[i] = pow(p[i],k[i]);

		c[i] = ex_lucas(n,m, p[i], k[i], mm[i]);

	}

	ll ans,kk;

	linear_congruence_equations(q,a,c,mm,ans,kk);

	return ans;

}

完整ac代码

// luogu p4720

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

// 求解不定方程ax+by=(a,b)的一组特解并返回a,b最大公约数

// x,y存储返回的一组特解。

ll ex_gcd(ll a, ll b, ll &x, ll &y) {

	if (b) {

		auto d = ex_gcd(b, a%b, y, x); // 注意x和y位置互换了。

		// x是后，无需赋值，y是 前-a/b*后 即 y -= a/b*x

		y -= a/b*x;

		return d;

	} else {

		x = 1; y = 0;

		return a;

	}

}

// 求解不定方程ax+by=c.

// 返回值表示是否有解

// d存储是(a,b)

// 当有解的情况下

// x,y存储一组特解,并且确保x是最小的非负整数。

// 通解是X=x+(b/d)*t,Y=y-(a/d)*t t是整数。

bool binary_linear_indefinite_equation(ll a, ll b, ll c, ll &x, ll &y, ll &d) {

	d = ex_gcd(a, b, x, y); // solve: ax+by=(a,b)

	if (c%d) return false;

	x *= c/d;

	y *= c/d;

	auto k = b/d;

	x = (x%k+k)%k; // 调为最小非负整数

	y = (d-a*x)/b;

	return true;

}

// 线性同余方程 Linear congruence equation

// ax = c (mod m) <===> ax+my=c

// x存储最小非负解，通解X=x+kt t为整数

// 有解的情况下，最小非负解x肯定在[0,m)范围内

bool linear_congruence_equation(ll a, ll c, ll m, ll &x, ll &k) {

	ll y, d;

	auto ans = binary_linear_indefinite_equation(a, m, c, x, y, d);

	k = m/d;

	return ans;

}

// 求a在Zm<+,*>中的乘法逆元x

// 返回逆元是否存在,x存储逆元

// ax = 1 (mod m)

bool multiplicative_inverse(ll a, ll m, ll &x) {

	ll k;

	return linear_congruence_equation(a, 1, m, x, k);

	// assert(k == m);

}

// 线性同余方程组 Linear congruence equations

// a_ix = c_i (mod m_i) 共n个

// 可能存在的问题，由于迭代过程中k一直在求最小公倍数，所以可能会爆long long，这个，最佳的方法是直接暴力把ll的定义改为__int128

// 但是要注意__int128的输入输出

// 如果还是爆，我没法子了

bool linear_congruence_equations(int n, ll a[], ll c[], ll m[], ll &x, ll &k) {

	ll x_i, k_i, t, t_i, d;

	x = 0; k = 1;

	for (int i = 0; i < n; ++i) {

		if (!linear_congruence_equation(a[i], c[i], m[i], x_i, k_i))

			return false;

		// kt+x

		// k_it_i+x_i

		if (!binary_linear_indefinite_equation(

			k, k_i, x_i-x, t, t_i, d))

			return false;

		x += k*t;

		k *= k_i/d;

	}

	return true;

}

ll pow(ll a, ll n)

{

	if (n == 0) return 1;

	// 始终维持要求的数可以表示为(a)^n*t

	ll t = 1;

	while (n > 1)

	{

		if (n&1) t = t*a;

		n >>= 1; a = a*a;

	}

	return a*t; // now n = 1

}

int factor(ll n, vector<ll>&p, vector<int>&k) {

	p.clear(); k.clear();

	if (n <= 1) return 0;

	int q = 0;

	for (ll i = 2; i*i <= n; ++i) { // 不必担心溢出，因为会溢出说明肯定会tle

		if (!(n%i)) {

			p.push_back(i);

			k.push_back(0);

			do {n /= i; ++k[q];} while (!(n%i));

			++q;

		}

	}

	if (n > 1) {

		p.push_back(n);

		k.push_back(1);

		++q;

	}

	return q;

}

// 求C(n,m)%(p^k)

const ll kMaxPk = 1000000;

// f[i]表示1..i中不是p的倍数的数的乘积(%pk) inv_f则是相应的逆元

ll f[kMaxPk],inv_f[kMaxPk];

ll ex_lucas(ll n, ll m, ll p, ll k, ll pk) {

	ll k1,k2,k3,r1,r2,r3,u1,u2,u3,v1,v2,v3;

	ll ans = 1;

	f[0] = 1; inv_f[0] = 1;

	for (ll j = 1; j < pk; ++j) {

		if (j%p) {

			f[j] = (f[j-1]*j)%pk;

			multiplicative_inverse(f[j],pk,inv_f[j]); // 肯定存在逆元

		} else {

			f[j] = f[j-1];

			inv_f[j] = inv_f[j-1];

		}

	}

	while(1) {

		if (m == 0) return ans;

		k1 = n/p, r1 = n%p;

		k2 = m/p, r2 = m%p;

		k3 = (n-m)/p, r3 = (n-m)%p;

		u1 = n/pk, v1 = n%pk;

		u2 = m/pk, v2 = m%pk;

		u3 = (n-m)/pk, v3 = (n-m)%pk;

		if (k1-k2-k3) { // == 1

			ans = (ans*p)%pk;

		} // else == 0

		ans = (ans*f[v1])%pk;

		ans = (ans*inv_f[v2])%pk;

		ans = (ans*inv_f[v3])%pk;

		if (u1-u2-u3) { // == 1

			ans = (ans*f[pk-1])%pk;

		} // else == 0

		if (r1 < r2) ans = ans*((k1-k2)%pk)%pk;

		n = k1; m = k2;

	}

}

const int kMaxQ = 30;

ll a[kMaxQ];

ll c[kMaxQ];

ll mm[kMaxQ];

// 返回C(n,m)%N N的分解因式后最大的x=p^k

// 可以开x大小的数组

ll ex_lucas_N(ll n, ll m, ll N) {

	vector<ll>p;

	vector<int>k;

	int q = factor(N,p,k);

	for (int i = 0; i < q; ++i) {

		a[i] = 1;

		mm[i] = pow(p[i],k[i]);

		c[i] = ex_lucas(n,m, p[i], k[i], mm[i]);

	}

	ll ans,kk;

	linear_congruence_equations(q,a,c,mm,ans,kk);

	return ans;

}

int main()

{

	ll n,m,N;

	scanf("%lld%lld%lld",&n,&m,&N);

	ll ans = ex_lucas_N(n,m,N);

	printf("%lld\n",ans);

	return 0;

}

巴特西

拓展lucas结论及模板

lucas及其拓展

lucas定理

拓展lucas的结论

(nm)&VeryThinSpace;mod&VeryThinSpace;N\tbinom{n}{m} \bmod N(mn)modN的求取

核心代码

完整ac代码

最新文章

热门文章

巴特西

拓展lucas结论及模板

lucas及其拓展

lucas定理

拓展lucas的结论

(nm)&VeryThinSpace;mod&VeryThinSpace;N\tbinom{n}{m} \bmod N(mn​)modN的求取

核心代码

完整ac代码

最新文章

热门文章

(nm)&VeryThinSpace;mod&VeryThinSpace;N\tbinom{n}{m} \bmod N(mn)modN的求取