分解质因数

问题引入

给定一个正整数 ,试快速找到它的一个因数。

考虑朴素算法,因数是成对分布的, 的所有因数可以被分成两块,即 。只需要把 里的数遍历一遍,再根据除法就可以找出至少两个因数了。这个方法的时间复杂度为

时,这个算法的运行时间我们是无法接受的,希望有更优秀的算法。一种想法是通过随机的方法,猜测一个数是不是 的因数,如果运气好可以在 的时间复杂度下求解答案,但是对于 的数据,成功猜测的概率是 , 期望猜测的次数是 。如果是在 里进行猜测,成功率会大一些。我们希望有方法来优化猜测。

朴素算法与 Pollard Rho 算法引入

最简单的算法即为从 进行遍历。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
// C++ Version
list<int> breakdown(int N) {
  list<int> result;
  for (int i = 2; i * i <= N; i++) {
    if (N % i == 0) {  // 如果 i 能够整除 N,说明 i 为 N 的一个质因子。
      while (N % i == 0) N /= i;
      result.push_back(i);
    }
  }
  if (N != 1) {  // 说明再经过操作之后 N 留下了一个素数
    result.push_back(N);
  }
  return result;
}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
# Python Version
def breakdown(N):
    result = []
    for i in range(2, int(sqrt(N)) + 1):
        if N % i == 0: # 如果 i 能够整除 N,说明 i 为 N 的一个质因子。
            while N % i == 0:
                N = N // i
                result.append(i)
    if N != 1: # 说明再经过操作之后 N 留下了一个素数
        result.append(N)
    return result

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

证明 result 中均为 的素因数

首先证明元素均为 的素因数:因为当且仅当 N % i == 0 满足时,result 发生变化:储存 ,说明此时 能整除 ,说明了存在一个数 使得 ,即 (其中, 自身发生变化后遇到 时所除的数。我们注意到 result 若在 push 之前就已经有数了,为 ,那么有 N ,被除的乘积即为 )。所以 的因子。

其次证明 result 中均为素数。我们假设存在一个在 result 中的合数 ,并根据整数基本定理,分解为一个素数序列 ,而因为 ,所以它一定会在 之前被遍历到,并令 while(N % k1 == 0) N /= k1,即让 N 没有了素因子 ,故遍历到 时,N 已经没有了整除关系了。

值得指出的是,如果开始已经打了一个素数表的话,时间复杂度将从 下降到 。去 筛法 处查阅更多打表的信息。

例题:CF 1445C

而下面复杂度复杂度更低的 Pollard-Rho 算法是一种用于快速分解非平凡因数的算法(注意!非平凡因子不是素因子)。而在此之前需要先引入生日悖论。

生日悖论

不考虑出生年份,问:一个房间中至少多少人,才能使其中两个人生日相同的概率达到 ?

解:假设一年有 天,房间中有 人,用整数 对这些人进行编号。假定每个人的生日均匀分布于 天之中,且两个人的生日相互独立。

设 k 个人生日互不相同为事件 , 则事件 的概率为

至少有两个人生日相同的概率为 。根据题意可知 , 那么就有

由不等式 可得

然而我们可以得到一个不等式方程,,其中 是一个概率。

代入,解得 。所以一个房间中至少 23 人,使其中两个人生日相同的概率达到 , 但这个数学事实十分反直觉,故称之为一个悖论。

时,出现两个人同一天生日的概率将大于 。那么在一年有 天的情况下,当房间中有 个人时,至少有两个人的生日相同的概率约为

考虑一个问题,设置一个数据 ,在 里随机选取 个数( 时就是它自己),使它们之间有两个数的差值为 。当 时成功的概率是 ,当 时成功的概率是 (考虑绝对值, 可以取 ),随着 的增大,这个概率也会增大最后趋向于 1。

构造伪随机函数

我们通过 来生成一个随机数序列 ,其中 ,是一个随机的常数。

随机取一个 ,令 ,在一定范围内可以认为这个数列是“随机”的。

举个例子,设 生成的数据为

可以发现数据在 3 以后都在 11,23,31 之间循环,这也是 被称为伪随机函数的原因。

如果将这些数如下图一样排列起来,会发现这个图像酷似一个 ,算法也因此得名 rho。

Pollard-rho1

优化随机算法

最大公约数一定是某个数的约数,即 ,只要选适当的 使得 ,就可以求得一个约数 。满足这样条件的 不少, 有若干个质因子,每个质因子及其倍数都是可行的。

将生日悖论应用到随机算法中,伪随机数序列中不同值的数量约为 个。设 的最小非平凡因子,显然有 。记 ,推导可得:

于是就得到了一个新序列 (当然也可以写作 ),并且根据生日悖论可以得知序列中不同值的个数约为

假设存在两个位置 ,使得 ,这意味着 , 因此我们可以通过 获得 的一个非平凡因子。

时间复杂度分析

我们期望枚举 来分解出 的一个非平凡因子 ,因此。Pollard-rho 算法能够在 的期望时间复杂度内分解出 的一个非平凡因子,通过上面的分析可知 ,那么 Pollard-rho 算法的总时间复杂度为 。下面介绍两种实现算法,两种算法都可以在 的时间复杂度内完成。

Floyd 判环

