模板 - Miller-Rabin 素性检验与 Pollard-Rho 算法

发表于 更新于 分类于 ， 本文字数： 159 阅读时长 ≈ 1 分钟

包含以 \(O(k\log^3n)\) 的时间复杂度对一个数进行素性检验的 Miller-Rabin 算法和期望以 \(O(n^\frac{1}{4})\) 的时间复杂度找出一个数最大素因子的 Pollard-Rho 算法

原理略

  • 仅在 GCC 下测试过
  • 求 gcd 的函数若在低于 C++17 下会调用非标准函数

https://cplib.tifa-233.com/src/code/nt/is_prime.hpp, https://cplib.tifa-233.com/src/code/nt/pfactors.hpp 存放了笔者对该算法 / 数据结构的最新实现，建议前往此处查看相关代码

Show code

Pollard_rho.hppview raw
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
89
90
91
92
93
94
95
96
97
98
99
100
101
namespace Pollard_rho {
using data_type = uint64_t;

using signed_data_t = make_signed_t<data_type>;
using unsigned_data_t = make_unsigned_t<data_type>;

constexpr signed_data_t
mul_mod(signed_data_t a, signed_data_t b, unsigned_data_t mod) {
  signed_data_t d = floor(1.0l * a * b / mod + 0.5l), ret = a * b - d * mod;
  return ret < 0 ? ret + mod : ret;
}

constexpr data_type
pow_mod(unsigned_data_t a, unsigned_data_t b, const unsigned_data_t &mod) {
  data_type res(1);
  a %= mod;
  for (; b; b >>= 1, a = mul_mod(a, a, mod))
    if (b & 1) res = mul_mod(res, a, mod);
  return res;
}

constexpr unsigned_data_t gcd(unsigned_data_t m, unsigned_data_t n) {
#if __cplusplus >= 201703L

#return std::gcd(m, n);

#else
  return std::__gcd(m, n);
#endif
}

namespace Primetest_miller_rabin {
// for int32
// constexpr unsigned_data_t bases[] = {2, 7, 61};

// for int64
constexpr unsigned_data_t bases[] = {
  2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};

constexpr bool is_prime(unsigned_data_t n) {
  if (n <= 1) return false;
  for (unsigned_data_t a : bases)
    if (n == a) return true;
  if (n % 2 == 0) return false;

  unsigned_data_t d = n - 1;
  while (d % 2 == 0) d /= 2;
  for (unsigned_data_t a : bases) {
    unsigned_data_t t = d, y = pow_mod(a, t, n);
    while (t != n - 1 && y != 1 && y != n - 1) {
      y = mul_mod(y, y, n);
      t <<= 1;
    }
    if (y != n - 1 && t % 2 == 0) return false;
  }
  return true;
}
}  // namespace Primetest_miller_rabin
using Primetest_miller_rabin::is_prime;

default_random_engine e(time(nullptr));
uniform_int_distribution<data_type> u;

#define Rand(x) (u(e) % (x) + 1)

data_type pollard_rho(const data_type x, const data_type y) {
  data_type t = 0, k = 1;
  signed_data_t v0 = Rand(x - 1), v = v0;
  data_type d, s = 1;
  while (true) {
    v = (mul_mod(v, v, x) + y) % x;
    s = mul_mod(s, abs(v - v0), x);
    if (!(v ^ v0) || !s) return x;
    if (++t == k) {
      if ((d = gcd(s, x)) ^ 1) return d;
      v0 = v, k <<= 1;
    }
  }
}

void resolve(data_type x, data_type &ans) {
  if (!(x ^ 1) || x <= ans) return;
  if (is_prime(x)) {
    if (ans < x) ans = x;
    return;
  }
  data_type y = x;
  while ((y = pollard_rho(x, Rand(x))) == x)
    ;
  while (x % y == 0) x /= y;
  resolve(x, ans);
  resolve(y, ans);
}
data_type get_max_prime_divisor(data_type x) {
  data_type ans = 0;
  resolve(x, ans);
  return ans;
}
}  // namespace Pollard_rho
using Pollard_rho::get_max_prime_divisor;
using Pollard_rho::is_prime;