组合数取模运算
参考博客：
大数量级组合数的快速计算方法,log转化
组合数取模
组合数取模2
组合公式
c(n,m)=p(n,m)/m!=n!/((n-m)!*m!)
c(n,m)=c(n,n-m)
c(n,m)=c(n-1,m)+c(n-1,m-1)
欧拉定理 （<---自家博客）
欧拉定理,（也称费马-欧拉定理）是一个关于同余的性质。欧拉定理表明，若n,a为正整数，且n,a互质，则:
φ(n)表示1~n中与n互质的数的个数
看一个基本的例子。令a = 3，n = 5，这两个数是互素的。比5小的正整数中与5互素的数有1、2、3和4，所以φ(5)=4（详情见[欧拉函数]）。计算:a^{φ(n)} = 3^4 =81，而81= 80 + 1 Ξ 1 （mod 5）
这个定理可以用来简化幂的模运算。比如计算7^{222}的个位数，实际是求7^{222}被10除的余数。7和10[[互素]]，且φ(10)=4。由欧拉定理知7^4Ξ1(mod 10)。所以7^{222}=(7^4)^55*(7^2)Ξ1^{55}*7^2Ξ49Ξ9 (mod 10)。
后话：欧拉定理其实是费马小定理的推广
费马小定理：
对于质数p，任意整数a，均满足：ap≡a（mod p）
#欧拉定理的推论：
若正整数a，n互质，那么对于任意正整数b，有ab≡ab mod φ（n）（mod n）
ps: 这个推论可以帮助我们在求幂运算的时候缩小数据范围和计算次数。具体的说：在求乘方运算时，可以先把底数对mod取模，再把指数对b mod φ（n）取模。
特别的，如果a，mod不互质，且b>φ（n）时，ab≡ab mod φ（n）+ φ（n）（mod n）
乘法逆元 (目前不会)
介绍：如果ax≡1 (mod p),且gcd(a,p)=1（a与p互质），则称a关于模p的乘法逆元为x。
(a/b) mod p=(a*b^(p-2)) mod p
条件：p是素数，gcd(b,p)=1，a%b=0
定义： 满足a*k≡1 (mod p)的k值就是a关于p的乘法逆元。
(PS：p一定是个素数才能对a有乘法逆元（除1），特别注意：当p是1时，对于任意a，k都为1）
为什么要有乘法逆元呢？
当我们要求（a/b） mod p的值（a/b一定是个整数），且a很大，无法直接求得a/b的值时，我们就要用到乘法逆元。
我们可以通过求b关于p的乘法逆元k，将a乘上k再模p，即(a*k) mod p。
其结果与(a/b) mod p等价。
接下来进入正题：
Lucas定理针对该取值范围较大又不太大的情况 （Lucas定理是用来求 c(n,m) mod p，p为素数的值。）
定理描述
这样将组合数的求解分解为小问题的乘积，下面考虑计算C(ni, mi) %p.
已知C(n, m) mod p = n!/(m!(n - m)!) mod p。当我们要求（a/b）mod p的值，且b很大，无法直接求得a/b的值时，我们可以转而使用乘法逆元k，将a乘上k再模p，即(a*k) mod p。 其结果与(a/b) mod p等价。
ACM数论之旅10---大组合数-卢卡斯定理（在下卢卡斯，你是我的Master吗？(。-`ω´-) ）
卢卡斯说：
C(n, m) % p = C(n / p, m / p) * C(n%p, m%p) % p
对于C(n / p, m / p)，如果n / p 还是很大，可以递归下去，一直到世界的尽头
上代码：
ll Lucas(ll n, ll m, ll p){
return m ? Lucas(n/p, m/p, p) * C(n%p, m%p, p) % p : 1;
}
....................算了直接给模板吧（Lucas组合数取模运算函数）
1 ll fact(ll n, ll p){//n的阶乘求余p 2 ll ret = 1; 3 for (ll i = 1; i <= n ; i ++) ret = ret * i % p ; 4 return ret ; 5 } 6 void ex_gcd(ll a, ll b, ll &x, ll &y, ll &d){ 7 if (!b) {d = a, x = 1, y = 0;} 8 else{ 9 ex_gcd(b, a % b, y, x, d); 10 y -= x * (a / b); 11 } 12 } 13 ll inv(ll t, ll p){//如果不存在，返回-1 //求逆元 14 ll d, x, y; 15 ex_gcd(t, p, x, y, d); 16 return d == 1 ? (x % p + p) % p : -1; 17 } 18 ll comb(ll n, ll m, ll p){//C(n, m) % p 19 if (m < 0 || m > n) return 0; 20 return fact(n, p) * inv(fact(m, p), p) % p * inv(fact(n-m, p), p) % p; 21 } 22 23 ll Lucas(ll n, ll m, ll p){//Lucas定理求C(n, m) % p 24 return m ? Lucas(n/p, m/p, p) * comb(n%p, m%p, p) % p : 1; 25 }
关于拓展卢卡斯，也就是卢卡斯的升级版，卢卡斯限定只能取余一个素数，而拓展卢卡斯则没有限定
其求解方法如下：
若不是素数，将p分解质因数，将C(n,m)分别按照Lucas定理中的方法求对p的质因数的模，然后用中国剩余定理合并。
上模板：
1 ll exgcd(ll a,ll b,ll &x,ll &y) 2 { 3 if(!b){x=1;y=0;return a;} 4 ll res=exgcd(b,a%b,x,y),t; 5 t=x;x=y;y=t-a/b*y; 6 return res; 7 } 8 9 inline ll power(ll a,ll b,ll mod) 10 { 11 ll sm; 12 for(sm=1;b;b>>=1,a=a*a%mod)if(b&1) 13 sm=sm*a%mod; 14 return sm; 15 } 16 17 ll fac(ll n,ll pi,ll pk) 18 { 19 if(!n)return 1; 20 ll res=1; 21 for(register ll i=2;i<=pk;++i) 22 if(i%pi)(res*=i)%=pk; 23 res=power(res,n/pk,pk); 24 for(register ll i=2;i<=n%pk;++i) 25 if(i%pi)(res*=i)%=pk; 26 return res*fac(n/pi,pi,pk)%pk; 27 } 28 29 inline ll inv(ll n,ll mod) 30 { 31 ll x,y; 32 exgcd(n,mod,x,y); 33 return (x+=mod)>mod?x-mod:x; 34 } 35 36 inline ll CRT(ll b,ll mod,ll p){return b*inv(p/mod,mod)%p*(p/mod)%p;} 37 38 const int MAXN=11; 39 40 static ll w[MAXN]; 41 42 inline ll C(ll n,ll m,ll pi,ll pk) 43 { 44 ll up=fac(n,pi,pk),d1=fac(m,pi,pk),d2=fac(n-m,pi,pk); 45 ll k=0; 46 for(register ll i=n;i;i/=pi)k+=i/pi; 47 for(register ll i=m;i;i/=pi)k-=i/pi; 48 for(register ll i=n-m;i;i/=pi)k-=i/pi; 49 return up*inv(d1,pk)%pk*inv(d2,pk)%pk*power(pi,k,pk)%pk; 50 } 51 52 inline ll exlucus(ll n,ll m,ll p)//求C(n, m) % p 53 { 54 ll res=0,tmp=p,pk; 55 static int lim=sqrt(p)+5; 56 for(register int i=2;i<=lim;++i)if(tmp%i==0) 57 { 58 pk=1;while(tmp%i==0)pk*=i,tmp/=i; 59 (res+=CRT(C(n,m,i,pk),pk,p))%=p; 60 } 61 if(tmp>1)(res+=CRT(C(n,m,tmp,tmp),tmp,p))%=p; 62 return res; 63 }
附加：快速乘（如果直接乘会爆ll）
ll mul(ll a, ll b, ll p){//快速乘，计算a*b%p
ll ret = 0;
while(b){
if(b & 1) ret = (ret + a) % p;
a = (a + a) % p;
b >>= 1;
}
return ret;
}
实战一下吧
http://acm.hdu.edu.cn/showproblem.php?pid=5446
题意：
给你三个数n, m, k
第二行是k个数，p1,p2,p3...pk
所有p的值不相同且p都是质数
求C(n, m) % (p1*p2*p3*...*pk)的值
范围：1≤m≤n≤1e18，1≤k≤10，pi≤1e5，保证p1*p2*p3*...*pk≤1e18
AC代码：
1 #include<cstdio> 2 typedef long long LL; 3 const int N = 100000 + 5; 4 LL mul(LL a, LL b, LL p){//快速乘，计算a*b%p 5 LL ret = 0; 6 while(b){ 7 if(b & 1) ret = (ret + a) % p; 8 a = (a + a) % p; 9 b >>= 1; 10 } 11 return ret; 12 } 13 LL fact(int n, LL p){//n的阶乘求余p 14 LL ret = 1; 15 for (int i = 1; i <= n ; i ++) ret = ret * i % p ; 16 return ret ; 17 } 18 void ex_gcd(LL a, LL b, LL &x, LL &y, LL &d){ 19 if (!b) {d = a, x = 1, y = 0;} 20 else{ 21 ex_gcd(b, a % b, y, x, d); 22 y -= x * (a / b); 23 } 24 } 25 LL inv(LL t, LL p){//如果不存在，返回-1 26 LL d, x, y; 27 ex_gcd(t, p, x, y, d); 28 return d == 1 ? (x % p + p) % p : -1; 29 } 30 LL comb(int n, int m, LL p){//C(n, m) % p 31 if (m < 0 || m > n) return 0; 32 return fact(n, p) * inv(fact(m, p), p) % p * inv(fact(n-m, p), p) % p; 33 } 34 LL Lucas(LL n, LL m, int p){ 35 return m ? Lucas(n/p, m/p, p) * comb(n%p, m%p, p) % p : 1; 36 } 37 LL china(int n, LL *a, LL *m){//中国剩余定理 38 LL M = 1, ret = 0; 39 for(int i = 0; i < n; i ++) M *= m[i]; 40 for(int i = 0; i < n; i ++){ 41 LL w = M / m[i]; 42 //ret = (ret + w * inv(w, m[i]) * a[i]) % M;//这句写了会WA，用下面那句 43 ret = (ret + mul(w * inv(w, m[i]), a[i], M)) % M; 44 //因为这里直接乘会爆long long ,所以我用快速乘(unsigned long long也是爆掉，除非用高精度) 45 } 46 return (ret + M) % M; 47 } 48 int main(){ 49 int T, k; 50 LL n, m, p[15], r[15]; 51 scanf("%d", &T); 52 while(T--){ 53 scanf("%I64d%I64d%d", &n, &m, &k); 54 for(int i = 0; i < k; i ++){ 55 scanf("%I64d", &p[i]); 56 r[i] = Lucas(n, m, p[i]); 57 } 58 printf("%I64d\n", china(k, r, p)); 59 } 60 }
我们知道题目要求C(n, m) % (p1*p2*p3*...*pk)的值
其实这个就是中国剩余定理最后算出结果后的最后一步求余
那C(n, m)相当于以前我们需要用中国剩余定理求的值
然而C(n, m)太大，我们只好先算出
C(n, m) % p1 = r1
C(n, m) % p2 = r2
C(n, m) % p3 = r3
.
.
.
C(n, m) % pk = rk
用Lucas，这些r1,r2,r3...rk可以算出来
然后又是用中国剩余定理求答案
附上大佬博客->>>>dalao
