Pollard's rho algorithm
(number-theory/pollard-rho.hpp)

• View this file on GitHub
• Last update: 2026-07-28 20:47:27+09:00
• Include: #include "number-theory/pollard-rho.hpp"

64 bit 整数の素因数分解と約数列挙を行う．

• PollardRho::factorize(N)：$N$ の素因数と指数の組を素因数の昇順に並べた列を返す．$1\leq N\leq 10^{18}$ とし，$N=1$ のとき空列を返す．
• PollardRho::divisors(N)：$N$ の正の約数を昇順に並べた列を返す．$1\leq N\leq 10^{18}$ とし，$N=1$ のとき {1} を返す．

Miller-Rabin 素数判定と Pollard’s rho algorithm を用いる．factorize の期待計算量は $O(N^{1/4}\log N)$，divisors は約数の個数を $\tau(N)$ として期待計算量 $O(N^{1/4}\log N+\tau(N)\log\tau(N))$ である．

Depends on

• Miller-Rabin 素数判定 (number-theory/miller-rabin.hpp)

Required by

• 原始根 (number-theory/primitive-root.hpp)

Verified with

• verify/number-theory/LC_factorize.test.cpp
• verify/number-theory/LC_primitive_root.test.cpp
• verify/number-theory/UNIT_pollard_rho_divisors.test.cpp

Code

#pragma once

#include "number-theory/miller-rabin.hpp"

namespace PollardRho {
using ll = long long;
using u64 = uint64_t;

namespace internal {
u64 random() {
  static u64 x = 0x243f6a8885a308d3ULL;
  x ^= x << 7;
  x ^= x >> 9;
  return x;
}

u64 find_factor(u64 n) {
  if (n % 2 == 0) return 2;
  if (n % 3 == 0) return 3;

  while (true) {
    u64 y = random() % (n - 1) + 1;
    u64 c = random() % (n - 1) + 1;
    u64 m = 128, g = 1, r = 1, q = 1, x = 0, z = 0;
    auto f = [&](u64 v) {
      return (MillerRabin::internal::multiply_mod(v, v, n) + c) % n;
    };
    while (g == 1) {
      x = y;
      for (u64 i = 0; i < r; i++) y = f(y);
      for (u64 k = 0; k < r && g == 1; k += m) {
        z = y;
        for (u64 i = 0; i < min(m, r - k); i++) {
          y = f(y);
          u64 d = x > y ? x - y : y - x;
          q = MillerRabin::internal::multiply_mod(q, d, n);
        }
        g = gcd(q, n);
      }
      r <<= 1;
    }
    if (g == n) {
      do {
        z = f(z);
        u64 d = x > z ? x - z : z - x;
        g = gcd(d, n);
      } while (g == 1);
    }
    if (g != n) return g;
  }
}

void factorize(u64 n, vector<u64>& factors) {
  if (n == 1) return;
  if (MillerRabin::is_prime(n)) {
    factors.push_back(n);
    return;
  }
  u64 d = find_factor(n);
  factorize(d, factors);
  factorize(n / d, factors);
}
};  // namespace internal

vector<pair<ll, int>> factorize(ll n) {
  assert(n >= 1);
  vector<u64> factors;
  internal::factorize(n, factors);
  sort(factors.begin(), factors.end());

  vector<pair<ll, int>> ret;
  for (u64 p : factors) {
    if (ret.empty() || ret.back().first != (ll)p)
      ret.emplace_back(p, 1);
    else
      ret.back().second++;
  }
  return ret;
}

vector<ll> divisors(ll n) {
  vector<ll> ret{1};
  for (auto [p, e] : factorize(n)) {
    size_t size = ret.size();
    ll q = 1;
    while (e--) {
      q *= p;
      for (size_t i = 0; i < size; i++) ret.push_back(ret[i] * q);
    }
  }
  sort(ret.begin(), ret.end());
  return ret;
}
};  // namespace PollardRho

/**
 * @brief Pollard's rho algorithm
 * @docs docs/number-theory/pollard-rho.md
 */

#line 2 "number-theory/pollard-rho.hpp"

