组合数学的一点应用
C(n,m)的含义是有n个小朋友妹子,要选出m个发奖品泡妹子。(无序)
一个常用的公式C(n,m)=(n!)/(m!*(n-m)!);
另一个常用公式C(n,m)=C(n-1,m-1)+C(n-1,m);
证明:考虑最后一个去不去,或代数证明。
这两个公式就对应了两种方法。
但它们并不是很高效,我们就要用Lucas定理。一个非常重要的东东就是要p一定要质数。
定理:
对所有mi,ni都小于p。
等等,p不是质数咋办?
这是就可以用扩展lucas了。
http://blog.csdn.net/clove_unique/article/details/54571216
模板
1 #include<algorithm> 2 #include<iostream> 3 #include<cstring> 4 #include<cstdio> 5 #include<cmath> 6 using namespace std; 7 #define LL long long 8 9 LL n,m,MOD,ans; 10 11 LL fast_pow(LL a,LL p,LL Mod)//快速幂 12 { 13 LL ans=1LL; 14 for (;p;p>>=1,a=a*a%Mod) 15 if (p&1) 16 ans=ans*a%Mod; 17 return ans; 18 } 19 void exgcd(LL a,LL b,LL &x,LL &y)//扩欧 20 { 21 if (!b) x=1LL,y=0LL; 22 else exgcd(b,a%b,y,x),y-=a/b*x; 23 } 24 LL inv(LL A,LL Mod)//A的逆元 25 { 26 if (!A) return 0LL; 27 LL a=A,b=Mod,x=0LL,y=0LL; 28 exgcd(a,b,x,y); 29 x=((x%b)+b)%b; 30 if (!x) x+=b; 31 return x; 32 } 33 LL Mul(LL n,LL pi,LL pk) 34 { 35 if (!n) return 1LL; 36 LL ans=1LL; 37 if (n/pk) 38 { 39 for (LL i=2;i<=pk;++i) 40 if (i%pi) ans=ans*i%pk; 41 ans=fast_pow(ans,n/pk,pk);//另一部分的值。 42 } 43 for (LL i=2;i<=n%pk;++i) 44 if (i%pi) ans=ans*i%pk;//一部分的值。 45 return ans*Mul(n/pi,pi,pk)%pk; 46 } 47 LL C(LL n,LL m,LL Mod,LL pi,LL pk) 48 { 49 if (m>n) return 0LL; 50 LL a=Mul(n,pi,pk),b=Mul(m,pi,pk),c=Mul(n-m,pi,pk);//把pi除去阶乘的摸值 51 LL k=0LL,ans;//k表示pi还有多少个。 52 for (LL i=n;i;i/=pi) k+=i/pi; 53 for (LL i=m;i;i/=pi) k-=i/pi; 54 for (LL i=n-m;i;i/=pi) k-=i/pi; 55 ans=a*inv(b,pk)%pk*inv(c,pk)%pk*fast_pow(pi,k,pk)%pk; 56 return ans*(Mod/pk)%Mod*inv(Mod/pk,pk)%Mod;//中国剩余定理。 57 } 58 int main() 59 { 60 scanf("%I64d%I64d%I64d",&n,&m,&MOD); 61 for (LL x=MOD,i=2;i<=MOD;++i) 62 if (x%i==0) 63 { 64 LL pk=1LL; 65 while (x%i==0) pk*=i,x/=i; 66 ans=(ans+C(n,m,MOD,i,pk))%MOD;//中国剩余定理。 67 } 68 printf("%I64d\n",ans); 69 }
我们再用最前面的公式解决mi,ni,只要预处理出这阶乘,和阶乘的逆就行了。
首先补充个线性求逆。
有:
k=(p/i)(向下取整),r=p%i,k*i+r=p=0 (mod p);
同乘i-1×r-1。有i-1=-k*r-1。
就有递推式:A[i] = -(p / i) * A[p % i];
其他就没什么问题了。
代码:
1 void Init() 2 { 3 int i; 4 fac[0]=1; 5 for(i=1;i<=n&&i<p;i++) 6 fac[i]=fac[i-1]*i%p; 7 inv[1]=1; 8 for(i=2;i<=n&&i<p;i++) 9 inv[i]=(p-p/i)*inv[p%i]%p; 10 inv[0]=1; 11 for(i=1;i<=n&&i<p;i++) 12 inv[i]=inv[i]*inv[i-1]%p; 13 } 14 ll C(ll x,ll y) 15 { 16 if(x<y) return 0; 17 if(x<p&&y<p) 18 return fac[x]*inv[y]%p*inv[x-y]%p; 19 return C(x/p,y/p)*C(x%p,y%p)%p; 20 }
一个公式: C(n,m)+C(n+1,m)+C(n+2,m)+...+C(n+p,m)=C(n+p+1,m+1)-C(n,m+1);
C(n,m)+C(n+1,m)+C(n+2,m)+...+C(n+p,m)
=-C(n,m+1)+C(n,m+1)+C(n,m)+C(n+1,m)+C(n+2,m)+...+C(n+p,m)
=-C(n,m+1)+C(n+1,m+1)+C(n+1,m)+C(n+2,m)+...+C(n+p,m)
=-C(n,m+1)+C(n+2,m+1)+C(n+2,m)+...+C(n+p,m)
=-C(n,m+1)+C(n+p,m+1)+C(n+p,m)
=C(n+p+1,m+1)-C(n,m+1)