假设两个人在赛跑,A 的速度快,B 的速度慢,经过一定时间后,A 一定会和 B 相遇,且相遇时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的 n 倍。

,每一次更新 ,只要检查在更新过程中 a、b 是否相等,如果相等了,那么就出现了环。

我们每次令 ,判断 d 是否满足 ,若满足则可直接返回 。由于 是一个伪随机数列,必定会形成环,在形成环时就不能再继续操作了,直接返回 n 本身,并且在后续操作里调整随机常数 ,重新分解。

基于 Floyd 判环的 Pollard-Rho 算法
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
// C++ Version
ll Pollard_Rho(ll N) {
  ll c = rand() % (N - 1) + 1;
  ll t = f(0, c, N);
  ll r = f(f(0, c, N), c, N);
  while (t != r) {
    ll d = gcd(abs(t - r), N);
    if (d > 1) return d;
    t = f(t, c, N);
    r = f(f(r, c, N), c, N);
  }
  return N;
}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# Python Version
def Pollard_Rho(N):
c = random.randint(0, 32767) % (N - 1) + 1
t = f(0, c, N)
r = f(f(0, c, N), c, N)
while t != r:
    d = gcd(abs(t - r), N)
    if d > 1:
        return d
    t = f(t, c, N)
    r = f(f(r, c, N), c, N)
return N

倍增优化

使用 求解的时间复杂度为 ,频繁地调用会使算法运行地很慢,可以通过乘法累积来减少求 的次数。如果 ,则有 ,并且有

我们每过一段时间将这些差值进行 运算,设 ,如果某一时刻得到 那么表示分解失败,退出并返回 本身。每隔 个数,计算是否满足 。此处取 ,可以根据实际情况进行调节。

参考实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
ll Pollard_Rho(ll x) {
  ll s = 0, t = 0;
  ll c = rand() % (x - 1) + 1;
  int step = 0, goal = 1;
  ll val = 1;
  for (goal = 1;; goal <<= 1, s = t, val = 1) {
    for (step = 1; step <= goal; ++step) {
      t = f(t, c, x);
      val = val * abs(t - s) % x;
      if ((step % 127) == 0) {
        ll d = gcd(val, x);
        if (d > 1) return d;
      }
    }
    ll d = gcd(val, x);
    if (d > 1) return d;
  }
}

例题:P4718【模板】Pollard-Rho 算法

对于一个数 ,用 Miller Rabin 算法 判断是否为素数,如果是就可以直接返回了,否则用 Pollard-Rho 算法找一个因子 ,将 除去因子 。再递归分解 ,用 Miller Rabin 判断是否出现质因子,并用 max_factor 更新就可以求出最大质因子了。由于这个题目的数据过于庞大,用 Floyd 判环的方法是不够的,这里采用倍增优化的方法。

参考实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

int t;
long long max_factor, n;

long long gcd(long long a, long long b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}

long long quick_pow(long long x, long long p, long long mod) {  // 快速幂
  long long ans = 1;
  while (p) {
    if (p & 1) ans = (__int128)ans * x % mod;
    x = (__int128)x * x % mod;
    p >>= 1;
  }
  return ans;
}

bool Miller_Rabin(long long p) {  // 判断素数
  if (p < 2) return 0;
  if (p == 2) return 1;
  if (p == 3) return 1;
  long long d = p - 1, r = 0;
  while (!(d & 1)) ++r, d >>= 1;  // 将d处理为奇数
  for (long long k = 0; k < 10; ++k) {
    long long a = rand() % (p - 2) + 2;
    long long x = quick_pow(a, d, p);
    if (x == 1 || x == p - 1) continue;
    for (int i = 0; i < r - 1; ++i) {
      x = (__int128)x * x % p;
      if (x == p - 1) break;
    }
    if (x != p - 1) return 0;
  }
  return 1;
}

long long Pollard_Rho(long long x) {
  long long s = 0, t = 0;
  long long c = (long long)rand() % (x - 1) + 1;
  int step = 0, goal = 1;
  long long val = 1;
  for (goal = 1;; goal *= 2, s = t, val = 1) {  // 倍增优化
    for (step = 1; step <= goal; ++step) {
      t = ((__int128)t * t + c) % x;
      val = (__int128)val * abs(t - s) % x;
      if ((step % 127) == 0) {
        long long d = gcd(val, x);
        if (d > 1) return d;
      }
    }
    long long d = gcd(val, x);
    if (d > 1) return d;
  }
}

void fac(long long x) {
  if (x <= max_factor || x < 2) return;
  if (Miller_Rabin(x)) {              // 如果x为质数
    max_factor = max(max_factor, x);  // 更新答案
    return;
  }
  long long p = x;
  while (p >= x) p = Pollard_Rho(x);  // 使用该算法
  while ((x % p) == 0) x /= p;
  fac(x), fac(p);  // 继续向下分解x和p
}

int main() {
  scanf("%d", &t);
  while (t--) {
    srand((unsigned)time(NULL));
    max_factor = 0;
    scanf("%lld", &n);
    fac(n);
    if (max_factor == n)  // 最大的质因数即自己
      printf("Prime\n");
    else
      printf("%lld\n", max_factor);
  }
  return 0;
}

评论