快速幂取模算法

问题：求解a^b%c

1.首先看最原始的做法：

int num=1;
for(int i=1;i<=b;i++){
    num*=a;
} 
num%=c;

这个算法的复杂度为O(b),而且容易知道，当a,b很大时，很容易溢出，改进吧！

2.有一个定理：积的取余等于取余的积的取余。(a*b)%c=( (a%c)*(b%c) )%c

所以可以以此优化：

int num=1;
a%=c;
for(int i=1;i<=b;i++){
    num*=a;
} 
num%=c;

3.更进一步，既然某个因子取余之后相乘再取余保持余数不变，那么新算得的num也可以进行取余，所以得到比较良好的改进版本。

int num=1;
a%=c;
for(int i=1;i<=b;i++){
    num=(num*a)%c;
} 
num=num%c;

这个算法在时间复杂度上没有改进，仍为O(b)，不过已经好很多的，但是在c过大的条件下，还是很有可能超时。怎么办呢。。。。那只有在b上下功夫啦

4.这时候想到讨论b的奇偶性。

（1）如果b是偶数，我们可以记k = a²mod c，那么求(k)^b/2 mod c就可以了。

（2）如果b是奇数，我们也可以记k =a² mod c，那么求

((k)^b/2 mod c × a ) mod c =((k)^b/2 mod c * a) mod c 就可以了。

于是有：

int num=1;
a%=c;
if(b%2==0){
    num=(num*2)%c;//如果b是奇数，要多求一步，可以提前算到num中 
}
int k=(a*a)%c;//取a2而不是a 

for(int i=1;i<=b/2;i++){
    num=(num*k)%c;
} 
num=num%c;

5.可以看到，我们把时间复杂度变成了O(b/2).当然，这样子治标不治本。但我们可以看到，当我们令k = (a * a)mod c时，状态已经发生了变化，我们所要求的最终结果即为(k)^b/2 mod c而不是原来的a^b mod c，所以我们发现这个过程是可以迭代下去的。当然，对于奇数的情形会多出一项a mod c，所以为了完成迭代，当b是奇数时，我们通过   num = (num * a) % c;来弥补多出来的这一项，此时剩余的部分就可以进行迭代了。                      形如上式的迭代下去后，当b=0时，所有的因子都已经相乘，算法结束。于是便可以在O（log b）的时间内完成了。于是，有了最终的算法：快速幂算法。

于是有：

int num=1;
a%=c;
while(b>0){
    if(b%2==1){
        num=(num*a)%c;    
    }
    b/=2;     //这一步将b->log2(b) 
    a=(a*a)%c;
}

函数化上述代码，解决a^b%c的问题：

unsigned fastmod(unsigned a,unsigned b,unsigned c){
    unsigned num=1;
    a%=c;//这里不进行初始化也是可以的
    
    while(b>0){
        if(b%2==1){
            num=(num*a)%c;    
        }
        b/=2;     //这一步将b->log2(b) 
        a=(a*a)%c;
    }
    return num;
}

现在时间复杂度为O（log b），经过中间的取模运算，也不必担心溢出的问题，问题基本解决。

这其中基本参数的解释参照  http://blog.csdn.net/chen77716/article/details/7093600   中反复平方法的介绍，这是目前见到的最为清楚的解释。

下面写一道ural  1528的运用

原题如下：

Sequence II

Time Limit:3000MS     Memory Limit:65536KB     64bit IO Format:%I64d & %I64u

Description

You are given a recurrent formula for a sequence f:

f( n) = 1 + f(1) g(1) + f(2) g(2) + … + f( n−1) g( n−1),

where g is also a recurrent sequence given by formula

g( n) = 1 + 2 g(1) + 2 g(2) + 2 g(3) + … + 2 g( n−1) − g( n−1) g( n−1).

It is known that f(1) = 1, g(1) = 1. Your task is to find f( n) mod  p.

Input

The input consists of several cases. Each case contains two numbers on a single line. These numbers are n (1 ≤  n ≤ 10000) and p (2 ≤ p ≤ 2·10 ⁹). The input is terminated by the case with n = p = 0 which should not be processed. The number of cases in the input does not exceed 5000.

Output

Output for each case the answer to the task on a separate line.

Sample Input

input	output
1 2 2 11 0 0	1 2

AC代码如下：

#include <iostream>
#include <cmath>
using namespace std;

int main(){
    int n,p;
    while(1){
        cin>>n>>p;
        if(n==0&&p==0){
            break;
        }
        
        long long temp=1;
        for(int i=1;i<=n;i++){
            temp=temp*i%p;//(a*b)%c=((a%c)*(b%c))%c  
        }
        
        cout<<temp<<endl;
        
    }
    return 0;   
} 

 --------------------------------------

新看到了一种计算乘法的方式，或许可以帮助理解该算法：

int mul(int x,int y){
    if(y==0) return 0;
    int z=mul(x,floor(y/2));
    if(y%2==0){
        return 2*z;
    } 
    else{
        return x+2*z;
    }
}

 根据这一思想，可以写出幂取模的一种递归实现：

int z=0;//即结果
//计算x^y mod n 
int modexp(int x,int y,int n){
    if(y==0) return 1;
    z=modexp(x,floor(y/2),n);
    if(z%2==0)
        return z*z%n;
    else
        return x*z*z%n;
}