快速幂取模算法
问题:求解ab%c
1.首先看最原始的做法:
1 int num=1; 2 for(int i=1;i<=b;i++){ 3 num*=a; 4 } 5 num%=c;
这个算法的复杂度为O(b),而且容易知道,当a,b很大时,很容易溢出,改进吧!
2.有一个定理:积的取余等于取余的积的取余。(a*b)%c=( (a%c)*(b%c) )%c
所以可以以此优化:
1 int num=1; 2 a%=c; 3 for(int i=1;i<=b;i++){ 4 num*=a; 5 } 6 num%=c;
3.更进一步,既然某个因子取余之后相乘再取余保持余数不变,那么新算得的num也可以进行取余,所以得到比较良好的改进版本。
1 int num=1; 2 a%=c; 3 for(int i=1;i<=b;i++){ 4 num=(num*a)%c; 5 } 6 num=num%c;
这个算法在时间复杂度上没有改进,仍为O(b),不过已经好很多的,但是在c过大的条件下,还是很有可能超时。怎么办呢。。。。那只有在b上下功夫啦
4.这时候想到讨论b的奇偶性。
(1)如果b是偶数,我们可以记k = a2mod c,那么求(k)b/2 mod c就可以了。
(2)如果b是奇数,我们也可以记k =a2 mod c,那么求
((k)b/2 mod c × a ) mod c =((k)b/2 mod c * a) mod c 就可以了。
于是有:
1 int num=1; 2 a%=c; 3 if(b%2==0){ 4 num=(num*2)%c;//如果b是奇数,要多求一步,可以提前算到num中 5 } 6 int k=(a*a)%c;//取a2而不是a 7 8 for(int i=1;i<=b/2;i++){ 9 num=(num*k)%c; 10 } 11 num=num%c;
5.可以看到,我们把时间复杂度变成了O(b/2).当然,这样子治标不治本。但我们可以看到,当我们令k = (a * a)mod c时,状态已经发生了变化,我们所要求的最终结果即为(k)b/2 mod c而不是原来的ab mod c,所以我们发现这个过程是可以迭代下去的。当然,对于奇数的情形会多出一项a mod c,所以为了完成迭代,当b是奇数时,我们通过 num = (num * a) % c;来弥补多出来的这一项,此时剩余的部分就可以进行迭代了。 形如上式的迭代下去后,当b=0时,所有的因子都已经相乘,算法结束。于是便可以在O(log b)的时间内完成了。于是,有了最终的算法:快速幂算法。
于是有:
1 int num=1; 2 a%=c; 3 while(b>0){ 4 if(b%2==1){ 5 num=(num*a)%c; 6 } 7 b/=2; //这一步将b->log2(b) 8 a=(a*a)%c; 9 }
函数化上述代码,解决ab%c的问题:
1 unsigned fastmod(unsigned a,unsigned b,unsigned c){ 2 unsigned num=1; 3 a%=c;//这里不进行初始化也是可以的 4 5 while(b>0){ 6 if(b%2==1){ 7 num=(num*a)%c; 8 } 9 b/=2; //这一步将b->log2(b) 10 a=(a*a)%c; 11 } 12 return num; 13 }
现在时间复杂度为O(log b),经过中间的取模运算,也不必担心溢出的问题,问题基本解决。
这其中基本参数的解释参照 http://blog.csdn.net/chen77716/article/details/7093600 中反复平方法的介绍,这是目前见到的最为清楚的解释。
下面写一道ural 1528的运用
原题如下:
Description
Input
Output
Sample Input
input | output |
---|---|
1 2 2 11 0 0 |
1 2 |
AC代码如下:
1 #include <iostream> 2 #include <cmath> 3 using namespace std; 4 5 int main(){ 6 int n,p; 7 while(1){ 8 cin>>n>>p; 9 if(n==0&&p==0){ 10 break; 11 } 12 13 long long temp=1; 14 for(int i=1;i<=n;i++){ 15 temp=temp*i%p;//(a*b)%c=((a%c)*(b%c))%c 16 } 17 18 cout<<temp<<endl; 19 20 } 21 return 0; 22 }
--------------------------------------
新看到了一种计算乘法的方式,或许可以帮助理解该算法:
1 int mul(int x,int y){ 2 if(y==0) return 0; 3 int z=mul(x,floor(y/2)); 4 if(y%2==0){ 5 return 2*z; 6 } 7 else{ 8 return x+2*z; 9 } 10 }
根据这一思想,可以写出幂取模的一种递归实现:
1 int z=0;//即结果 2 //计算x^y mod n 3 int modexp(int x,int y,int n){ 4 if(y==0) return 1; 5 z=modexp(x,floor(y/2),n); 6 if(z%2==0) 7 return z*z%n; 8 else 9 return x*z*z%n; 10 }