POJ 1811 Prime Test (大素数判断和素因子分解）

给你一个数N(2 <= N < 2^54) ，若N是素数 ， 输出Prime， 否则输出最小的素因子。

由于N很大，所以只能先用Miller_Rabin算法进行素数判断，然后用Pollard_rho分解因子。

Miller-Rabin算法本质上是一种概率算法，存在误判的可能性，但是出错的概率非常小。出错的概率到底是多少，存在严格的理论推导。

费马小定理

• 如果p是质数且(a,p)=1，则有a^p−1≡1(mod p)。
  当然反过来不一定成立。即当(a^p−1)%p=1时，p未必是质数。但是这个概率比较小。所以利用费尔马小定理来检测素数，不能保证时刻都对，只能保证出错的概率比较小。
  给定正整数n，问n是否为质数（显然只需判断正奇数），最基本的做法就是计算2^(n−1)%n是否为1。如果不是1，n肯定为合数；否则，n可能为质数。

二次探测定理

• 如果p是一个素数且x^2≡1(mod p) , 则 x=1或者x=−1，注意到在模p的意义下，x=−1等价于x=p−1。
• 证明：方程x^2≡1(mod p)可以改写为(x + 1)(x - 1) % p = 0.所以x为1或-1.但在mod p的情况下，x=-1 等价于 x=p-1

Miller-Rabin算法

利用上面两个定理，就可以构造出Miller-Rabin算法。考虑到n肯定是奇数,则n一定可以表示为n−1=2^s∗d，其中s≥1且d是奇数。则

    a^(n−1)=a^(2^s∗d)=(((a^d)^2)^...^)2

也就是说，a^(n−1)相当于a^d平方若干次。例如当n=7时，a^(n−1)就是a^6，就是a^3的平方。当n=13时，a^(n−1)就是a^12，就是a^3的平方的平方。
以n=13的情况进行说明（所有运算都是在模n的意义下，以下的文字说明省略了这一点），任取一个a , 1< a <13，计算a^3，再将其平方一次得到a^6，注意到a^3是a^6的平方根，根据平方根定理的推论，如果a^6=1且a^3≠±1，则n肯定是合数。将a^6平方一次得到a^12，同样，如果a^12=1且a^6≠±1，则n肯定是合数。最后，根据费尔马小定理，如果a^12≠1，则n肯定是合数。否则，n有极大概率为质数。
为了增加得到正确判断的概率，可以将a重复取不同的值，对每一个a验证一次a^d到a^(n−1)的过程。不过，考虑到ACM的特殊性，测试数据应该不会选择伪素数特别是强伪素数。所以很多题目的AC程序实际上只来一次即可。

代码

typedef long long ll;//利用二进制计算a*b%modll multiMod(ll a,ll b,ll mod){    a %= mod , b %= mod;    ll ret = 0;    while(b){        if (b & 1) {            ret += a;            if(ret >= mod) ret -= mod;        }        a <<= 1;        if(a >= mod) a -= mod;        b >>= 1;    }    return ret;}//计算a^b%modll powerMod(ll a,ll b,ll mod){    ll ret = 1;    a %= mod;    while(b){        if (b & 1) ret = multiMod(ret , a , mod);        b >>= 1;        a = multiMod(a , a , mod);    }    return ret;}//Miller-Rabin测试,测试n是否为素数bool Miller_Rabin(ll n,int repeat){    if (2 == n || 3 == n) return true;    if (!(n & 1)) return false;    //将n分解为2^s*d    ll d = n - 1;    int s = 0;    while((d&1)== 0) ++s, d>>=1LL;    //srand((unsigned)time(0)); G++不能使用srand(time(NULL)) ， 会RE    for(int i = 0; i < repeat; ++i){//重复repeat次        ll a = rand() % (n - 3) + 2;//取一个随机数,[2,n-1)        ll x = powerMod(a , d , n);        ll y = 0;        for(int j = 0; j < s; ++j){            y = multiMod(x , x , n);            if ( 1 == y && 1 != x && n-1 != x ) return false;            x = y;        }        if ( 1 != y ) return false;    }    return true;}

参考：http://blog.csdn.net/u012061345/article/details/48241581#

pollard_rho 算法

• 通过某种方法得到两个整数a和b，而待分解的大整数为n，计算p=gcd(a-b,n)，直到p不为1，或者a，b出现循环为止。然后再判断p是否为n，如果p=n成立，那么返回n是一个质数，否则返回p是n的一个因子，那么我们又可以递归的计算Pollard(p)和Pollard(n/p)，这样，我们就可以求出n的所有质因子。
• 具体操作中，我们通常使用函数x2=x1*x1+c来计算逐步迭代计算a和b的值，实践中，通常取c为1，即b=a*a+1，在下一次计算中，将b的值赋给a，再次使用上式来计算新的b的值，当a，b出现循环时，即可退出进行判断。
• 在实际计算中，a和b的值最终肯定一出现一个循环，而将这些值用光滑的曲线连接起来的话，可以近似的看成是一个ρ型的。
  对于Pollard rho，它可以在O(sqrt(p))的时间复杂度内找到n的一个小因子p，可见效率还是可以的，但是对于一个因子很少、因子值很大的大整数n来说，Pollard rho算法的效率仍然不是很好，那么，我们还得寻找更加的方法了。

代码

ll factor[100]; //素因子的结果，返回时是无序的int tol; //素因子的个数ll gcd(ll a , ll b){    ll t;    while(b) {t = a; a = b; b = t%b;}    if(a >= 0) return a;    return -a;}ll pollard_rho(ll n , ll c){    ll i = 1 , k = 2;    //srand(time(NULL));    ll a = rand() % (n - 1) + 1 , b = a;    while(1)    {        i++;        a = (multiMod(a , a , n) + c) % n;        ll p = gcd(b - a , n);        if(p != 1 && p != n) return p;        if(b == a) return n;        if(i == k) {b = a; k += k;}    }}void findfac(ll n , int k)//k为自己设定的值，一般取1{    if(1 == n) return ;    if(Miller_Rabin(n , 8))    {        factor[tol++] = n;        return ;    }    ll p = n;    int c = k;    while(p >= n) p = pollard_rho(p , c--);    findfac(p , k);    findfac(n/p , k);}

参考：http://blog.sina.com.cn/s/blog_86a9d97201015cj7.html