数论 · exLucas 定理
前言
终于看懂了!!!连夜丢掉卡特兰来做了两道 exLucas 的题目。
定理
作为一个没什么用的铺垫:数论 · Lucas 定理
求解 C n m % p C_n^m \%p Cnm%p(不保证 p p p 为质数)。
证明
1
先把 p p p 拆分一下: p = p 1 a 1 ∗ p 2 a 2 ∗ ⋯ p k a k p = p_1^{a_1} * p_2^{a_2} * \cdots p_k^{a_k} p=p1a1∗p2a2∗⋯pkak。
然后逐一求解 C n m % p i a i C_n^m \% p_i^{a_i} Cnm%piai,
最后用 CRT 求解
{
x
≡
b
1
(
m
o
d
p
1
a
1
)
x
≡
b
2
(
m
o
d
p
2
a
2
)
⋯
x
≡
b
k
(
m
o
d
p
k
a
k
)
2
假设求解 C n m % p i k C_n^m \% p_i^{k} Cnm%pik,设 p k = p i k p_k=p_i^k pk=pik。
有 C n m = n ! m ! ( n − m ) ! C_n^m=\dfrac{n!}{m!(n - m)!} Cnm=m!(n−m)!n!,
所以就变成求解 n ! % p k m ! % p k ( n − m ) ! % p k \dfrac{n!\%p_k}{m!\%p_k(n - m)!\%p_k} m!%pk(n−m)!%pkn!%pk。
但是因为分母部分需要逆元,有可能 m ! % p k m!\%p_k m!%pk 或 ( n − m ) ! % p k (n - m)!\%p_k (n−m)!%pk 与 p k p_k pk 不互质,所以我们要先除去他们中所有的 p i p_i pi,然后取模完之后再乘回去即可。
假设 x x x 为 n ! n! n! 中包含的 p i p_i pi 的个数, y y y、 z z z 同理。
则变为求解 ( n ! p i x % p k m ! p i y % p k ( n − m ) ! p i z % p k ∗ p i x − y − z ) % p k (\dfrac{\dfrac{n!}{p_i^x} \% p_k}{\dfrac{m!}{p_i^y} \% p_k\dfrac{(n-m)!}{p_i^z} \% p_k} * p_i^{x-y-z}) \% p_k (piym!%pkpiz(n−m)!%pkpixn!%pk∗pix−y−z)%pk
3
考虑求解 n ! p i x % p k \dfrac{n!}{p_i^x}\% p_k pixn!%pk。
先举个例子:
设 n = 19 , p i = 3 , k = 2 n=19,\ p_i = 3,\ k = 2 n=19, pi=3, k=2。
先把 n ! n! n! 中 p i p_i pi 的倍数先提取出来再合并一下:
n = ( 1 ∗ 2 ∗ 4 ∗ 5 ∗ 7 ∗ 8 ∗ 10 ∗ 11 ∗ 13 ∗ 14 ∗ 16 ∗ 17 ∗ 19 ) ∗ 3 6 ∗ ( 1 ∗ 2 ∗ 3 ∗ 4 ∗ 5 ∗ 6 ) n=(1* 2* 4* 5* 7* 8* 10* 11* 13* 14* 16* 17* 19) * 3^6 * (1* 2* 3* 4* 5* 6) n=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗36∗(1∗2∗3∗4∗5∗6)
按照 2 的后半部分所述, 3 6 3^6 36 可以先忽略掉。
而柿子的后半部分 ( 1 ∗ 2 ∗ 3 ∗ 4 ∗ 5 ∗ 6 ) (1* 2* 3* 4* 5* 6) (1∗2∗3∗4∗5∗6),也就是 ⌊ n p i ⌋ ! \left\lfloor\dfrac{n}{p_i}\right\rfloor! ⌊pin⌋! 我们可以用递归再次求解。
对于柿子的前半部分,会发现它对 p i k p_i^k pik 取余会有一个周期:
( 1 ∗ 2 ∗ 4 ∗ 5 ∗ 7 ∗ 8 ) ≡ ( 10 ∗ 11 ∗ 13 ∗ 14 ∗ 16 ∗ 17 ) ( m o d p i k ) (1* 2* 4* 5* 7* 8) \equiv (10 * 11 * 13 * 14 * 16 * 17) \pmod {p_i^k} (1∗2∗4∗5∗7∗8)≡(10∗11∗13∗14∗16∗17)(modpik)
所以只用求出第一个周期: r e s = ∏ i = 1 p k i ( p i ∤ i ) res=\prod_{i=1}^{p_k} i\ ( p_i \nmid i) res=∏i=1pki (pi∤i),然后再用快速幂求出 r e s ⌊ n p k ⌋ res^{\left\lfloor\frac{n}{p_k}\right\rfloor} res⌊pkn⌋ 即可。
然后就会发先剩下个 19 19 19,这些数直接暴力乘进 r e s res res(进行完快速幂之后)即可。
4
关于 2 中 x , y , z x,\ y,\ z x, y, z 的计算。
小奥中学过试除法(忘了叫啥),可以求解类似“100 的阶乘中有多少个因子 5”的问题。
这里直接除然后求出 x − y − z x-y-z x−y−z 即可。
代码如下:
int al = n;
while(al)
k += al / pi, al /= pi;
al = m;
while(al)
k -= al / pi, al /= pi;
al = n - m;
while(al)
k -= al / pi, al /= pi;
最后用 扩展欧几里得 求出 m ! p i y % p k ( n − m ) ! p i z % p k \dfrac{m!}{p_i^y} \% p_k\dfrac{(n-m)!}{p_i^z} \% p_k piym!%pkpiz(n−m)!%pk 的逆元即可。
5
关于使用 CRT 合并。
每次求完 C n m % p i k = b i C_n^m \% p_i^{k}=b_i Cnm%pik=bi,就有了 x ≡ b i ( m o d p i a i ) x\equiv{b_i}\pmod{p_i^{a_i}} x≡bi(modpiai)。
先给代码:ans += bi * inv(p / pi, pi) % p * (p / pi) % p
bi * inv(p / pi, pi) % p
是为了让它在取模
p
i
p_i
pi 的情况下余数为
b
i
b_i
bi;
* (p / pi) % p
是为了让它在取模其他质数的情况下余数为 0。
详见 CRT。
代码
#include<algorithm>
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
#define int long long
#define rep(i, a, b) for(int i = a; i <= b; ++i)
int n, m, p;
inline int power(int x, int t, int mod)
{
int res = 1;
while(t >= 1)
{
if(t & 1)
res = res * x % mod;
t /= 2, x = x * x % mod;
}
return res;
}
inline int fac(int nw, int pi, int pk)
{
if(!nw) return 1;
int res = 1;
rep(i, 2, pk)
if(i % pi) res = res * i % pk;
res = power(res, nw / pk, pk);
rep(i, 2, nw % pk)
if(i % pi) res = res * i % pk;
return res * fac(nw / pi, pi, pk) % pk;
}
int x, y;
inline void exGcd(int a, int b)
{
if(!b)
{
x = 1, y = 0;
return;
}
exGcd(b, a % b);
int tmp = x;
x = y, y = tmp - a / b * y;
return;
}
inline int inv(int nw, int mod)
{
exGcd(nw, mod);
return (x + mod) > mod ? x : x + mod;
}
inline int C(int pi, int pk)
{
int fn = fac(n, pi, pk), fm = fac(m, pi, pk), fnm = fac(n - m, pi, pk);
int k = 0;
int al = n;
while(al)
k += al / pi, al /= pi;
al = m;
while(al)
k -= al / pi, al /= pi;
al = n - m;
while(al)
k -= al / pi, al /= pi;
return fn * inv(fm, pk) % pk * inv(fnm, pk) % pk * power(pi, k, pk) % pk;
}
inline int Crt(int bi, int mi)
{
return bi * inv(p / mi, mi) % p * (p / mi) % p;
}
inline int exLucas()
{
int pk, res = 0;
int lim = ceil(sqrt(p)), tmp = p;
rep(i, 2, lim)
{
if(tmp % i) continue;
pk = 1;
while(!(tmp % i))
tmp /= i, pk *= i;
res = (res + Crt(C(i, pk), pk)) % p;
}
if(tmp > 1)
res = (res + Crt(C(tmp, tmp), tmp)) % p;
return res;
}
signed main()
{
scanf("%lld %lld %lld", &n, &m, &p);
printf("%lld\n", exLucas());
return 0;
}
—— E n d End End——
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 地球OL攻略 —— 某应届生求职总结
· 周边上新:园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源!
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】