求 \(n\) 的所有质因子（\(n\ge 2\)）。

Notation

\(\mathbb{N}:=\{0,1,2,3\ldots \}\) 表示自然数。

\(\mathbb{P}:=\{2,3,5,7\ldots \}\) 表示素数集。

\(\mathbb{P}^\infty:=\{2,2^2,2^3\ldots 3,3^2,3^3\ldots \}\) 表示素数的幂次构成的集合。

\(\mathbb{Z}\) 表示整数环。

\(\pi(x)\) 表示 \([1,x]\cap\mathbb{N}\) 中的素数个数。

我会暴力

注意到，如果 \(n\not\in\mathbb{P}\)，那么其至少有一个因子 \(2\le d\le\lfloor\sqrt{n}\rfloor\)，于是枚举 \(d\) 即可。显然是 \(O(\sqrt{n})\) 的。

注意力惊人

为了方便，我们首先排除 \(n\in\mathbb{P}^\infty\) 的情况，因为 \((\log n)^{\log n}>n\)，所以这一步也就是 \(O(\log ^2n)\) 的（这里一个 \(\log\) 来自快速幂，某些更严格的 \(\text{unit cost}\) 会给出不同的代价，但是至多就是若干个 \(\log\)，根据常识，分解因数目前没有 \(P\) 的算法所以这里一定不是瓶颈）。

此时我们就可以说 \(n=\prod\limits_{i=1}^kp_i^{\alpha_i}\)，并且 \(k\ge 2,\alpha_i>0\)。

我们在 \(\mathbb{Z}/n\mathbb{Z}\) 上讨论，注意到此时 \(\sqrt1\) 不只是 \(1,-1\)，还存在别的解。具体地，假设 \(n=ab,\gcd(a,b)=1\)，那么

分别有 \(\pm1\) 两个解，用中国剩余定理可以合并出四个 \(\bmod n\) 下本质不同的解，所以一定存在 \(u\in \mathbb{Z}/n\mathbb{Z}\)，使得 \(u\not\equiv\pm 1\pmod n,u^2\equiv 1\pmod n\)。

不难注意到，\(1^2\equiv u^2\pmod n\Rightarrow n|(1-u)(1+u)\)，并且因为 \(u\mp1\not\equiv0\pmod n\)，从而只能是 \(\gcd(n,1-u),\gcd(n,1+u)>1\)，于是求解 \(\gcd\) 即可求得因数。

所以我们的目标就变成了求解 \(u\)。这里因为没有别的信息所以一般都选择随机化，从 \((1,n)\cap\mathbb{N}\) 里随机选一个数 \(x\)，如果 \(\gcd(x,n)>1\)，那么已经找到一个非平凡因子；否则 \(x\) 在 \(\bmod n\) 下存在逆元，用同样的办法找一个 \(y\)，若有 \(x^2\equiv y^2\pmod n\)，就有 \(y\in\{x,-x,ux,-ux\}\)（如果 \(n\) 不止两个素因子，这里可能的情况会更多），其中只有 \(x,-x\) 是平凡的情况，于是至少有 \(\frac 12\) 的概率可以得到 \(u\)。

Pollard Rho

非常 \(\text{naive}\) 的，我们随机取上述的 \(x\)，得到序列 \(\{x_i\}\)，重复直到存在 \(x_i^2\equiv x_j^2\pmod n\)。根据之前的分析，这个流程期望最多只需要重复 \(2\) 次。

然后我们考虑这个序列的期望长度，首先可能的平方值个数 \(\le n\)。假设已经有了 \(i\) 个这样的 \(x\) 并且还没有发生碰撞，那么此时碰撞的概率 \(\ge\frac in\)。

用 \(L\) 表示序列 \(\{x_i\}\) 的长度，就有 \(P(L>i)\le \prod\limits_{j=1}^i(1-\frac jn)\)，运用不等式 \(1+x\le e^x\) 就有

取 \(i=\sqrt n\) 就可以保证 \(O(1)\) 的成功率，这么实现和暴力就差不多了，也是 \(O(\sqrt n)\) 的。

