POJ 3641 Pseudoprime numbers (快速幂模)

题目链接：http://poj.org/problem?id=3641 题目大意：给定a,p，判断一个数是不是伪素数（伪素数是指p---p满足a^p = a (mod p)但p本身不是素数），2<=a,p<=1000000000. 思路： 我们发现指数p很大，不过因为结果mod p并且由定理(a * a) % p = (a % p) * (a % p)可知结果在int范围内，但是暴力解a^p显然还是不可能的. 这里就要用到一个著名的算法：蒙哥马利模算法（快速幂模）来快速计算a^p % m: 这种算法利用了二分的思想，可以达到O(logn)。 可以把b按二进制展开为：b = p(n)*2^n  +  p(n-1)*2^(n-1)  +…+   p(1)*2  +  p(0)（其中p(i) (0<=i<=n)为 0 或 1） 这样 a^b =   a^(p(n)*2^n)  *  a^(p(n-1)*2^(n-1))  *...*  a^(p(1)*2)  *  a^p(0) 对于p(i)=0的情况, a^（p(i) * 2^(i-1) ） = 1不用处理.我们要考虑的仅仅是p(i)=1的情况： 化简：a^(2^i)  = (  a^(  2^(i-1)  )  )^2 （这里很重要！！具体请参阅秦九韶算法：http://zh.wikipedia.org/wiki/%E7%A7%A6%E4%B9%9D%E9%9F%B6%E7%AE%97%E6%B3%95） 利用这一点，我们可以递推地算出所有的a^(2^i) 当然由算法1的结论，我们加上取模运算：a^(2^i)%c = ( (a^(2^(i-1))%c) * a^(2^(i-1)))  %c 于是再把所有满足p(i)=1的a^(2^i)%c按照算法1乘起来再%c就是结果， 即二进制扫描从最高位一直扫描到最低位. 扩展：当a为一个矩阵时，算法就是矩阵快速幂.   递归代码（我一直觉得递归真的很优美~~~）：

int exp_mod(int a, int p, int m){   //a^p (mod m)
    if (p == 0)
        return 1;
    if (p == 1)
        return a % m;
    
	long long t = exp_mod(a, p/2, m);	//long long防止下面t*t时溢出
    
	t = t * t % m;
    
	if (p & 1)
        t = t * a % m;
    
	return t;
}

  非递归代码：  

//计算a^bmodn   
int modexp(int a, int b, int n)   
{   
    int ret = 1;   
    long long tmp = a;   
    while(b)   
    {   
       //基数存在   
       if (b & 1) ret = ret * tmp % n;   
       tmp = tmp * tmp % n;   
       b >>= 1;   
    }   
    return ret;   
}     

  POJ 3641代码：  

#include 
using namespace std;
bool prime(int n){
    for (int i = 2; i * i <= n; i ++){
        if (n % i == 0){
            return false;
        }
    }
    return true;
}
int exp_mod(int a, int p, int m){   //a^p (mod m)
    if (p == 0)
        return 1;
    if (p == 1)
        return a % m;
    long long t = exp_mod(a, p/2, m);
    t = t * t % m;
    if (p & 1)
        t = t * a % m;
    return t;
}
int main(){
    int p, a;
    while(cin >> p >> a){
        if (p + a == 0){
            return 0;
        }
        if (exp_mod(a, p, p) == a){
            if (!prime(p)){
                cout << "yes\n";
            }
            else{
                cout << "no\n";
            }
        }
        else{
            cout << "no\n";
        }
    }
    return 0;
}