$素性检测问题$

给定一个数 $n$，判断其是否为素数。

$费马素性检测$

我们有费马小定理

当 $p$ 为素数时，对于任意整数 $a \in [1,p)$ 满足 $a^{p-1}\equiv 1 \pmod p$。

这启发我们可以随机一个 $a\in[2,p)$，若 $a^{p-1}\equiv 1 \pmod p$ 不成立，那么 $p$ 一定不是素数。

当然也有伪素数，即使不是素数，但也满足 $\forall a \in [1, p), a^{p-1}\equiv 1 \pmod p$，例如 $561$ 就是一个伪素数。

$\textrm{Miller-Rabin} 素性检测$

列出一个同余方程，$b^2\equiv 1\pmod p$，若 $p$ 为素数，那么就一定满足，$b-1$ 或 $b+1$ 为 $p$ 的倍数，于是解出 $b\equiv \pm 1\pmod p$

现在开始判断 $p$ 是否为素数：

• 首先用费马小定理，得出 $a^{p-1} \equiv 1 \pmod p$。

• $t\gets p - 1$

• 然后若 $t$ 为偶数，则使用二次检测，由 $(a^{\frac{t}{2}})^2\equiv 1\pmod p$ 得到 $a^{\frac{t}{2}} \equiv \pm 1 \pmod p$，若不成立则 $p$ 为合数。

• 若 $a^{\frac{t}{2}} \equiv -1\pmod p$，那么判断 $p$ 为素数，结束。

• 否则 $t \gets \frac{t}{2}$，且任然满足 $a^t\equiv 1\pmod p$，回到步骤 $3$。

以上为执行一次 Miller-Rabin 素性检测，实践中一般带入多个 $a$ 判断。

选取 $a = \{2, 7, 61\}$ 可以判断 $2^{32}$ 以内的数；

选取 $a = \{2,3,5,7,11,13,17,19,23,29,31,37\}$ 可以判断 $2^{78}$ 以内的数。

il ll Abs(ll x) {
    return x < 0 ? -x : x;
}

il ll gcd(ll x, ll y) {
    return y == 0 ? x : gcd(y, x % y);
}

il ll mul(ull x, ull y, ll p) {  //快速乘
    return (x * y - (ull)((long double)x / p * y) * p + p) % p;
}

il ll qpow(ll x, ll b, ll c) {
    ll res = 1;
    while (b) {
        if (b & 1) res = mul(res, x, c);
        x = mul(x, x, c);
        b >>= 1;
    }
    return res;
}

il bool mr(int x, long long p) {
    long long t = p - 1;
    if (qpow(x, t, p) != 1) return false;
    while (t % 2 == 0) {
        t >>= 1;
        ll h = qpow(x, t, p);
        if (h == p - 1) return true;
        else if (h != 1) return false;
    }
    return true;
}

il bool isprime(ll p) {
    if (p < 2) return 0;
    if (p == 2 || p == 3 || p == 5 || p == 7 || p == 11 || p == 13
    || p == 17 || p == 19 || p == 23 || p == 29 || p == 31 || p == 37) return 1;
    return mr(2, p) && mr(3, p) && mr(5, p) && mr(7, p) && mr(11, p) && mr(13, p) &&
    mr(17, p) && mr(19, p) && mr(23, p) && mr(29, p) && mr(31, p) && mr(37, p);
}

$\textrm{Pollard-Rho} 算法$

该算法用于大数的分解质因数。

$a,c$ 为随机数，数列 $<x>$ 定义为，$x_0 = a \bmod p，x_i = x_{i-1}^2 + c \bmod p$。

我们把每一个数想成一个点，由 $x_{i-1}$ 向 $x_i$ 连边，最终一定形成一个环，根据生日悖论，环长期望为 $O(\sqrt{p})$

