快速幂取模算法

问题：求解ab%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 = 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 就可以了。

于是有：

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而不是原来的ab 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;
}

函数化上述代码，解决ab%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 9). 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;
}