This documentation is automatically generated by competitive-verifier/competitive-verifier
#include "number/pollard_rho.hpp"
#include "miller_rabin.hpp"
#include "../hash/splitmix64.hpp"
// n must not be prime
long long find_prime_factor(long long n) {
constexpr int m = 200;
for (long long c_i = 1;; ++c_i) {
long long c = splitmix64(c_i) % n;
if (c == 0) continue;
auto f = [=](long long x) { return (__int128_t(x) * x + c) % n; };
long long x = 0, y = 0, y_b = 0, q = 1;
long long g = 1;
int r = 1, k = 0;
do {
x = y;
while (k < 3 * r / 4) {
y = f(y);
++k;
}
while (k < r && g == 1) {
y_b = y;
for (int i = 0; i < min(m, r - k); ++i) {
y = f(y);
q = __int128_t(q) * abs(x - y) % n;
}
g = gcd(q, n);
k += m;
}
k = r;
r <<= 1;
} while (g == 1);
if (g == n) {
g = 1;
y = y_b;
do {
y = f(y);
g = gcd(abs(x - y), n);
} while (g == 1);
}
if (g == n) continue;
if (is_prime(g)) return g;
if (is_prime(n / g)) return n / g;
return find_prime_factor(min(n / g, g));
}
return n;
}
vector<long long> factorize(long long n) {
vector<long long> ps;
if (n >= 2) {
while ((n & 1) == 0) {
ps.push_back(2);
n >>= 1;
}
if (is_prime(n)) {
ps.push_back(n);
} else {
while (n > 1) {
long long p = find_prime_factor(n);
do {
ps.push_back(p);
n /= p;
} while (n % p == 0);
if (is_prime(n)) {
ps.push_back(n);
break;
}
}
}
}
sort(ps.begin(), ps.end());
return ps;
}
#line 1 "number/miller_rabin.hpp"
bool miller_rabin(long long n, const vector<long long> &ts) {
const long long s2 = (n - 1) & (1 - n), d = (n - 1) / s2;
const auto pow_mod_n = [=](long long b, long long p) {
long long res = 1;
while (p > 0) {
if (p & 1) res = __int128_t(res) * b % n;
b = __int128_t(b) * b % n;
p >>= 1;
}
return res;
};
for (const long long a : ts) {
if (a >= n) break;
auto ad = pow_mod_n(a, d);
if (ad == 1) continue;
long long t = 1;
for (; t < s2; t <<= 1) {
if (ad == n - 1) break;
ad = __int128_t(ad) * ad % n;
}
if (t == s2) return false;
}
return true;
}
bool is_prime(long long n) {
if (n == 1) return false;
if (n == 2) return true;
if ((n & 1) == 0) return false;
if (n < 4759123141) return miller_rabin(n, {2, 7, 61});
return miller_rabin(n, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}
#line 1 "hash/splitmix64.hpp"
constexpr uint64_t splitmix64(uint64_t x) {
// http://xorshift.di.unimi.it/splitmix64.c
x += 0x9e3779b97f4a7c15;
x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
return x ^ (x >> 31);
}
#line 3 "number/pollard_rho.hpp"
// n must not be prime
long long find_prime_factor(long long n) {
constexpr int m = 200;
for (long long c_i = 1;; ++c_i) {
long long c = splitmix64(c_i) % n;
if (c == 0) continue;
auto f = [=](long long x) { return (__int128_t(x) * x + c) % n; };
long long x = 0, y = 0, y_b = 0, q = 1;
long long g = 1;
int r = 1, k = 0;
do {
x = y;
while (k < 3 * r / 4) {
y = f(y);
++k;
}
while (k < r && g == 1) {
y_b = y;
for (int i = 0; i < min(m, r - k); ++i) {
y = f(y);
q = __int128_t(q) * abs(x - y) % n;
}
g = gcd(q, n);
k += m;
}
k = r;
r <<= 1;
} while (g == 1);
if (g == n) {
g = 1;
y = y_b;
do {
y = f(y);
g = gcd(abs(x - y), n);
} while (g == 1);
}
if (g == n) continue;
if (is_prime(g)) return g;
if (is_prime(n / g)) return n / g;
return find_prime_factor(min(n / g, g));
}
return n;
}
vector<long long> factorize(long long n) {
vector<long long> ps;
if (n >= 2) {
while ((n & 1) == 0) {
ps.push_back(2);
n >>= 1;
}
if (is_prime(n)) {
ps.push_back(n);
} else {
while (n > 1) {
long long p = find_prime_factor(n);
do {
ps.push_back(p);
n /= p;
} while (n % p == 0);
if (is_prime(n)) {
ps.push_back(n);
break;
}
}
}
}
sort(ps.begin(), ps.end());
return ps;
}