POJ 1811 Prime Test （Pollard rho 大整数分解）

#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <string.h>
#include <cstdlib>
#include <cmath>

using namespace std;
long long n;
long long minf;
long long multi(long long a,long long b,long long mod){
    a=a%n;
    long long ret=0;
    while(b){
        if(b&1){
            //ret=(ret+a)%mod;
            ret=ret+a;
            if(ret>=mod)
                ret-=mod; //原来用取模1500ms，现在547ms。速度变为原来的1/3，没想到取模那么慢。。。
        }
        a=a<<1;
        if(a>=mod)
            a-=mod;
        //a=(a+a)%mod;
        b=b>>1;
    }
    return ret;
}
long long quickPow(long long a,long long b,long long mod) {
    long long ret=1;
    while(b) {
        if(b&1)
            ret=multi(ret,a,mod);
        a=multi(a,a,mod);
        b=b>>1;
    }
    return ret;
}

bool witness(long long a,long long n) {
    long long m=n-1;
    int j=0;
    while(!(m&1)) {
        j++;
        m=m>>1;
    }
    long long x=quickPow(a,m,n);
    if(x==1||x==n-1)
        return false;
    while(j--) {
        x=(x*x)%n;
        if(x==n-1)
            return false;
    }
    return true;
}
//用这个547ms
bool Miller_Rabin(long long n) {
    if(n<2)
        return false;
    if(n==2)
        return true;
    if(!(n&1))
        return false;
    long long m=n-1;
    int j=0;
    while(!(m&1)) {
        j++;
        m=m>>1;
    }
    for(int i=1; i<=50; i++) {
        long long a=rand()%(n-1)+1;
        long long x=quickPow(a,m,n);
        if(x==1||x==n-1)
            continue; //n可能为素数
        int k=j;
        bool flag=false;
        while(k--) {
            x=(x*x)%n;
            if(x==n-1){
                flag=true; //n可能为素数
                break;
            }
        }
        if(flag)
            continue;
        return false;  //n为合数
    }
    return true;
}
//用这个454ms
bool Miller_Rabin2(long long n)
{
    if(n==2)return 1;
    if(n<2||!(n&1))return 0;

    long long a, u=n-1, x, y;
    int t=0;
    while(u%2==0){
        t++;
        u/=2;
    }
    srand(100);
    for(int i=0;i<6;i++)
    {
        a = rand() % (n-1) + 1;
        x = quickPow(a, u, n);
        for(int j=0;j<t;j++)
        {
            y = multi(x,x,n);
            //如果n是素数那么x^2mod n = 1,仅有两个根（不同余），+1和-1(n-1)。如果当x^2 mod n=1，但是x不为1或n-1，那么n是合数
            if ( y == 1 && x != 1 && x != n-1 )
                return false; //为合数
            x = y;
        }
        if( y!=1) return false;//最后y的值为a^(n-1)。由费马小定理 a^(n-1) mod n 恒等于 1 则可以认为 n 是素数。若有不为1的情况，则n为合数
    }
    return true;
}

long long gcd(long long a,long long b) {
    return b==0?a:gcd(b,a%b);
}
long long pollard_rho(long long n,long long c) {
    long long x,y,d,i=1,k=2;
    x=rand()%(n-1)+1;
    y=x;
    while(1) {
        i++;
        x=(multi(x,x,n)+c)%n;
        d=gcd(y-x,n);
        if(1<d && d<n)
            return d;
        if(y==x)
            return n;
        if(i==k) {
            y=x;
            k=k<<1;
        }
    }
}
//对n进行素数分解
void factorFind(long long n,int k) {
    if(n==1)
        return;
    if(Miller_Rabin(n)) {
        if(n<minf)
            minf=n;
        return;
    }
    long long p=n;
    while(p>=n)
        p=pollard_rho(p,k--);
    factorFind(p,k);
    factorFind(n/p,k);
}
int main() {
    int t;
    scanf("%d",&t);
    while(t--) {
        scanf("%I64d",&n);
        minf=n+1;
        if(Miller_Rabin(n))
            printf("Prime\n");
        else {
            factorFind(n,150);
            printf("%I64d\n",minf);
        }
    }
    return 0;
}