This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub maspypy/library
#include "seq/famous/stirling_number_query.hpp"
#include "nt/primetest.hpp" // O(p^2) 時間の前計算のもと、O(log n) 時間 struct Stirling_Number_Query { const int p; vvc<int> MEMO_C; vvc<int> MEMO_S1; vvc<int> MEMO_S2; Stirling_Number_Query(int p, bool first_kind = true, bool second_kind = true) : p(p) { assert(primetest(p)); assert(p <= (1 << 15)); build_C(); if (first_kind) build_S1(); if (second_kind) build_S2(); } int C(ll n, ll k) { if (k < 0 || k > n) return 0; int res = 1; while (n) { int i = n % p, j = k % p; if (j > i) return 0; res = res * MEMO_C[i][j] % p; n /= p; k /= p; } return res; } int S1(ll n, ll k) { if (k < 0 || k > n) return 0; ll i = n / p; int j = n % p; if (i > k) return 0; ll a = (k - i) / (p - 1); int b = (k - i) % (p - 1); if (b == 0 && j > 0) { b += (p - 1); a -= 1; } if (a < 0 || i < a || b > j) return 0; int x = C(i, a); int y = MEMO_S1[j][b]; int res = x * y % p; if ((i + a) % 2 == 1 && res) { res = p - res; } return res; } int S2(ll n, ll k) { if (k < 0 || k > n) return 0; if (n == 0) return 1; ll i = k / p; int j = k % p; if (n < i) return 0; ll a = (n - i) / (p - 1); int b = (n - i) - (p - 1) * a; if (b == 0) { b += p - 1; a -= 1; } if (a < 0 || j > b) return 0; if (b < p - 1) { return C(a, i) * MEMO_S2[b][j] % p; } if (j == 0) return C(a, i - 1); return C(a, i) * MEMO_S2[p - 1][j] % p; } private: void build_C() { auto& A = MEMO_C; A.resize(p); A[0] = {1}; FOR(i, 1, p) { A[i] = A[i - 1]; A[i].emplace_back(0); FOR(j, 1, i + 1) { A[i][j] += A[i - 1][j - 1]; if (A[i][j] >= p) A[i][j] -= p; } } } void build_S1() { auto& A = MEMO_S1; A.resize(p); A[0] = {1}; FOR(i, 1, p) { A[i].assign(i + 1, 0); FOR(j, i + 1) { if (j) A[i][j] += A[i - 1][j - 1]; if (j < i) A[i][j] += A[i - 1][j] * (p - i + 1); A[i][j] %= p; } } } void build_S2() { auto& A = MEMO_S2; A.resize(p); A[0] = {1}; FOR(i, 1, p) { A[i].assign(i + 1, 0); FOR(j, i + 1) { if (j) A[i][j] += A[i - 1][j - 1]; if (j < i) A[i][j] += A[i - 1][j] * j; A[i][j] %= p; } } } };
#line 2 "mod/mongomery_modint.hpp" // odd mod. // x の代わりに rx を持つ template <int id, typename U1, typename U2> struct Mongomery_modint { using mint = Mongomery_modint; inline static U1 m, r, n2; static constexpr int W = numeric_limits<U1>::digits; static void set_mod(U1 mod) { assert(mod & 1 && mod <= U1(1) << (W - 2)); m = mod, n2 = -U2(m) % m, r = m; FOR(5) r *= 2 - m * r; r = -r; assert(r * m == U1(-1)); } static U1 reduce(U2 b) { return (b + U2(U1(b) * r) * m) >> W; } U1 x; Mongomery_modint() : x(0) {} Mongomery_modint(U1 x) : x(reduce(U2(x) * n2)){}; U1 val() const { U1 y = reduce(x); return y >= m ? y - m : y; } mint &operator+=(mint y) { x = ((x += y.x) >= m ? x - m : x); return *this; } mint &operator-=(mint y) { x -= (x >= y.x ? y.x : y.x - m); return *this; } mint &operator*=(mint y) { x = reduce(U2(x) * y.x); return *this; } mint operator+(mint y) const { return mint(*this) += y; } mint operator-(mint y) const { return mint(*this) -= y; } mint operator*(mint y) const { return mint(*this) *= y; } bool operator==(mint y) const { return (x >= m ? x - m : x) == (y.x >= m ? y.x - m : y.x); } bool operator!=(mint y) const { return not operator==(y); } mint pow(ll n) const { assert(n >= 0); mint y = 1, z = *this; for (; n; n >>= 1, z *= z) if (n & 1) y *= z; return y; } }; template <int id> using Mongomery_modint_32 = Mongomery_modint<id, u32, u64>; template <int id> using Mongomery_modint_64 = Mongomery_modint<id, u64, u128>; #line 3 "nt/primetest.hpp" bool primetest(const u64 x) { assert(x < u64(1) << 62); if (x == 2 or x == 3 or x == 5 or x == 7) return true; if (x % 2 == 0 or x % 3 == 0 or x % 5 == 0 or x % 7 == 0) return false; if (x < 121) return x > 1; const u64 d = (x - 1) >> lowbit(x - 1); using mint = Mongomery_modint_64<202311020>; mint::set_mod(x); const mint one(u64(1)), minus_one(x - 1); auto ok = [&](u64 a) -> bool { auto y = mint(a).pow(d); u64 t = d; while (y != one && y != minus_one && t != x - 1) y *= y, t <<= 1; if (y != minus_one && t % 2 == 0) return false; return true; }; if (x < (u64(1) << 32)) { for (u64 a: {2, 7, 61}) if (!ok(a)) return false; } else { for (u64 a: {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) { if (!ok(a)) return false; } } return true; } #line 2 "seq/famous/stirling_number_query.hpp" // O(p^2) 時間の前計算のもと、O(log n) 時間 struct Stirling_Number_Query { const int p; vvc<int> MEMO_C; vvc<int> MEMO_S1; vvc<int> MEMO_S2; Stirling_Number_Query(int p, bool first_kind = true, bool second_kind = true) : p(p) { assert(primetest(p)); assert(p <= (1 << 15)); build_C(); if (first_kind) build_S1(); if (second_kind) build_S2(); } int C(ll n, ll k) { if (k < 0 || k > n) return 0; int res = 1; while (n) { int i = n % p, j = k % p; if (j > i) return 0; res = res * MEMO_C[i][j] % p; n /= p; k /= p; } return res; } int S1(ll n, ll k) { if (k < 0 || k > n) return 0; ll i = n / p; int j = n % p; if (i > k) return 0; ll a = (k - i) / (p - 1); int b = (k - i) % (p - 1); if (b == 0 && j > 0) { b += (p - 1); a -= 1; } if (a < 0 || i < a || b > j) return 0; int x = C(i, a); int y = MEMO_S1[j][b]; int res = x * y % p; if ((i + a) % 2 == 1 && res) { res = p - res; } return res; } int S2(ll n, ll k) { if (k < 0 || k > n) return 0; if (n == 0) return 1; ll i = k / p; int j = k % p; if (n < i) return 0; ll a = (n - i) / (p - 1); int b = (n - i) - (p - 1) * a; if (b == 0) { b += p - 1; a -= 1; } if (a < 0 || j > b) return 0; if (b < p - 1) { return C(a, i) * MEMO_S2[b][j] % p; } if (j == 0) return C(a, i - 1); return C(a, i) * MEMO_S2[p - 1][j] % p; } private: void build_C() { auto& A = MEMO_C; A.resize(p); A[0] = {1}; FOR(i, 1, p) { A[i] = A[i - 1]; A[i].emplace_back(0); FOR(j, 1, i + 1) { A[i][j] += A[i - 1][j - 1]; if (A[i][j] >= p) A[i][j] -= p; } } } void build_S1() { auto& A = MEMO_S1; A.resize(p); A[0] = {1}; FOR(i, 1, p) { A[i].assign(i + 1, 0); FOR(j, i + 1) { if (j) A[i][j] += A[i - 1][j - 1]; if (j < i) A[i][j] += A[i - 1][j] * (p - i + 1); A[i][j] %= p; } } } void build_S2() { auto& A = MEMO_S2; A.resize(p); A[0] = {1}; FOR(i, 1, p) { A[i].assign(i + 1, 0); FOR(j, i + 1) { if (j) A[i][j] += A[i - 1][j - 1]; if (j < i) A[i][j] += A[i - 1][j] * j; A[i][j] %= p; } } } };