高速素因数分解 (Math/fast_factorizer.hpp)
- View this file on GitHub
- View document part on GitHub
- Last update: 2025-04-13 19:22:19+09:00
- Include:
#include "Math/fast_factorizer.hpp"
FastFactorizer
Pollard’s rho アルゴリズムを用いた素因数分解を行います。 素数判定にはMiller-Rabin素数判定法を使用しています。
FastFactorizer::factorize
static vector<long long> FastFactorizer::factorize(long long N)
整数 N
を素因数分解し、因数を昇順に並べたリストを返します。
制約
- $0 \lt N$
計算量
- $O(\sqrt[4]{N})$
Depends on
Verified with
Code
#include <bits/stdc++.h>
using namespace std;
#include "Math/primality_test.hpp"
struct FastFactorizer
{
private:
static unsigned long long abs_diff(unsigned long long a, unsigned long long b)
{
return (a > b) ? (a - b) : (b - a);
}
// pollards_rho法によりNの約数を返す
static unsigned long long pollards_rho(unsigned long long N)
{
if (PrimalityTest::is_prime(N))
{
return N;
}
unsigned long long m = 1ll << (bit_width(N) / 8);
// 100回試して見つからなかったら失敗
for (unsigned long long c = 1; c < 100; c++)
{
auto f = [&](unsigned long long x) -> unsigned long long
{ return ((__uint128_t(x) * x) + c) % N; };
unsigned long long x0 = c % (N + 1);
unsigned long long y = x0, q = 1, g = 1, ys = 0, x = 0;
for (int r = 1; g <= 1; r <<= 1)
{
x = y;
for (int i = 0; i < r; i++)
{
y = f(y);
}
unsigned long long k = 0;
while (k < r && g <= 1)
{
ys = y;
for (int i = 0; i < min(m, r - k); i++)
{
y = f(y);
q = __uint128_t(q) * abs_diff(x, y) % N;
}
g = gcd(q, N);
k += m;
}
}
if (g == N)
{
while (g <= 1)
{
ys = f(ys);
g = gcd(abs_diff(x, ys), N);
}
}
if (g != N)
{
return g;
}
}
assert(false && "fail to factorize by pollards rho algorithm");
}
public:
static vector<long long> factorize(long long N)
{
assert(0 < N);
if (N == 1)
{
return {};
}
long long p = pollards_rho(N);
if (p == N)
{
return {p};
}
vector<long long> left = factorize(p);
vector<long long> right = factorize(N / p);
vector<long long> result;
merge(left.begin(), left.end(), right.begin(), right.end(), back_inserter(result));
return result;
}
};
#line 1 "Math/fast_factorizer.hpp"
#include <bits/stdc++.h>
using namespace std;
#line 2 "Math/primality_test.hpp"
using namespace std;
#line 2 "Math/mod_power.hpp"
using namespace std;
#line 2 "Other/fast_power.hpp"
using namespace std;
template <class S>
S fast_pow(S x, long long n, function<S(S, S)> mul, function<S()> e)
{
assert(0 <= n);
S ans = e();
while (n)
{
if (n & 1)
{
ans = mul(ans, x);
}
x = mul(x, x);
n >>= 1;
}
return ans;
}
#line 5 "Math/mod_power.hpp"
template <typename T>
enable_if_t<is_integral_v<T> || is_same_v<T, __int128_t>, T>
mod_pow(T x, T n, T mod)
{
assert(0 <= n);
assert(0 < mod);
assert(mod <= numeric_limits<T>::max() / mod);
x %= mod;
if (x < 0)
{
x += mod;
}
auto mul = [&](T a, T b) -> T
{
return (a * b) % mod;
};
auto e = [&]() -> T
{
return 1;
};
return fast_pow<T>(x, n, mul, e);
}
#line 5 "Math/primality_test.hpp"
struct PrimalityTest
{
private:
static bool miller_rabin(long long N, const vector<long long> &A)
{
const int s = countr_zero((unsigned long long)N - 1);
const long long d = (N - 1) >> s;
auto is_composite = [&](long long a) -> bool
{
// Fermat の小定理が成り立たない
if (a % N == 0)
{
return false;
}
long long x = mod_pow<__int128_t>(a, d, N);
if (x == 1 || x == N - 1)
{
return false;
}
for (int t = 0; t < s - 1; ++t)
{
x = __int128_t(x) * x % N;
if (x == N - 1)
{
return false;
}
}
return true;
};
return none_of(A.begin(), A.end(), is_composite);
}
public:
static bool is_prime(long long N)
{
assert(0 <= N);
if (N <= 2)
{
return N == 2;
}
if (N % 2 == 0)
{
return false;
}
return miller_rabin(N, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}
};
#line 5 "Math/fast_factorizer.hpp"
struct FastFactorizer
{
private:
static unsigned long long abs_diff(unsigned long long a, unsigned long long b)
{
return (a > b) ? (a - b) : (b - a);
}
// pollards_rho法によりNの約数を返す
static unsigned long long pollards_rho(unsigned long long N)
{
if (PrimalityTest::is_prime(N))
{
return N;
}
unsigned long long m = 1ll << (bit_width(N) / 8);
// 100回試して見つからなかったら失敗
for (unsigned long long c = 1; c < 100; c++)
{
auto f = [&](unsigned long long x) -> unsigned long long
{ return ((__uint128_t(x) * x) + c) % N; };
unsigned long long x0 = c % (N + 1);
unsigned long long y = x0, q = 1, g = 1, ys = 0, x = 0;
for (int r = 1; g <= 1; r <<= 1)
{
x = y;
for (int i = 0; i < r; i++)
{
y = f(y);
}
unsigned long long k = 0;
while (k < r && g <= 1)
{
ys = y;
for (int i = 0; i < min(m, r - k); i++)
{
y = f(y);
q = __uint128_t(q) * abs_diff(x, y) % N;
}
g = gcd(q, N);
k += m;
}
}
if (g == N)
{
while (g <= 1)
{
ys = f(ys);
g = gcd(abs_diff(x, ys), N);
}
}
if (g != N)
{
return g;
}
}
assert(false && "fail to factorize by pollards rho algorithm");
}
public:
static vector<long long> factorize(long long N)
{
assert(0 < N);
if (N == 1)
{
return {};
}
long long p = pollards_rho(N);
if (p == N)
{
return {p};
}
vector<long long> left = factorize(p);
vector<long long> right = factorize(N / p);
vector<long long> result;
merge(left.begin(), left.end(), right.begin(), right.end(), back_inserter(result));
return result;
}
};