poj-3641Pseudoprime numbers

poj-3641Pseudoprime numbers

Pseudoprime numbers
Time Limit: 1000MS
Memory Limit: 65536K
Total Submissions: 10836
Accepted: 4686
Description
Fermat’s theorem states that for any prime number p and for any integer a > 1, ap = a (mod p). That is, if we raise a to the pth power and divide by p, the remainder is a. Some (but not very many) non-prime values of p, known as base-a pseudoprimes, have this property for some a. (And some, known as Carmichael Numbers, are base-a pseudoprimes for all a.)
Given 2 < p ≤ 1000000000 and 1 < a < p, determine whether or not p 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 containing p and a.
Output
For each test case, output “yes” if p is a base-a pseudoprime; otherwise output “no”.
Sample Input
3 2
10 3
341 2
341 3
1105 2
1105 3
0 0
Sample Output
no
no
yes
no
yes
yes
Source
Waterloo Local Contest, 2007.9.23

Pollard’s rho + 快速幂
由于poj不能用rand（），所以随便给了个数。
正式的Pollard’s rho的代码：

#include <iostream>#include <cstdio> #include <cstring>#include <cmath>#include <algorithm>#include <ctime>using namespace std;int gcd(int a, int b) {    return b == 0 ? a : gcd(b, a%b);}int Pollard_rho(int n){    int x = rand()%n;    int y = x, k = 2, i = 1;    while(1) {      i++;      x = ((2*x*x-4*x)%n+n)%n;      int d = gcd(abs(y-x), n);      if (d != 1 && d != n) return d;      if (i == k) {        y = x;        k = 2 * k;      }      if(i >= 100) break;    }     x = rand()%n;    y = x, k = 2, i = 1;    while(1) {      i++;      x = ((x*x-1)%n+n)%n;      int d = gcd(abs(y-x), n);      if (d != 1 && d != n) return d;      if (i == k) {        y = x;        k = 2 * k;      }      if(i >= 100) break;    }     return 1;}int main(){    time(NULL);    int n;    cin>>n;    cout<<Pollard_rho(n);    return 0;}

pollard’s rho的原理：
（1）xi = f(xi-1)%n,{xi}一定会形成一个环
（2）如果q是n的一个约数且q|（y-x），那么，1 < gcd(abs(y-x), n)<=n。
（3）生日悖论：如果随机选取（x， y）这样的数组sqrt（n）次，那么有50%的概率（x，y）模n同余。
（4）以伪随机的方法取（x，y），1 < gcd（y-x，n）< n即为所求。
（5）时间复杂度：期望O(logq)时间内找到n的一个较小的约数

答案：

#include <iostream>#include <cstdio>#include <cstring>#include <cmath>#include <cstring>#include <ctime>#include <algorithm>using namespace std;typedef long long LL;LL gcd(LL a, LL b){    return b == 0 ? a : gcd(b, a%b);}LL abs_(LL n){    return n >= 0 ? n : -n;}int pollard_rho(int n){    LL x = 12362 % n;    LL y = x, k = 2, i = 1;    while (1) {      i++;      x = ((x*x-3*x+451)%n+n)%n;      LL d = gcd(abs_(y-x), n);      if (d != 1 && d != n) return d;      if (i == k) {        y = x;        k *= 2;      }      if (i >= 100) break;    }    x = 1798237 % n;    y = x, k = 2, i = 1;    while (1) {      i++;      x = ((-3*x*x+4*x+12321)%n+n)%n;      LL d = gcd(abs_(x-y), n);      if (d != 1 && d != n) return d;      if (i == k) {        y = x;        k *= 2;      }      if (i >= 100) break;    }    return 1;}LL Fast(LL x, LL s){    LL res = 1, mod = s;    while (s) {      if (s & 1) res = (res * x) % mod;      x = (x * x) % mod;      s = s >> 1;    }     return (res+mod)%mod;}int main(){    LL n, p;    while (scanf("%lld%lld", &p, &n) && n) {      if (pollard_rho(p) != 1 && Fast(n, p)%p == n%p) printf("yes\n");       else printf("no\n");    }    return 0;}