这时，我们需要使用一点技巧来提高效率，这里我们选择重新分析复杂度（并且进行细节上的修改）。上面的说法可能比较类似把所有 \(x\) 存到一个 \(\text{hash table}\) 里，一直存直到碰撞，这样的话的确是 \(O(\sqrt n)\) 的。

但是还有一种做法类似贪心，先生成 \(k\) 个 \(x\)，这里的 \(k\) 可以任取（但是太小很显然不能保证有碰撞）。然后从其中选出 \(m\) 组 \((x_i,x_j)\)，直接判断 \(\gcd(x_i-x_j,n)\) 和 \(\gcd(x_i+x_j, n)\)。

设 \(p\) 为 \(n\) 最小的素因子，那么上面的复杂度应当是 \(O(p)\) 的，因为随机取数期望 \(O(\sqrt p)\) 次就可以得到数对 \((x_i,x_j)\) 满足 \(\gcd\) 不为 \(1\)。不过这样依然有 \(O(\sqrt n)\) 的复杂度（考虑一对孪生素数即可）。

我们的目标是在 \(\mathbb{Z}/p\mathbb{Z}\) 中找到一对碰撞的 \((x_i,x_j)\)，也就是使得 \(x_i-x_j\equiv 0\pmod p\)。这个时候我们又有了一个 \(\text{naive}\) 的想法，我们找一个有周期的随机函数，这样生成出来的随机数就是形如 \(x_1\rightarrow x_2\rightarrow\ldots x_a \rightarrow x_{a+1}\rightarrow\ldots x_{a+b}\rightarrow x_a\rightarrow\ldots\) 这样的 "\(\rho\)" 形函数。

根据大量的尝试，高于一次的多项式函数几乎都可以符合这个性质。这里真的没有证明了，因为人类目前还不会证（虽然有一些不够完美的结论，比如这里）。最常用的例子就是 \(g(x)=x^2+c\)，这里的 \(c\) 是 \((1,n)\) 中的随机数，也就是 \(x_{i+1}=g(x_i)\)，\(x_1\) 也是 \((1,n)\) 中的随机数。

如果函数是真的随机的话，证明其中循环的长度是 \(O(\sqrt p)\) 就是之前分析的碰撞概率，碰撞了显然是循环。

其实从 \([1,n]\) 中随机选择 \(x_i\) 也是对的。

Quadratic Sieve

思路

我们考虑 \(B\)-\(\text{smooth}\) 的那些数，也就是所有素因子都不超过 \(B\) 的数（\(B\) 是一个任选的常数）。

