library
This documentation is automatically generated by online-judge-tools/verification-helper
Project maintained by tko919 Hosted on GitHub Pages — Theme by mattgraham
Mod Fast
(Math/modfast.hpp)
- View this file on GitHub
- Last update: 2026-11-06 05:58:50+09:00
- Include:
#include "Math/modfast.hpp"
Depends on
Fast Division
(Math/fastdiv.hpp)
- lpf table
(Math/lpftable.hpp)
- Miller-Rabin
(Math/miller.hpp)
- Pollard-Rho
(Math/pollard.hpp)
- Primitive Function
(Math/primitive.hpp)
- Prime Sieve
(Math/sieve.hpp)
- Random
(Utility/random.hpp)
Code
#pragma once
#include "Math/primitive.hpp"
#include "Math/lpftable.hpp"
struct ModFast {
const int p;
array<pair<int, int>, (1 << 20) + 1> frac;
static constexpr int K = 1 << 21;
array<int, K * 2 + 1> LOG, INV;
array<ll, (1 << 15) + 1> phi, plo;
ModFast(int p_, int rt = -1) : p(p_) {
if (rt == -1) {
if (p == 998244353)
rt = 3;
else if (p == 1000000007)
rt = 5;
else
rt = getPrimitiveRoot(p);
}
// build POW,INV
{
plo[0] = phi[0] = 1;
rep(i, 0, 1 << 15) plo[i + 1] = (plo[i] * rt) % p;
rep(i, 0, 1 << 15) phi[i + 1] = (phi[i] * plo[1 << 15]) % p;
INV[K + 1] = 1;
rep(i, 2, K + 1) {
ll q = (p + i - 1) / i;
INV[K + i] = ll(INV[K + i * q - p]) * q % p;
}
rep(i, 1, K + 1) INV[K - i] = p - INV[K + i];
}
// build LOG
{
const int S = 1 << 17;
unordered_map<int, int> mp;
ll pw = 1;
rep(k, 0, S) {
mp[pw] = k;
pw = (pw * rt) % p;
}
auto lpf = lpftable(K);
rep(i, 2, K + 1) {
if (lpf[i] < i) {
int q = lpf[i];
LOG[K + i] = (LOG[K + q] + LOG[K + i / q]) % (p - 1);
continue;
}
if (ll(i) * i > p) {
int q = p / i, r = p - i * q;
LOG[K + i] =
(ll(LOG[K + r] + (p - 1) / 2) + (p - 1 - LOG[K + q])) %
(p - 1);
continue;
}
{
ll tar = i, add = 0;
for (;;) {
if (mp.count(tar)) {
LOG[K + i] = (add + mp[tar]) % (p - 1);
break;
}
add += p - 1 - S, tar = (tar * pw) % p;
}
continue;
}
for (;;) {
int k = Random::get(0, p - 1);
ll ret = p - 1 - k, x = (ll(i) * pow_r(k)) % p;
for (int q : {2, 3, 5, 7, 11, 13, 17, 19}) {
while (x % q == 0) {
x /= q;
ret += LOG[K + q];
}
}
if (x >= K)
continue;
while (i < x && x < K && lpf[x] < i) {
int q = lpf[x];
ret += LOG[K + q];
x /= q;
}
if (1 < x && x < i)
ret += LOG[K + x];
if (x == 1) {
LOG[K + i] = ret % (p - 1);
break;
}
}
}
rep(i, 1, K + 1) LOG[K - i] = (LOG[K + i] + (p - 1) / 2) % (p - 1);
}
// build frac
{
stack<tuple<int, int, int, int>> st;
st.push({0, 1, 1, 1});
while (!st.empty()) {
auto [a, b, c, d] = st.top();
st.pop();
if (b + d < 2048) {
st.emplace(a, b, a + c, b + d);
st.emplace(a + c, b + d, c, d);
continue;
}
int L = (ll(a) * p) / b / 1024, R = (ll(c) * p) / d / 1024;
frac[L] = {a, b}, frac[R] = {c, d};
rep(i, L + 1, R) frac[i] = {min(a, c), min(b, d)};
}
}
}
pair<int, int> get(int x) {
auto [a, b] = frac[x >> 10];
int t = x * b - a * p;
return {t % p, b % p};
}
int inv(int x) {
auto [t, b] = get(x);
return ll(INV[K + t]) * b % p;
}
int pow(int x, ll t) {
return pow_r(t * log(x) % (p - 1));
}
int pow_r(ll t) {
return phi[t >> 15] * plo[t & 32767] % p;
}
int log(int x) {
auto [t, b] = get(x);
return ((LOG[K + t]) + (p - 1) - LOG[K + b]) % (p - 1);
}
};
/**
* @brief Mod Fast
*/#line 2 "Math/fastdiv.hpp"
struct FastDiv {
using u64 = uint64_t;
using u128 = __uint128_t;
constexpr FastDiv() : m(), s(), x() {}
constexpr FastDiv(int _m)
: m(_m), s(__lg(m - 1)), x(((u128(1) << (s + 64)) + m - 1) / m) {}
constexpr int get() {
return m;
}
constexpr friend u64 operator/(u64 n, const FastDiv &d) {
return (u128(n) * d.x >> d.s) >> 64;
}
constexpr friend int operator%(u64 n, const FastDiv &d) {
return n - n / d * d.m;
}
constexpr pair<u64, int> divmod(u64 n) const {
u64 q = n / (*this);
return {q, n - q * m};
}
int m, s;
u64 x;
};
struct FastDiv64 {
using u64 = uint64_t;
using u128 = __uint128_t;
u128 mod, mh, ml;
explicit FastDiv64(u64 mod = 1) : mod(mod) {
u128 m = u128(-1) / mod;
if (m * mod + mod == u128(0))
++m;
mh = m >> 64;
ml = m & u64(-1);
}
u64 umod() const {
return mod;
}
u64 modulo(u128 x) {
u128 z = (x & u64(-1)) * ml;
z = (x & u64(-1)) * mh + (x >> 64) * ml + (z >> 64);
z = (x >> 64) * mh + (z >> 64);
x -= z * mod;
return x < mod ? x : x - mod;
}
u64 mul(u64 a, u64 b) {
return modulo(u128(a) * b);
}
};
/**
* @brief Fast Division
*/
#line 2 "Math/miller.hpp"
struct m64 {
using i64 = int64_t;
using u64 = uint64_t;
using u128 = __uint128_t;
static u64 mod;
static u64 r;
static u64 n2;
static u64 get_r() {
u64 ret = mod;
rep(_,0,5) ret *= 2 - mod * ret;
return ret;
}
static void set_mod(u64 m) {
assert(m < (1LL << 62));
assert((m & 1) == 1);
mod = m;
n2 = -u128(m) % m;
r = get_r();
assert(r * mod == 1);
}
static u64 get_mod() { return mod; }
u64 a;
m64() : a(0) {}
m64(const int64_t &b) : a(reduce((u128(b) + mod) * n2)){};
static u64 reduce(const u128 &b) {
return (b + u128(u64(b) * u64(-r)) * mod) >> 64;
}
u64 get() const {
u64 ret = reduce(a);
return ret >= mod ? ret - mod : ret;
}
m64 &operator*=(const m64 &b) {
a = reduce(u128(a) * b.a);
return *this;
}
m64 operator*(const m64 &b) const { return m64(*this) *= b; }
bool operator==(const m64 &b) const {
return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
}
bool operator!=(const m64 &b) const {
return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
}
m64 pow(u128 n) const {
m64 ret(1), mul(*this);
while (n > 0) {
if (n & 1) ret *= mul;
mul *= mul;
n >>= 1;
}
return ret;
}
};
typename m64::u64 m64::mod, m64::r, m64::n2;
bool Miller(ll n){
if(n<2 or (n&1)==0)return (n==2);
m64::set_mod(n);
ll d=n-1; while((d&1)==0)d>>=1;
vector<ll> seeds;
if(n<(1<<30))seeds={2, 7, 61};
else seeds={2, 325, 9375, 28178, 450775, 9780504};
for(auto& x:seeds){
if(n<=x)break;
ll t=d;
m64 y=m64(x).pow(t);
while(t!=n-1 and y!=1 and y!=n-1){
y*=y;
t<<=1;
}
if(y!=n-1 and (t&1)==0)return 0;
} return 1;
}
/**
* @brief Miller-Rabin
*/
#line 2 "Utility/random.hpp"
namespace Random {
mt19937_64 randgen(chrono::steady_clock::now().time_since_epoch().count());
using u64 = unsigned long long;
u64 get() {
return randgen();
}
template <typename T> T get(T L) { // [0,L]
return get() % (L + 1);
}
template <typename T> T get(T L, T R) { // [L,R]
return get(R - L) + L;
}
double uniform() {
return double(get(1000000000)) / 1000000000;
}
string str(int n) {
string ret;
rep(i, 0, n) ret += get('a', 'z');
return ret;
}
template <typename Iter> void shuffle(Iter first, Iter last) {
if (first == last)
return;
int len = 1;
for (auto it = first + 1; it != last; it++) {
len++;
int j = get(0, len - 1);
if (j != len - 1)
iter_swap(it, first + j);
}
}
template <typename T> vector<T> select(int n, T L, T R) { // [L,R]
if (n * 2 >= R - L + 1) {
vector<T> ret(R - L + 1);
iota(ALL(ret), L);
shuffle(ALL(ret));
ret.resize(n);
return ret;
} else {
unordered_set<T> used;
vector<T> ret;
while (SZ(used) < n) {
T x = get(L, R);
if (!used.count(x)) {
used.insert(x);
ret.push_back(x);
}
}
return ret;
}
}
void relabel(int n, vector<pair<int, int>> &es) {
shuffle(ALL(es));
vector<int> ord(n);
iota(ALL(ord), 0);
shuffle(ALL(ord));
for (auto &[u, v] : es)
u = ord[u], v = ord[v];
}
template <bool directed, bool multi, bool self>
vector<pair<int, int>> genGraph(int n, int m) {
vector<pair<int, int>> cand, es;
rep(u, 0, n) rep(v, 0, n) {
if (!self and u == v)
continue;
if (!directed and u > v)
continue;
cand.push_back({u, v});
}
if (m == -1)
m = get(SZ(cand));
// chmin(m, SZ(cand));
vector<int> ord;
if (multi)
rep(_, 0, m) ord.push_back(get(SZ(cand) - 1));
else {
ord = select(m, 0, SZ(cand) - 1);
}
for (auto &i : ord)
es.push_back(cand[i]);
relabel(n, es);
return es;
}
vector<pair<int, int>> genTree(int n) {
vector<pair<int, int>> es;
rep(i, 1, n) es.push_back({get(i - 1), i});
relabel(n, es);
return es;
}
}; // namespace Random
/**
* @brief Random
*/
#line 4 "Math/pollard.hpp"
vector<ll> Pollard(ll n) {
if (n <= 1)
return {};
if (Miller(n))
return {n};
if ((n & 1) == 0) {
vector<ll> v = Pollard(n >> 1);
v.push_back(2);
return v;
}
for (ll x = 2, y = 2, d;;) {
ll c = Random::get(2LL, n - 1);
do {
x = (__int128_t(x) * x + c) % n;
y = (__int128_t(y) * y + c) % n;
y = (__int128_t(y) * y + c) % n;
d = __gcd(x - y + n, n);
} while (d == 1);
if (d < n) {
vector<ll> lb = Pollard(d), rb = Pollard(n / d);
lb.insert(lb.end(), ALL(rb));
return lb;
}
}
}
vector<pair<ll, int>> Pollard2(ll n) {
auto ps = Pollard(n);
sort(ALL(ps));
using P = pair<ll, int>;
vector<P> pes;
for (auto &p : ps) {
if (pes.empty() or pes.back().first != p) {
pes.push_back({p, 1});
} else {
pes.back().second++;
}
}
return pes;
}
vector<ll> EnumDivisors(ll n) {
auto pes = Pollard2(n);
vector<ll> ret;
auto rec = [&](auto &rec, int id, ll d) -> void {
if (id == SZ(pes)) {
ret.push_back(d);
return;
}
rec(rec, id + 1, d);
rep(e, 0, pes[id].second) {
d *= pes[id].first;
rec(rec, id + 1, d);
}
};
rec(rec, 0, 1);
sort(ALL(ret));
return ret;
}
/**
* @brief Pollard-Rho
*/
#line 4 "Math/primitive.hpp"
ll mpow(ll a, ll t, ll m) {
ll res = 1;
FastDiv64 im(m);
while (t) {
if (t & 1)
res = im.modulo(__int128_t(res) * a);
a = im.modulo(__int128_t(a) * a);
t >>= 1;
}
return res % m;
}
ll minv(ll a, ll m) {
ll b = m, u = 1, v = 0;
while (b) {
ll t = a / b;
a -= t * b;
swap(a, b);
u -= t * v;
swap(u, v);
}
u = (u % m + m) % m;
return u;
}
ll getPrimitiveRoot(ll p) {
vector<ll> ps = Pollard(p - 1);
sort(ALL(ps));
rep(x, 1, inf) {
for (auto &q : ps) {
if (mpow(x, (p - 1) / q, p) == 1)
goto fail;
}
return x;
fail:;
}
assert(0);
}
template <typename T> T extgcd(T a, T b, T &p, T &q) {
if (b == 0) {
p = 1;
q = 0;
return a;
}
T d = extgcd(b, a % b, q, p);
q -= a / b * p;
return d;
}
template <typename T> pair<T, T> crt(vector<T> vs, vector<T> ms) {
T V = vs[0], M = ms[0];
rep(i, 1, vs.size()) {
T p, q, v = vs[i], m = ms[i];
if (M < m)
swap(M, m), swap(V, v);
T d = extgcd(M, m, p, q);
if ((v - V) % d != 0)
return {0, -1};
T md = m / d, tmp = (v - V) / d % md * p % md;
V += M * tmp;
M *= md;
}
V = (V % M + M) % M;
return {V, M};
}
ll garner(vector<ll> vs, vector<ll> p, int mod) {
int sz = SZ(vs);
vector<ll> kp(sz + 1), rmul(sz + 1, 1);
p.push_back(mod);
rep(i, 0, sz) {
ll x = (vs[i] - kp[i]) * minv(rmul[i], p[i]) % p[i];
if (x < 0)
x += p[i];
rep(j, i + 1, sz + 1) {
kp[j] += rmul[j] * x;
kp[j] %= p[j];
rmul[j] *= p[i];
rmul[j] %= p[j];
}
}
return kp.back();
}
ll ModLog(ll a, ll b, ll p) {
ll g = 1;
for (ll t = p; t; t >>= 1)
g = g * a % p;
g = __gcd(g, p);
ll t = 1, c = 0;
for (; t % g; c++) {
if (t == b)
return c;
t = t * a % p;
}
if (b % g)
return -1;
t /= g, b /= g;
ll n = p / g, h = 0, gs = 1;
for (; h * h < n; h++)
gs = gs * a % n;
unordered_map<ll, ll> bs;
for (ll s = 0, e = b; s < h; bs[e] = ++s)
e = e * a % n;
for (ll s = 0, e = t; s < n;) {
e = e * gs % n, s += h;
if (bs.count(e)) {
return c + s - bs[e];
}
}
return -1;
}
ll TonelliShanks(ll a, ll p) {
a %= p;
if (a == 0)
return 0;
if (p == 2)
return a;
if (mpow(a, (p - 1) >> 1, p) != 1)
return -1;
ll b = 1;
while (mpow(b, (p - 1) >> 1, p) == 1)
b = Random::get(1LL, p - 1);
ll q = p - 1, k = 0;
while (q % 2 == 0) {
q >>= 1;
k++;
}
ll x = mpow(a, (q + 1) >> 1, p);
b = mpow(b, q, p);
k -= 2;
while (mpow(x, 2, p) != a) {
ll err = minv(a, p) * mpow(x, 2, p) % p;
if (mpow(err, 1 << k, p) != 1)
x = x * b % p;
b = mpow(b, 2, p);
k--;
}
return x;
}
ll mod_root(ll k, ll a, ll m) {
if (a == 0)
return k ? 0 : -1;
if (m == 2)
return a & 1;
k %= m - 1;
ll g = gcd(k, m - 1);
if (mpow(a, (m - 1) / g, m) != 1)
return -1;
a = mpow(a, minv(k / g, (m - 1) / g), m);
FastDiv64 im(m);
auto _subroot = [&](ll p, int e, ll a) -> ll { // x^(p^e)==a(mod m)
ll q = m - 1;
int s = 0;
while (q % p == 0) {
q /= p;
s++;
}
int d = s - e;
ll pe = mpow(p, e, m),
res = mpow(a, ((pe - 1) * minv(q, pe) % pe * q + 1) / pe, m), c = 1;
while (mpow(c, (m - 1) / p, m) == 1)
c++;
c = mpow(c, q, m);
map<ll, ll> mp;
ll v = 1, block = sqrt(d * p) + 1,
bs = mpow(c, mpow(p, s - 1, m - 1) * block % (m - 1), m);
rep(i, 0, block + 1) mp[v] = i, v = im.modulo(i128(v) * bs);
ll gs = minv(mpow(c, mpow(p, s - 1, m - 1), m), m);
rep(i, 0, d) {
ll err = im.modulo(i128(a) * minv(mpow(res, pe, m), m));
ll pos = mpow(err, mpow(p, d - 1 - i, m - 1), m);
rep(j, 0, block + 1) {
if (mp.count(pos)) {
res = im.modulo(i128(res) *
mpow(c,
(block * mp[pos] + j) *
mpow(p, i, m - 1) % (m - 1),
m));
break;
}
pos = im.modulo(i128(pos) * gs);
}
}
return res;
};
for (ll d = 2; d * d <= g; d++)
if (g % d == 0) {
int sz = 0;
while (g % d == 0) {
g /= d;
sz++;
}
a = _subroot(d, sz, a);
}
if (g > 1)
a = _subroot(g, 1, a);
return a;
}
ull floor_root(ull a, ull k) {
if (a <= 1 or k == 1)
return a;
if (k >= 64)
return 1;
if (k == 2)
return sqrtl(a);
constexpr ull LIM = -1;
if (a == LIM)
a--;
auto mul = [&](ull &x, const ull &y) {
if (x <= LIM / y)
x *= y;
else
x = LIM;
};
auto pw = [&](ull x, ull t) -> ull {
ull y = 1;
while (t) {
if (t & 1)
mul(y, x);
mul(x, x);
t >>= 1;
}
return y;
};
ull ret = (k == 3 ? cbrt(a) - 1 : pow(a, nextafter(1 / double(k), 0)));
while (pw(ret + 1, k) <= a)
ret++;
return ret;
}
/**
* @brief Primitive Function
*/
#line 2 "Math/sieve.hpp"
template <int L = 101010101> vector<int> sieve(int N) {
bitset<L> isp;
int n, sq = ceil(sqrt(N));
for (int z = 1; z <= 5; z += 4) {
for (int y = z; y <= sq; y += 6) {
for (int x = 1; x <= sq and (n = 4 * x * x + y * y) <= N; ++x) {
isp[n].flip();
}
for (int x = y + 1; x <= sq and (n = 3 * x * x - y * y) <= N;
x += 2) {
isp[n].flip();
}
}
}
for (int z = 2; z <= 4; z += 2) {
for (int y = z; y <= sq; y += 6) {
for (int x = 1; x <= sq and (n = 3 * x * x + y * y) <= N; x += 2) {
isp[n].flip();
}
for (int x = y + 1; x <= sq and (n = 3 * x * x - y * y) <= N;
x += 2) {
isp[n].flip();
}
}
}
for (int y = 3; y <= sq; y += 6) {
for (int z = 1; z <= 2; ++z) {
for (int x = z; x <= sq and (n = 4 * x * x + y * y) <= N; x += 3) {
isp[n].flip();
}
}
}
for (int n = 5; n <= sq; ++n)
if (isp[n]) {
for (int k = n * n; k <= N; k += n * n) {
isp[k] = false;
}
}
isp[2] = isp[3] = true;
vector<int> ret;
for (int i = 2; i <= N; i++)
if (isp[i]) {
ret.push_back(i);
}
return ret;
}
/**
* @brief Prime Sieve
*/
#line 3 "Math/lpftable.hpp"
vector<int> lpftable(ll n) {
auto ps = sieve(n);
vector<int> ret(n + 1, -1);
reverse(ALL(ps));
for (auto &p : ps) {
for (ll q = p; q <= n; q += p)
ret[q] = p;
}
return ret;
}
/**
* @brief lpf table
*/
#line 4 "Math/modfast.hpp"
struct ModFast {
const int p;
array<pair<int, int>, (1 << 20) + 1> frac;
static constexpr int K = 1 << 21;
array<int, K * 2 + 1> LOG, INV;
array<ll, (1 << 15) + 1> phi, plo;
ModFast(int p_, int rt = -1) : p(p_) {
if (rt == -1) {
if (p == 998244353)
rt = 3;
else if (p == 1000000007)
rt = 5;
else
rt = getPrimitiveRoot(p);
}
// build POW,INV
{
plo[0] = phi[0] = 1;
rep(i, 0, 1 << 15) plo[i + 1] = (plo[i] * rt) % p;
rep(i, 0, 1 << 15) phi[i + 1] = (phi[i] * plo[1 << 15]) % p;
INV[K + 1] = 1;
rep(i, 2, K + 1) {
ll q = (p + i - 1) / i;
INV[K + i] = ll(INV[K + i * q - p]) * q % p;
}
rep(i, 1, K + 1) INV[K - i] = p - INV[K + i];
}
// build LOG
{
const int S = 1 << 17;
unordered_map<int, int> mp;
ll pw = 1;
rep(k, 0, S) {
mp[pw] = k;
pw = (pw * rt) % p;
}
auto lpf = lpftable(K);
rep(i, 2, K + 1) {
if (lpf[i] < i) {
int q = lpf[i];
LOG[K + i] = (LOG[K + q] + LOG[K + i / q]) % (p - 1);
continue;
}
if (ll(i) * i > p) {
int q = p / i, r = p - i * q;
LOG[K + i] =
(ll(LOG[K + r] + (p - 1) / 2) + (p - 1 - LOG[K + q])) %
(p - 1);
continue;
}
{
ll tar = i, add = 0;
for (;;) {
if (mp.count(tar)) {
LOG[K + i] = (add + mp[tar]) % (p - 1);
break;
}
add += p - 1 - S, tar = (tar * pw) % p;
}
continue;
}
for (;;) {
int k = Random::get(0, p - 1);
ll ret = p - 1 - k, x = (ll(i) * pow_r(k)) % p;
for (int q : {2, 3, 5, 7, 11, 13, 17, 19}) {
while (x % q == 0) {
x /= q;
ret += LOG[K + q];
}
}
if (x >= K)
continue;
while (i < x && x < K && lpf[x] < i) {
int q = lpf[x];
ret += LOG[K + q];
x /= q;
}
if (1 < x && x < i)
ret += LOG[K + x];
if (x == 1) {
LOG[K + i] = ret % (p - 1);
break;
}
}
}
rep(i, 1, K + 1) LOG[K - i] = (LOG[K + i] + (p - 1) / 2) % (p - 1);
}
// build frac
{
stack<tuple<int, int, int, int>> st;
st.push({0, 1, 1, 1});
while (!st.empty()) {
auto [a, b, c, d] = st.top();
st.pop();
if (b + d < 2048) {
st.emplace(a, b, a + c, b + d);
st.emplace(a + c, b + d, c, d);
continue;
}
int L = (ll(a) * p) / b / 1024, R = (ll(c) * p) / d / 1024;
frac[L] = {a, b}, frac[R] = {c, d};
rep(i, L + 1, R) frac[i] = {min(a, c), min(b, d)};
}
}
}
pair<int, int> get(int x) {
auto [a, b] = frac[x >> 10];
int t = x * b - a * p;
return {t % p, b % p};
}
int inv(int x) {
auto [t, b] = get(x);
return ll(INV[K + t]) * b % p;
}
int pow(int x, ll t) {
return pow_r(t * log(x) % (p - 1));
}
int pow_r(ll t) {
return phi[t >> 15] * plo[t & 32767] % p;
}
int log(int x) {
auto [t, b] = get(x);
return ((LOG[K + t]) + (p - 1) - LOG[K + b]) % (p - 1);
}
};
/**
* @brief Mod Fast
*/