This documentation is automatically generated by online-judge-tools/verification-helper
#include "Math/stirlingquery.hpp"
#pragma once
#include "Math/fastdiv.hpp"
class StirlingNumberQuery {
const int p;
FastDiv ip;
vector<vector<int>> binom, F, S;
ll nCr(ll n, ll k) {
if (n < 0 or k < 0 or n < k)
return 0;
ll res = 1;
while (n) {
res = (res * binom[n % ip][k % ip]) % ip;
n /= p;
k /= p;
}
return res;
}
public:
StirlingNumberQuery(int _p) : p(_p), ip(p) {
binom.resize(p, vector<int>(p));
F.resize(p, vector<int>(p));
S.resize(p, vector<int>(p));
binom[0][0] = F[0][0] = S[0][0] = 1;
rep(n, 1, p) rep(k, 0, n + 1) {
if (k)
binom[n][k] = binom[n - 1][k - 1];
binom[n][k] = (binom[n][k] + binom[n - 1][k]) % ip;
if (k)
F[n][k] = F[n - 1][k - 1];
F[n][k] = (F[n][k] + ll(p - n + 1) * F[n - 1][k]) % ip;
if (k)
S[n][k] = S[n - 1][k - 1];
S[n][k] = (S[n][k] + ll(k) * S[n - 1][k]) % ip;
}
}
int FirstKind(ll n, ll k) {
if (n < 0 or k < 0 or k > n)
return 0;
ll i = n / ip, j = n % ip;
if (k < i)
return 0;
ll a = (k - i) / (p - 1), b = (k - i) % (p - 1);
if (b == 0 and j)
b += p - 1, a--;
if (a < 0 or a > i or b > j)
return 0;
int res = (nCr(i, a) * F[j][b]) % ip;
if ((i + a) & 1)
res = (p - res) % ip;
return res;
}
int SecondKind(ll n, ll k) {
if (n < 0 or k < 0 or k > n)
return 0;
if (n == 0)
return 1;
ll i = k / ip, j = k % ip;
if (n < i)
return 0;
ll a = (n - i) / (p - 1), b = (n - i) % (p - 1);
if (b == 0)
b += p - 1, a--;
if (a < 0 or b < j)
return 0;
if (b == p - 1 and j == 0)
return nCr(a, i - 1);
else
return (nCr(a, i) * S[b][j]) % ip;
}
};
/**
* @brief Stirling Number for query
*/
#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 3 "Math/stirlingquery.hpp"
class StirlingNumberQuery {
const int p;
FastDiv ip;
vector<vector<int>> binom, F, S;
ll nCr(ll n, ll k) {
if (n < 0 or k < 0 or n < k)
return 0;
ll res = 1;
while (n) {
res = (res * binom[n % ip][k % ip]) % ip;
n /= p;
k /= p;
}
return res;
}
public:
StirlingNumberQuery(int _p) : p(_p), ip(p) {
binom.resize(p, vector<int>(p));
F.resize(p, vector<int>(p));
S.resize(p, vector<int>(p));
binom[0][0] = F[0][0] = S[0][0] = 1;
rep(n, 1, p) rep(k, 0, n + 1) {
if (k)
binom[n][k] = binom[n - 1][k - 1];
binom[n][k] = (binom[n][k] + binom[n - 1][k]) % ip;
if (k)
F[n][k] = F[n - 1][k - 1];
F[n][k] = (F[n][k] + ll(p - n + 1) * F[n - 1][k]) % ip;
if (k)
S[n][k] = S[n - 1][k - 1];
S[n][k] = (S[n][k] + ll(k) * S[n - 1][k]) % ip;
}
}
int FirstKind(ll n, ll k) {
if (n < 0 or k < 0 or k > n)
return 0;
ll i = n / ip, j = n % ip;
if (k < i)
return 0;
ll a = (k - i) / (p - 1), b = (k - i) % (p - 1);
if (b == 0 and j)
b += p - 1, a--;
if (a < 0 or a > i or b > j)
return 0;
int res = (nCr(i, a) * F[j][b]) % ip;
if ((i + a) & 1)
res = (p - res) % ip;
return res;
}
int SecondKind(ll n, ll k) {
if (n < 0 or k < 0 or k > n)
return 0;
if (n == 0)
return 1;
ll i = k / ip, j = k % ip;
if (n < i)
return 0;
ll a = (n - i) / (p - 1), b = (n - i) % (p - 1);
if (b == 0)
b += p - 1, a--;
if (a < 0 or b < j)
return 0;
if (b == p - 1 and j == 0)
return nCr(a, i - 1);
else
return (nCr(a, i) * S[b][j]) % ip;
}
};
/**
* @brief Stirling Number for query
*/