最近做题遇到一个问题,具体题目不表,我把关键问题抽象一下:
求组合数
C(n1,m1)*C(n2,m2)*...*C(nk,mk) % T
(T是一个非常大的素数)
那么C(n,m)
怎么求呢?回忆一下高中数学,其实就是n!/(m!*(n-m)!)
。所以怎么算呢?天真的我先把所有的分子累乘起来,up = n! % T * n1! % T * ... nk! % T
,再把所有分母累乘起来down=m1!*(n1-m1)! % T * ... mk! * (nk-mk)! % T
,然后ans = up / down
令人绝望的是,N次提交全部出错,而且提交结果显示78/100
,啥意思呢?就是100组测试数据,我的程序过了前77个第78个出错了。这种情况一般是哪个边界case没有处理到,然后进入了疯狂debug模式,丝毫没有意识到其实是我公式完全用错了,有除法的取模根本不能这么算!最后无可奈何去看答案,才发现了除法这个坑!不过从中也确实学习到了很多。
乘法和加法是满足以下公式的:
- (a+b)%c = (a%c + b%c) % c
- (a*b)%c = (a%c * b%c) % c
减法比较特殊,不同语言的实现方式导致不同结果,不过在python中以下公式也是成立的:
- (a-b)%c = (a%c-b%c)%c
为什么说减法特殊呢?这涉及到具体取模运算的过程和程序实现时的选择。先想想我们自己怎么计算 18 % 5,读作18除以5的余数。小学就学过了,18÷5=3...3,先用18对5进行整除得到3,然后18-5*3=3,这样就得到余数。所以编程语言也帮我们做了这个转换:a%b ==> a-(a/b)*b
。ab同号还好,如果ab异号就有点麻烦了。简单说:
- 如果
-7/5=-1
,那么(a-b)%c
≠(a%c-b%c)%c
- 如果
-7/5=-2
,那么式子就是满足的
你可以自己把数字带进去验证一下。
在python中,-7//5==-2
,你可能觉得有点匪夷所思,这里有篇文章详细介绍了python为什么要这么设计。其实你可以仔细想想,负数整除一个正数是否有什么实际含义?我确实没有想到有什么真实的场景需要用-7/5=-1
来表达,换句话说,-7/5
如果等于-1将没有任何物理意义和数学意义,同时它会使得减法的同余公式不成立,因此-7/5=-2
将是一个更好的选择。当然这只是题外话,由于各种历史原因和backward compatible的考虑,C系语言中还是使用的是-7/5=-1
某种意义上,你也可以认为这是为什么那么多人喜欢用python来实现科学计算库的原因——对数学支持好
说完了减法我们这才进入到今天的主题——除法
简单地举个例子你就能很容易看出(a/b)%c
≠((a%c) / (b%c))%c
,比如(50/20)%5=2
,但是(50%5)/(20%5)
根本都没法算,因为分母为0!
那计算C(n,m)时我们除了用高精度计算就毫无办法了呢?不是的,下面就来讲讲如何计算。
在这之前我们先来学习一下非常重要的费马小定理。即,当 p为质数时,对于任何整数a都存在 (a^p)-a = kp
,其中k为整数。这个式子读作
当p为质数时,a的p次方减去a一定是p的整数倍
我们可以验证一下,比如a=2 p=7,2^7-2=126=18*7
再推广一下,如果a p互质,那么存在a^(p-1) = kp+1
,读作:
当p为质数且ap互质时,a的p-1次方等于整数倍p加上1
依然可以验证下,a=2 p=7,ap也是互质的,2^(7-1)=64=9*7+1。如果我们把式子a^(p-1)=kp+1
两边同时对p取模,很容易发现结果等于1。也就是说:
a的p-1次方对p取模结果恒等于1(ap互质且p为质数)
以上便是费马小定理。
接下来引入一个逆元的概念,当方程 ax % p ≡ 1有唯一解,那么我们称这个解x是a的逆元(≡读作恒等于)。根据费马小定理,我们知道a^(p-1) % p ≡ 1
,把 a^(p-1)
写作 a*a^(p-2)
,我们很容易就看出,a的逆元是a^(p-2)%p
。
有了逆元的概念,我们再来看最开始的那个问题 (a/b)%c
,我们可以把式子等价变换一下:
-
(a/b) % c = a*1*b^(-1) % c
---> (1) - 由费马小定理知
b^(c-1) % c ≡ 1
- 把(1)中的1替换成
b^(c-1)%c
,则a*1*b^(-1) %c = a*b^(c-1)%c *b^(-1) %c
- 由于乘法是满足同余的,因此上式又可以写成
(a*(b^(c-2))%c)%c
到这里我们发现(b^(c-2))%c
其实就是b关于c的逆元,换句话说,(a/b)%c = a乘以(b关于c的逆元)再 %c。
到这里我们已经能够通过转换计算除法的取模运算了,那么怎么快速计算出b关于c的逆元呢?也就是怎么快速计算b^(c-2) % c
。BTW,你可能会问,为什么要强调快速计算,累乘起来不就好了吗,比如:
func pow(a,c,m int) int {
ans := 1
for i:=0; i<c; i++ {
ans *= a
ans %= m
}
return ans
}
这么算当然没问题,但是考虑到实际情况,各种题目为了保证你答案的准确性,会尽量找一个大的质数来让你取模,一般都是10^9+7,也就是c,这种数据量线性累乘是不行的,严重超时。怎么办呢?
快速幂
快速幂利用分治的思想:
a^2n=a^n*a^n=(a^n)^2
a^(2n+1)=a*a^n*a^n=a*(a^n)^2
时间复杂度为logN,显著地减少了我们的计算量。
递归的实现:
func quickPow(a,n,m int) int {
if n < 0 {
return 0
}
a %= m
if n == 1 {
return a
}
ans := 1
if n % 2 == 1 {
ans = a
}
t := quickPow(a,n/2,m)
return ans*t%m*t%m
}
非递归的实现:
func quickPow(a,n,m int) int {
if n < 0 {
return 0
}
ans := 1
a %= m
for n>0 {
if n%2==1 {
ans *= a
ans %= m
}
n /=2
a *= a
a %= m
}
return ans
}
写在最后
从这个里面,我发现最有意思的就是各种语言当ab异号时对 a/b 的处理。负数整除正整数,或者正整数整除负整数,到底有没有什么实际物理意义。如果不能整除,结果向下取整,那么-7/5=-1.4,向下取整,也就是取一个刚好比-1.4小的整数,应该是-2才对,-1是不是一个设计错误?
平时的代码中以后遇到除法、减法和取模,也一定要特别小心