【数学】PollardRho算法

这个是整数分解PollardRho算法。不是离散对数PollardRho算法，在竞赛里面离散对数一般使用BSGS算法解决，虽然消耗一些空间，但是复杂度稳定。

快速乘可以替换成其他正确形式的乘法。
MR可以适当去掉一些，例如只使用2,7,61。

每次调用PR返回其中一个非1的因子，所以在外面套一层getFactor返回其中一个非平凡因子。

对速度有要求时应该只对超过线性筛范围的数字进行getFactor。

不要直接调用PR，对质数调用PR会直接死循环。（MR算法只会把强伪质数误判为质数，不会导致PR算法TLE）

ll qmul(ll a, ll b, ll mod) {
    return (__int128)a * b % mod;
}

ll qpow(ll x, ll n, ll mod) {
    if(x >= mod)
        x %= mod;
    ll res = 1;
    while(n) {
        if(n & 1)
            res = qmul(res, x, mod);
        x = qmul(x, x, mod);
        n >>= 1;
    }
    return res;
}

namespace MillerRabin {

    // private
    bool MR(ll n, ll p) {
        for(ll k = n - 1; k; k >>= 1) {
            ll t = qpow(p, k, n);
            if(t != 1 && t != n - 1)
                return false;
            if((k & 1) == 1 || t == n - 1)
                // probably true
                return true;
        }
        // probably true
        return true;
    }

    // public
    bool isPrime(ll n) {
        if(n <= 1)
            return false;
        if(n <= 3)
            return true;
        if(!(n & 1))
            return false;
        static int basePrime[5] = {2, 3, 7, 61, 24251};
        for(int i = 1; i <= 5; ++i) {
            if(n == basePrime[i - 1])
                return true;
            if(!MR(n, basePrime[i - 1]))
                return false;
        }
        if(n == 46856248255981LL)
            return false;
        // probably true
        return true;
    }

}

using namespace MillerRabin;

namespace PollardRho {

    // private
    ll PR(ll n) {
        ll s = 0, t = 0, c = rand() % (n - 1) + 1;
        for(int i = 1;; i <<= 1, s = t) {
            ll p = 1;
            for(int j = 1; j <= i; ++j) {
                t = (qmul(t, t, n) + c) % n;
                p = qmul(p, abs(t - s), n);
                if(p == 0)
                    return n;
                if(j % 127 == 0) {
                    ll d = __gcd(p, n);
                    if(d > 1)
                        return d;
                }
            }
            ll d = __gcd(p, n);
            if(d > 1)
                return d;
        }
    }

    // public
    ll getFactor(ll n) {
        if(n <= 3)
            return n;
        if(!(n & 1))
            return 2;
        if(isPrime(n))
            return n;
        ll d = PR(n);
        while(d == n)
            d = PR(n);
        return d;
    }

}

using namespace PollardRho;