组合数:
在N个数中选取M个数,问选的方式有几种?
直接递归暴力简单
#include<cstdio>
const int N = 2000 + 5;
const int MOD = (int)1e9 + 7;
int comb[N][N];//comb[n][m]就是C(n,m)
void init(){
for(int i = 0; i < N; i ++){
comb[i][0] = comb[i][i] = 1;
for(int j = 1; j < i; j ++){
comb[i][j] = comb[i-1][j] + comb[i-1][j-1];
comb[i][j] %= MOD;
}
}
}
int main(){
init();
}
卢卡斯定理:
求
C
m
n
m
o
d
p
C_m^n~mod~p
Cmn mod p
设
m
=
a
0
p
0
+
a
1
p
1
+
⋯
+
a
k
p
k
m={a_0}^{p_0}+{a_1}^{p_1}+\cdots+{a_k}^{p_k}
m=a0p0+a1p1+⋯+akpk
设
n
=
b
0
p
0
+
b
1
p
1
+
⋯
+
b
k
p
k
n={b_0}^{p_0}+{b_1}^{p_1}+\cdots+{b_k}^{p_k}
n=b0p0+b1p1+⋯+bkpk
则
C
m
n
≡
∏
C
a
i
b
i
(
m
o
d
p
)
C_m^n\equiv\prod{C_{a_i}^{b_i}}(mod~p)
Cmn≡∏Caibi(mod p)
证明我就不打了,百度百科上有,数学符号太多了~ _ ~
前置知识:
1.欧几里得或者费马小定理求逆元
2.快速幂
实现代码:
#include<iostream>
using namespace std;
typedef long long LL;
const LL N=1e5+2;
LL a[N];
void init(LL p)
{
a[1]=1;
for(int i=2;i<=p;++i)a[i]=a[i-1]*i%p;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
if(!b){
x=1;
y=0;
return;
}
exgcd(b,a%b,y,x);
y-=a/b*x;
}
LL ksm(LL x,LL n,LL mod)
{
LL ans=1;
while(n){
if(n&1)ans=ans*x%mod;
n>>=1;
x=x*x%mod;
}
return ans;
}
LL C(LL n,LL m,LL p)
{
if(n==m||m==0)return 1;
if(n<m)return 0;
if(m*2>n)m=n-m; /*C(n,m)=c(n,n-m)*/
return a[n]*ksm(a[m]*a[n-m],p-2,p)%p; /*求(a[m]*a[n-m])在(mod p)意义下的乘法逆元*/
/*拓展欧几里得与费马小定理均可*/
/*LL x,y;
exgcd(a[m]*a[n-m],p,x,y);
return (a[n]*x%p+p)%p;*/
}
LL lucas(LL n,LL m,LL p)
{
if(!m)return 1;
return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int main()
{
ios::sync_with_stdio(false);
LL T,n,m,p;
cin>>T;
while(T--){
cin>>n>>m>>p;
init(p);
cout<<lucas(n+m,m,p)<<endl;
}
return 0;
}
为了快,放了一个数组,限制了P的大小,然后再发一个没有限制的卢卡斯。
#define MOD 1000000007
typedef long long LL;
LL quickPower(LL a, LL b)
{
LL ans = 1;
a %= MOD;
while (b)
{
if (b & 1)
{
ans = ans * a % MOD;
}
b >>= 1;
a = a * a % MOD;
}
return ans;
}
LL c(LL n, LL m)
{
if (m > n)
{
return 0;
}
LL ans = 1;
for (int i = 1; i <= m; i++)
{
LL a = (n + i - m) % MOD;
LL b = i % MOD;
ans = ans * (a * quickPower(b, MOD - 2) % MOD) % MOD;
}
return ans;
}
LL lucas(LL n, LL m)
{
if (m == 0)
{
return 1;
}
return c(n % MOD, m % MOD) * lucas(n / MOD, m / MOD) % MOD;
}
int main(int argc, const char * argv[])
{
LL n, m;
while (~scanf("%lld %lld", &n, &m))
{
LL max, min;
max = n + m - 3;
min = m - 1;
printf("%lld\n", lucas(max - 1, min - 1));
}
return 0;
}
ExLucas扩展卢卡斯定理:
求解
C
n
m
%
P
C_{n}^{m}\% P
Cnm%P,其中m,n较大,P较小且不一定为素数
转化为CRT(中国剩余定理)
好像这也不是什么定理,只是一个计算方法
计
算
C
n
m
%
P
,
其
中
p
=
p
1
q
1
×
p
2
q
2
×
⋯
p
k
q
k
时
,
我
们
可
以
先
求
出
C
m
n
m
o
d
p
i
q
i
,
然
后
用
C
R
T
合
并
。
那
么
怎
么
计
算
C
m
n
m
o
d
p
i
q
i
呢
?
C
m
n
=
m
!
n
!
(
m
−
n
)
!
,
我
们
只
需
要
算
出
m
!
,
n
!
−
1
,
(
m
−
n
)
!
−
1
,
然
后
乘
在
一
起
。
注
:
n
!
可
能
在
模
p
i
q
i
的
意
义
下
没
有
逆
元
啊
,
那
这
就
是
错
的
了
其
实
这
里
求
得
不
是
逆
元
(
可
能
没
有
逆
元
)
,
求
出
来
的
是
a
×
p
i
b
(
g
c
d
(
a
,
p
)
=
1
)
,
前
面
的
a
a
用
逆
元
,
后
面
的
次
数
加
加
减
减
一
下
就
好
了
问
题
转
换
成
求
n
!
m
o
d
p
q
例
如
n
=
19
,
p
=
3
,
q
=
2
n
=
19
,
p
=
3
,
q
=
2
:
计算C_{n}^{m}\% P,其中p=p1^{q1}×p2^{q2}×⋯pk^{qk}时,我们可以先求出C_m^n~mod~{p_i}^{q_i},然后用CRT合并。\\ 那么怎么计算C_m^n~mod~{p_i}^{q_i}呢?\\ C_m^n=\frac{m!}{n!(m-n)!},我们只需要算出m!,n!^{−1},(m−n)!^{−1} ,然后乘在一起。\\ 注:n! 可能在模{p_i}^{q_i}的意义下没有逆元啊,那这就是错的了\\ 其实这里求得不是逆元(可能没有逆元),求出来的是a\times {p_i}^b(gcd(a,p)=1),\\前面的aa用逆元,后面的次数加加减减一下就好了\\ 问题转换成求n!~mod~p^q\\ 例如n=19,p=3,q=2n=19,p=3,q=2:
计算Cnm%P,其中p=p1q1×p2q2×⋯pkqk时,我们可以先求出Cmn mod piqi,然后用CRT合并。那么怎么计算Cmn mod piqi呢?Cmn=n!(m−n)!m!,我们只需要算出m!,n!−1,(m−n)!−1,然后乘在一起。注:n!可能在模piqi的意义下没有逆元啊,那这就是错的了其实这里求得不是逆元(可能没有逆元),求出来的是a×pib(gcd(a,p)=1),前面的aa用逆元,后面的次数加加减减一下就好了问题转换成求n! mod pq例如n=19,p=3,q=2n=19,p=3,q=2:
参考
19
!
=
1
×
2
×
3
×
⋯
×
19
=
(
1
×
2
×
4
×
5
×
7
×
8
⋯
×
16
×
17
×
19
)
×
(
3
×
6
×
9
×
12
×
15
×
18
)
=
(
1
×
2
×
4
×
5
×
7
×
8
⋯
×
16
×
17
)
×
19
×
36
×
(
1
×
2
×
3
×
4
×
5
×
6
)
=
(
1
×
2
×
4
×
5
×
7
×
8
)
2
×
19
×
36
×
(
1
×
2
×
3
×
4
×
5
×
6
)
19!\\ =1×2×3×⋯×19\\ =(1×2×4×5×7×8⋯×16×17×19)×(3×6×9×12×15×18)\\ =(1×2×4×5×7×8⋯×16×17)×19×36×(1×2×3×4×5×6)\\ =(1×2×4×5×7×8)2×19×36×(1×2×3×4×5×6)
19!=1×2×3×⋯×19=(1×2×4×5×7×8⋯×16×17×19)×(3×6×9×12×15×18)=(1×2×4×5×7×8⋯×16×17)×19×36×(1×2×3×4×5×6)=(1×2×4×5×7×8)2×19×36×(1×2×3×4×5×6)
上面这个式子分为四部分:
第
一
部
分
:
(
1
×
2
×
4
×
5
×
7
×
8
)
2
(
1
×
2
×
4
×
5
×
7
×
8
)
2
。
这
部
分
的
数
不
超
过
p
q
p
q
个
,
可
以
暴
力
算
第
二
部
分
:
1919
。
这
部
分
的
数
不
超
过
p
q
p
q
个
,
可
以
暴
力
算
第
三
部
分
:
3636
。
这
个
在
最
后
处
理
时
求
出
m
!
,
n
!
,
(
m
−
n
)
!
m
!
,
n
!
,
(
m
−
n
)
!
分
别
有
多
少
个
p
p
(
设
为
x
,
y
,
z
x
,
y
,
z
)
,
则
答
案
要
乘
上
p
x
−
y
−
z
p
x
−
y
−
z
第
四
部
分
:
1
×
2
×
3
×
4
×
5
×
61
×
2
×
3
×
4
×
5
×
6
。
这
个
是
⌊
n
p
⌋
!
⌊
n
p
⌋
!
,
可
以
递
归
处
理
第一部分:(1×2×4×5×7×8)2(1×2×4×5×7×8)2。这部分的数不超过pqpq个,可以暴力算\\ 第二部分:1919。这部分的数不超过pqpq个,可以暴力算\\ 第三部分:3636。这个在最后处理时求出m!,n!,(m−n)!m!,n!,(m−n)!\\ \qquad\ \ \ \ \ \ \ \ \ \ \ \ 分别有多少个pp(设为x,y,zx,y,z),则答案要乘上px−y−zpx−y−z\\ 第四部分:1×2×3×4×5×61×2×3×4×5×6。\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 这个是⌊np⌋!⌊np⌋!,可以递归处理
第一部分:(1×2×4×5×7×8)2(1×2×4×5×7×8)2。这部分的数不超过pqpq个,可以暴力算第二部分:1919。这部分的数不超过pqpq个,可以暴力算第三部分:3636。这个在最后处理时求出m!,n!,(m−n)!m!,n!,(m−n)! 分别有多少个pp(设为x,y,zx,y,z),则答案要乘上px−y−zpx−y−z第四部分:1×2×3×4×5×61×2×3×4×5×6。 这个是⌊np⌋!⌊np⌋!,可以递归处理
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N=1e5+9;
LL A[N],M[N];
LL ksm(LL x,LL n,LL mod)
{
LL ans=1;
while(n){
if(n&1)ans=ans*x%mod;
n>>=1,x=x*x%mod;
}
return ans;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
if(!b)x=1,y=0;
else exgcd(b,a%b,y,x),y-=a/b*x;
}
LL inv(LL a,LL p)
{
LL x,y;
exgcd(a,p,x,y);
return (x+p)%p?x:x+p;
}
LL get(LL n,LL pi,LL p) /*求(与pi互素后的n!)%M[i]*/
{
if(!n)return 1;
LL ans=1;
if(n/p){ /*判断有无循环节 */
for(LL i=2;i<=p;++i)if(i%pi)ans=ans*i%p;
ans=ksm(ans,n/p,p);
}
for(LL i=2;i<=n%p;++i)if(i%pi)ans=ans*i%p; /*循环节剩余部分*/
return ans*get(n/pi,pi,p)%p;
}
LL exlucas(LL n,LL m,LL pi,LL p) /*求A[i]*/
{
LL nn=get(n,pi,p); /*求(与pi互素后的n)%M[i]*/
LL mm=get(m,pi,p); /*求(m!与pi互素后的m!)%M[i]*/
LL nm=get(n-m,pi,p); /*求(与pi互素后的(n-m)!)%M[i]*/
LL k=0; /*含质因数pi的数量*/
for(LL i=n;i;i/=pi)k+=i/pi;
for(LL i=m;i;i/=pi)k-=i/pi;
for(LL i=n-m;i;i/=pi)k-=i/pi;
return nn*inv(mm,p)*inv(nm,p)*ksm(pi,k,p)%p;
}
LL crt(LL len,LL Lcm)
{
LL ans=0;
for(LL i=1;i<=len;++i){
LL Mi=Lcm/M[i];
ans=((ans+A[i]*inv(Mi,M[i])*Mi)%Lcm+Lcm)%Lcm;
}
return ans;
}
int main()
{
ios::sync_with_stdio(false);
LL n,m,P,num;
while(cin>>n>>m>>P){
if(n<m){
cout<<0<<endl;
continue;
}
num=0;
memset(A,0,sizeof(A));
memset(M,0,sizeof(M));
for(LL x=P,i=2;i<=P;++i)
if(x%i==0){
M[++num]=1;
while(x%i==0){
M[num]*=i;
x/=i;
}
A[num]=exlucas(n,m,i,M[num])%P;
}
cout<<crt(num,P)<<endl;
}
return 0;
}