Lucas定理
求解 C n m m o d    p C_n^m\mod p Cnmmodp (p为素数)
蒟蒻不会证明,记着递推公式就好=_=
L
u
c
a
s
(
n
,
m
)
=
C
(
n
,
m
)
 
m
o
d
 
p
Lucas(n,m)=C(n,m)\bmod p
Lucas(n,m)=C(n,m)modp
=
C
(
n
 
m
o
d
 
p
=C(n\bmod p
=C(nmodp
,
m
 
m
o
d
 
p
,m\bmod p
,mmodp
)
∗
L
u
c
a
s
(
n
/
p
,
m
/
p
)
 
m
o
d
 
p
)*Lucas(n/p,m/p)\bmod p
)∗Lucas(n/p,m/p)modp
边界为
L
u
c
a
s
(
i
,
0
)
=
1
Lucas(i,0)=1
Lucas(i,0)=1
int C(int n,int m)
{
int ans=1;
for(int i=1;i<=m;++i)
ans*=(n-m+i)*inv[i]%p;
return ans;
}
int lucas(int n,int m)
{
if(m==0) return 1;
else return C(n%p,m%p)*lucas(n/p,m/p)%p;
}
扩展Lucas
求解 C n m m o d    p C_n^m\mod p Cnmmodp (p不一定为素数)
先导知识 中国剩余定理、扩欧、逆元、组合数、质因数分解
首先讲p分解质因数
p
=
p
1
k
1
+
p
2
k
2
+
.
.
.
+
p
q
k
q
p=p_1^{k_1}+p_2^{k_2}+...+p_q^{k_q}
p=p1k1+p2k2+...+pqkq
于是可以列出同余方程组
{
a
n
s
≡
C
n
m
m
o
d
  
p
1
k
1
(
m
o
d
  
p
1
k
1
)
a
n
s
≡
C
n
m
m
o
d
  
p
2
k
2
(
m
o
d
  
p
2
k
2
)
.
.
.
a
n
s
≡
C
n
m
m
o
d
  
p
q
k
q
(
m
o
d
  
p
q
k
q
)
\left\{\begin{aligned} ans\equiv C_n^m\mod p_1^{k_1}(\mod p_1^{k_1})\\ ans\equiv C_n^m\mod p_2^{k_2}(\mod p_2^{k_2})\\ ...\\ ans\equiv C_n^m\mod p_q^{kq}(\mod p_q^{kq}) \end{aligned}\right.
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧ans≡Cnmmodp1k1(modp1k1)ans≡Cnmmodp2k2(modp2k2)...ans≡Cnmmodpqkq(modpqkq)
显然ans就是我们所要的答案
因为
p
i
k
i
p_i^{k_i}
piki是分解质因数得到的,所以必定两两互质
只要用中国剩余定理合并就好
我们要讨论的是怎样求出
C
n
m
m
o
d
  
p
i
k
i
C_n^m\mod p_i^{k_i}
Cnmmodpiki
直接考虑从定义式入手
C
n
m
m
o
d
  
p
i
k
i
=
n
!
m
o
d
  
p
i
k
i
m
!
m
o
d
  
p
i
k
i
∗
(
n
−
m
)
!
m
o
d
  
p
i
k
i
C_n^m\mod p_i^{k_i}=\frac{n!\mod p_i^{k_i}}{m!\mod p_i^{k_i}*(n-m)!\mod p_i^{k_i}}
Cnmmodpiki=m!modpiki∗(n−m)!modpikin!modpiki
直接看例子会比较容易理解,看这个例子 n = 19 , p i = 3 , k i = 2 n=19,pi=3,ki=2 n=19,pi=3,ki=2
先试着把
n
!
n!
n!中所有
p
i
p_i
pi的倍数分离出来
得
3
∗
6
∗
9
∗
12
∗
15
∗
18
3*6*9*12*15*18
3∗6∗9∗12∗15∗18
再把这里面的
p
i
p_i
pi全部提出来
得
3
6
∗
(
1
∗
2
∗
3
∗
4
∗
5
∗
6
)
=
3
6
∗
6
!
3^6*(1*2*3*4*5*6)=3^6*6!
36∗(1∗2∗3∗4∗5∗6)=36∗6!
形式化的描述这部分就是
p
i
⌊
n
/
p
i
⌋
∗
⌊
n
/
p
i
⌋
!
p_i^{\lfloor n/p_i \rfloor}*\lfloor n/p_i \rfloor!
pi⌊n/pi⌋∗⌊n/pi⌋!
p
i
p_i
pi得幂次方部分可以直接快速幂,阶乘部分可以继续递归
处理完上面来看剩下一部分
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*19
1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19
仔细观察可以发现这一串在膜
p
i
k
i
pi^{ki}
piki意义下是有周期的
(
1
∗
2
∗
4
∗
5
∗
7
∗
8
)
≡
(
10
∗
11
∗
13
∗
14
∗
16
∗
17
)
(
m
o
d
  
p
i
k
i
)
(1*2*4*5*7*8)\equiv(10*11*13*14*16*17)(\mod pi^{ki})
(1∗2∗4∗5∗7∗8)≡(10∗11∗13∗14∗16∗17)(modpiki)
每一个循环节长度为
p
i
k
i
pi^{ki}
piki
只要求出第一个循环节,然后在快速幂
⌊
n
/
p
i
k
i
⌋
\lfloor n/pi^{k_i} \rfloor
⌊n/piki⌋次方就可以了
最后发现还剩下一个19
也就是说我们算完上面两部分还有剩余
n
−
⌊
n
/
p
i
k
i
⌋
∗
p
i
k
i
n-\lfloor n/pi^{k_i} \rfloor*pi^{k_i}
n−⌊n/piki⌋∗piki个数字
不过这些数字长度显然不超过
p
i
k
i
pi^{ki}
piki,直接暴力就好
最后再次总结一下分解的几部分
1.一个幂次方以及一个阶乘, p i ⌊ n / p i ⌋ ∗ ⌊ n / p i ⌋ ! p_i^{\lfloor n/p_i \rfloor}*\lfloor n/p_i \rfloor! pi⌊n/pi⌋∗⌊n/pi⌋!
2.一个循环节的 ⌊ n / p k ⌋ \lfloor n/p_k \rfloor ⌊n/pk⌋次方
3.剩下的 n − ⌊ n / p i k i ⌋ ∗ p i k i n-\lfloor n/pi^{k_i} \rfloor*pi^{k_i} n−⌊n/piki⌋∗piki个数字暴力乘
洛谷P4720 【模板】扩展卢卡斯
#include<iostream>
#include<cmath>
#include<algorithm>
#include<queue>
#include<cstring>
#include<cstdio>
using namespace std;
typedef long long lt;
lt read()
{
lt f=1,x=0;
char ss=getchar();
while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
return f*x;
}
const int maxn=50001;
lt a[maxn],b[maxn],cnt;
lt qpow(lt ai,lt k,lt mod)
{
lt mul=1;
while(k>0)
{
if(k&1)mul=(mul*ai)%mod;
ai=(ai*ai)%mod;
k>>=1;
}
return mul;
}
lt fac(lt n,lt pi,lt pk)
{
if(!n) return 1;
lt mul=1;
for(lt i=2;i<=pk;++i)//分解阶乘第二部分,循环节
if(i%pi)mul=(mul*i)%pk;
mul=qpow(mul,n/pk,pk);
for(lt i=2;i<=n%pk;++i)//分解阶乘第三部分,求剩余数字
if(i%pi)mul=(mul*i)%pk;
return mul*fac(n/pi,pi,pk)%pk;//分解阶乘第一部分的另一个阶乘递归
}
void exgcd(lt a,lt b,lt &x,lt &y)
{
if(b==0){ x=1; y=0; return;}
exgcd(b,a%b,x,y);
lt tp=x;
x=y; y=tp-a/b*y;
}
lt inv(lt a,lt b)
{
lt x,y;
exgcd(a,b,x,y);
return (x%b+b)%b;
}
lt C(lt n,lt m,lt pi,lt pk)
{
lt facn=fac(n,pi,pk);//分别求n,m,n-m膜pi^ki的阶乘
lt facm=fac(m,pi,pk);
lt facnm=fac(n-m,pi,pk);
lt kk=0;
for(lt i=n;i;i/=pi)kk+=i/pi;//上述分解阶乘第一部分的pi幂次方
for(lt i=m;i;i/=pi)kk-=i/pi;
for(lt i=n-m;i;i/=pi)kk-=i/pi;
return facn*inv(facm,pk)%pk*inv(facnm,pk)%pk*qpow(pi,kk,pk)%pk;//注意求逆元
}
void div(lt n,lt m,lt x)
{
for(lt i=2;i<=sqrt(x);++i)
{
if(x%i==0)
{
lt pi=i,ki=0;
while(x%i==0)x/=i,ki++;
b[++cnt]=qpow(pi,ki,1e7);
a[cnt]=C(n,m,pi,b[cnt]);//C(n,m)%pi^ki
}
}
if(x>1)b[++cnt]=x,a[cnt]=C(n,m,x,b[cnt]);
}
lt china()
{
lt ans=0,M=1,x,y;
for(int i=1;i<=cnt;++i) M*=b[i];
for(int i=1;i<=cnt;++i)
{
lt tp=M/b[i];
exgcd(tp,b[i],x,y);
x=(x%b[i]+b[i])%b[i];
ans=(ans+tp*x*a[i])%M;
}
return (ans+M)%M;
}
int main()
{
lt n=read(),m=read(),p=read();
div(n,m,p);
printf("%lld",china());
return 0;
}