题目链接：
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;
}