分解质因数#

引入#

给定一个正整数 $N∈N^+$ ，试快速找到它的一个非平凡因数 ( 除了1和数本身以外的因数 )。

朴素算法#

考虑朴素算法——试除法，因数是成对分布的， $N$ 的所有因数可以被分成两块，即 $[2,\sqrt{N}]$ 和 $[\sqrt{N}+1,N]$ 。只需要把 $[2, \sqrt {N}]$ ) 里的数遍历一遍，再根据除法就可以找出至少两个因数了。这个方法的时间复杂度为 $O(\sqrt {N})$ 。

最简单的算法即从 $[2,\sqrt{N}]$ 进行遍历。

1
vector<int> breakdown(int N) {
2
 vector<int> result;
3
 for (int i = 2; i * i <= N; i++) {
4
  if (N % i == 0) {
5
  // 如果 i 能够整除 N，说明 i 为 N 的一个质因子。
6
   while (N % i == 0) N /= i;
7
   result.push_back(i);
8
  }
9
 }
10
 if (N != 1) {
11
  // 说明再经过操作之后 N 留下了一个素数
12
  result.push_back(N);
13
 }
14
 return result;
15
}

我们能够证明 result 中的所有元素均为 N 的素因数。

首先证明元素均为 $N$ 的素因数：因为当且仅当 N % i == 0 满足时，result 发生变化： $i$ 储存，说明此时 $i$ 能整除 $\frac{N}{A}$ ，说明了存在一个数 $p$ 使得 $pi = \frac{N}{A}$ ，即（其中， $A$ 为 $N$ 自身发生变化后遇到 $i$ 时所除的数。我们注意到 result 若在 $push\ i$ 之前就已经有数了，为 $R_1,R_2,...,R_n$ ，那么有 $N = \frac{N}{R_1^{q_1}·R_2^{q_2}...R_n^{q_n}}$ ，被除的乘积即为 $A$ ）。所以 i 为 N 的因子。

其次证明 result 中均为素数。我们假设存在一个在 result 中的合数 $K$ ，并根据算术基本定理，分解为一个素数序列 $K = K_1^{e_1}·K_2^{e_2}·...·K_n^{e_n}$ ，而因为 $K_1 < K$ ，所以它一定会在 $K$ 之前被遍历到，并令 while(N % k1 == 0) N /= k1，即让 N 没有了素因子 $K_1$ ，故遍历到 $K$ 时，N 和 K 已经没有了整除关系了。

值得指出的是，如果开始已经打了一个素数表的话，时间复杂度将从 $O(\sqrt{N})$ 下降到 $O(\sqrt \frac{N}{ln N})$ 。具体可以去了解筛法的相关知识。【埃筛欧筛】

但是，当 $N$ 是个大整数时( $\ge 10^{18}$ )，这个算法的运行时间是不优秀的。我们期望一个更优秀的算法。

这里给出一种想法，即：通过随机的方法，猜测一个数是不是 $N$ 的因数，如果运气好可以在 $O(1)$ 的时间复杂度下求解答案，但是对于的数据 $N$ 是个大整数时( $\ge 10^{18}$ )，成功猜测的概率是 $\frac{1}{10^{18}}$ , , 期望猜测的次数是 $10^{18}$ 。如果是在 $[2, \sqrt {N}]$ ) 里进行猜测，成功率会大一些。

1
template <class T>
2
T randint(T l, T r = 0) // 生成随机数建议用<random>里的引擎和分布，而不是rand()模数，那样保证是均匀分布
3
{
4
   static mt19937 eng(time(0));
5
   if (l > r)
6
       swap(l, r);
7
   uniform_int_distribution<T> dis(l, r);
8
   return dis(eng);
9
}
10
int find_factor(int n)
11
{
12
   if (is_prime(n)) // 特判素数
13
       return n;
14
   int x, d;
15
   do
16
   {
17
       x = randint(2, n - 1); // 生成2和n-1之间的随机数
18
       d = gcd(n, x);       // 求gcd
19
   } while (d == 1);
20
   return d;
21
}

在最差的情况下， $n = p^2$ (p是质数)，这时 $[1,p^2-1]$ 里只有 1 个 p 是 n 的因数；然而，当 $x$ 取 $p,2p,3p,...,(p-1)p$ 时，都有 $gcd(x,n) = p > 1$ ，而公因数自然是 n 的因数。所以此时最差期望时间复杂度可以降到 $O(\sqrt{n}logn)$ (gcd的复杂度log，根号来自 $x$ 有约 $p - 1 ≈ \sqrt{n}$ 种满足条件的选择) ，当然，它连朴素的试除法都比不过。

因此，我们仍希望有方法来优化猜测。

Pollard Rho 算法#

生日悖论#

为了尝试优化，我们或许需要了解生日悖论。

生日悖论不算严格意义上的悖论，只是反直觉罢了。

