引入
有卢卡斯定理, ( n m ) ≡ ( n m o d P m m o d P ) × ( ⌊ n P ⌋ ⌊ m P ⌋ ) ( m o d P ) \begin{pmatrix}n\\m\end{pmatrix}\equiv\begin{pmatrix}n\bmod P\\m\bmod P\end{pmatrix}\times\begin{pmatrix}\left\lfloor\frac{n}{P}\right\rfloor\\\left\lfloor\frac{m}{P}\right\rfloor\end{pmatrix}\pmod P (nm)≡(nmodPmmodP)×(⌊Pn⌋⌊Pm⌋)(modP)。卢卡斯定理的适用条件是: P P P 为质数。对于 P P P 是一般正整数的情况,可以使用扩展卢卡斯定理计算 C n m ( m o d P ) C_n^m\pmod P Cnm(modP)。
算法分析
原问题的转化
要求
C
n
m
m
o
d
P
C_n^m\bmod P
CnmmodP,首先对
P
P
P 进行质因数分解,即
P
=
p
1
k
1
p
2
k
2
⋯
p
c
n
t
k
c
n
t
P=p_1^{k_1}p_2^{k_2}\cdots p_{cnt}^{k_{cnt}}
P=p1k1p2k2⋯pcntkcnt,记
a
i
=
C
n
m
m
o
d
p
i
k
i
a_i=C_n^m\bmod p_i^{k_i}
ai=Cnmmodpiki,可以获得一个一元线性同余方程组。
(
S
)
:
{
C
n
m
≡
a
1
(
m
o
d
p
1
k
1
)
C
n
m
≡
a
2
(
m
o
d
p
2
k
2
)
⋮
C
n
m
≡
a
c
n
t
(
m
o
d
p
c
n
t
k
c
n
t
)
(S):\left\{\begin{array}{c} C_n^m\equiv a_{1}\pmod {p_{1}^{k_1}} \\ C_n^m\equiv a_{2}\pmod {p_{2}^{k_2}} \\ \vdots \\ C_n^m\equiv a_{cnt}\pmod {p_{cnt}^{k_{cnt}}} \end{array}\right.
(S):⎩⎪⎪⎪⎨⎪⎪⎪⎧Cnm≡a1(modp1k1)Cnm≡a2(modp2k2)⋮Cnm≡acnt(modpcntkcnt)
由于
p
1
,
p
2
,
⋯
,
p
c
n
t
p_1,p_2,\cdots,p_{cnt}
p1,p2,⋯,pcnt 是互不相同的质数,因此
p
1
k
1
,
p
2
k
2
,
⋯
,
p
c
n
t
k
c
n
t
p_1^{k_1},p_2^{k_2},\cdots,p_{cnt}^{k_{cnt}}
p1k1,p2k2,⋯,pcntkcnt 两两互质,由中国剩余定理,其解集为
A
=
{
x
∈
Z
∣
k
P
+
∑
i
=
1
c
n
t
a
i
t
i
P
i
,
k
∈
Z
}
A=\left\{x\in\mathbb Z\mid kP+\sum\limits_{i=1}^{cnt} a_it_iP_i,k\in \mathbb Z\right\}
A={x∈Z∣kP+i=1∑cntaitiPi,k∈Z};其中,
P
i
=
P
p
i
k
i
P_i=\frac{P}{p_{i}^{k_i}}
Pi=pikiP,
P
i
t
i
≡
1
(
m
o
d
P
)
P_it_i\equiv 1\pmod P
Piti≡1(modP)。尽管求出的
C
n
m
C_n^m
Cnm 并不是准确值,仅仅确定了
C
n
m
∈
A
C_n^m\in A
Cnm∈A,但是可以注意到,此集合中的元素在模
P
P
P 意义下的值都为
(
∑
i
=
1
c
n
t
a
i
t
i
P
i
)
m
o
d
P
\left(\sum\limits_{i=1}^{cnt} a_it_iP_i\right)\bmod P
(i=1∑cntaitiPi)modP;根据逻辑三段论,
C
n
m
≡
∑
i
=
1
c
n
t
a
i
t
i
P
i
(
m
o
d
P
)
C_n^m\equiv \sum\limits_{i=1}^{cnt} a_it_iP_i\pmod P
Cnm≡i=1∑cntaitiPi(modP)。所以,一元线性同余方程组
{
x
≡
a
1
(
m
o
d
p
1
k
1
)
x
≡
a
2
(
m
o
d
p
2
k
2
)
⋮
x
≡
a
c
n
t
(
m
o
d
p
c
n
t
k
c
n
t
)
\left\{\begin{array}{c} x\equiv a_{1}\pmod {p_{1}^{k_1}} \\ x\equiv a_{2}\pmod {p_{2}^{k_2}} \\ \vdots \\ x\equiv a_{cnt}\pmod {p_{cnt}^{k_{cnt}}} \end{array}\right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧x≡a1(modp1k1)x≡a2(modp2k2)⋮x≡acnt(modpcntkcnt)的最小正整数解即为
C
n
m
m
o
d
P
C_n^m\bmod P
CnmmodP。
问题在于,要如何先求出
a
i
=
C
n
m
m
o
d
p
i
k
i
a_i=C_n^m\bmod p_i^{k_i}
ai=Cnmmodpiki。
计算 C n m m o d p k C_n^m\bmod p^k Cnmmodpk
下面讨论如何计算
C
n
m
m
o
d
p
k
C_n^m\bmod p^k
Cnmmodpk,其中
p
p
p 为质数。
首先由组合数的公式
C
n
m
=
n
!
m
!
(
n
−
m
)
!
C_{n}^{m}=\frac{n!}{m!(n-m)!}
Cnm=m!(n−m)!n!,分别计算出
n
!
,
m
!
,
(
n
−
m
)
!
n!,m!,(n-m)!
n!,m!,(n−m)! 在模
p
k
p^k
pk 意义下的值,再求
m
!
,
(
n
−
m
)
!
m!,(n-m)!
m!,(n−m)! 模
p
k
p^k
pk 意义下的逆元,就可以得到答案。但是
m
!
m!
m! 和
(
n
−
m
)
!
(n-m)!
(n−m)! 可能包含质因子
p
p
p,根据裴蜀定理,此时它们在模
p
k
p^k
pk 意义下的逆元不存在。不妨将
n
!
,
m
!
,
(
n
−
m
)
!
n!,m!,(n-m)!
n!,m!,(n−m)! 中的质因子
p
p
p 除尽,对
C
n
m
C_n^m
Cnm 做如下变换。
C
n
m
=
n
!
p
x
m
!
p
y
×
(
n
−
m
)
!
p
z
×
p
x
−
y
−
z
\text{C}_{n}^{m}=\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\times\frac{(n-m)!}{p^z}}\times p^{x-y-z}
Cnm=pym!×pz(n−m)!pxn!×px−y−z
上式中,
x
x
x 为
n
!
n!
n! 中质因子
p
p
p 的个数,
y
y
y 为
m
!
m!
m! 中质因子
p
p
p 的个数,
z
z
z 为
(
n
−
m
)
!
(n-m)!
(n−m)! 中质因子
p
p
p 的个数。
n
!
p
x
\frac{n!}{p^x}
pxn!,
m
!
p
y
\frac{m!}{p^y}
pym!,
(
n
−
m
)
!
p
z
\frac{(n-m)!}{p^z}
pz(n−m)! 均与
p
k
p^k
pk 互质,都存在模
p
k
p^k
pk 意义下的逆元。
以上,我们已经解决了阶乘求逆元的问题。只要求得
n
!
p
x
,
m
!
p
y
,
(
n
−
m
)
!
p
z
\frac{n!}{p^x},\frac{m!}{p^y},\frac{(n-m)!}{p^z}
pxn!,pym!,pz(n−m)! 在模
p
k
p^k
pk 意义下的值,再求两次逆元,即可得到
C
n
m
m
o
d
p
k
C_n^m\bmod p^k
Cnmmodpk。但是,问题在于:如何求形如
n
!
p
x
m
o
d
p
k
\frac{n!}{p^x}\bmod p^k
pxn!modpk 的式子。
计算 n ! p x m o d p k \frac{n!}{p^x}\bmod p^k pxn!modpk
n
!
=
1
⋅
2
⋅
3
⋅
⋯
⋅
n
n!=1⋅2⋅3⋅\ \cdots\ ⋅n
n!=1⋅2⋅3⋅ ⋯ ⋅n=
(
p
⋅
2
p
⋅
3
p
⋯
)
(
1
⋅
2
⋅
⋯
)
(p\cdot 2p\cdot 3p\cdots)(1\cdot 2\cdot \cdots)
(p⋅2p⋅3p⋯)(1⋅2⋅⋯),左括号内为
1
∼
n
1\sim n
1∼n 中
p
p
p 的倍数的乘绩,右括号内反之。
1
∼
n
1\sim n
1∼n 中
p
p
p 的倍数有
⌊
n
p
⌋
\left\lfloor\frac{n}{p}\right\rfloor
⌊pn⌋ 个,故
n
!
=
p
⌊
n
p
⌋
⌊
n
p
⌋
!
∏
i
=
1
,
p
∤
i
n
i
n!=p^{\left\lfloor\frac{n}{p}\right\rfloor}\left\lfloor\frac{n}{p}\right\rfloor!\prod\limits_{i=1,p\nmid i}^n i
n!=p⌊pn⌋⌊pn⌋!i=1,p∤i∏ni。
计算
n
!
p
x
\frac{n!}{p^x}
pxn! 时,可以直接忽略
p
⌊
n
p
⌋
p^{\left\lfloor\frac{n}{p}\right\rfloor}
p⌊pn⌋ ,但是
⌊
n
p
⌋
!
\left\lfloor\frac{n}{p}\right\rfloor !
⌊pn⌋! 还可能包含质因子
p
p
p,难以直接忽略其中的质因子;不妨定义函数
f
a
c
(
n
,
p
,
p
k
)
=
n
!
p
x
m
o
d
p
k
fac(n,p,pk)=\frac{n!}{p^x}\bmod pk
fac(n,p,pk)=pxn!modpk,则
f
a
c
(
n
,
p
,
p
k
)
=
f
a
c
(
⌊
n
p
⌋
,
p
,
p
k
)
∏
i
=
1
,
p
∤
i
n
i
fac(n,p,pk)=fac\left(\left\lfloor\frac{n}{p}\right\rfloor,p,pk\right)\prod\limits_{i=1,p\nmid i}^n i
fac(n,p,pk)=fac(⌊pn⌋,p,pk)i=1,p∤i∏ni,如此一来可以进行递归求解。
(
∏
i
=
1
,
p
∤
i
n
i
)
m
o
d
p
k
\left(\prod\limits_{i=1,p\nmid i}^n i\right)\bmod p^k
(i=1,p∤i∏ni)modpk 是有循环节的,比如当
p
=
3
,
t
=
2
,
n
=
19
p=3,t=2,n=19
p=3,t=2,n=19 时,有:
19
!
=
3
6
×
(
1
×
2
×
3
×
4
×
5
×
6
)
×
(
1
×
2
×
4
×
5
×
7
×
8
×
10
×
11
×
13
×
14
×
16
×
17
×
19
)
19!=3^6\times(1×2×3×4×5×6)×(1×2×4×5×7×8×10×11×13×14×16×17×19)
19!=36×(1×2×3×4×5×6)×(1×2×4×5×7×8×10×11×13×14×16×17×19),右括号中
1
×
2
×
4
×
5
×
7
×
8
≡
10
×
11
×
13
×
14
×
16
×
17
(
m
o
d
3
2
)
1×2×4×5×7×8\equiv 10×11×13×14×16×17\pmod {3^2}
1×2×4×5×7×8≡10×11×13×14×16×17(mod32),可以看到,
(
∏
i
=
1
,
p
∤
i
n
i
)
m
o
d
p
k
\left(\prod\limits_{i=1,p\nmid i}^n i\right)\bmod p^k
(i=1,p∤i∏ni)modpk 以
p
k
p^k
pk 为一个周期。
(
∏
i
=
1
,
p
∤
i
n
i
)
m
o
d
p
k
\left(\prod\limits_{i=1,p\nmid i}^n i\right)\bmod p^k
(i=1,p∤i∏ni)modpk 是以
p
k
p^k
pk 为周期的,每个循环节长度为 $p^k $,所以只需要计算最后不满足一个周期的数的乘积(例子中只需要计算
19
19
19)。显然,完整的周期数为
⌊
n
p
k
⌋
\left\lfloor\frac{n}{p^k}\right\rfloor
⌊pkn⌋, 每个周期模
p
k
p^k
pk 意义下的成绩为
∏
i
=
1
,
p
∤
i
p
k
i
m
o
d
p
k
\prod\limits_{i=1,p\nmid i}^{p^k} i\bmod p^k
i=1,p∤i∏pkimodpk ;不满足一个周期的数的个数为
n
m
o
d
p
k
n\bmod p^k
nmodpk,模
p
k
p^k
pk 意义下的乘积为
∏
i
=
1
,
p
∤
i
n
m
o
d
p
k
i
m
o
d
p
k
\prod\limits_{i=1,p\nmid i}^{n\bmod p^k} i\bmod p^k
i=1,p∤i∏nmodpkimodpk。
洛谷P4720 【模板】扩展卢卡斯
题目背景
这是一道模板题。
题目描述
求 C n m m o d P C_n^m \bmod{P} CnmmodP,其中 C C C 为组合数。
输入格式
一行三个整数 n , m , P n,m,P n,m,P,含义由题所述。
输出格式
一行一个整数,表示答案。
输入输出样例
输入 #1
5 3 3
输出 #1
1
输入 #2
666 233 123456
输出 #2
61728
说明/提示
1 ≤ m ≤ n ≤ 1 0 18 1\leq m\leq n\leq 10^{18} 1≤m≤n≤1018, 2 ≤ P ≤ 1 0 6 2\leq P\leq 10^6 2≤P≤106,不保证 P P P 是质数。
代码
#include<iostream>
#define N 12
#define ll long long
using namespace std;
ll b[N],a[N];
int cnt;
//-------------------------------------------------------------------------------------------------------
//扩展欧几里得算法
//用于求逆元
void ex_gcd(ll a,ll b,ll &x,ll &y)
{
if(!b)
{
x=1;
y=0;
return;
}
ex_gcd(b,a%b,x,y);
ll temp=x;
x=y;
y=temp-a/b*y;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//运用扩展欧几里得算法求逆元
//ax≡1(mod b)
ll inv(ll a,ll b)
{
ll x,y;
ex_gcd(a,b,x,y);
//通解 x=x0+kb
return (x%b+b)%b;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//快速幂并取模
ll qpow_mod(ll a,ll b,ll mod)
{
ll res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod;
b>>=1;
}
return res%mod;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//阶乘除去质因子模质数幂
//pk表示p的k次幂
ll fac_mod(ll n,ll p,ll pk)
{
if(!n) return 1;
ll res=1;
int i;
//----------------------------------------
//一个循环节内的乘积
//周期为pk
for(i=1;i<=pk;i++) if(i%p) res=res*i%pk;
//n/pk个完整周期
res=qpow_mod(res,n/pk,pk);
//不满足一个周期的数的乘积
for(i=1;i<=n%pk;i++) if(i%p) res=res*i%pk;
//-----------------------------------------
//-----------------------------------------
//递归求 (n/p)!除去质因子p并对p的k次幂取模
return res*fac_mod(n/p,p,pk)%pk;
//-----------------------------------------
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//组合数对p的k次幂取模
ll C_mod(ll n,ll m,ll p,ll pk)
{
ll x,y,z;
x=y=z=0;
ll i;
//------------------------------
//阶乘除尽质因子
//最后乘p的(x-y-z)次幂
//n!中p质因子的数目
for(i=p;i<=n;i*=p) x+=n/i;
//m!中p质因子的数目
for(i=p;i<=m;i*=p) y+=m/i;
//(n-m)!中p质因子的数目
for(i=p;i<=n-m;i*=p) z+=(n-m)/i;
//-------------------------------
//n!/m!/(n-m)!
//求两次逆元
//三次阶乘除以质因子模质数幂
return fac_mod(n,p,pk)*inv(fac_mod(m,p,pk),pk)%pk*inv(fac_mod(n-m,p,pk),pk)%pk*qpow_mod(p,x-y-z,pk);
//--------------------------------------------------------------------------------------------------
}
//-------------------------------------------------------------------------------------------------------
//中国剩余定理
ll CRT()
{
int i;
ll mul=1;
ll res=0;
for(i=1;i<=cnt;i++) mul*=a[i];
for(i=1;i<=cnt;i++)
{
ll Mi=mul/a[i];
ll x,y;
res+=b[i]*inv(Mi,a[i])*Mi;
}
return res%mul;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//扩展卢卡斯定理
ll ex_Lucas(ll n,ll m,ll P)
{
ll i;
//-------------------------------------
//分解质因数
for(i=2;i*i<=P;i++)
{
if(P%i==0)
{
ll pk=1;
while(P%i==0)
{
pk*=i;
P/=i;
}
cnt++;//不同的质因子个数
//----------------------------
//质因子i
//pk表示i的k次幂
a[cnt]=pk;//模数
b[cnt]=C_mod(n,m,i,pk);//余数
//----------------------------
}
}
if(P>1)
{
cnt++;
//质因子P
//pk=P
a[cnt]=P;
b[cnt]=C_mod(n,m,P,P);
}
return CRT();
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
int main()
{
ll n,m,P;
cin>>n>>m>>P;
cout<<ex_Lucas(n,m,P)<<endl;
return 0;
}