这玩意我其实搞了很久。

Upd on 2024.7.26：这篇文章是我在 2023.4.29 编写的，而我从 2023.2.1 就开始搞分解，一直到现在，详见 该专栏。

Miller-Rabin && Pollard $\rho$

见 这篇文章。

Pollard P-1 算法

设要分解的数为 $n$。

如果 $n-1$ 的所有素因子都是小于等于 $B$ 的，则直接计算 $\gcd(2^{B!}-1,n)$ 就可以求出 $n$ 的一个素因子。

原理：将费马小定理改成同余形式：$a^{t(p-1)}-1=kp$。也就是说，如果指数 $exp$ 是 $p-1$ 的倍数，那么 $p|a^{exp}-1$。

对于合数 $n=pq$，如果我们能找到一个指数 $B$，使得对于某一个底数 $a$，有 $p\mid a^B-1$ 且 $q\nmid a^B-1$，那么计算 $\gcd(a^B-1,n)$ 就分解了 $n$。

这里 $B$ 就是 $p$ 的最大素因子。也即如果 $p$ 是 $B-\text{smooth}$ 的且 $q$ 不是，那么就可以选择这样的 $B$。反之，找到的就是 $q$。

复杂度证明：不会。

__uint128_t pollard_pm1(__uint128_t value) {
    if (__uint128_t ret = sqrtl(1.0 * value) + 1e-6; ret * ret == value) return ret;
    mint::set_mod(value, false);
    mint a(2);
    for (uint64_t j = 2;; j++) {
        a = a.pow(j);
        __uint128_t ret = std::__gcd((a - 1).val(), value);
        if (ret > 1) return ret;
    }
}

SQUFOF

即 Shanks 的平方型分解算法。当要分解的数的两个素因数比较近时这个算法跑得很快。

不知道如果两个素因子差为 $m$，那么时间复杂度是否为 $O(\sqrt{m})$。反正这个算法运算量有一个精确上界 $3\sqrt{2\sqrt{n}}$（还要再乘以一个常数，视实现决定）。

__uint128_t squfof(__uint128_t value) {
    __uint128_t s = sqrtl(1.0l * value) + 1e-6;
    while (s * s > value) s--;
    if (s * s == value) return s;
    __int128 d, po, p, p_prev, q, q_prev, q_, b, r;
    long long l, b_, i;
    int k = 0;
    l = 2 * sqrtl(2.0l * s);
    b_ = 3 * l;
    while (++k) {
        d = k * value;
        po = p_prev = p = sqrtl(1.0l * d);
        q_prev = 1;
        q = d - po * po;
        while (q < 0) {
            po--;
            p_prev--;
            p--;
            q = d - po * po;
        }
        if (!q) return gcd(value, k);
        for (i = 2; i < b_; i++) {
            b = (po + p) / q;
            p = b * q - p;
            q_ = q;
            q = q_prev + b * (p_prev - p);
            r = sqrtl(1.0 * q) + 1e-6;
            if (!(i & 1) && r * r == q) break;
            q_prev = q_;
            p_prev = p;
        }
        if (i >= b_) continue;
        b = (po - p) / r;
        p_prev = p = b * r + p;
        q_prev = r;
        q = (d - p_prev * p_prev) / q_prev;
        i = 0;
        do {
            b = (po + p) / q;
            p_prev = p;
            p = b * q - p;
            q_ = q;
            q = q_prev + b * (p_prev - p);
            q_prev = q_;
            i++;
        } while (p != p_prev && i < b_);
        if (i >= b_) continue;
        r = gcd(value, q_prev);
        if (r != 1) return r;
    }
    return 0;
}

Dixon 的随机平方算法

如果我们手上有一系列（不同的）$x^2\equiv y^2\pmod n$ 的式子，那么每一个都有 $50\%$ 的概率分解 $n$。

随机一个 $\lbrack 1,n\rbrack$ 的 $r$，令

$$g(r)\equiv r^2\pmod n$$

然后用 $B$ 内的素数 $g(r)$。如果 $g(r)$ 不是 $B-\text{smooth}$ 的，那么就重新选一个 $r$。

对于每一个 $r$ 可以最终得到这样的式子：

$$\begin{array}{l} z&exp_{B_1}&exp_{B_2}&\cdots&exp_{B_{|B|}} \end{array} $$

其中 $exp_{B_i}$ 是 $g(r)$ 除以 $B_i$ 的余数（$\bmod 2$）。

将所有的这些关系组合，就可以得到一个 $0-1$ 矩阵。对这个矩阵进行高斯消元。即选出一些关系，使它们异或起来为

$$\begin{array}{l} \prod\limits_{r\in A}r&0&0&\cdots&0 \end{array} $$

通过这些式子就可以分解 $n$。

例：分解 $998244359987710471$。

选取 $B=\{2,3,5,7,11,13\}$。

矩阵如下：

$$\begin{array}{l} &2&3&5&7&11&13\\ 895200474254967936&0&1&0&1&1&0\\ 470344441344285771&1&0&1&0&0&1\\ 741511338075501596&1&1&1&1&1&1\\ \end{array} $$

将这些关系的左右乘起来得到 $258641729810462066^2\equiv 30030^2$。计算 $\gcd(258641729810462066+30030,998244359987710471)$ 得到 $998244359987710471$ 的一个非平凡因子 $1000000007$。

二次筛

前置知识：二次剩余

上面这个算法的问题在于如何正确选取 $r$ 和如何构造 $g(r)$。幸运的是，我们可以取 $g(r)=(r+\lfloor x\rfloor)^2-x\sim 2r\lfloor x\rfloor$。

通过实践可以发现，$g(r)$ 是 $B-\text{smooth}$ 的。

然后就是要在模意义下求解 $p\mid g(x)$，即 $x^2\equiv a\pmod p$。这可以使用 Cipolla 算法。

这就是单多项式二次筛（SPQS，Simple Polynomial Quadratic Sieve）。

优化们：

1. 模数只取 $p\bmod 4=3$ 的情况。这样可以避免使用 Cipolla 算法。
2. 我们不仅可以选用正值的 $r$，也可以选用负值的 $r$，只需选取偶数个负值的 $r$ 符号即可抵消，所以我们在矩阵中使用一列记录是否为负即可。
3. 可能会有部分数不 $B-\text{smooth}$，但是分解后剩余的数在 $(B_{|B|},B_{|B|}^2)$，这必然剩的是一个大质数 $p$。对于这些数，我们并不能直接在表格中加入 $p$ 列，但是如果有相同 $p$ 的多个数 $\{r_1,r_2,\cdots,r_m\}$，我们就可以将 $r_1r_2,r_1r_3,\cdots,r_1r_m$ 加入，这样 $p$ 就被平方了。
4. 使用多个多项式。

接下来重点讲第 4 个优化。设 $g(r)=ar^2+2br+c$，则 $ag(r)=a^2r^2+2abr+ac$，若 $b^2-ac=n$，我们就可以有 $ag(r)=(ar+b)^2-n$，则 $(ar+b)^2\equiv ag(r)\pmod n$。若 $a$ 为平方数，我们就可以考虑筛这些新的 $g$。我们可以取小质数 $p$，令 $a=p^2,b^2\equiv n\pmod p$ 并解出 $c$。

这就是重多项式二次筛（MPQS，Multiple Polynomial Quadratic Sieve）。

时间复杂度 $\exp(\sqrt{(1+o(1))\ln n\ln\ln n})$。

---

一份集成实现见 提交记录 #1720023 - LibreOJ。（没有集成 ECM）

参考资料：

1. 整数分解
2. SQUFOF
3. 二次筛