【科技】扩展Lucas随想
2018-07-24 15:10
211 查看
扩展Lucas解决的还是一个很Simple的问题:
求:$C_{n}^{m} \; mod \; p$。
其中$n,m$都会比较大,而$p$不是很大,而且不一定是质数。
扩展Lucas可以说和Lucas本身并没有什么关系,重要的是中国剩余定理。扩展Lucas这个算法中教会我们的除了算组合数,还有在模数不是质数的时候,往往可以用$CRT$来合并答案。
将原模数质因数分解:$P = \prod_{i = 1}^{m} p_{i}^{k_{i}}$。
列出$m$个同余方程,第$i$个形如:$C_{n}^{m} \; \equiv a_{i} (mod \; p_{i}^{k_{i}})$。
由于$m$个方程中模数互质,则$CRT$后就是原答案。
现在来对于某个方程求解$a_{i}$是多少,即$C_{n}^{m} \; mod \; p^{k}$的答案。
把组合数转化成阶乘:$\frac{n!}{m!(n - m)!}$,我们先求一个阶乘在$mod \; p^{k}$下的值,设这个函数为$Fac(n)$。
常规的对于式子下方的阶乘我们需要求逆元,而阶乘中存在$p$的倍数,这意味可能不与$p^{k}$互质。为了解决这个问题,我们将有关$p$单独考虑,于是一个算阶乘的函数将包括两部分:
- 首先考虑所有$p$的倍数,总共有$\lfloor \frac{n}{p} \rfloor$个,将$p$提出来,这$\lfloor \frac{n}{p} \rfloor$个数又成为一个阶乘的形式,递归即可,总层数不会超过$log$。这部分的答案就是$p^{\lfloor \frac{n}{p} \rfloor} * Fac(\lfloor \frac{n}{p} \rfloor)$。
- 剩下的数都将与$p^{k}$互质。我们考虑以$p^{k}$分块,我们可以证明每段$p^{k}$中所有不是$p$的倍数的数的乘积在模$p^{k}$意义下是相同的。具体原因在于$i + p^{k} \equiv i (mod \; p^{k})$。通过暴力计算,这部分的复杂度就是$O(p^{k})$的。
接下来就没有什么问题了,用扩展欧几里得求逆元,有关$p$的幂次在除法时指数相减就行了。
#include <cstdio> typedef long long LL; int P; int Pow(int x, LL b, int p) { static int re; for (re = 1; b; b >>= 1, x = (LL) x * x % p) if (b & 1) re = (LL) re * x % p; return re; } int Ex_gcd(int a, int b, int &x, int &y) { if (b == 0) return x = 1, y = 0, a; int gcd = Ex_gcd(b, a % b, y, x); y -= a / b * x; return gcd; } int Inv(int a, int p) { static int x, y; int gcd = Ex_gcd(a, p, x, y); if (gcd != 1) throw; return (x % p + p) % p; } int Fac(LL n, int p, int pk) { if (n == 0) return 1; int re = 1; for (int i = 1; i <= pk; ++i) if (i % p != 0) re = (LL) re * i % pk; re = Pow(re, n / pk, pk); for (int i = 1; i <= n % pk; ++i) if (i % p != 0) re = (LL) re * i % pk; return (LL) re * Fac(n / p, p, pk) % pk; } int Crt(LL n, LL m, int p, int pk) { int fn = Fac(n, p, pk), fm = Fac(m, p, pk), fnm = Fac(n - m, p, pk); int cp = 0; for (LL i = n; i; i /= p) cp += i / p; for (LL i = m; i; i /= p) cp -= i / p; for (LL i = n - m; i; i /= p) cp -= i / p; int a = (LL) fn * Inv(fm, pk) % pk * Inv(fnm, pk) % pk * Pow(p, cp, pk) % pk; return (LL) a * (P / pk) % P * Inv(P / pk, pk) % P; } int Lucas(LL n, LL m, int p) { int re = 0, x = p; for (int i = 2; i <= p; ++i) { if (x % i != 0) continue; int pk = 1; while (x % i == 0) pk *= i, x /= i; re = (re + Crt(n, m, i, pk)) % p; } return re; } int main() { LL n, m; scanf("%lld%lld%d", &n, &m, &P); printf("%d\n", Lucas(n, m, P)); return 0; }View Code
相关文章推荐
- 【数学/扩展欧几里得/Lucas定理】BZOJ 1951 :[Sdoi 2010]古代猪文
- BZOJ_P1951&Codevs_P1830 [SDOI2010]古代猪文(Lucas定理+扩展欧几里得+中国剩余定理)
- NKOJ 2640 (SDOI 2013)方程(扩展Lucas+容斥原理)
- 扩展lucas定理bzoj2142待修改(这是我自己yy的。。。所以看看就好23333)
- [BZOJ2142]礼物(扩展Lucas定理+中国剩余定理)
- BZOJ3738 [Ontak2013]Kapitał 【扩展Lucas】
- BZOJ 1951 [Sdoi2010]古代猪文 欧拉定理+(扩展)lucas定理+逆元+快速幂
- [BZOJ 2142]礼物:扩展Lucas定理
- 乘法逆元+扩展欧几里得+Lucas
- BZOJ4830 [Hnoi2017]抛硬币 【扩展Lucas】
- lucas定理及扩展lucas定理
- Lucas定理与扩展Lucas
- 扩展CRT&&扩展lucas
- 【SDOI2010】古代猪文 扩展Lucas+中国剩余定理
- hdu3037 Saving Beans(Lucas定理+费马小定理or扩展欧几里德算发)
- [BZOJ3129][Sdoi2013]方程(容斥原理+扩展lucas)
- [BZOJ3129][SDOI2013]方程(扩展Lucas定理+容斥)
- [BZOJ4830][HNOI2017]抛硬币-扩展Lucas定理
- GYM100633J. Ceizenpok’s formula 扩展lucas模板
- [BZOJ 1951] 古代猪文【Lucas定理/费马小定理/中国剩余定理/扩展欧几里得】