CP-library
This documentation is automatically generated by online-judge-tools/verification-helper
NT/prime/factorize.hpp
- View this file on GitHub
- Last update: 2024-06-05 18:50:12+07:00
- Include:
#include "NT/prime/factorize.hpp"
Depends on
Required by
Verified with
Code
#pragma once
#include "Modint/Barrett_2.hpp"
#include "NT/prime/prime_test.hpp"
u64 pollard(u64 n) {
Barrett_u64 br;
br.set(n);
u64 x = 0, y = 0, t = 30, prd = 2, i = 1, q;
auto f = [&](u64 x) { return mul_b64(&br,x, x) + i; };
while (t++ % 40 || __gcd(prd, n) == 1) {
if (x == y) x = ++i, y = f(x);
if ((q = mul_b64(&br,prd, max(x,y) - min(x,y)))) prd = q;
x = f(x), y = f(f(y));
}
return __gcd(prd, n);
}
vector<u64> factor(u64 n) {
if (n == 1) return {};
if (isPrime(n)) return {n};
u64 x = pollard(n);
auto l = factor(x), r = factor(n / x);
l.insert(l.end(), all(r));
return l;
}#line 2 "Mod/mod_inv.hpp"
#include <cassert>
#include <type_traits>
// gcd(a, m) != 1 return -1
template <typename T>
T inv_mod(T a, T m) {
if (m == 1) return 0;
if (a >= m) a %= m;
if (a < 0) a += m;
if(__gcd(a,m)!=1) return -1;
T b = m, s = 1, t = 0;
while (true) {
if (a == 1) return s;
t -= b / a * s;
b %= a;
if (b == 1) return t + m;
s -= a / b * t;
a %= b;
}
}
#line 3 "Modint/Barrett_2.hpp"
struct Barrett_u32 {
u64 m, m2, im;
u64 div, rem;
void set(u64 mod) {
m = mod;
m2 = mod << 1;
im = ((((u128)(1ull) << 64)) + mod - 1) / mod;
div = 0;
rem = 0;
}
void reduce(u64 a) {
u64 x = (u64)(((u128)(a) * im) >> 64);
u64 y = x * m;
unsigned long long z;
u32 w = __builtin_usubll_overflow(a, y, &z) ? m : 0;
div = x;
rem = z + w;
}
u32 add(u32 a, u32 b) {
a += b;
a -= (a >= (u32)m ? (u32)m : 0);
return a;
}
u32 sub(u32 a, u32 b) {
a += (a < b ? (u32)m : 0);
a -= b;
return a;
}
u32 min(u32 a) {
return sub(0, a);
}
u32 shl(u32 a) {
return (a <<= 1) >= m ? a - m : a;
}
u32 shr(u32 a) {
return (a & 1) ? ((a >> 1) + (m >> 1) + 1) : (a >> 1);
}
};
struct Barrett_u64 {
u128 m, m2, im;
u128 div, rem;
void set(u128 mod) {
m = mod;
m2 = mod << 1;
im = (~((u128)0ull)) / mod;
div = 0;
rem = 0;
}
void reduce(u128 x) {
if (m == 1) {
div = x;
rem = 0;
return;
}
uint8_t f;
u128 a = x >> 64;
u128 b = x & 0xffffffffffffffffull;
u128 c = im >> 64;
u128 d = im & 0xffffffffffffffffull;
u128 ac = a * c;
u128 bd = (b * d) >> 64;
u128 ad = a * d;
u128 bc = b * c;
f = (ad > ((u128)((i128)(-1L)) - bd));
bd += ad;
ac += f;
f = (bc > ((u128)((i128)(-1L)) - bd));
bd += bc;
ac += f;
u128 q = ac + (bd >> 64);
u128 r = x - q * m;
if (m <= r) {
r -= m;
q += 1;
}
div = q;
rem = r;
}
u64 add(u64 a, u64 b) {
a += b;
a -= (a >= (u64)m ? (u64)m : 0);
return a;
}
u64 sub(u64 a, u64 b) {
a += (a < b ? (u64)m : 0);
a -= b;
return a;
}
u64 min(u64 a) {
return sub(0, a);
}
u64 shl(u64 a) {
return (a <<= 1) >= m ? a - m : a;
}
u64 shr(u64 a) {
return (a & 1) ? ((a >> 1) + (m >> 1) + 1) : (a >> 1);
}
};
// Barrett reduction - 32-bit
u32 mul_b32(struct Barrett_u32 *b, u32 a, u32 x) {
b->reduce((u64)a * x);
return (u32)b->rem;
}
u32 div_b32(struct Barrett_u32 *b, u32 a, u32 x) {
b->reduce((u64)a << 32);
return (u32)(b->div * inv_mod(x,(u32)b->m));
}
u32 pow_b32(struct Barrett_u32 *b, u32 a, u64 k) {
u32 ret = 1u;
u64 deg = k;
while (deg > 0) {
if (deg & 1) {
ret = mul_b32(b, ret, a);
}
a = mul_b32(b, a, a);
deg >>= 1;
}
return ret;
}
// Barrett reduction - 64-bit
u64 mul_b64(struct Barrett_u64 *b, u64 a, u64 x) {
b->reduce((u128)a * x);
return (u64)b->rem;
}
u64 div_b64(struct Barrett_u64 *b, u64 a, u64 x) {
b->reduce((u128)a << 64);
return (u64)(b->div * inv_mod(x,(u64)b->m));
}
u64 pow_b64(struct Barrett_u64 *b, u64 a, u64 k) {
u64 ret = 1ull, deg = k;
while (deg > 0) {
if (deg & 1) {
ret = mul_b64(b, ret, a);
}
a = mul_b64(b, a, a);
deg >>= 1;
}
return ret;
}
#line 2 "Mod/mod_mul.hpp"
u64 get_nr(u64 M) {
u64 IV = 2 - M;
for (int i = 0; i < 5; ++i) {
IV *= 2 - M * IV;
}
return IV;
}
u64 modmul(u64 x, u64 y, u64 M) {
u64 IV=get_nr(M);
auto t = u128(x) * y;
u64 lo = t, hi = t >> 64;
return (hi + M) - u64((u128(lo * IV) * M) >> 64);
}
u64 modpow(u64 a, u64 b, u64 M) {
u64 res = 1;
u64 IV = get_nr(M);
while (b) {
if (b & 1) {
res = modmul(res, a, M);
}
a = modmul(a, a, M);
b >>= 1;
}
return res;
}
//or
//only good for long long or int64_t
long long modmul2(long long a,long long b,long long mod){
return (a*b)%mod;
}
#line 3 "NT/prime/prime_test.hpp"
bool isPrime(u64 x) {
if (x < 64) {
return (u64(1) << x) & 0x28208a20a08a28ac;
}
if (x % 2 == 0) {
return false;
}
const int k = __builtin_ctzll(x - 1);
const u64 d = (x - 1) >> k, IV = get_nr(x), R = (-x) % x, R2 = (-u128(x)) % x, nR = x - R;
auto mr7 = [&](u64 t1, u64 t2, u64 t3, u64 t4, u64 t5, u64 t6, u64 t7) {
u64 r1 = R, r2 = R, r3 = R, r4 = R, r5 = R, r6 = R, r7 = R;
t1 = modmul(t1, R2, x), t2 = modmul(t2, R2, x), t3 = modmul(t3, R2, x);
t4 = modmul(t4, R2, x), t5 = modmul(t5, R2, x), t6 = modmul(t6, R2, x), t7 = modmul(t7, R2, x);
for (u64 b = d; b; b >>= 1) {
if (b & 1) {
r1 = modmul(r1, t1, x), r2 = modmul(r2, t2, x), r3 = modmul(r3, t3, x);
r4 = modmul(r4, t4, x), r5 = modmul(r5, t5, x), r6 = modmul(r6, t6, x), r7 = modmul(r7, t7, x);
}
t1 = modmul(t1, t1, x), t2 = modmul(t2, t2, x), t3 = modmul(t3, t3, x);
t4 = modmul(t4, t4, x), t5 = modmul(t5, t5, x), t6 = modmul(t6, t6, x), t7 = modmul(t7, t7, x);
}
r1 = min(r1, r1 - x), r2 = min(r2, r2 - x), r3 = min(r3, r3 - x);
r4 = min(r4, r4 - x), r5 = min(r5, r5 - x), r6 = min(r6, r6 - x), r7 = min(r7, r7 - x);
int res1 = (r1 == R) | (r1 == nR), res2 = (r2 == R) | (r2 == nR), res3 = (r3 == R) | (r3 == nR);
int res4 = (r4 == R) | (r4 == nR), res5 = (r5 == R) | (r5 == nR), res6 = (r6 == R) | (r6 == nR), res7 = (r7 == R) | (r7 == nR);
for (int j = 0; j < k - 1; ++j) {
r1 = modmul(r1, r1, x), r2 = modmul(r2, r2, x), r3 = modmul(r3, r3, x);
r4 = modmul(r4, r4, x), r5 = modmul(r5, r5, x), r6 = modmul(r6, r6, x), r7 = modmul(r7, r7, x);
res1 |= (min(r1, r1 - x) == nR), res2 |= (min(r2, r2 - x) == nR), res3 |= (min(r3, r3 - x) == nR);
res4 |= (min(r4, r4 - x) == nR), res5 |= (min(r5, r5 - x) == nR), res6 |= (min(r6, r6 - x) == nR), res7 |= (min(r7, r7 - x) == nR);
}
return res1 & res2 & res3 & res4 & res5 & res6 & res7;
};
if (x == 2 || x == 3 || x == 5 || x == 13 || x == 19 || x == 73 || x == 193 || x == 407521 || x == 299210837) {
return true;
}
return mr7(2, 325, 9375, 28178, 450775, 9780504, 1795265022);
}
#line 4 "NT/prime/factorize.hpp"
u64 pollard(u64 n) {
Barrett_u64 br;
br.set(n);
u64 x = 0, y = 0, t = 30, prd = 2, i = 1, q;
auto f = [&](u64 x) { return mul_b64(&br,x, x) + i; };
while (t++ % 40 || __gcd(prd, n) == 1) {
if (x == y) x = ++i, y = f(x);
if ((q = mul_b64(&br,prd, max(x,y) - min(x,y)))) prd = q;
x = f(x), y = f(f(y));
}
return __gcd(prd, n);
}
vector<u64> factor(u64 n) {
if (n == 1) return {};
if (isPrime(n)) return {n};
u64 x = pollard(n);
auto l = factor(x), r = factor(n / x);
l.insert(l.end(), all(r));
return l;
}