#include "fps/stirling-matrix.hpp"
#pragma once #include "../modulo/binomial.hpp" #include "composite-exp.hpp" #include "formal-power-series.hpp" #include "multipoint-evaluation.hpp" #include "pascal-matrix.hpp" #include "polynomial-interpolation.hpp" // S_{i, j} = stirling{i, j} を満たす行列 S を縦ベクトルに作用 template <typename mint> FormalPowerSeries<mint> stirling_matrix(FormalPowerSeries<mint> a, int rev = false) { using fps = FormalPowerSeries<mint>; if (a.empty()) return {}; int N = a.size(); Binomial<mint> binom(N + 10); if (rev == false) { for (int i = 0; i < N; i++) a[i] *= binom.finv(i); fps f = pascal_matrix_trans(a, true); fps b = composite_exp<mint>(f, 1, N); for (int i = 0; i < N; i++) b[i] *= binom.fac(i); return b; } else { for (int i = 0; i < N; i++) a[i] *= binom.finv(i); fps f = inverse_of_composite_exp<mint>(a, 1); fps b = pascal_matrix_trans(f, false); for (int i = 0; i < N; i++) b[i] *= binom.fac(i); return b; } } // S_{i, j} = stirling{j, i} を満たす行列 S を縦ベクトルに作用 template <typename mint> FormalPowerSeries<mint> stirling_matrix_trans(FormalPowerSeries<mint> a, int rev = false) { using fps = FormalPowerSeries<mint>; if (a.empty()) return {}; int N = a.size(); Binomial<mint> binom(N + 10); vector<mint> xs(N); for (int i = 0; i < N; i++) xs[i] = i; if (rev == false) { auto _f = MultipointEvaluation(a, xs); fps f{begin(_f), end(_f)}; fps g = pascal_matrix(f, true); for (int i = 0; i < N; i++) g[i] *= binom.finv(i); return g; } else { for (int i = 0; i < N; i++) a[i] *= binom.fac(i); auto g = pascal_matrix(a, false); return PolynomialInterpolation(xs, g); } }
#line 2 "fps/stirling-matrix.hpp" #line 2 "modulo/binomial.hpp" #include <cassert> #include <type_traits> #include <vector> using namespace std; // コンストラクタの MAX に 「C(n, r) や fac(n) でクエリを投げる最大の n 」 // を入れると倍速くらいになる // mod を超えて前計算して 0 割りを踏むバグは対策済み template <typename T> struct Binomial { vector<T> f, g, h; Binomial(int MAX = 0) { assert(T::get_mod() != 0 && "Binomial<mint>()"); f.resize(1, T{1}); g.resize(1, T{1}); h.resize(1, T{1}); if (MAX > 0) extend(MAX + 1); } void extend(int m = -1) { int n = f.size(); if (m == -1) m = n * 2; m = min<int>(m, T::get_mod()); if (n >= m) return; f.resize(m); g.resize(m); h.resize(m); for (int i = n; i < m; i++) f[i] = f[i - 1] * T(i); g[m - 1] = f[m - 1].inverse(); h[m - 1] = g[m - 1] * f[m - 2]; for (int i = m - 2; i >= n; i--) { g[i] = g[i + 1] * T(i + 1); h[i] = g[i] * f[i - 1]; } } T fac(int i) { if (i < 0) return T(0); while (i >= (int)f.size()) extend(); return f[i]; } T finv(int i) { if (i < 0) return T(0); while (i >= (int)g.size()) extend(); return g[i]; } T inv(int i) { if (i < 0) return -inv(-i); while (i >= (int)h.size()) extend(); return h[i]; } T C(int n, int r) { if (n < 0 || n < r || r < 0) return T(0); return fac(n) * finv(n - r) * finv(r); } inline T operator()(int n, int r) { return C(n, r); } template <typename I> T multinomial(const vector<I>& r) { static_assert(is_integral<I>::value == true); int n = 0; for (auto& x : r) { if (x < 0) return T(0); n += x; } T res = fac(n); for (auto& x : r) res *= finv(x); return res; } template <typename I> T operator()(const vector<I>& r) { return multinomial(r); } T C_naive(int n, int r) { if (n < 0 || n < r || r < 0) return T(0); T ret = T(1); r = min(r, n - r); for (int i = 1; i <= r; ++i) ret *= inv(i) * (n--); return ret; } T P(int n, int r) { if (n < 0 || n < r || r < 0) return T(0); return fac(n) * finv(n - r); } // [x^r] 1 / (1-x)^n T H(int n, int r) { if (n < 0 || r < 0) return T(0); return r == 0 ? 1 : C(n + r - 1, r); } }; #line 2 "fps/composite-exp.hpp" #line 4 "fps/composite-exp.hpp" #include <utility> #line 6 "fps/composite-exp.hpp" using namespace std; #line 2 "fps/formal-power-series.hpp" template <typename mint> struct FormalPowerSeries : vector<mint> { using vector<mint>::vector; using FPS = FormalPowerSeries; FPS &operator+=(const FPS &r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i]; return *this; } FPS &operator+=(const mint &r) { if (this->empty()) this->resize(1); (*this)[0] += r; return *this; } FPS &operator-=(const FPS &r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i]; return *this; } FPS &operator-=(const mint &r) { if (this->empty()) this->resize(1); (*this)[0] -= r; return *this; } FPS &operator*=(const mint &v) { for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v; return *this; } FPS &operator/=(const FPS &r) { if (this->size() < r.size()) { this->clear(); return *this; } int n = this->size() - r.size() + 1; if ((int)r.size() <= 64) { FPS f(*this), g(r); g.shrink(); mint coeff = g.back().inverse(); for (auto &x : g) x *= coeff; int deg = (int)f.size() - (int)g.size() + 1; int gs = g.size(); FPS quo(deg); for (int i = deg - 1; i >= 0; i--) { quo[i] = f[i + gs - 1]; for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j]; } *this = quo * coeff; this->resize(n, mint(0)); return *this; } return *this = ((*this).rev().pre(n) * r.rev().inv(n)).pre(n).rev(); } FPS &operator%=(const FPS &r) { *this -= *this / r * r; shrink(); return *this; } FPS operator+(const FPS &r) const { return FPS(*this) += r; } FPS operator+(const mint &v) const { return FPS(*this) += v; } FPS operator-(const FPS &r) const { return FPS(*this) -= r; } FPS operator-(const mint &v) const { return FPS(*this) -= v; } FPS operator*(const FPS &r) const { return FPS(*this) *= r; } FPS operator*(const mint &v) const { return FPS(*this) *= v; } FPS operator/(const FPS &r) const { return FPS(*this) /= r; } FPS operator%(const FPS &r) const { return FPS(*this) %= r; } FPS operator-() const { FPS ret(this->size()); for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i]; return ret; } void shrink() { while (this->size() && this->back() == mint(0)) this->pop_back(); } FPS rev() const { FPS ret(*this); reverse(begin(ret), end(ret)); return ret; } FPS dot(FPS r) const { FPS ret(min(this->size(), r.size())); for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i]; return ret; } // 前 sz 項を取ってくる。sz に足りない項は 0 埋めする FPS pre(int sz) const { FPS ret(begin(*this), begin(*this) + min((int)this->size(), sz)); if ((int)ret.size() < sz) ret.resize(sz); return ret; } FPS operator>>(int sz) const { if ((int)this->size() <= sz) return {}; FPS ret(*this); ret.erase(ret.begin(), ret.begin() + sz); return ret; } FPS operator<<(int sz) const { FPS ret(*this); ret.insert(ret.begin(), sz, mint(0)); return ret; } FPS diff() const { const int n = (int)this->size(); FPS ret(max(0, n - 1)); mint one(1), coeff(1); for (int i = 1; i < n; i++) { ret[i - 1] = (*this)[i] * coeff; coeff += one; } return ret; } FPS integral() const { const int n = (int)this->size(); FPS ret(n + 1); ret[0] = mint(0); if (n > 0) ret[1] = mint(1); auto mod = mint::get_mod(); for (int i = 2; i <= n; i++) ret[i] = (-ret[mod % i]) * (mod / i); for (int i = 0; i < n; i++) ret[i + 1] *= (*this)[i]; return ret; } mint eval(mint x) const { mint r = 0, w = 1; for (auto &v : *this) r += w * v, w *= x; return r; } FPS log(int deg = -1) const { assert(!(*this).empty() && (*this)[0] == mint(1)); if (deg == -1) deg = (int)this->size(); return (this->diff() * this->inv(deg)).pre(deg - 1).integral(); } FPS pow(int64_t k, int deg = -1) const { const int n = (int)this->size(); if (deg == -1) deg = n; if (k == 0) { FPS ret(deg); if (deg) ret[0] = 1; return ret; } for (int i = 0; i < n; i++) { if ((*this)[i] != mint(0)) { mint rev = mint(1) / (*this)[i]; FPS ret = (((*this * rev) >> i).log(deg) * k).exp(deg); ret *= (*this)[i].pow(k); ret = (ret << (i * k)).pre(deg); if ((int)ret.size() < deg) ret.resize(deg, mint(0)); return ret; } if (__int128_t(i + 1) * k >= deg) return FPS(deg, mint(0)); } return FPS(deg, mint(0)); } static void *ntt_ptr; static void set_fft(); FPS &operator*=(const FPS &r); void ntt(); void intt(); void ntt_doubling(); static int ntt_pr(); FPS inv(int deg = -1) const; FPS exp(int deg = -1) const; }; template <typename mint> void *FormalPowerSeries<mint>::ntt_ptr = nullptr; /** * @brief 多項式/形式的冪級数ライブラリ * @docs docs/fps/formal-power-series.md */ #line 9 "fps/composite-exp.hpp" // 多項式 f に exp(cx) 代入 // 次数 : mod x^{deg} まで計算, 指定がない場合 f と同じ長さ計算 template <typename mint> FormalPowerSeries<mint> composite_exp(FormalPowerSeries<mint> f, mint c = 1, int deg = -1) { using fps = FormalPowerSeries<mint>; assert(c != 0); if (deg == -1) deg = f.size(); if (f.empty()) return {}; int N = f.size(); vector<pair<fps, fps>> fs; for (int i = 0; i < N; i++) fs.emplace_back(fps{f[i]}, fps{1, -c * i}); while (fs.size() > 1u) { vector<pair<fps, fps>> nx; for (int i = 0; i + 1 < (int)fs.size(); i += 2) { pair<fps, fps>& f0 = fs[i]; pair<fps, fps>& f1 = fs[i + 1]; fps s = f0.first * f1.second + f1.first * f0.second; fps t = f0.second * f1.second; nx.emplace_back(s, t); } if (fs.size() % 2) nx.push_back(fs.back()); fs = nx; } fps g = (fs[0].first * fs[0].second.inv(deg)).pre(deg); mint b = 1; for (int i = 0; i < deg; i++) g[i] *= b, b /= i + 1; return g; } // 入力 f(x) = sum_{0 <= k < N} a_i exp(ckx) を満たす g(x) (mod x^N) // 出力 a(x) = sum_{0 <= k < N} a_i x^i template <typename mint> FormalPowerSeries<mint> inverse_of_composite_exp(FormalPowerSeries<mint> f, mint c = 1) { using fps = FormalPowerSeries<mint>; if (f.empty()) return {}; int N = f.size(); mint b = 1; for (int i = 0; i < N; i++) f[i] *= b, b *= i + 1; int B = 1; while (B < N) B *= 2; vector<fps> mod(2 * B, fps{1}); for (int i = 0; i < N; i++) mod[B + i] = fps{-c * i, 1}; for (int i = B - 1; i; i--) mod[i] = mod[2 * i] * mod[2 * i + 1]; fps denom = mod[1].rev(); fps numer = (f * denom).pre(N); vector<mint> a(N); auto dfs = [&](auto rc, int i, int l, int r, fps g) -> void { if (N <= l) return; if (l + 1 == r) { a[l] = g.eval(0); return; } int m = (l + r) / 2; rc(rc, i * 2 + 0, l, m, g % mod[i * 2 + 0]); rc(rc, i * 2 + 1, m, r, g % mod[i * 2 + 1]); }; dfs(dfs, 1, 0, B, numer.rev()); vector<mint> fac(N); fac[0] = 1; for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * c * i; for (int i = 0; i < N; i++) { a[i] /= fac[N - 1 - i] * fac[i] * ((N - 1 - i) % 2 ? -1 : 1); } return fps{begin(a), end(a)}; } /** * @brief $f(exp(cx))$ の計算 */ #line 2 "fps/multipoint-evaluation.hpp" #line 4 "fps/multipoint-evaluation.hpp" template <typename mint> struct ProductTree { using fps = FormalPowerSeries<mint>; const vector<mint> &xs; vector<fps> buf; int N, xsz; vector<int> l, r; ProductTree(const vector<mint> &xs_) : xs(xs_), xsz(xs.size()) { N = 1; while (N < (int)xs.size()) N *= 2; buf.resize(2 * N); l.resize(2 * N, xs.size()); r.resize(2 * N, xs.size()); fps::set_fft(); if (fps::ntt_ptr == nullptr) build(); else build_ntt(); } void build() { for (int i = 0; i < xsz; i++) { l[i + N] = i; r[i + N] = i + 1; buf[i + N] = {-xs[i], 1}; } for (int i = N - 1; i > 0; i--) { l[i] = l[(i << 1) | 0]; r[i] = r[(i << 1) | 1]; if (buf[(i << 1) | 0].empty()) continue; else if (buf[(i << 1) | 1].empty()) buf[i] = buf[(i << 1) | 0]; else buf[i] = buf[(i << 1) | 0] * buf[(i << 1) | 1]; } } void build_ntt() { fps f; f.reserve(N * 2); for (int i = 0; i < xsz; i++) { l[i + N] = i; r[i + N] = i + 1; buf[i + N] = {-xs[i] + 1, -xs[i] - 1}; } for (int i = N - 1; i > 0; i--) { l[i] = l[(i << 1) | 0]; r[i] = r[(i << 1) | 1]; if (buf[(i << 1) | 0].empty()) continue; else if (buf[(i << 1) | 1].empty()) buf[i] = buf[(i << 1) | 0]; else if (buf[(i << 1) | 0].size() == buf[(i << 1) | 1].size()) { buf[i] = buf[(i << 1) | 0]; f.clear(); copy(begin(buf[(i << 1) | 1]), end(buf[(i << 1) | 1]), back_inserter(f)); buf[i].ntt_doubling(); f.ntt_doubling(); for (int j = 0; j < (int)buf[i].size(); j++) buf[i][j] *= f[j]; } else { buf[i] = buf[(i << 1) | 0]; f.clear(); copy(begin(buf[(i << 1) | 1]), end(buf[(i << 1) | 1]), back_inserter(f)); buf[i].ntt_doubling(); f.intt(); f.resize(buf[i].size(), mint(0)); f.ntt(); for (int j = 0; j < (int)buf[i].size(); j++) buf[i][j] *= f[j]; } } for (int i = 0; i < 2 * N; i++) { buf[i].intt(); buf[i].shrink(); } } }; template <typename mint> vector<mint> InnerMultipointEvaluation(const FormalPowerSeries<mint> &f, const vector<mint> &xs, const ProductTree<mint> &ptree) { using fps = FormalPowerSeries<mint>; vector<mint> ret; ret.reserve(xs.size()); auto rec = [&](auto self, fps a, int idx) { if (ptree.l[idx] == ptree.r[idx]) return; a %= ptree.buf[idx]; if ((int)a.size() <= 64) { for (int i = ptree.l[idx]; i < ptree.r[idx]; i++) ret.push_back(a.eval(xs[i])); return; } self(self, a, (idx << 1) | 0); self(self, a, (idx << 1) | 1); }; rec(rec, f, 1); return ret; } template <typename mint> vector<mint> MultipointEvaluation(const FormalPowerSeries<mint> &f, const vector<mint> &xs) { if(f.empty() || xs.empty()) return vector<mint>(xs.size(), mint(0)); return InnerMultipointEvaluation(f, xs, ProductTree<mint>(xs)); } /** * @brief Multipoint Evaluation */ #line 2 "fps/pascal-matrix.hpp" #line 5 "fps/pascal-matrix.hpp" // P_{i, j} = binom(i, j) を満たす行列 P を縦ベクトルに作用 template <typename mint> FormalPowerSeries<mint> pascal_matrix(FormalPowerSeries<mint> a, int rev = false) { using fps = FormalPowerSeries<mint>; if (a.empty()) return {}; int N = a.size(); Binomial<mint> binom(N + 10); if (rev == false) { fps e(N); for (int i = 0; i < N; i++) { a[i] *= binom.finv(i); e[i] = binom.finv(i); } fps b = (a * e).pre(N); for (int i = 0; i < N; i++) b[i] *= binom.fac(i); return b; } else { fps ie(N); for (int i = 0; i < N; i++) { a[i] *= binom.finv(i); ie[i] = binom.finv(i) * (i % 2 ? -1 : 1); } fps b = (a * ie).pre(N); for (int i = 0; i < N; i++) b[i] *= binom.fac(i); return b; } } // P_{i, j} = binom(j, i) を満たす行列 P を縦ベクトルに作用 template <typename mint> FormalPowerSeries<mint> pascal_matrix_trans(FormalPowerSeries<mint> a, int rev = false) { using fps = FormalPowerSeries<mint>; if (a.empty()) return {}; int N = a.size(); Binomial<mint> binom(N + 10); if (rev == false) { fps e(N); for (int i = 0; i < N; i++) { a[i] *= binom.fac(i); e[i] = binom.finv(i); } fps b = (a.rev() * e).pre(N).rev(); for (int i = 0; i < N; i++) b[i] *= binom.finv(i); return b; } else { fps ie(N); for (int i = 0; i < N; i++) { a[i] *= binom.fac(i); ie[i] = binom.finv(i) * (i % 2 ? -1 : 1); } fps b = (a.rev() * ie).pre(N).rev(); for (int i = 0; i < N; i++) b[i] *= binom.finv(i); return b; } } #line 2 "fps/polynomial-interpolation.hpp" #line 5 "fps/polynomial-interpolation.hpp" template <class mint> FormalPowerSeries<mint> PolynomialInterpolation(const vector<mint> &xs, const vector<mint> &ys) { using fps = FormalPowerSeries<mint>; assert(xs.size() == ys.size()); ProductTree<mint> ptree(xs); fps w = ptree.buf[1].diff(); vector<mint> vs = InnerMultipointEvaluation<mint>(w, xs, ptree); auto rec = [&](auto self, int idx) -> fps { if (idx >= ptree.N) { if (idx - ptree.N < (int)xs.size()) return {ys[idx - ptree.N] / vs[idx - ptree.N]}; else return {mint(1)}; } if (ptree.buf[idx << 1 | 0].empty()) return {}; else if (ptree.buf[idx << 1 | 1].empty()) return self(self, idx << 1 | 0); return self(self, idx << 1 | 0) * ptree.buf[idx << 1 | 1] + self(self, idx << 1 | 1) * ptree.buf[idx << 1 | 0]; }; return rec(rec, 1); } #line 9 "fps/stirling-matrix.hpp" // S_{i, j} = stirling{i, j} を満たす行列 S を縦ベクトルに作用 template <typename mint> FormalPowerSeries<mint> stirling_matrix(FormalPowerSeries<mint> a, int rev = false) { using fps = FormalPowerSeries<mint>; if (a.empty()) return {}; int N = a.size(); Binomial<mint> binom(N + 10); if (rev == false) { for (int i = 0; i < N; i++) a[i] *= binom.finv(i); fps f = pascal_matrix_trans(a, true); fps b = composite_exp<mint>(f, 1, N); for (int i = 0; i < N; i++) b[i] *= binom.fac(i); return b; } else { for (int i = 0; i < N; i++) a[i] *= binom.finv(i); fps f = inverse_of_composite_exp<mint>(a, 1); fps b = pascal_matrix_trans(f, false); for (int i = 0; i < N; i++) b[i] *= binom.fac(i); return b; } } // S_{i, j} = stirling{j, i} を満たす行列 S を縦ベクトルに作用 template <typename mint> FormalPowerSeries<mint> stirling_matrix_trans(FormalPowerSeries<mint> a, int rev = false) { using fps = FormalPowerSeries<mint>; if (a.empty()) return {}; int N = a.size(); Binomial<mint> binom(N + 10); vector<mint> xs(N); for (int i = 0; i < N; i++) xs[i] = i; if (rev == false) { auto _f = MultipointEvaluation(a, xs); fps f{begin(_f), end(_f)}; fps g = pascal_matrix(f, true); for (int i = 0; i < N; i++) g[i] *= binom.finv(i); return g; } else { for (int i = 0; i < N; i++) a[i] *= binom.fac(i); auto g = pascal_matrix(a, false); return PolynomialInterpolation(xs, g); } }