Prime factorization (Pollard's rho algorithm)
(Mylib/Number/Prime/pollard_rho.cpp)
Operations
Requirements
Notes
Problems
References
Depends on
Verified with
Code
#pragma once
#include <algorithm>
#include <cassert>
#include <numeric>
#include <optional>
#include <utility>
#include <vector>
#include "Mylib/Misc/int128.cpp"
#include "Mylib/Number/Prime/miller_rabin.cpp"
namespace haar_lib {
namespace pollard_rho_impl {
int128_t f(int128_t x) {
return x * x + 1;
}
std::optional<int64_t> rho(int64_t n) {
int64_t x = 2, y = 2, d = 1;
while (d == 1) {
x = f(x) % n;
y = f(f(y) % n) % n;
d = std::gcd(std::abs(x - y), n);
if (d == n) return {};
}
return {d};
}
} // namespace pollard_rho_impl
auto pollard_rho(int64_t n) -> std::vector<std::pair<int64_t, int64_t>> {
std::vector<std::pair<int64_t, int64_t>> ret;
for (int i = 2; i <= 1000000; ++i) {
if (n % i == 0) {
int c = 0;
while (n % i == 0) {
n /= i;
++c;
}
ret.emplace_back(i, c);
}
if (i > n) break;
}
while (n > 1) {
if (miller_rabin(n)) {
ret.emplace_back(n, 1);
break;
}
auto res = pollard_rho_impl::rho(n);
if (not res) {
assert(false);
}
auto r = *res;
if (r == 1) break;
int c = 0;
while (n % r == 0) {
n /= r;
++c;
}
ret.emplace_back(r, c);
}
std::sort(ret.begin(), ret.end());
return ret;
}
} // namespace haar_lib
#line 2 "Mylib/Number/Prime/pollard_rho.cpp"
#include <algorithm>
#include <cassert>
#include <numeric>
#include <optional>
#include <utility>
#include <vector>
#line 2 "Mylib/Misc/int128.cpp"
namespace haar_lib {
#ifdef __SIZEOF_INT128__
using uint128_t = __uint128_t;
using int128_t = __int128_t;
#else
#include <boost/multiprecision/cpp_int.hpp>
using uint128_t = boost::multiprecision::uint128_t;
using int128_t = boost::multiprecision::int128_t;
#endif
} // namespace haar_lib
#line 2 "Mylib/Number/Prime/miller_rabin.cpp"
#include <cstdint>
#include <initializer_list>
#line 5 "Mylib/Number/Prime/miller_rabin.cpp"
namespace haar_lib {
namespace miller_rabin_impl {
uint128_t power(uint128_t a, uint128_t b, uint128_t p) {
uint128_t ret = 1;
while (b > 0) {
if (b & 1) ret = ret * a % p;
a = a * a % p;
b >>= 1;
}
return ret;
}
bool is_composite(uint64_t a, uint64_t p, int s, uint64_t d) {
uint128_t x = power(a, d, p);
if (x == 1) return false;
for (int i = 0; i < s; ++i) {
if (x == p - 1) return false;
x = x * x % p;
}
return true;
}
} // namespace miller_rabin_impl
bool miller_rabin(uint64_t n) {
if (n <= 1) return false;
if (n == 2) return true;
if (n % 2 == 0) return false;
int s = 0;
uint64_t d = n - 1;
while ((d & 1) == 0) {
s += 1;
d >>= 1;
}
if (n < 4759123141) {
for (uint64_t x : {2, 7, 61}) {
if (x < n and miller_rabin_impl::is_composite(x, n, s, d)) return false;
}
return true;
}
for (uint64_t x : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
if (x < n and miller_rabin_impl::is_composite(x, n, s, d)) return false;
}
return true;
}
} // namespace haar_lib
#line 10 "Mylib/Number/Prime/pollard_rho.cpp"
namespace haar_lib {
namespace pollard_rho_impl {
int128_t f(int128_t x) {
return x * x + 1;
}
std::optional<int64_t> rho(int64_t n) {
int64_t x = 2, y = 2, d = 1;
while (d == 1) {
x = f(x) % n;
y = f(f(y) % n) % n;
d = std::gcd(std::abs(x - y), n);
if (d == n) return {};
}
return {d};
}
} // namespace pollard_rho_impl
auto pollard_rho(int64_t n) -> std::vector<std::pair<int64_t, int64_t>> {
std::vector<std::pair<int64_t, int64_t>> ret;
for (int i = 2; i <= 1000000; ++i) {
if (n % i == 0) {
int c = 0;
while (n % i == 0) {
n /= i;
++c;
}
ret.emplace_back(i, c);
}
if (i > n) break;
}
while (n > 1) {
if (miller_rabin(n)) {
ret.emplace_back(n, 1);
break;
}
auto res = pollard_rho_impl::rho(n);
if (not res) {
assert(false);
}
auto r = *res;
if (r == 1) break;
int c = 0;
while (n % r == 0) {
n /= r;
++c;
}
ret.emplace_back(r, c);
}
std::sort(ret.begin(), ret.end());
return ret;
}
} // namespace haar_lib
Back to top page