定义 \(\Psi(x;B)=\#\{1<n\le x:\forall p|n\land p\in \mathbb{P},p\le B\}\)，其中 \(\#\) 表示对集合大小计数（这大概是解析数论中比较常用的记号，反正和 \(||\) 也没区别）。

我们考虑选取一个表达式 \(Q(s)=s^2+cn\equiv s^2\pmod n\)，然后通过某种办法选出 \(k\) 个数 \(\{s_i\}\) 使得 \(Q(s_i)\) 是 \(B\)-\(\text{smooth}\) 的。

此时我们可以列出 \([1,B]\cap\mathbb{N}\) 中的 \(\pi(B)\) 个素数，每个 \(Q(s_i)\) 都可以写成

的形式，而我们已经保证 \(Q(s_i)\equiv s^2\pmod n\) 了，我们想要的就是找到一个指标集 \(S\subset\{1,2,3\ldots k\}\) 使得 \(\prod\limits_{i\in S}Q(s_i)\equiv m^2\pmod n\)，这样 \((\prod\limits_{i\in S}s_i, m)\) 就是一个 \(\bmod n\) 下平方相等的数对。

不难想到我们只关心那些不是平方的因子，也就是我们只关心 \(\alpha_{i,j}\bmod 2\) 的值。不妨把一个 \(Q(s_i)\) 的素因子指数在 \(\bmod 2\) 意义下看成一个 \(\pi(B)\) 维向量，那么显然只要有 \(\pi(B)+t\) 个 \(Q(s_i)\)，我们就可以找到至少 \(t\) 组本质不同的线性组合使得其为全 \(0\) 向量，也就是平方。（这里其实不用说本质不同，因为在 \(\bmod 2\) 的意义下系数只可能是 \(\{0,1\}\)）

根据之前的分析，这样就有至少 \(1-\frac{1}{2^t}\) 的概率成功得到 \(n\) 的一个因子。实际上一般取一个 \(10\sim20\) 的 \(t\)，然后 \(k=\pi(B)+t\) 即可。

细节与复杂度

很显然，上面的算法主要有两步，找 \(B\text{-smooth number}\) 和进行高斯消元。

高斯消元的时候因为一般可以认为 \(t=O(1)\)，所以是 \(O(\pi^3)\) 的（这里的 \(\pi\) 是 \(\pi(B)\) 的缩写）。当然如果你能注意到一个数 \(x\) 最多有 \(\log x\) 个不同的素因子，那么可以用 Berlekamp–Massey 算法做到 \(O(\pi^2\min\{\log n,\pi(B)\})\) 的复杂度。结合后面比较繁杂的分析可以说明认为 \(\log n\ll\pi(B)\) 不太会导致问题。

找 \(B\text{-smooth}\) 数主要有两步，第一个是检验，暴力做的话是试除 \(\pi(B)\) 个小素数，就是 \(O(\pi(B))\) 的。但是假如我们知道要筛的区间（这也是为什么这个算法叫筛法），比如是 \(s\in [L,R]\)，因为 \(Q(s)=s^2+cn\)，我们可以预先求出 \(s^2+cn\equiv 0\pmod p\) 的根 \(r_1,r_2\)，只有形如 \(r_1+kp,r_2+kp\) 的 \(s\) 才是 \(p\) 的倍数。

但是这样筛了我们还是要暴力试除，于是我们考虑启发式地选取，对于一个 \(Q(s)\)，我们已经知道其有若干因子 \(q_1,q_2\ldots q_l\)，那么就有

一个质朴的想法是考虑到 \(\sum \frac 1{p^2}\) 收敛，于是钦定 \(e_i=1\) 也只是需要 \(O(1)\) 倍次检验。当然实际上处于常数考虑一般是只检测那些 \(\log Q(s)\ge \alpha\sum\limits_{i=1}^l{q_i}\) 的 \(s\)，\(\alpha\) 一般取 \(1.2\sim 1.5\) 的常数。

求解 \(Q(s)\equiv 0\pmod p\) 则就是开方而已，在模意义下可以用这里提供的算法，用其中复杂度比较优秀的就可以做到 \(O((R-L)\log\pi(B))\)。（朴素的暴力显然是 \(O((R-L)\pi(B))\) 的）

第二步就是每次找有多大的概率成功，这个和 \(B\text{-smooth}\) 数的分布是有关的，设 \(m=\max\{Q(s_i)\}\)，那么成功的概率就是 \(\frac{\Psi(m;B)}{m}\)（希望你还记得 \(\Psi(m;B)\) 表示 \(\le m\) 的 \(B\text{-smooth}\) 数的数量）。上面的 \(R-L\) 期望上就是 \(\frac{m}{\Psi(m;B)}\)。

这里发挥一下注意力，不难想到，取 \(s=\sqrt{-cn}+O(1)\) 可以使得 \(Q(s)=s^2=O(\sqrt n)\)。这样比直接取 \(m=O(n)\) 要好不少。（实际上中间这些步骤我也没有全部试过，有空的人可以告诉我那些取法能比 \(\text{pollard rho}\) 优秀）

整体的复杂度合起来就好了。得到

其中 \(\pi(B)=\frac{B}{\log B}\) 是大家熟知的结果，难点在于 \(\Psi(n;B)\) 的求解。

省流：\(u=\log_B n,\Psi(n;B)=u^{-u(1+o(1))}\)。

根据省流的知识，我们再经过一番辛苦的初等代数运算便可以知道 \(B=\exp(\sqrt{2\log n\log\log n})\) 时取得最小值，总复杂度最优值为 \(O(\exp(\sqrt{(2+o(1))\log n\log\log n}))\)。

至于这个定理的证明，那下次再说，口古口古……

给出一个非常 \(\text{naive}\) 的 python 实现，常数非常巨大，即使在 python 中也喜提最劣解。

from math import floor, sqrt, log, gcd
from random import randint

def mod_num(x : int, mod : int) -> int:
  if x >= mod: x %= mod
  elif x < 0: x -= (x // mod) * mod
  return x

def mod_add(x : int, y : int, mod : int) -> int:
  z = x + y
  return z if z < mod else z - mod

def mod_mul(x : int, y : int, mod : int) -> int:
  z = x * y
  if z >= mod: z %= mod
  return z

def mod_pow(base : int, power : int, mod : int) -> int:
  base = mod_num(base, mod)
  if base == 0 or base == 1: return base
  ret = 1
  while power > 0:
    if power % 2 == 1: ret = mod_mul(ret, base, mod)
    base = mod_mul(base, base, mod)
    power //= 2
  return ret

class prime_checker:
  def __init__(self, small_prime_bound : int = 37):
    self.small_prime = [
      x for x in range(2, small_prime_bound + 1) if self.bf_checker(x)
    ]

  def bf_checker(self, x : int) -> bool:
    if x == 1: return False
    if x % 2 == 0: return x == 2
    bound = floor(sqrt(x)) + 1
    for i in range(3, bound, 2):
      if x % i == 0:
        return False
    return True

  def fermat_test(self, x : int) -> bool:
    if x == 1: return False
    if x % 2 == 0: return x == 2
    a = randint(2, x - 1)
    return mod_pow(a, x - 1, x) == 1
 
  def quadratic_test(self, x : int, a : int | None = None) -> bool:
    if x == 1: return False
    if x % 2 == 0: return x == 2
    if a is None: a = randint(2, x - 2)
    if mod_pow(a, x - 1, x) != 1: return False # fermat ctest
    d, s = x - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    t = mod_pow(a, d, x)
    if t == 1 or t == x - 1:
        return True
    for _ in range(1, s):
        t = mod_mul(t, t, x)
        if t == x - 1:
            return True
    return False

  def miller_rabin_test(self, x : int, times : int = 0) -> bool:
    if x == 1: return False
    if x in self.small_prime: return True
    for p in self.small_prime:
      if not self.quadratic_test(x, p):
        return False
    for i in range(times):
      if not self.quadratic_test(x):
        return False
    return True

def generate_random_prime(l : int, r : int) -> int:
  if r < l: raise ValueError(f"The range is illegal.")
  checker = prime_checker()
  for _ in range(10000):
    x = randint(l, r)
    if checker.miller_rabin_test(x): return x
  raise RuntimeError(f"Cannot find prime in [{l}, {r}], maybe there's no prime or need more test")

def generate_factor_hard_number(r : list[tuple[int, int]]) -> tuple[int, list[int]]:
  ret = 1
  lst = []
  for L, R in r:
    x = generate_random_prime(L, R)
    lst.append(x)
    ret *= x
  return (ret, lst)

class image_number:
  def __init__(
    self,
    mod : int,
    non_quadratic_residue : int,
    num : tuple[int, int] = (0, 0)):
    self.mod = mod
    self.non_quadratic_residue = non_quadratic_residue
    self.number = num
    # A+Bi, i^2 = non quadratic residue

  def __mul__(self, other : "image_number") -> "image_number":
    mod = self.mod
    square = self.non_quadratic_residue
    a, b = self.number
    c, d = other.number
    e = mod_num(a * c + b * d * square, mod)
    f = mod_num(b * c + a * d, mod)
    return image_number(mod, square, (e, f))

  def __add__(self, other : "image_number") -> "image_number":
    mod = self.mod
    square = self.non_quadratic_residue
    a, b = self.number
    c, d = other.number
    e, f = mod_num(a + c, mod), mod_num(b + d, mod)
    return image_number(mod, square, (e, f))

  def __pow__(self, other : int) -> "image_number":
    mod = self.mod
    square = self.non_quadratic_residue
    ret = image_number(mod, square, (1, 0))
    base = image_number(mod, square, self.number)
    while other > 0:
      if other % 2 == 1: ret = ret * base
      other //= 2
      base = base * base
    return ret

class quadratic_residue:    
  def __init__(self, safety : bool = True):
    self.safety = safety

  def is_quadratic_residue(self, x : int, p : int) -> bool:
    """判断 x 是否是 mod p 下的二次剩余"""
    return mod_pow(x, (p - 1) // 2, p) == 1  # 勒让德符号

  def _generate_non_quadratic_residue(self, a : int, p : int) -> tuple[int, int]:
    """
    找到 n 使得 n^2-a 是 mod p 的非二次剩余，返回 (n, n^2-a)
    """
    # 随机选取的成功率是 1/2
    while True:
      n = randint(1, p - 1)
      num = mod_num(n * n - a, p)
      if not self.is_quadratic_residue(num, p):
        return (n, num)

  def solve_equation(self, a : int, p : int) -> list[int]:
    """
    使用 Cipolla 算法，求出 x^2\equiv a\pmod p 在 \mathbb{Z}_p 下的所有解
    无解时返回 [] 而不是 None。
    保证返回的解有序。
    """
    if self.safety and not prime_checker().miller_rabin_test(p):
      raise TypeError(f"{p} is not a prime number!")
    if a == 0: return [0]
    if not self.is_quadratic_residue(a, p): return []
    # corner case
    r, i2 = self._generate_non_quadratic_residue(a, p)
    # i2 = i ^2, x = (r+i)^\frac{p+1}{2}
    base = image_number(p, i2, (r, 1))
    x = (base ** ((p + 1) // 2)).number[0]  # 可以保证虚部为 0
    return sorted([x, p - x])  # 还是保证有序好用一点

class quadratic_sieve:
  # 存在把 (B, B^2) 间素数存下来的优化，但是似乎没有更快
  def __init__(self, B_smooth, alpha : float = 0.7):
    self.B = B_smooth  # 找 B smooth 的数
    self.alpha = alpha  # 要求 \sum log p >= alpha * log |Q|
    checker = prime_checker()
    self.small_prime = [2] + [
      x for x in range(3, B_smooth + 1, 2) if checker.miller_rabin_test(x)
    ]
    self.smooth_num : list[tuple[int, int]] = []
    # 同时存放了 s 和 Q(s)
    self.smooth_num_p : list[int] = []
    # list[int]，只用 0/1 表示是否不是平方数，1 表示次数为奇数，用 int 压位处理以提速
    self.smooth_num_exp = []
    # 完整的素因数指数列表
    self.qr_solver = quadratic_residue(safety=False)
    self.prime_test = checker
    self.quadratic_param = (0, 0, 0)
    self.quadratic_extra = (0, 0, 0)

  def _make_quadratic_param(self, n : int, M : int) -> int:
    """
    生成 Ax^2+Bx+C 使得 B^2/4-AC=n，这样 x\sim\sqrt n 时，Q(s) 更小
    因为此时 Ax^2+Bx+C=\frac{(Ax+B)^2-n}{A}
    额外保存: D, B0, E, Q(s) ≡ (E + D*s)^2 (mod n)
    """
    """basic method"""
    # sn = floor(sqrt(n))
    # A = 1
    # B = 2 * sn
    # C = sn * sn - n
    # self.quadratic_param = (A, B, C)
    """MPQS method"""
    target_D = max(3, int(sqrt(sqrt(2.0 * n) / M)))
    while True:
      L = max(3, target_D // 2)
      R = max(L + 10, target_D * 2)
      D = generate_random_prime(L, R)
      if gcd(n, D) != 1:
        if n % D == 0:
          raise RuntimeError(f"Found factor directly: {D}")
        continue
      roots = self.qr_solver.solve_equation(n % D, D)
      if roots: break
    A = D * D
    cand = []
    for y0 in roots:
      if y0 % D == 0: continue
      t = ((n - y0 * y0) // D) % D
      inv = mod_pow((2 * y0) % D, D - 2, D)
      k = (t * inv) % D
      B0 = y0 + k * D
      for b in (B0, A - B0):
        C = (b * b - n) // A
        cand.append((abs(C), b, C))
    _, B0, C = min(cand)
    E = mod_num(B0 * pow(D, -1, n), n)
    self.quadratic_param = (A, 2 * B0, C)
    self.quadratic_extra = (D, B0, E)
    # Hensel lifting

  def _calc_quadratic_formula(self, s : int) -> int:
    A, B, C = self.quadratic_param
    return A * s * s + B * s + C
  def _msb_pos(self, x: int) -> int:
    return x.bit_length() - 1
  def _is_smooth(self, x : int) -> bool:
    for p in self.small_prime:
      while x % p == 0:
        x //= p
    return x == 1

  def _factor_exponents(self, x: int) -> tuple[bool, dict[int, int], int]:
    """
    返回 (是否 B-smooth, 指数表, 符号) sign = 0/1 表示是否含有 -1
    """
    exps = {}
    sign = 0
    if x < 0:
      sign = 1
      x = -x
    for p in self.small_prime:
      cnt = 0
      while x % p == 0:
        x //= p
        cnt += 1
      if cnt:
        exps[p] = cnt
    return (x == 1, exps, sign)

  def _roots_mod_prime(self, n: int, p: int) -> list[int]:
    A, _, C = self.quadratic_param
    D, B0, _ = self.quadratic_extra
    if p == 2: return [C & 1]
    a = A % p
    if a == 0: return []
    roots = self.qr_solver.solve_equation(mod_num(n, p), p)
    if not roots: return []
    inv_a = mod_pow(a, p - 2, p)
    r1 = mod_num((roots[0] - B0) * inv_a, p)
    r2 = mod_num((roots[-1] - B0) * inv_a, p)
    if r1 == r2: return [r1]
    return [r1, r2]

  def _find_smooth(self, n : int) -> None:
    """针对指定的 n 生成一定数量的 B-smooth number"""
    self.smooth_num.clear()
    self.smooth_num_p.clear()
    self.smooth_num_exp.clear()
    # 清空历史数据
    need_num = len(self.small_prime) + 20  # 冗余 20 个，如此可能导致对小的情况需要暴力，解不出来的概率只有 1e-6
    search_num = need_num * 50  # 经验公式，乘的是 smooth number 出现概率
    self._make_quadratic_param(n, search_num)  # 参数和搜索长度有关
    A, B, C = self.quadratic_param
    D, B0, E = self.quadratic_extra
    # MPQS 生成表达式
    # 这些是 basic method 用的
    sieve_score = [0.0] * search_num
    for p in self.small_prime:
      roots = self._roots_mod_prime(n, p)
      if not roots: continue
      lp = log(p)
      for r in roots:
        for s in range(r, search_num, p):
          sieve_score[s] += lp
    for s in range(search_num):
      q = self._calc_quadratic_formula(s)
      if q == 0: continue
      if sieve_score[s] < self.alpha * log(abs(q)): continue
      ok, exps, sign = self._factor_exponents(q)
      if not ok: continue  # 暴力检测发现不是 smooth 的
      mask = sign
      for i, p in enumerate(self.small_prime, start=1):
        if exps.get(p, 0) & 1:
          mask |= (1 << i)
      self.smooth_num.append((s, q))
      self.smooth_num_p.append(mask)
      self.smooth_num_exp.append((exps, sign))
      if len(self.smooth_num) >= need_num: return
    # 判断 Q(s) 是否生成了一个 B-smooth num，如果是，则加入 self.smooth_num 并建立对应的 smooth_num_p
    # 可能不够，此时自动填充
    offset = search_num
    while len(self.smooth_num) < need_num:
      extra = search_num
      sieve_score = [0.0] * extra
      for p in self.small_prime:
        roots = self._roots_mod_prime(n, p)
        if not roots: continue
        lp = log(p)
        for r in roots:
          start = (r - offset) % p
          for t in range(start, extra, p):
            sieve_score[t] += lp
      for t in range(extra):
        s = offset + t
        q = self._calc_quadratic_formula(s)
        if q == 0: continue
        if sieve_score[t] < self.alpha * log(abs(q)): continue
        ok, exps, sign = self._factor_exponents(q)
        if not ok: continue
        mask = sign
        for i, p in enumerate(self.small_prime, start=1):
          if exps.get(p, 0) & 1: mask |= (1 << i)
        self.smooth_num.append((s, q))
        self.smooth_num_p.append(mask)
        self.smooth_num_exp.append((exps, sign))
        if len(self.smooth_num) >= need_num: return
      offset += extra

  def _find_dependencies(self) -> list[list[int]]:
    """
    返回所有线性依赖，得到 X^2 - Y^2 \equiv 0 \pmod p
    """
    basis = {}
    deps = []
    for idx, vec0 in enumerate(self.smooth_num_p):
      vec = vec0
      comb = 1 << idx
      while vec:
        pivot = self._msb_pos(vec)
        if pivot not in basis:
          basis[pivot] = (vec, comb)
          break
        bvec, bcomb = basis[pivot]
        vec ^= bvec
        comb ^= bcomb
      if vec == 0:
        dep, j, tmp = [], 0, comb
        while tmp:
          if tmp & 1: dep.append(j)
          tmp >>= 1
          j += 1
        deps.append(dep)
    return deps

  def _build_xy(self, n: int, dep: list[int]) -> tuple[int, int]:
    X = 1
    total_exp = {p: 0 for p in self.small_prime}
    total_sign = 0
    for idx in dep:
      s, q = self.smooth_num[idx]
      D, B0, E = self.quadratic_extra
      X = mod_num(X * mod_num(E + D * s, n), n)
      exps, sign = self.smooth_num_exp[idx]
      total_sign ^= sign
      for p, e in exps.items():
        total_exp[p] += e
    # 理论上 total_sign 必须是 0
    if total_sign != 0: raise ValueError("Invalid dependence")
    Y = 1
    for p, e in total_exp.items():
      if e: Y = mod_mul(Y, mod_pow(p, e // 2, n), n)
    return X, Y

  def _bf_factorize(self, n : int) -> list[int]:
    """对 <= B^2 的数进行暴力分解"""
    ret = []
    for p in self.small_prime:
      while n % p == 0:
        ret.append(p)
        n //= p
    if n != 1: ret.append(n)  # 保证此时的 n 是 (B, B^2) 中的素数
    return ret

  def factorize(self, n : int) -> list[int]:
    """从小到大返回 n 的所有因数"""
    if n <= 0: raise ValueError(f"Cannot factorize {n} less than 1")
    if n == 1: return [1]
    elif self.prime_test.miller_rabin_test(n): return [n]
    elif n <= self.B ** 2: return sorted(self._bf_factorize(n))
    else:
      r = floor(sqrt(n))
      if r * r == n: return sorted(self.factorize(r) * 2)
    self._find_smooth(n)
    deps = self._find_dependencies()
    if not deps: raise RuntimeError("Failed to find a linear dependence.")
    for dep in deps:
      X, Y = self._build_xy(n, dep)
      g1 = gcd(mod_num(X - Y, n), n)
      g2 = gcd(mod_num(X + Y, n), n)
      for g in [g1, g2]:
        if g <= 1 or g >= n: continue
        return sorted(self.factorize(g) + self.factorize(n // g))
    raise RuntimeError("Dependency found but gcd was trivial; need another dependency.")

还有不少可以优化的地方，不过我常数太大了，测试结果可能不正常（）。想知道的可以看这里，不过其中的 \(3\) 我本地没有跑出优化效果，我的代码里面实现了其中的第 \(1,4\) 个优化。