Pseudoprime numbers POJ - 3641

Fermat's theorem states that for any prime number p and for any integera > 1, a^p = a (mod p). That is, if we raise a to the p^th power and divide byp, the remainder is a. Some (but not very many) non-prime values ofp, known as base-a pseudoprimes, have this property for some a. (And some, known as Carmichael Numbers, are base-a pseudoprimes for alla.)

Given 2 < p ≤ 1000000000 and 1 < a < p, determine whether or notp is a base-a pseudoprime.

Input

Input contains several test cases followed by a line containing "0 0". Each test case consists of a line containingp and a.

Output

For each test case, output "yes" if p is a base-a pseudoprime; otherwise output "no".

Sample Input

3 210 3341 2341 31105 21105 30 0

Sample Output

nonoyesnoyesyes

题意：给你两个数p和a，如果a的p次方对p取余等于a并且p不是素数，输出yes，否则输出no。
思路：因为数据过大,直接写for循环求次方会超时.这里就要用到快速幂这个技巧了.简单说就是降幂,比如说a的b次方,如果b是偶数,该数就可以表示为(a^2)^(b/2),如果b是奇数,可以 另写一个数保存此时的底数然后幂数减1,又可以继续降幂了,知道最后。优化代码。注意__int64输入防止爆int

#include <iostream>#include <algorithm>#include <cstdio>#include <cstring>#include <cmath>using namespace std;bool is_prime(long long n){    for(long long i=2;i<=sqrt(n);i++)    {        if(n%i==0)        {            return 0;        }    }    return 1;}long long mod_exp(long long p,long long a){    long long ans=1;    long long b=p;//b为变化的指数    while(b)    {        if(b%2==1)        {            ans=(ans*a)%p;        }        a=(a*a)%p;//此处为ans才是快速幂(更新底数才对        b=b/2;    }    return ans%p;}int main(){    long long a,p;    while(~scanf("%I64d %I64d",&p,&a))    {        if(p==0&&a==0)            break;        if(is_prime(p)==0&&mod_exp(p,a)==a)//不是质数           printf("yes\n");        else            printf("no\n");    }    return 0;}

附上经典解析

2.6 数学问题的解题窍门
快速幂运算
用反复平方法做快速幂运算判断条件②，条件①嘛，朴素的素性测试就行了。
有读者问我为什么特意写个函数LL mod_mult(LL a, LL b, LL m)来求(a * b) % m，因为这样不容易溢出。如果直接用运算符先*后%的话，哪怕是unsigned long long也可能溢出。
这个求余乘法的思想是，先将一个数用2进制表示：
bn表示b的二进制的第n个bit，当然，首个比特是从0开始算的。将a乘入括号中，得到：
由于bn要么是0要么是1，所以只需计算为1的部分就可以了，比如3*5：
每加一次就求一次余，这样每次加上去的都是小于m的余数，这样就不怕溢出了。由于每个bit都需要计算一次，所以复杂度是O(log(N))。
#ifndef ONLINE_JUDGE
#pragma warning(disable : 4996)
#endif
#include <iostream>
using namespace std;
 
typedef long long LL;
 
// return (a * b) % m
LL mod_mult(LL a, LL b, LL m)
{
LL res = 0;
LL exp = a % m;
while (b)
{
if (b & 1)
{
res += exp;
if (res > m) res -= m;
}
exp <<= 1;
if (exp > m) exp -= m;
b >>= 1;
}
return res;
}
 
// return (a ^ b) % m
LL mod_exp(LL a, LL b, LL m) 
{
LL res = 1;
LL exp = a % m;
while (b)
{
if (b & 1) res = mod_mult(res, exp, m);
exp = mod_mult(exp, exp, m);
b >>= 1;
}
return res;
}
 
//************************************
// Method:    is_prime
// FullName:  is_prime
// Access:    public 
// Returns:   bool
// Qualifier: 素性测试
// Parameter: const int & n
//************************************
bool is_prime(const int& n)
{
for (int i = 2; i * i <= n; ++i)
{
if (n % i == 0)
{
return false;
}
}
 
return n != 1;// 1是例外
}
 
///////////////////////////SubMain//////////////////////////////////
int main(int argc, char *argv[])
{
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
freopen("out.txt", "w", stdout);
#endif
int p, a;
while (cin >> p >> a && p && a)
{
if (!is_prime(p) && (mod_exp(a, p, p) == a))
{
cout << "yes" << endl;
}
else
{
cout << "no" << endl;
}
}
#ifndef ONLINE_JUDGE
fclose(stdin);
fclose(stdout);
system("out.txt");
#endif
return 0;
}
///////////////////////////End Sub//////////////////////////////////