POJ 1811 大整数素数判断 Miller_Rabin

#include <cstdio>
#include <cstring>
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <iostream>
using namespace std;

#define ll long long
const int S = 10;
ll ans;

ll gcd(ll a,ll b){
    if(a < 0) return gcd(-a , b);
    if (b==0) return a;
    return gcd(b , a%b);
}

ll multi_mod(ll a , ll b , ll k) //a*b%k
{
    a %= k;
    b %= k;
    ll ret = 0;
    //这里因为是长整数，for(i=0->b)循环一个个加太耗时，需要采用这种logn复杂度的方法
    while(b){
        if(b & 1){
            ret += a;
            ret %= k;
        }
        a <<= 1;
        a %= k;
        b >>= 1;
    }
    return ret;
}

ll pow_mod(ll a, ll b, ll mod) //求(a^b)%mod
{
    if(b == 1) return a%mod;
    int bit[1000] , t = 0;
    while(b){
        bit[t++] = b&1;
        b>>=1;
    }
    //类似矩阵快速幂
    ll ret = 1;
    for(int i = t-1 ; i>=0 ; i--){
        ret = multi_mod(ret , ret , mod);
        if(bit[i]) ret = multi_mod(ret , a , mod);
    }
    return ret;
}
//n = x*2^t , 这里以a为底 , 检验n是否为合数
bool check(ll a , ll n , ll x , ll t)
{
    ll ret = pow_mod(a , x , n);
    ll last = ret; //last记录上一次的数据，保证最后结果模1时，检查到last不为n-1或1，就表示为合数
    for(int i = 1 ; i<=t ; i++){
        ret = multi_mod(ret,ret,n);
        if(ret == 1 && last!=1 && last!=n-1) return true;
        last = ret;
    }
    //始终无法找到余数为1的结果，则表示n为合数
    if(ret!=1) return true;
    return false;
}
//n为合数返回true
bool miller_rabin(ll n)
{
    ll x = n-1 , t = 0;
    while((x&1) == 0){
        x>>=1;
        t++;
    }
    bool flag = true;

    //保证最后得到的 x^2 mod p = 1 只有x=1或x=p-1两个解，才满足二次探测，否则有其他解就直接证明不是素数
    if(t >= 1 && (x&1)){
        //miller_robin是一种取随机测试数据的算法，这里S=20,给定20组测试数据，当然数据组数越多，正确率越大
        for(int i = 0 ; i<S ; i++){
            ll a = rand()%(n-1)+1;
            //二次探测在check中检测
            if(check(a,n,x,t)){
                flag = true;
                break;
            }
            flag = false;
        }
    }
    if(!flag || n==2)
        return false;
    else return true;
}

//长整数利用随机数找到其中一个因子
ll Pollard_rho(ll x,ll c){
    ll i=1,x0=rand()%x,y=x0,k=2;
    while (1){
        i++;
        x0=(multi_mod(x0,x0,x)+c)%x;
        ll d=gcd(y-x0,x);
        if (d!=1&& d!=x){
            return d;
        }
        if (y==x0) return x;
        if (i==k){
            y=x0;
            k+=k;
        }
    }
}

//不断利用Pollard_rho递归来查找到n的所有素数因子
void find_factor(ll n)
{
    if(!miller_rabin(n)){
        ans = min(ans , n);
        return;
    }
    ll p = n;
    while(p >= n)
        p = Pollard_rho(p , rand() % (n-1) + 1);
    find_factor(p);
    find_factor(n / p);
}

int main()
{
    srand(time(NULL));
    int T;
    scanf("%d" , &T);
    while(T--){
        ll n;
        scanf("%lld" , &n);
        bool flag = miller_rabin(n);
        if(!flag) printf("Prime
");
        else{
            ans = n;
            find_factor(n);

            printf("%lld
" , ans);
        }
    }
    return 0;
}