若 $p = p1^{c1} \times h$，数列 $<x’>$ 定义为 $x’_0 = a \bmod p1, x’_i = x_{i-1}’^{2} + c \bmod p1$。

那么显然 $x’_i \equiv x_i \pmod p1$

设 $i,j$ 为数列上两个位置，那么若 $x_i \equiv x_j \pmod p1$， 则我们可以用 $\gcd(|x_i - x_j|, p)$ 求出 $p1$。

问题是如何试出 $i,j$，若 $i$ 必定在 $p1$ 的环上，那么 $j$ 只要走期望 $p^{\frac{1}{4}}$ 步，每走一步求一下 $\gcd(|x_i - x_j|, p)$，若大于 $1$，则求出 $p1$。

现在重点在于 $i$ 不一定在环上，我们使用倍增法，先让 $i=0,j=0$，$j$ 走 $2^0$ 步后，$i$ 直接向后跳到 $j$ 的位置，让 $j$ 走 $2^1$ 步，以此类推。

$进一步优化$

由于 $j$ 每走一步，都要求一次 $\gcd$，所以我们干脆记录一个 $M$，始终不清空，$M = \prod |x_i-x_j| \bmod p$，当 $j$ 每次走完 $2^k$ 步后，判断 $\gcd(M, p)$ 的值。

这样复杂度从 $O(p^{\frac{1}{4}}\log p)$ 降为 $O(p^{\frac{1}{4}})$。

#include <bits/stdc++.h>
#define ull unsigned long long
#define ll long long
#define il inline
using namespace std;

il ll Abs(ll x) {
    return x < 0 ? -x : x;
}

il ll gcd(ll x, ll y) {
    return y == 0 ? x : gcd(y, x % y);
}

il ll mul(ull x, ull y, ll p) {
    return (x * y - (ull)((long double)x / p * y) * p + p) % p;
}

il ll qpow(ll x, ll b, ll c) {
    ll res = 1;
    while (b) {
        if (b & 1) res = mul(res, x, c);
        x = mul(x, x, c);
        b >>= 1;
    }
    return res;
}

il bool mr(int x, long long p) {
    long long t = p - 1;
    if (qpow(x, t, p) != 1) return false;
    while (t % 2 == 0) {
        t >>= 1;
        ll h = qpow(x, t, p);
        if (h == p - 1) return true;
        else if (h != 1) return false;
    }
    return true;
}

il bool isprime(ll p) {
    if (p < 2) return 0;
    if (p == 2 || p == 3 || p == 5 || p == 7 || p == 11 || p == 13
    || p == 17 || p == 19 || p == 23 || p == 29 || p == 31 || p == 37) return 1;
    return mr(2, p) && mr(3, p) && mr(5, p) && mr(7, p) && mr(11, p) && mr(13, p) &&
    mr(17, p) && mr(19, p) && mr(23, p) && mr(29, p) && mr(31, p) && mr(37, p);
}

ll ans;

il ll rho(ll p) {
    while (1) {
        ll x, y;
        x = y = rand() % p;
        ll c = rand() % p;
        ll i = 0, j = 1, e = 1;
        while (++ i) {
            x = (mul(x, x, p) + c) % p;
            e = mul(e, Abs(y - x), p);
            if (x == y || !e) break;        
            if ((i % 127 == 0) || i == j) {
                ll d = gcd(e, p);
                if (d > 1) return d;
                if (i == j) y = x, j <<= 1;
            }
        }
    }
}

il void prho(ll p) {
    if (p <= ans) return;
    if (isprime(p)) return ans = p, void();
    ll pi = rho(p);
    while (p % pi == 0) p /= pi;
    prho(pi), prho(p);
}

int main() {
    int T;
    cin >> T;
    srand(time(0));
    while (T -- ) {
        long long p;
        scanf("%lld", &p);
        ans = 1;
        prho(p);
        if (ans == p) puts("Prime");
        else printf("%lld\n", ans);
    }
    return 0;
}