算法伪代码#

1
x ← 2; y ← 2; d ← 1
2
while d = 1:
3
    x ← g(x)
4
    y ← g(g(y))
5
    d ← gcd(|x - y|, n)
6
if d = n:
7
    return failure
8
else:
9
    return d

核心思想#

假设要分解的整数 $n = pq$，其中 $p, q$ 都是质数，不妨设 $p \le q$。我们使用 $g(x) = (x^2 + 1) \ (\textrm{mod} \ n)$ 来生成一个伪随机数序列。

不妨设初始的 x 为 2，那么我们有一个序列为 $\{x_1 = g(2), x_2 = g(g(2)), ... x_k = g^k(2),...\}$，因为序列中的值一定是有穷的，并且序列的每一个数只依赖于前一个数，所以该序列一定会循环。

定义 $\{y_k = x_k\ (\textrm{mod} \ p), k = 1, 2, ...\}$ 序列，那么这个序列也一定会循环，并且由于 $p$ 比 $n$ 小得多，所以序列循环也一定会比 $\{x_k\}$ 序列循环得早许多。

那么我们只需要检测到这个环，假设环长为 r，他一定有$y_{r+k} - y_{k} \equiv 0 \ (\textrm{mod} \ p)$，所以使用 Floyd Cycle Detection ，并使用最小公约数来检测环的出现 (最小公约数不为 1)，此时如果最小公约数为不等于 n，那么它一定是 p 或者 q。

唯一需要注意的是环出现的时候可能最小公约数为 n，也就是上述算法中 x 等于 y。此时我们随机更换序列种子 x，并进行下一轮分解。

时间复杂度#

期望运行时间正比于 n 的最小素数的开根，此处大约为 O(3.2e4)。

递归深度#

递归深度顶多为 min(128, b - 1)，证明如下：

递归结束条件有两个，一个是 b 等于 1，或者欧拉函数为 1。

我们证明 $m \ge 2$ 时，$\phi(\phi(m)) \le m/2$ 即可，由 $\phi(m)$ 定义, $1 \le \phi(m) \le m$ 恒成立。

分以下两种情况：

1. m 为偶数，则由定义 $\phi(m) \le m / 2$，$\phi(\phi(m)) \le \phi(m) \le m / 2$
2. m 为奇数，则由定义$\phi(m)$一定为偶数，$\phi(\phi(m)) \le \phi(m) / 2 \le m / 2$

证毕。

所以如果 $\phi^{s}(m) = 1, m < 2^{64}$，那么 $s \le 64 \times 2 = 128$ 恒成立，所以递归深度最大不会超过 128。

扩展欧拉定理的应用#

扩展形式由于有一个大小比较的条件，所以需要估算当前栈迭代次幂的大小，而迭代次幂 (tetration) 增长十分快，所以值在 (2^63 - 1) 以内的 a 和 b 甚至可以枚举出来：

这些数我们可以轻易的计算出来，与模数进行比较，而其余的一定大于模数。

大整数分解#

由于 $m$ 的范围，我先打表 1e6 内所有的素数，并先将 m 中这样的质因子全部分解，假设剩余的数为 $m'$，那么 $m'$ 如果是合数，一定是两个质数 $p, q \gt 1000000$ 的乘积。

如果 $m' \le 1e12$，那么 $m'$ 一定是素数，否则我们先对 $m'$ 进行素数测试 (Miller-Rabin)，如果为合数再使用 Polland-Rho 进行分解，一旦分解一定为两个素数。

64 位模乘 trick#

Gcc/Clang 均提供内置类型 __int128_t，可以支持 128 位的运算，而如果运行在 intel x86_64 架构上，可以直接使用如下汇编

1
long mulmod(long a, long b, long m) {
2
    long res;
3
    asm("mulq %2; divq %3" : "=d"(res), "+a"(a) : "S"(b), "c"(m));
4
    return res;
5
}

算法流程#

$T(a, b, m)$:

1. 如果 a == 1 或者 b == 1，返回 a % m
2. 如果 m == 1，返回 0
3. 对 m 进行质因数分解，并计算欧拉函数 $\phi(m)$
4. 递归计算 $e = T(a, b - 1, \phi(m))$
5. 估算 $T(a, b - 1)$ 并与 $\phi(m)$ 比较大小，如果 $T(a, b - 1) \ge \phi(m)$，$e = e + \phi(m)$
6. 计算 $a ^ e \% m$

整个流程保证了所有操作数都在 64 位以内，递归深度最大为 min(128, b - 1)。

算法代码托管在 https://github.com/arkbriar/hackerrank-projecteuler/blob/master/cpp/188.cc 。