#line 2 "number-theory/miller-rabin.hpp"

namespace MillerRabin {
using u64 = uint64_t;
using u128 = __uint128_t;

namespace internal {
u64 multiply_mod(u64 a, u64 b, u64 mod) { return u128(a) * b % mod; }

u64 power_mod(u64 a, u64 n, u64 mod) {
  u64 ret = 1;
  while (n) {
    if (n & 1) ret = multiply_mod(ret, a, mod);
    a = multiply_mod(a, a, mod);
    n >>= 1;
  }
  return ret;
}
};  // namespace internal

bool is_prime(long long n) {
  if (n < 2) return false;
  u64 x = n;
  for (u64 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
    if (x % p == 0) return x == p;
  }

  int s = __builtin_ctzll(x - 1);
  u64 d = (x - 1) >> s;
  for (u64 a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
    if (a % x == 0) continue;
    u64 y = internal::power_mod(a % x, d, x);
    if (y == 1 || y == x - 1) continue;
    bool composite = true;
    for (int r = 1; r < s; r++) {
      y = internal::multiply_mod(y, y, x);
      if (y == x - 1) {
        composite = false;
        break;
      }
    }
    if (composite) return false;
  }
  return true;
}
};  // namespace MillerRabin

/**
 * @brief Miller-Rabin 素数判定
 * @docs docs/number-theory/miller-rabin.md
 */
#line 4 "number-theory/pollard-rho.hpp"

namespace PollardRho {
using ll = long long;
using u64 = uint64_t;

namespace internal {
u64 random() {
  static u64 x = 0x243f6a8885a308d3ULL;
  x ^= x << 7;
  x ^= x >> 9;
  return x;
}

u64 find_factor(u64 n) {
  if (n % 2 == 0) return 2;
  if (n % 3 == 0) return 3;

  while (true) {
    u64 y = random() % (n - 1) + 1;
    u64 c = random() % (n - 1) + 1;
    u64 m = 128, g = 1, r = 1, q = 1, x = 0, z = 0;
    auto f = [&](u64 v) {
      return (MillerRabin::internal::multiply_mod(v, v, n) + c) % n;
    };
    while (g == 1) {
      x = y;
      for (u64 i = 0; i < r; i++) y = f(y);
      for (u64 k = 0; k < r && g == 1; k += m) {
        z = y;
        for (u64 i = 0; i < min(m, r - k); i++) {
          y = f(y);
          u64 d = x > y ? x - y : y - x;
          q = MillerRabin::internal::multiply_mod(q, d, n);
        }
        g = gcd(q, n);
      }
      r <<= 1;
    }
    if (g == n) {
      do {
        z = f(z);
        u64 d = x > z ? x - z : z - x;
        g = gcd(d, n);
      } while (g == 1);
    }
    if (g != n) return g;
  }
}

void factorize(u64 n, vector<u64>& factors) {
  if (n == 1) return;
  if (MillerRabin::is_prime(n)) {
    factors.push_back(n);
    return;
  }
  u64 d = find_factor(n);
  factorize(d, factors);
  factorize(n / d, factors);
}
};  // namespace internal

vector<pair<ll, int>> factorize(ll n) {
  assert(n >= 1);
  vector<u64> factors;
  internal::factorize(n, factors);
  sort(factors.begin(), factors.end());

  vector<pair<ll, int>> ret;
  for (u64 p : factors) {
    if (ret.empty() || ret.back().first != (ll)p)
      ret.emplace_back(p, 1);
    else
      ret.back().second++;
  }
  return ret;
}

vector<ll> divisors(ll n) {
  vector<ll> ret{1};
  for (auto [p, e] : factorize(n)) {
    size_t size = ret.size();
    ll q = 1;
    while (e--) {
      q *= p;
      for (size_t i = 0; i < size; i++) ret.push_back(ret[i] * q);
    }
  }
  sort(ret.begin(), ret.end());
  return ret;
}
};  // namespace PollardRho

/**
 * @brief Pollard's rho algorithm
 * @docs docs/number-theory/pollard-rho.md
 */

Back to top page