对于组合数C(n,m),意为在n个球中任取m个的方案数,
当n < m时C(n,m)=0,
当n>=m时,
C(n,m)=n!m!(n−m)!
组合数的递推式,通常用的一种:
C(n,m)=C(n−1,m−1)+C(n−1,m)
讨论第n个球取不取即可得出递推式。
通常计算组合数C(n,m)%P,有以下几种方法,
可以根据数据特性和范围来选择:
1.n,m较小时,直接用递推式/定义式计算
2.P为素数时,对于定义式分母里的m!(n−m)!
求它对于P的逆元,就可以直接计算了
3.P为素数且P较小时,用Lucas定理计算
Lucas定理(当P为质数):C(n,m)=C(n%P,m%P)∗C(n/P,m/P)
那么可以递归计算,当n<P且m<P,直接预处理出P内的阶乘逆元等
就可以求了,非常快……似乎是log级别的
4.P不是素数,那只能用扩展Lucas定理了
假设P=pa11∗pa22∗……∗pakk,pi为素数且pi互不相同
扩展Lucas定理求出C(n,m)%paii=ansi,由于paii两两互质,所以用中国剩余定理(孙子定理)合并即可。
这里引出孙子定理:
中国剩余定理对于w个同余方程n=xi(mod pi)当所有pi两两互质的时候,令M=p1∗p2∗……∗pw,mi=M/pi则分别求出mi对于pi的逆元ai,答案就是∑xi∗ai∗mi (%M)
于是问题变成C(n,m)%paii,pi是一个质数
这可以根据定义式来,计算阶乘%paii
阶乘%paii可以分组,比如19!%32(好像都是这个例子……)
19!=1∗2∗3∗4∗5∗6∗7∗8∗9∗10∗11∗12∗13∗14∗15∗16∗17∗18∗19
=(1∗2∗4∗5∗7∗8)∗36∗(1∗2∗4∗5∗7∗8)∗19
(%32)
发现中间括号内的部分是重复的,而且长度不超过pcii
19这种多余项个数不超过pcii
36是直接计算来的,即⌊npi⌋
当然先要把所有pi的因子提取出来,最后再乘上去
∗∗∗注意!如果ci=1,直接Lucas即可,扩lucas会慢到飞!∗∗∗
bzoj2142,扩展Lucas的模板题
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int
MAX=100000;
ll n,mod;
int pcnt,len;
ll prime[6000],p[10000],pk[10000],c[10000];
bool notprime[MAX];
void Get_Prime(){
notprime[1]=1,pcnt=0;
for (int i=2;i<=MAX;i++){
if (!notprime[i]) prime[++pcnt]=i;
for (int j=1;j<=pcnt;j++){
if (prime[j]*i>MAX) break;
notprime[prime[j]*i]=1;
if (!(i%prime[j])) break;
}
}
}
void FenJie(ll x){
len=0;int j=1;
while (x!=1){
if (x%prime[j]==0){
p[++len]=prime[j],c[len]=0;
pk[len]=1;
while (x%prime[j]==0)
c[len]++,x/=prime[j],pk[len]*=prime[j];
}
j++;
}
}
ll ksm(ll x,ll y,ll mod){
ll z=1LL;
while (y){
if (y&1LL) z=z*x%mod;
y>>=1LL;
x=x*x%mod;
}
return z;
}
ll exgcd(ll a,ll b,ll &x,ll &y){
if (!b){
x=1LL,y=0LL;
return a;
}
ll tt=exgcd(b,a%b,x,y),tmp=x;
x=y,y=tmp-a/b*y;
return tt;
}
ll inv(ll a,ll b){
ll t1,t2;
exgcd(a,b,t1,t2);
t1=(t1+b)%b;
return t1;
}
ll Fac(ll n,ll pi,ll pk){
if (!n) return 1LL;
ll t=1LL;
if (n/pk){
for (ll i=2;i<=pk;i++)
if (i%pi) t=t*i%pk;
t=ksm(t,n/pk,pk);
}
for (ll i=2;i<=n%pk;i++)
if (i%pi) t=t*i%pk;
return t*Fac(n/pi,pi,pk)%pk;
}
ll Lucas(ll n,ll m,ll pi,ll pk){
if (n<m) return 0LL;
ll x=Fac(n,pi,pk),y=Fac(m,pi,pk),z=Fac(n-m,pi,pk);
ll t=0LL;
for (ll a=n;a;a/=pi) t+=a/pi;
for (ll a=m;a;a/=pi) t-=a/pi;
for (ll a=n-m;a;a/=pi) t-=a/pi;
return x*ksm(pi,t,pk)%pk*inv(y,pk)%pk*inv(z,pk)%pk;
}
ll C(ll n,ll m){
ll t=0LL,t1;
for (int i=1;i<=len;i++){
t1=Lucas(n,m,p[i],pk[i]);
t1=t1*(mod/pk[i])%mod*inv(mod/pk[i],pk[i])%mod;
t=(t+t1)%mod;
}
return t;
}
int main(){
Get_Prime();
scanf("%lld",&mod);
FenJie(mod);int m;
scanf("%lld%d",&n,&m);
ll x,ans=1LL;
while (m--){
scanf("%lld",&x);
ans=(ans*C(n,x))%mod;
n-=x;
if (n<0) return puts("Impossible"),0;
}
printf("%lld\n",ans);
return 0;
}
那么组合数的各种问题就大致解决了,有2个小tricks值得记记:
1.n!中x因子的个数=⌊nx⌋+⌊nx2⌋+⌊nx3⌋……
2.卡特兰数:(n个数的出栈顺序数列统计)C(2n,n)n+1
或者
C(2n,n)−C(2n,n−1)
接下来是欧几里得:gcd(a,b)=gcd(b,a%b),
也就是“辗转相除法”
当a,b很大比如高精度,可以用“更相减损术”,
即如果都能/2就/2,不然辗转相除。
扩展欧几里得exgcd(a,b,x,y)当a和b互质的时候,exgcd(a,b,x,y),更改了x和y的值能够求出a关于b的数论逆元(x)也就是ax=1(%b)扩欧还可以用来解不定方程ax+by=c方法:首先求出w=gcd(a,b),如果c%w≠0,则此方程无解,不然a,b,c均除以w然后a和b就互质了,求出ax=1(%b)的最小x,再将x∗c就得到了x的一个基础解。整个解集是一个周期,x以b为周期,y以a为周期。
扩欧用处多多啊!得会得会。。
//用了long long
ll exgcd(ll a,ll b,ll &x,ll &y){
if (!b){
x=1LL,y=0LL;
return a;
}
ll tt=exgcd(b,a%b,x,y),tmp=x;
x=y,y=tmp-a/b*y;
return tt;
}
费马小定理:当p为素数的时候,ap−1=1 (mod p)欧拉定理:当gcd(a,p)=1,aϕ(p)=1 (mod p)扩展欧拉定理:也就是当gcd(a,p)≠1的时候也成立的定理当b>ϕ(p),ab=ab%phi(p)+phi(p) (mod p)注意p不停变成phi(p)是只用log(p)次的,这说明了x不停变成ax,(a必须固定)那么到了log(p)次左右就x=ax了。
扩展欧拉定理一个应用是bzoj4869
……应该不会考什么这种题了。。。
数论可能真没什么题了……快速幂啥的虚无。。