【模板】扩展卢卡斯定理/exLucas
题目链接:luogu P4720
题目大意
要你求一个组合数在模一个数的结果。
模数不一定是质数。
思路
我们考虑质因数分解模数。
考虑分成了 p = p 1 a 1 , p 2 a 2 , . . . p=p_1^{a_1},p_2^{a_2},... p=p1a1,p2a2,...
然后你考虑求出:
{
C
n
m
m
o
d
p
1
a
1
C
n
m
m
o
d
p
2
a
2
…
…
\begin{cases}C_{n}^m\mod p_1^{a_1}\\ C_{n}^m\mod p_2^{a_2}\\ ……\end{cases}
⎩⎪⎨⎪⎧Cnmmodp1a1Cnmmodp2a2……
不难看出求出来的结果我们可以直接用 CRT(中国剩余定理)合并。
那接着就变成如何求
C
n
m
m
o
d
p
k
C_n^m\mod p^{k}
Cnmmodpk。
首先我们用普通的式子:
n
!
m
!
(
n
−
m
)
!
m
o
d
p
k
\dfrac{n!}{m!(n-m)!}\mod p^k
m!(n−m)!n!modpk
显然是不能直接搞,因为我们不能求出下面的两个的逆元。
然后又逆元的条件是这个数和模数互质,所以我们可以换一下式子:
n
!
p
x
m
!
p
y
(
n
−
m
)
!
p
z
p
x
−
y
−
z
m
o
d
p
k
\dfrac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z}\mod p^k
pym!pz(n−m)!pxn!px−y−zmodpk
然后 x x x 是 n ! n! n! 中包含 p p p 因子的个数。( y , z y,z y,z 同理)
那问题就是求
n
!
p
x
m
o
d
p
k
\dfrac{n!}{p^x}\mod p^k
pxn!modpk
考虑处理一下
n
!
n!
n!:
=
(
p
⋅
2
p
⋅
3
p
⋅
.
.
.
)
(
1
⋅
2
⋅
3
⋅
.
.
.
)
=(p\cdot 2p\cdot3p\cdot...)(1\cdot 2\cdot 3\cdot...)
=(p⋅2p⋅3p⋅...)(1⋅2⋅3⋅...)
(左边是
1
∼
n
1\sim n
1∼n 中
p
p
p 的倍数,右边就是不是倍数的)
然后提出左边的
p
p
p:
=
p
⌊
n
p
⌋
(
⌊
n
p
⌋
)
!
∏
i
=
1
,
i
≢
0
(
m
o
d
p
)
n
i
=p^{\left\lfloor\frac{n}{p}\right\rfloor}(\left\lfloor\dfrac{n}{p}\right\rfloor)!\prod\limits_{i=1,i\not\equiv0\pmod p}^ni
=p⌊pn⌋(⌊pn⌋)!i=1,i≡0(modp)∏ni
然后后面的显然是有循环节:
=
p
⌊
n
p
⌋
(
⌊
n
p
⌋
)
!
(
∏
i
=
1
,
i
≢
0
(
m
o
d
p
)
p
k
i
)
⌊
n
p
k
⌋
(
∏
i
=
p
k
⌊
n
p
k
⌋
,
i
≢
0
(
m
o
d
p
)
n
i
)
=p^{\left\lfloor\frac{n}{p}\right\rfloor}(\left\lfloor\dfrac{n}{p}\right\rfloor)!(\prod\limits_{i=1,i\not\equiv0\pmod p}^{p^k}i)^{\left\lfloor\frac{n}{p^k}\right\rfloor}(\prod\limits_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor,i\not\equiv0\pmod p}^{n}i)
=p⌊pn⌋(⌊pn⌋)!(i=1,i≡0(modp)∏pki)⌊pkn⌋(i=pk⌊pkn⌋,i≡0(modp)∏ni)
两个分别代表循环节和余项。
然后因为取模,
p
⌊
n
p
⌋
p^{\left\lfloor\frac{n}{p}\right\rfloor}
p⌊pn⌋ 会没掉,但是
(
⌊
n
p
⌋
)
!
(\left\lfloor\dfrac{n}{p}\right\rfloor)!
(⌊pn⌋)! 中可能还有
p
p
p。
那就是一个递归下去的过程,设
F
(
n
)
=
n
!
p
x
F(n)=\dfrac{n!}{p^x}
F(n)=pxn!
然后递推就是:
F
(
n
)
=
F
(
⌊
n
p
⌋
)
(
∏
i
=
1
,
i
≢
0
(
m
o
d
p
)
p
k
i
)
⌊
n
p
k
⌋
(
∏
i
=
p
k
⌊
n
p
k
⌋
,
i
≢
0
(
m
o
d
p
)
n
i
)
F(n)=F(\left\lfloor\dfrac{n}{p}\right\rfloor)(\prod\limits_{i=1,i\not\equiv0\pmod p}^{p^k}i)^{\left\lfloor\frac{n}{p^k}\right\rfloor}(\prod\limits_{i=p^k\left\lfloor\frac{n}{p^k}\right\rfloor,i\not\equiv0\pmod p}^{n}i)
F(n)=F(⌊pn⌋)(i=1,i≡0(modp)∏pki)⌊pkn⌋(i=pk⌊pkn⌋,i≡0(modp)∏ni)
(这个的复杂度是
log
p
n
\log_pn
logpn,边界是
f
(
0
)
=
1
f(0)=1
f(0)=1)
然后在原来的式子,就变成了
F
(
n
)
F
(
m
)
F
(
n
−
m
)
p
x
−
y
−
z
m
o
d
p
k
\dfrac{F(n)}{F(m)F(n-m)}p^{x-y-z}\mod p^k
F(m)F(n−m)F(n)px−y−zmodpk
这下分母的两个就是一定跟
p
k
p^k
pk 互质,就可以用扩欧求啦。
然后你会发现还有一个
p
x
−
y
−
z
p^{x-y-z}
px−y−z。
那我们就是要求
F
(
n
)
=
n
!
p
x
F(n)=\dfrac{n!}{p^x}
F(n)=pxn! 中的
x
x
x。
想想上面的式子,我们要的就是被除掉的部分,就是
p
⌊
n
p
⌋
p^{\left\lfloor\frac{n}{p}\right\rfloor}
p⌊pn⌋ 中的
⌊
n
p
⌋
\left\lfloor\dfrac{n}{p}\right\rfloor
⌊pn⌋。
每个
⌊
n
p
⌋
\left\lfloor\dfrac{n}{p}\right\rfloor
⌊pn⌋ 加起来就是结果,所以:
g
(
n
)
=
⌊
n
p
⌋
+
g
(
⌊
n
p
⌋
)
g(n)=\left\lfloor\dfrac{n}{p}\right\rfloor+g(\left\lfloor\dfrac{n}{p}\right\rfloor)
g(n)=⌊pn⌋+g(⌊pn⌋)
(复杂度也是
log
p
n
\log_pn
logpn,边界是
g
(
n
)
=
0
(
n
<
p
)
g(n)=0(n<p)
g(n)=0(n<p))
然后再带入到原来的式子就是:
F
(
n
)
F
(
m
)
F
(
n
−
m
)
p
G
(
n
)
−
G
(
m
)
−
G
(
n
−
m
)
m
o
d
p
k
\dfrac{F(n)}{F(m)F(n-m)}p^{G(n)-G(m)-G(n-m)}\mod p^k
F(m)F(n−m)F(n)pG(n)−G(m)−G(n−m)modpk
然后就好啦!!!
代码
#include<cstdio>
#define ll long long
using namespace std;
ll n, m, p;
ll a[1001], b[1001];
ll ksm(ll x, ll y, ll p) {
ll re = 1;
while (y) {
if (y & 1) re = re * x % p;
x = x * x % p;
y >>= 1;
}
return re;
}
ll exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) {
x = 1; y = 0;
return a;
}
ll re = exgcd(b, a % b, y, x);
y -= x * (a / b);
return re;
}
ll inv(ll a, ll p) {//别用快速幂求,要用扩欧
ll x, y;
exgcd(a, p, x, y);
return (x % p + p) % p;
}
ll F(ll n, ll p, ll pk) {
if (!n) return 1;
ll rou = 1, el = 1;
for (ll i = 1; i <= pk; i++) {//循环节部分
if (i % p) rou = rou * i % pk;
}
rou = ksm(rou, n / pk, pk);
for (ll i = pk * (n / pk); i <= n; i++) {//余项部分
if (i % p) el = el * (i % pk) % pk;
}
return F(n / p, p, pk) * rou % pk * el % pk;
}
ll G(ll n, ll p) {
if (n < p) return 0;
return n / p + G(n / p, p);
}
ll work_C(ll n, ll m, ll p, ll pk) {
ll fz = F(n, p, pk), fm1 = inv(F(m, p, pk), pk), fm2 = inv(F(n - m, p, pk), pk);
ll mi = ksm(p, G(n, p) - G(m, p) - G(n - m, p), pk);
return fz * fm1 % pk * fm2 % pk * mi % pk;
}
ll exLucas(ll n, ll m, ll p) {
ll tmpp = p, tot = 0, di;
for (ll i = 2; i * i <= tmpp; i++) {
if (tmpp % i == 0) {
tot++; di = 1;
while (tmpp % i == 0) {
di *= i;
tmpp /= i;
}
a[tot] = di;
b[tot] = work_C(n, m, i, di);
}
}
if (tmpp > 1) a[++tot] = tmpp, b[tot] = work_C(n, m, tmpp, tmpp);
ll ans = 0;//CRT
for (int i = 1; i <= tot; i++) {
ll M = p / a[i], V = inv(M, a[i]);
ans = (ans + b[i] * M % p * V % p) % p;
}
return ans;
}
int main() {
scanf("%lld %lld %lld", &n, &m, &p);
printf("%lld", exLucas(n, m, p));
return 0;
}