bzoj3667 Rabin-Miller算法

传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=3667

【题解】

PollardRho，讲解见http://www.cnblogs.com/galaxies/p/bzoj4802.html

# include <stdio.h>
# include <string.h>
# include <iostream>
# include <algorithm>
// # include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
const int M = 3e4 + 10;

# define RG register
# define ST static

ll n;

inline ll mul(ll a, ll b, ll P) {
    a %= P, b %= P;
    ll ret = 0;
    while(b) {
        if(b&1) {
            ret = ret + a;
            if(ret >= P) ret -= P;
        }
        a <<= 1;
        if(a >= P) a -= P;
        b >>= 1;
    }
    return ret;
}

inline ll pwr(ll a, ll b, ll P) {
    a %= P;
    ll ret = 1;
    while(b) {
        if(b&1) ret = mul(ret, a, P);
        a = mul(a, a, P);
        b >>= 1;
    }
    return ret;
}

namespace Millar_Rabin {
    const int Prime[] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
    const int PN = 10;
    inline bool witness(int pr, ll res, int times, ll n) {
        ll p = pwr((ll)pr, res, n);
        for (int i=0; i<times; ++i) {
            if(p == 1) return 0;
            if(p == n-1) return 0; 
            p = mul(p, p, n); 
        }
        return 1;
    }
    inline bool main(ll n) {
        if(n <= 1) return 0;
        for (int i=1; i<=PN; ++i) {
            if(n == Prime[i]) return 1;
            if(n % Prime[i] == 0) return 0;
        }
        ll p = n-1;
        int times = 0;
        while(!(p&1)) {
            times ++;
            p >>= 1;
        }
        for (int i=1; i<=PN; ++i)
            if(witness(Prime[i], p, times, n)) return false;
        return true;
    }
}

namespace Pollard_Rho {
    int p[M], pn; 
    
    inline void PollardRho(ll n) {
         if(Millar_Rabin::main(n)) {
             p[++pn] = n;
            return ;
        }
        ll a, b, c, del;
        while(1) {
            c = rand() % n;
            a = b = rand() % n;
            b = (mul(b, b, n) + c) % n;
            while(a != b) {
                del = __gcd(abs(a-b), n);
                if(del > 1 && del < n) {
                    PollardRho(del);
                    PollardRho(n/del);
                    return ;
                }
                a = (mul(a, a, n) + c) % n;
                b = (mul(b, b, n) + c) % n;
                b = (mul(b, b, n) + c) % n;
            }         
        }
    }
    
    inline void main(ll n) {
        if(Millar_Rabin::main(n)) {
            cout << "Prime" << endl;  
            return ;
        }
        pn = 0;
        PollardRho(n); 
        sort(p+1, p+pn+1);
        cout << p[pn] << endl; 
    }
}

int main() {
    srand(20001130); 
    int T; cin >> T;
    while(T--) {
        cin >> n;
        Pollard_Rho::main(n);
    }
    return 0;
}