题解SP4948 Integer Factorization (29 digits)

## Miller-Rabin 素性检验算法

### 欧拉函数

欧拉函数定义为 $\varphi(x)$，其意义为 $1\sim x$ 中与 $x$ 互素的正整数的个数（两端都可以取到）。

不难得到，当 $x$ 为质数时 $\varphi(x)=x-1$。

例如 $\varphi(12)=4$，因为 $\lbrack 1,12\rbrack$ 中 $\{1,5,7,11\}$ 和 $12$ 互素。

### 缩系

在上述定义中，我们称 $12$ 的缩系为 $\{1,5,7,11\}$。

不难证明 $x$ 的缩系是一个阶数为 $\varphi(x)$ 的群。

### 欧拉定理

定理内容：$\forall a,n\in\mathbb{Z}^+,a\perp n$，有 $a^{\varphi(n)}\equiv 1\pmod n$。

定理证明：

我们设集合 $R$ 为 $n$ 的缩系以及集合 $R',R''$，满足 $\forall i\in\lbrack 1,\varphi(n)\rbrack$，有 $R'_ i=aR_i\bmod n,R''_ i=aR_i$（其中 $R_i$ 代表集合 $R$ 的第 $i$ 个元素，这里认为下标从 $1$ 开始）。

如果 $R\cong R'$，则
$$
\begin{aligned}
\prod_{x\in R}x&=\prod_{x\in R'}x\\
&\equiv\prod_{x\in R''}x\\
&=\prod_{x\in R}x\times a^{\varphi(n)}\pmod n
\end{aligned}
$$
两边同时乘以 $\prod\limits_{x\in R}x$ 的逆元即可得到 $a^{\varphi(n)}\equiv 1\pmod n$。

下面是 $R\cong R'$ 的证明：

$$
\begin{array}{l}
\forall i\in\lbrack 1,\varphi(n)\rbrack\\
\because \gcd(a,n)=1,\gcd(R_i,n)=1\\
\therefore \gcd(aR_i,n)=1\\
\therefore \gcd(R'_ i,n)=1\\
\because \forall j\in\lbrack 1,\varphi(n)\rbrack\setminus\{i\},R'_ i\neq R'_ j\\
\end{array}
$$
于是我们假设 $R'_ i=R'_ j$，有 $aR_i\equiv aR_j\pmod n$。
$$
\begin{array}{l}
\because \gcd(a,n)=1\\
\therefore R_i=R_j\\
\therefore i=j
\end{array}
$$
与 $i\neq j$ 矛盾，故 $R\cong R'$。定理得证。

### 费马小定理

有了上面的欧拉定理，我们就有了费马小定理。

定理内容：$\forall a,p\in\mathbb{Z}^+,p\in\text{prime}$，有 $a^{p-1}\equiv 1\pmod p$。

定理证明：由于 $p\in\text{prime}$，$\varphi(p)=p-1$。将这个式子带入欧拉定理中即可。

### 费马素性检验法

很遗憾，费马小定理的逆定理并不一定成立。但是它仍然可以筛掉很多合数。

所以我们可以随机选取多次 $a$，极大程度上降低错误概率。

### 二次探测定理

如果 $p$ 是奇素数且 $x^2\equiv 1\pmod p$，则 $x\equiv\pm 1\pmod p$。

证明：移项然后平方差公式分解。

### Miller 的研究

但是有一类数，你无论如何选取 $a$ 费马素性检验法都会判断成素数，我们称这类数为[卡迈克尔数（OEIS-A006931）](https://oeis.org/A006931)。

卡迈克尔数有这么一条性质：如果 $n$ 为卡迈克尔数，则 $2^n-1$ 也是卡迈克尔数。于是卡迈克尔数的个数是无穷的。

于是有了 Miller-Rabin 素性检验法。

#### 算法内部逻辑

将 $a^{p-1}\equiv 1\pmod p$ 中的 $p-1$ 分解成 $u\times 2^t$，其中 $t$ 应该尽量大。

然后 $\forall s\in\lbrack 0,t\rbrack$ 进行检验：

- 如果答案是 $p-1$ 则直接判定为素数。
- 如果答案不是 $\pm 1$ 则直接判定为合数。

基于广义黎曼猜想，Miller 取的底数为 $\lbrack 2,\text{min}(n-2,\lfloor 2\ln^2 n\rfloor)\rbrack$ 中的所有正整数。

如果广义黎曼猜想成立则 Miller-Rabin 为目前时间复杂度最低的确定性素数判定算法（即能准确判断一个数是否是素数），时间复杂度为 $\Theta(\log^5 n)$。

### Rabin 的研究

但是如果广义黎曼猜想不成立怎么办？（当然，检验这么多底数也没用）

于是 Rabin 想到了一种随机化算法：随机底数然后判定。这样的出错几率是 $\leq 4^{-k}$ 的，其中 $k$ 是不同的底数个数。

但是每次随机感觉太慢了，而且很容易就会出错。

有一天，数学家们发现了这样两个集合：$\{2, 7, 61\}$ 和 $\{2, 325, 9375, 28178, 450775, 9780504, 1795265022\}$，分别可以准确判断 $2^{32}$ 和 $2^{64}$ 内的所有数的素性。

当然，$\{2,3,5,7,11,13,17,19,23,29,31,37\}$ 可以准确判断 $2^{78}$ 内所有数的素性。

我更习惯用 $2$ 初判然后底数用 $\{3,7,11,41,43,61\}$。

代码：
```cpp
// mint 改编自 hly1204 的 LongMontgomeryModInt
inline bool mr(mint x) {
    const mint one(1), mone(-1);
    if (x.pow(mone.val()) != one) return 0;
    __uint128_t y = mone.val();
    mint z(0);
    while (!(y & 1)) {
        y >>= 1;
        z = x.pow(y);
        if (z != one && z != mone) return false;
        if (z == mone) return true;
    }
    return true;
}
inline bool prime(__uint128_t x) { 
    if (x < 2) return false;
    if (x == 2 || x == 3 || x == 7 || x == 11 || x == 41 || x == 43 || x == 61) return true;
    if (x % 2 == 0 || x % 3 == 0 || x % 7 == 0 || x % 11 == 0 || x % 41 == 0 || x % 43 == 0 || x % 61 == 0) return false;
    mint::set_mod(x);
    return mr(3) && mr(7) && mr(11) && mr(41) && mr(43) && mr(61);
}
```
## Pollard $\rho$ 算法

### 生日悖论

根据组合数可以得到只要一个房间里有 $\color{red}23$ 个人就有 $50\%$ 的概率有两个人的生日会相同。

当 $k>\color{red}56$，$n=365$ 时，出现两个人同一天生日的概率将大于 $\color{red}99\%$。

考虑一个问题，设置一个数据 $n$，在 $[1,1000]$ 里随机选取 $i$ 个数（$i=1$ 时就是它自己），使它们之间有两个数的差值为 $k$。当 $i=1$ 时成功的概率是 $\frac{1}{1000}$，当 $i=2$ 时成功的概率是 $\frac{1}{500}$（考虑绝对值，$k_2$ 可以取 $k_1-k$ 或 $k_1+k$），随着 $i$ 的增大，这个概率也会增大，最后趋向于 $100\%$。

### 随机数

我们设 $f(x)=(x^2+c)\bmod n$，其中 $n$ 是待分解的整数；

再设数列 $x$ 的递推式为 $x_i=f(x_{i-1})$。初值为 $x_0=1$。

这样可以获得足够的随机性，并且根据生日悖论，只要生成 $\sqrt{n}$ 个数就可以找到环。

设 $m$ 为 $n$ 的最小非平凡因子，则一定有 $1<m<n$。

于是我们再设数列 $y$，满足 $\forall i\in\mathbb{Z}^+,y_i=x_i\bmod m$。如果 $\exists i,j$ 使得 $i\neq j,x_i\neq x_j,y_i=y_j$，我们就找到了 $n$ 的最小非平凡因子 $m$。

再应用一次生日悖论可以得到时间复杂度为 $O(\sqrt{m})\leq\Theta(n^\frac{1}{4})$。

### Floyd 判环算法

假设两个人同时同向同地开始赛跑，A 的速度快，B 的速度慢，经过一定时间后，A 一定会追上 B，且 B 第一次被追上时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的 $1$ 倍。

于是我们可以写出这样的代码：
```cpp
#ifdef __USE_MY_RAND__
static unsigned long long seed;
unsigned long long u64_hasher(unsigned long long x) {
    return x *= 0x1952BB648B74B19Eull, x ^= x >> 47,
           x *= 0x1952BB648B74B19Eull, x ^= x >> 47,
           x * 0x1952BB648B74B19Eull;
}
unsigned long long rng() { return seed += u64_hasher(seed); }
#else
static std::mt19937_64 rng;
#endif
inline uint64_t rho(uint64_t value) {
    if (value % 2 == 0) return 2;
    if (value % 3 == 0) return 3;
    if (value % 5 == 0) return 5;
    if (value % 7 == 0) return 7;
    mint::set_mod(value);
    mint x, y;
    mint c;
    while (true) {
        c = rng();
        y = c;
        x = c * c + c;
        while (x != y) {
            mint z((uint64_t)abs((__int128)y.val() - (__int128)x.val()));
            if (uint64_t g = std::__gcd(z.val(), value); g > 1 && g < value) return g;
            y = y * y + c;
            x = x * x + c;
            x = x * x + c;
        }
    }
}
```

### 倍增优化

上面的做法要调用很多次 $\gcd$，导致算法乘上了一个 $\log$。

其实这个 $\log$ 是可以避免的（和常数抵消）。

修改上述结论：假设两个人同时同向同地开始赛跑，A 的速度快，B 的速度慢，经过一定时间后，A 一定会追上 B，且此时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的 $k(k\in\mathbb{Z}^+)$ 倍。（无论之前 B 有没有被 A 追上过）

我们每过一段时间将这些差值进行 $\gcd$ 运算。设 $z=\prod|x_i-y_i|\bmod n$，由欧几里得算法得到 $\gcd(ab\bmod n,n)=\gcd(ab,n)$，所以 $\gcd(z,n)=\gcd(\prod|x_i-y_i|,n)$。

如果某一时刻得到 $z=0$ 那么表示分解失败，退出并返回 $n$ 本身。每隔 $k$ 个数，计算是否满足 $1<\gcd(z,n)<n$。此处取 $k=127$，可以根据实际情况进行调节。

更正：由拉梅定理得到 $\gcd$ 的辗转相除次数 $\leq 5\log_{10}\text{min}(n,m)=5\log_{10}2\times\log\text{min}(n,m)$，所以对于 $64$ 位整型可取 $k=100$。

### 代码

用几个小素数判一下效率会有显著提升（？）。

附带 `readu128()`、`putu128()`。这份代码中采用倍增优化，$k=32768$。

这里有一个小技巧：只改变 $x$ 在环上的位置，不改变 $y$ 在环上的位置。可以证明这和一个前进两格、另一个前进一格是等价的。
```cpp
#ifdef __USE_MY_RAND__
static unsigned long long seed;
unsigned long long u64_hasher(unsigned long long x) {
    return x *= 0x1952BB648B74B19Eull, x ^= x >> 47,
           x *= 0x1952BB648B74B19Eull, x ^= x >> 47,
           x * 0x1952BB648B74B19Eull;
}
unsigned long long rng() { return seed += u64_hasher(seed); }
#else
static std::mt19937_64 rng;
#endif
inline __uint128_t rho(__uint128_t value) {
    if (value % 2 == 0) return 2;
    if (value % 3 == 0) return 3;
    if (value % 5 == 0) return 5;
    if (value % 7 == 0) return 7;
    mint::set_mod(value);
    mint x, y, z, c, tmp;
    uint64_t i, j;
    __uint128_t g;
    while (true) {
        y = x = rng();
        z = 1;
        c = rng();
        i = 0;
        j = 1;
        while (++i) {
            x = x * x + c;
            tmp = (x.val() > y.val()) ? (x - y) : (y - x);
            z *= tmp;
            if (x == y) break;
            if (!z.val()) return tmp.val();
            if (!(i & 32767) || i == j) {
                g = std::__gcd(z.val(), value);
                if (g > 1ull) return g;
                if (i == j) {
                    y = x;
                    j <<= 1;
                }
            }
        }
    }
}
inline auto crack(__uint128_t value) {
    struct _RetType { __uint128_t base; uint64_t power; };
    std::vector<_RetType> ret;
    if (prime(value)) {
        ret.push_back({value, 1ul});
        return ret;
    }
    _RetType node{2, 0};
    while (value % 2 == 0) ++(node.power), value >>= 1;
    if (node.power) ret.push_back(node);
    auto dfs = [&](auto self, __uint128_t v) {
        if (v <= 1) return;
        if (!prime(v)) {
            __uint128_t tmp = rho(v);
            self(self, tmp);
            self(self, v / tmp);
        }
        else {
            auto pos = std::find_if(ret.begin(), ret.end(), [=](const auto& x) { return x.base == v; });
            if (pos == ret.end()) ret.push_back({v, 1ul});
            else ++(pos->power);
        }
    };
    dfs(dfs, value);
    std::sort(ret.begin(), ret.end(), [=](const auto& x, const auto& y) { return x.base < y.base; });
    return ret;
}
inline __uint128_t readu128() {
    __uint128_t ret = 0;
    char c = getchar();
    while (!isdigit(c)) c = getchar();
    while (isdigit(c)) ret = ret * 10 + (c & 15), c = getchar();
    return ret;
}
inline void putu128(__uint128_t x) {
    static char buf[45];
    static int curpos;
    curpos = 0;
    do {
        buf[curpos++] = x % 10;
        x /= 10;
    } while (x);
    while (curpos--) putchar(buf[curpos] | 48);
}
```

## 题外话

这个 Pollard $\rho$ 模板是我有一次打模拟赛的时候想在暴力基础上加速分解素因数然后写的。那个是挂掉的板子。

我也是被这个 `random_device` 气疯了。在 [LOJ](https://loj.ac/p/6466) 上一会 3000ms，一会 5600ms。

[![](https://cdn.luogu.com.cn/upload/image_hosting/nmvfpawo.png)](http:// "原 LOJ#6466，现在没了")

## 参考文献

1. [数论四大定理之欧拉定理](https://www.jianshu.com/p/4fadf3e45764)
2. [素数 - OI-Wiki](https://oi-wiki.org/math/number-theory/prime/)
2. [分解素因数 - OI-Wiki](https://oi-wiki.org/math/number-theory/pollard-rho/)