数学学习小记（二） 组合数求模：Lucas 定理 LightOJ 1067 + Hdu 3037

最近被数论搞得生活不能自理。。。。

电路终于考完了，一道非线性原件的题我算出电阻是负的。。。。。生死未卜。。。

五一假期还要忙数学建模。。。。

从以下博文中学习的，代码中有许多借鉴的地方也来自这里：

hdu - 4349 - Xiao Ming's Hope - 大大的Lucas定理 && 小小的乘法逆元 - Julyana_Lin_夜 - 博客频道 - CSDN.NET
http://blog.csdn.net/julyana_lin/article/details/7840491

插板法公式。_百度知道
http://zhidao.baidu.com/question/399163465.html

hdu3037 大组合数取模(Lucas定理) - tju_virus的专栏 - 博客频道 - CSDN.NET
http://blog.csdn.net/tju_virus/article/details/7843248

欧几里德辗转相除法 费马小定理  欧拉定理  扩展欧几里德算法简介_xiaxia_新浪博客
http://blog.sina.com.cn/s/blog_48e3f9cd010002uc.html

#include <cmath>#include <cstdio>#include <cstring>#include <iostream>#include <algorithm>using namespace std;class Combinations_Mod{private:__int64 factorial[1000010];public:__int64 Montgomery (__int64 s,__int64 index,__int64 mod)//蒙哥马利幂模算法  //快速幂//返回值(s^index)%mod{__int64 ans=1;  s%=mod;  while (index>=1){if ((index&1)==1)   //奇数ans=(ans*s)%mod;  index>>=1;  s=s*s%mod;  }return ans;}__int64 exp_mod (__int64 s,__int64 index,__int64 p){   // 快速幂递归if(index==0)return 1;if(index==1)return s%p;__int64 t=exp_mod(s,index/2,p);t=t*t%p;if (index&1)t=(t*s)%p;return t;}__int64 exGcd (__int64 a,__int64 b,__int64 &x,__int64 &y){//扩展欧几里得求乘法逆元 ax+by=Gcd(a,b)__int64 r,t;if (b==0){x=1;y=0;return a;}r=exGcd(b,a%b,x,y);t=x;x=y;y=t-a/b*y;return r;}__int64 cm (__int64 n,__int64 m,__int64 p)  //计算C(n,m)%p{/*C(a,b) = a! / (b! * (a-b)!) mod p其实就是求 (a!/(a-b)!) * (b!)^(p-2) mod p(上面这一步变换是根据费马小定理：假如p是质数，且a,p互质，那么a的（p-1）次方除以p的余数恒为1,那么a和a^(p-2)互为乘法逆元，则 b/a = (b * a^(p-2) ) mod p)*/if (m>n)return 0;__int64 a=1,b=1;while (m){a=(a*n)%p;b=(b*m)%p;m--;n--;}return a*Montgomery(b,p-2,p)%p;}__int64 cal (__int64 n,__int64 mod){__int64 x,y;exGcd(n,mod,x,y);return (x%mod+mod)%mod;} __int64 getC (__int64 n,__int64 m,__int64 mod)  //计算C(n,m)%p,适用于p较大且不变{return factorial[n]*cal(factorial[m],mod)*cal(factorial[n-m],mod)%mod;}__int64 Lucas (__int64 n,__int64 m,__int64 p)   //该算法不需要预处理{//组合数求摸:求C(n,m)%p, p为素数,最大为10^5if (m==0)return 1;return (Lucas(n/p,m/p,p)*cm(n%p,m%p,p))%p;}void Init (__int64 p){factorial[0]=1;for (int i=1;i<=p;i++)factorial[i]=factorial[i-1]*i%p;}__int64 Lucas_init (__int64 a,__int64 m,__int64 p)  //Lucas预处理版实现{//适用于p较小且变化__int64 ans = 1;Init (p);while (a && m){  __int64 aa = a%p;__int64 bb = m%p;if (aa < bb)return 0; //C(aa,bb) 表示 在aa里面取bb个,取法为0;//由于p是素数，所以 a / b % p ,b 对于 mod a 肯定存在逆元//这儿的求逆不可先处理ans = ans*factorial[aa]*exp_mod(factorial[bb]*factorial[aa-bb]%p,p-2,p)%p;a/=p;m/=p;}return ans;}}ob;/*m个相同的球放到n个盒子里 C(m+n-1,m)现在就需要求不大于m的，相当于对i = 0,1...,m对C(n+i-1,i)求和,根据公式C(n,k) = C(n-1,k)+C(n-1,k-1)得C(n-1,0)+C(n,1)+...+C(n+m-1,m)= C(n,0)+C(n,1)+C(n+1,2)+...+C(n+m-1,m)= C(n+m,m)现在就是要求C(n+m,m) % p,其中p是素数*/int main ()  //Hdu 3037{int T;__int64 n,m,p;scanf("%d",&T);for (int Cas=1;Cas<=T;Cas++){scanf("%I64d %I64d %I64d",&n,&m,&p);printf("%I64d\n",ob.Lucas (n+m,m,p));   //较快//printf("%I64d\n",ob.Lucas_init (n+m,m,p));}return 0;}*/int main ()  //LightOJ 1067{int T;__int64 n,m;scanf("%d",&T);ob.Init(1000003);for (int Cas=1;Cas<=T;Cas++){scanf("%lld%lld",&n,&m);printf("Case %d: %lld\n",Cas,ob.getC(n,m,1000003));}return 0;}/*LightOJ 1067的测试数据1100 50output :Case 1: 44000461000000 4999981000000 4999991000000 5000001000000 5000011000000 500002999983 99991OutputCase 1: 375003Case 2: 125000Case 3: 375001Case 4: 125000Case 5: 375003Case 6: 759191*/