【学习笔记】Miller-Rabin素性测试与Pollard-Rho大数分解

【BZOJ 3667】
第一行：CAS,代表数据组数（不大于350），以下CAS行，每行一个数字，保证在64位长整形范围内，并且没有负数。你需要对于
每个数字：第一，检验是否是质数，是质数就输出Prime
第二，如果不是质数，输出它最大的质因子是哪个。

---

检验一个数是否是质数的朴素方法是$O(\sqrt{n})$的，并不能满足这道题的需求。
其实直到现在，人们对一个数的素性判定依旧没有正确率为100%的方法。但素数简直是太重要了，我们生活中方方面面都要用到它（如著名的RSA加密算法）。所以，为了应付日常的生产需求，人们发明了一个能较快判断素数的算法。虽然正确率达不到100%，但如果我们的参数选取适当，就可以做到十分的准确。这就是今天要介绍的Miller-Rabin素性测试。

---

这个算法的准确性保证来自于这两个定理：
$$\begin{gathered} 费马小定理：对于任意素数p与非负整数a，满足： \\ {a}^{p-1}\equiv 1modp \\ 二次探测定理：若a^2\equiv 1modp且a\in [1,p-1]，则a=1或p-1 \end{gathered}$$
第二个定理的证明比较简单：
$$\begin{gathered} \because a^2\equiv 1modp \\ \therefore p|a^2-1,即p|(a+1)(a-1) \\ \because a\in [1,p-1]且p为质数 \\ \therefore a=1或p-1 \end{gathered}$$

---

终于可以开始正题了。
但是在开始前，还是得先讲一个东西。

·Miller-Rabin素性测试的前身-费马测试。
还记得费马小定理吗？$$a^{p-1}\equiv 1modp$$
虽然满足上式的数并不一定是一个素数，但在很多情况下是成立的。所以，我们可以随便选一个$a$，然后代入上式验证即可。

---

但可惜的是，上面的方法的容错率还达不到我们的需求。怎么办呢？，接下来，该正式介绍Miller-Rabin素性测试了。

首先我们对偶数特判（反正它绝对不是质数）。对于待测数$n$，我们随机选择一个底$a\in [1,n-1]$，再尝试把$n-1$写成形如下形式
$$a^{r\ast 2^s}$$
写成这样有什么用呢？我们先算$a^{2r}$，然后尝试用费马小定理。如果此时它通过了本次测试，根据二次探测定理，此时$a^r$一定为1或$n-1$。如果不是，则$n$一定不是质数。否则，我们将这个东西再平方，重复上述过程，直到这s次平方全部用完，最后再判断是否符合费马小定理即可。
实践中，选择十个随机数进行判定，其准确度已经相当高了。如果还不满意，则可以尝试使用随机的素数。
这个算法的时间复杂度仅为$O(lo{g}^{2}n)$，已经可以应付大部分情况了。
现将代码贴在下面：

inline bool check(ll a, ll n, ll r, ll s)
{
    ll ans = ksm(a,r,n), p = ans;
    for(int i = 1; i <= s; i++)
    {
        ans = mul(ans, ans, n);
        if(ans == 1 && p != 1 && p != n - 1) return 1;
        p = ans;
    }
    if(ans != 1) return 1;
    return 0;
}
inline bool miller_rabbin(ll n)
{
    if(n <= 1) return 0;
    if(n == 2) return 1;
    if(n % 2 == 0) return 0;
    ll r = n - 1, s = 0;
    while(r % 2 == 0) r >>= 1, ++s;
    for(int i = 0; i < 10; i++)
        if(check(rand() % (n-1) + 1, n, r, s))
            return 0;
    return 1;
}

---

讲了这么多，最开头那个题的第一问，不就是一道板题吗？
好了接下来我们讨论如何分解出最大的那个质因数。

根据唯一分解定理，一个数只有一种分解方法。又由于素数的定义，我们可以得出这样一个方法：先随便拆一个因数出来，然后对这两个数分别分解。在分解过程中记录最大质数即可。
那怎么做到“随便拆一个因数出来”呢？这个任务就交给Pollard-Rho算法吧。
这也是个随机算法。信仰随机神教吧qwq。每次随机生成一个数$x_1$，尝试去构造另一个数$x_2$，使$gcd(n,|x_1-x_2|)>1$，然后将$gcd(n,|x_1-x_2|)$作为$n$的一个因子，递归地处理这个问题。
那么这个 { x i } \{x_i\} {xi​}到底是怎么生成的呢？这个算法的发明人提出按照下列公式生成这个数列，直到出现重复元素为止。
$$x_i=x_{i-1}^2+c，c为随机生成的一个数。$$
接下来遇到了一个问题：如果找遍这个数列，都不满足$gcd(n,|x_i-x_2|)>1$怎么办？
换一个c和$x_1$再试一下即可。
……
我这么非，一直都找不出来那我不凉了？

