Time Limit: 1000MS | Memory Limit: 32768KB | 64bit IO Format: %I64d & %I64u |
Description
Input
Output
Sample Input
2 5 2 3 5 2 61
Sample Output
1 10
组合数取模就是求的值,根据,和的取值范围不同,采取的方法也不一样。
下面,我们来看常见的两种取值情况(m、n在64位整数型范围内)
(1) ,
此时较简单,在O(n2)可承受的情况下组合数的计算可以直接用杨辉三角递推,边做加法边取模。
(2) , ,并且是素数
本文针对该取值范围较大又不太大的情况(2)进行讨论。
这个问题可以使用Lucas定理,定理描述:
其中
这样将组合数的求解分解为小问题的乘积,下面考虑计算C(ni, mi) %p.
已知C(n, m) mod p = n!/(m!(n - m)!) mod p。当我们要求(a/b)mod p的值,且a很大,无法直接求得a/b的值时,我们可以转而使用乘法逆元k,将a乘上k再模p,即(a*k) mod p。 其结果与(a/b) mod p等价。
那么逆元是什么?
<span style="font-size: 18px;">定义:满足a*k≡1 (mod p)的k值就是a关于p的乘法<strong>逆元</strong>(当p是1时,对于任意a,k都为1)</span>
除法取模,这里要用到m!(n-m)!的逆元。
根据费马小定理:
已知gcd(a, p) = 1,则 ap-1 ≡ 1 (mod p), 所以 a*ap-2 ≡ 1 (mod p)。
也就是 (m!(n-m)!)的逆元为 (m!(n-m)!)p-2 ;
上面的代码中用到了求幂取模操作来计算(m!(n-m)!)p-2 % p.下面解释幂取模算法:
反复平方法 求ab%m
通过研究指数b的二进制表示发现,对任意的整数b都可表示为:
- n表示b的实际二进制位数
- bi表示该位是0或1
因此,ab可表示为:
即用b的每一位表示a的每一项,而对任意相邻的两项存在平方关系,即:
因此我们构造下面的算法:
- 把b转换为二进制表示,并从右至左扫描其每一位(从低到高)
- 当扫描到第i位时,根据同余性质(2)计算a的第i项的模:
base变量表示第i-1位时计算出的模,通过递归能很容易地确定所有位的模。 - 如果第i位为1,即bi=1,则表示该位需要参与模运算,计算结果 result = (result*base) mod m;其中result为前i-1次的计算结果;若bi=0,则表示a的第i项为1,不必参与模运算
int mod(int a,int b,int m){
int result = 1;
int base = a;
while(b>0){
if(b & 1==1){
result = (result*base) % m;
}
base = (base*base) %m;
b >>=1;
}
return result;
}
其中运用了两个同余性质:
同余性质1:ab≡bc (mod m)
同余性质2: a≡c (mod m) => a2≡c2 (mod m)
理解要点:
- base记录了a的每项的模,无论b在该位是0还是1,该结果都记录,目的是给后续位为1的项使用,计算方式是前一结果的平方取模,这也是反复平方法的由来
- result只记录了位为1的项的模结果,该计算方式使用了同余性质1
- 通过地把a使用二进制表示,并结合同余性质1,2,巧妙地化解了大数取模的运算。对1024位这样的大数,也最多进行1024次循环便可计算模值,性能非常快。
该方法是许多西方数学家努力的结果,通常也称为Montgomery算法。
对于C(n, m) mod p。这里的n,m,p(p为素数)都很大的情况。就不能再用C(n, m) = C(n - 1,m) + C(n - 1, m - 1)的公式递推了。
这里用到Lusac定理
For non-negative integers m and n and a prime p, the following congruence relation holds:
where
and
are the base p expansions of m and n respectively.
对于单独的C(ni, mi) mod p,已知C(n, m) mod p = n!/(m!(n - m)!) mod p。显然除法取模,这里要用到m!(n-m)!的逆元。
根据费马小定理:
已知(a, p) = 1,则 ap-1 ≡ 1 (mod p), 所以 a*ap-2 ≡ 1 (mod p)。
也就是 (m!(n-m)!)的逆元为 (m!(n-m)!)p-2 ;
#include <iostream> #include <cstdio> #include <cstring> #include <cmath> using namespace std; const int N = 2000100; typedef long long LL; LL quick(LL n,LL m,LL p); LL lucas(LL n,LL m,LL p); LL power(LL n,LL m,LL p); int main() { int t; scanf("%d", &t); while(t--) { LL n, m, p; scanf("%lld %lld %lld",&n, &m, &p); printf("%lld\n",lucas(n, m, p)); } return 0; } LL lucas(LL n,LL m,LL p) { LL ans=1; while(m&&n&&ans) { ans=((ans%p)*(power(n,m,p)%p))%p; m/=p; n/=p; } return ans; } LL power(LL n,LL m,LL p) { LL sum1=1, sum2=1; for(int i=1;i<=m;i++) { sum1=((sum1%p)*(n-i+1)%p)%p; sum2=((sum2%p)*(i%p))%p; } LL ans=((sum1%p)*(quick(sum2,p-2,p)))%p; return ans; } LL quick(LL n,LL m,LL p) { LL r=1; while(m!=0) { if(m&1) { r=((r%p)*(n%p))%p; } n=n*n%p; m>>=1; } return r%p; }