Miller_rabin & pollard_rho 快速素数检和合数因式分解

关于 POJ 1811 题解

http://poj.org/problem?id=1811

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <string>

using namespace std;

class Solution {
public:
    long long factors[100];
    long long f_index;
    const static int rabin_check_times = 5;

    long long mult_mod(long long m1, long long m2, long long mod) {
        // return m1 * m2 % mod
        long long ret = 0;
        m1 %= mod;
        m2 %= mod;
        while (m2 > 0) {
            if (m2 & 1) { ret = (ret + m1) % mod; }
            m1 = (m1 << 1) % mod;
            m2 = m2 >> 1;
        }
        return ret;
    }

    long long pow_mod(long long a, long long pow, long long mod) {
        a %= mod;
        // return a^pow % mod
        if (pow == 0) return 1L;
        else if (pow == 1) return a;
        long long t = a, ret = 1;
        while (pow > 0) {
            if (pow & 1) ret = mult_mod(ret, t, mod);
            pow = pow >> 1;
            t = mult_mod(t, t, mod);
        }
        return ret;
    }

    bool miller_rabin_check(long long a, long long p) {
        // need to check a^(p-1) mod p = 1?
        // also this indicates to check whether a ^ (u * (2^t)) mod p = 1?
        // Then need to check that squaring (a^u) for t times, during which 
        // that the case where b1 * b1 mod p = 1 but b1 != 1 && b1 != p-1 cannot exist.
        long long u = p-1, t=0;
        while (u % 2 == 0) {
            t++;
            u >>= 1;
        }
        long long b1 = pow_mod(a, u, p);
        long long b2 = b1;
        for (int i = 0; i < t; i++) {
            b2 = mult_mod(b1, b1, p);
            if (b2 == 1 && b1 != 1 && b1 != p-1) return false; // is not a prime
            b1 = b2;
        }
        if (b1 != 1) return false; // is not a prime
        return true;    // is a prime (BE careful, this is only possible for a great certainty)
    }

    bool miller_rabin(long long p, long long check_times) {
        long long a;
        for (int i = 0; i < check_times; i++) {
            a = rand() % (p-1) + 1; // a's range: [1, p-1]
            if (!(miller_rabin_check(a, p))) return false;
        }
        return true;
    }

    // greatest common divisor
    long long gcd(long long a, long long b) {
        if (a < 0) a = -a;
        if (b < 0) b = -b;
        if (a < b) {
            a = a^b;
            b = a^b;
            a = a^b;
        }
        if (b == 0) return a;
        while (b > 0) {
            long long t = a % b;
            a = b;
            b = t;
        }
        return a;
    }

    long long pollard_rho(long long x, long long c) {
        long long i = 1, k = 2, x0 = rand() % x, y0 = x0;
        long long d;
        while (1) {
            i++;
            x0 = (mult_mod(x0, x0, x) + c) % x;
            d = gcd(y0 - x0, x);
            if (d > 1 && d < x) return d;
            if (y0 == x0) return x;
            if (i == k) {
                y0 = x0;
                k <<= 1;
            }
        }
    }

    // using pollard_rho reduces the time complexity to n^0.25
    void find_factors(long long x) {
        //printf("# %lld\n", x);
        if (miller_rabin(x, rabin_check_times)) {
            factors[f_index++] = x;
            return; 
        }
        long long p = x;
        while (p >= x) p = pollard_rho(x, rand() % (x-1) + 1);
        find_factors(p);
        find_factors(x/p);
    }

    void solve() {
        int kase; 
        long long n;
        scanf("%d", &kase);
        for (int i = 0; i < kase; i++) {
            scanf("%lld", &n);
            f_index = 0;
            find_factors(n);
            if (f_index == 1) printf("Prime\n");
            else {
                sort(factors, factors+f_index);
                printf("%lld\n", factors[0]);
            }
        }
    }
};

int main() {
    Solution solution;
    solution.solve();
    return 0;
}