---

Pollard-Rho的效率保证：生日悖论
问一个问题：在一个班中，随机选k（$k\ge 2$）个人，他们的生日相同的概率。
不难列出式子：
$$ans=1-\prod _{i=1}^k\frac{365-i+1}{365}$$
经计算发现，当$k\ge 23$时，$ans$超过了50%；当k为59时，ans接近100%。然而n多年的读书经历告诉我们并非如此。
这个东西告诉我们：如果我们采取某种组合方式，从一个样本空间中随机取样，抽中答案的概率比单次抽取的概率高。
再看回上面的算法流程。我们从“寻找$x$使$gcd(x,n)>1$”改为“寻找$x_1,x_2$使$gcd(|x_1-x_2|,n)>1$”。根据上面的结论，这种方式将降低Pollard-Rho算法的重选次数，从而将时间开销控制在可以接受的范围内。
（网上的很多博客中写的复杂度是$O(n^{\frac{1}{4}})$。）

inline ll pollard_rho(ll n, ll c)
{
    ll k = 2, x = rand() % n, y = x, p = 1;
    for(ll i = 1; p == 1; i++)
    {
        x = (mul(x,x,n) + c) % n;
        p = y > x ? y - x : x - y;
        p = gcd(n, p);
        if(i == k) y = x, k += k;
    }
    return p;
}
void solve(ll n)
{
    if(n == 1) return;
    if(miller_rabbin(n)) {ans = max(ans, n); return;}
    ll t = n;
    while(t == n) t = pollard_rho(n, rand() % (n - 1) + 1);
    solve(t), solve(n / t);
}

至此，我们便完美地解决了这个问题，现将完整代码贴在下面：

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int mn = 100005;
ll ans;
ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
inline ll mul(ll a, ll b, ll p)
{
    a = (a % p + p) % p, b = (b % p + p) % p;
    return (((a * b) - (ll)((long double)a * b / p) * p) % p + p) % p;
}
inline ll ksm(ll a, ll b, ll p)
{
    ll ret = 1;
    while(b)
    {
        if(b & 1) ret = mul(ret, a, p);
        a = mul(a, a, p), b >>= 1;
    }
    return ret;
}
inline bool check(ll a, ll n, ll r, ll s)
{
    ll ans = ksm(a,r,n), p = ans;
    for(int i = 1; i <= s; i++)
    {
        ans = mul(ans, ans, n);
        if(ans == 1 && p != 1 && p != n - 1) return 1;
        p = ans;
    }
    if(ans != 1) return 1;
    return 0;
}
inline bool miller_rabbin(ll n)
{
    if(n <= 1) return 0;
    if(n == 2) return 1;
    if(n % 2 == 0) return 0;
    ll r = n - 1, s = 0;
    while(r % 2 == 0) r >>= 1, ++s;
    for(int i = 0; i < 10; i++)
        if(check(rand() % (n-1) + 1, n, r, s))
            return 0;
    return 1;
}
inline ll pollard_rho(ll n, ll c)
{
    ll k = 2, x = rand() % n, y = x, p = 1;
    for(ll i = 1; p == 1; i++)
    {
        x = (mul(x,x,n) + c) % n;
        p = y > x ? y - x : x - y;
        p = gcd(n, p);
        if(i == k) y = x, k += k;
    }
    return p;
}
void solve(ll n)
{
    if(n == 1) return;
    if(miller_rabbin(n)) {ans = max(ans, n); return;}
    ll t = n;
    while(t == n) t = pollard_rho(n, rand() % (n - 1) + 1);
    solve(t), solve(n / t);
}
int main()
{
    ll n; int T;
    scanf("%d", &T);
    while(T--)
    {
        scanf("%lld", &n);
        ans = 0, solve(n);
        if(ans == n) puts("Prime");
        else printf("%lld\n", ans);
    }
}