这里给出一个简易的表述：一个房间里有23个人，则他们中有两人生日相同的概率超过一半（不考虑闰年）。关于证明：即 $\frac{365}{365} × \frac{364}{365} × \frac{363}{365} × ... \frac{365 - 22}{365} < \frac{1}{2}$ 。

为了破除旧的生日悖论的直觉，建立新的认识，我们尝试去解决生日悖论的提问：一个房间中至少多少人，才能使其中两个人生日相同的概率达到 $50\%$ ?

解：假设一年有 n 天，房间中有 k 人，用整数 1，2，…，k 对这些人进行编码。假定每个人的生日均匀分布于 n 天之中，且两个人的生日相互独立。设 k 个人生日互不相同为事件 A，则事件 A 的概率为： $P(A) = \prod \limits_{i=0}^{k-1} \frac{n-i}{n}$ ，至少有两个人生日相同的概率为 $P(\bar{A})=1-P(A)$ 。根据题意得 $P(\bar{A})\ge \frac{1}{2}$ ，那么有 $P(A)= \prod \limits_{i=0}^{k-1} \frac{n-i}{n}\le \frac{1}{2}$ 由不等式 $1+x\le e^x$ 可得： $P(A)=\prod \limits_{i=0}^{k-1}exp(-\frac{i}{n})=exp(-\frac{k(k-1)}{2n})$
因此， $exp(-\frac{k(k-1)}{2n})\le \frac{1}{2} => P(A)\le\frac{1}{2}$ 将 $n=365$ 带入得： $k\ge23$ 。所以一个房间中至少23人，使其中两个人生日相同的概率达到 $50\%$ 。

而当 $k>56,n=365$ 时，出现两个人同一天生日的概率将大于一 $99\%$ 。那么在一年有 $n$ 天的情况下，当房间中有 $\frac{1}{2}(\sqrt{8nln2 + 1}+1)≈\sqrt{2nln2}$ 个人时，至少有两个人的生日相同的概率约为 $50\%$ 。

这时，我们可以建立这样一个直觉：如果我们不断在某个范围内生成随机整数，很快便会生成到重复的数，期望大约在根号级别。精确地说，对于一个 $[1,𝑁]$ 内整数的理想随机数生成器，生成序列中第一个重复数字前期望有 $\sqrt{\frac{π𝑁}{2}}$ 个数。证明参见知乎问答

这里我们回归到原题，对于最坏的情形 $n=p^2$ ，如果我们不断在 $[1,n-1]$ 间生成随机数，那么期望在生成大约 $\sqrt{p}=n^{\frac{1}{4}}$ 个数后，可以出现 $2$ 个在模 $p$ 下相同的数（注意 $[1,n-1]$ 间的随机数模 $p$ 大致是 $[0,p-1]$ 间的随机数）。那么这 $2$ 个数的差的绝对值 $d$ ，就一定满足 $d \equiv 0\ (mod\ p)$ ，也就满足 $gcd(d,n)>1$ 。

但是知晓这件事的意义并不是那么大。正如生日悖论虽然正确，但你不一定等在班上遇到和自己生日相同的人，因为这个高概率是在两两比较下才成立的。对于这 $n^{\frac{1}{4}}$ 个数两两验证，复杂度还是立刻回到 $O(\sqrt{n}logn)$ ，这并没有什么进步，所以我们需要一些技巧。

利用最大公约数求出一个因数 || 伪随机数序列#

我们再次回顾一下这个性质： $n$ 和某个数的最大公约数一定是 $n$ 的因数，即 $\forall k ∈ N_+,gcd(k,n)|n$ ，只要选择恰当的 $k$ 使得 $1 <gcd(k,n)<n$ ，就可以求得 $n$ 的一个约数 $gcd(k,n)$ 。满足这样条件的 $k$ 不少， $n$ 有若干个质因子，每个质因子及其大部分倍数都是可行的。

我们通过 $f(x)=(x^2+c)\ mod\ n$ 来生成一个序列 ${x_i}$ ：随机取一个 $x_1$ ，令 $x_2 = f(x_1),x_3 = f(x_2),...,x_i=f(x_{i-1})$ 。其中 $c$ 是一个随机选取的常数。

举个例子，取 $n=9400,c=24,x_1=0$ ， $f(x)$ 生成数据为： $0,24,600,2824,3800,1624,5400,1224,3600,6824,8800,2824,3800,1624,...$

可以发现数据在 $x_4$ 以后都在 $2824,3800,1624,5400,1224,3600,6824,8800$ 之间循环。如果将这些数如上图一样排列起来，会发现这个图像酷似一个𝜌，算法也因此得名 $rho$ 。不难发现这是显然的，因为每个数都是由前一个数决定的，可生成的数又是有限的，那么迟早会进入循环。当然，这个循环极大可能是混循环。