#include "fps/multivariate-fps.hpp"
#pragma once #include <algorithm> #include <cassert> #include <vector> using namespace std; #include "../modulo/binomial.hpp" #include "../ntt/multivariate-multiplication.hpp" #include "formal-power-series.hpp" // FFT mod でないと使えない // f.size() != (base の総積) でも動くが, 二項演算は長さが等しい列同士しかダメ // 添え字アクセスは operator() を使う template <typename mint> struct MultivariateFormalPowerSeries { using mfps = MultivariateFormalPowerSeries<mint>; using fps = FormalPowerSeries<mint>; fps f; vector<int> base; MultivariateFormalPowerSeries() = default; MultivariateFormalPowerSeries(const vector<int>& _base) : base(_base) { int n = 1; for (auto& x : base) n *= x; f.resize(n); } MultivariateFormalPowerSeries(const fps& _f, const vector<int>& _base) : f(_f), base(_base) {} friend ostream& operator<<(ostream& os, const mfps& rhs) { os << "[ "; for (int i = 0; i < (int)rhs.f.size(); i++) { os << rhs.f[i] << (i + 1 == (int)rhs.f.size() ? "" : ", "); } return os << " ]"; } long long _id(int) { return 0; } template <typename Head, typename... Tail> long long _id(int i, Head&& head, Tail&&... tail) { assert(i < (int)base.size() && (int)head < base[i]); return head + _id(i + 1, forward<Tail>(tail)...) * base[i]; } template <typename... Args> long long id(Args&&... args) { return _id(0, forward<Args>(args)...); } template <typename... Args> mint& operator()(Args&&... args) { return f[id(forward<Args>(args)...)]; } mfps& operator+=(const mfps& rhs) { assert(base == rhs.base && f.size() == rhs.f.size()); for (int i = 0; i < (int)f.size(); i++) f[i] += rhs.f[i]; return *this; } mfps& operator-=(const mfps& rhs) { assert(base == rhs.base && f.size() == rhs.f.size()); for (int i = 0; i < (int)f.size(); i++) f[i] -= rhs.f[i]; return *this; } mfps& operator*=(const mfps& rhs) { assert(base == rhs.base && f.size() == rhs.f.size()); f = multivariate_multiplication(f, rhs.f, base); return *this; } mfps& operator/=(const mfps& rhs) { assert(base == rhs.base); assert(f.size() == rhs.f.size()); return (*this) *= rhs.inv(); } mfps& operator+=(const mint& rhs) { assert(!f.empty()); f[0] += rhs; return *this; } mfps& operator-=(const mint& rhs) { assert(!f.empty()); f[0] -= rhs; return *this; } mfps& operator*=(const mint& rhs) { for (auto& x : f) x *= rhs; return *this; } mfps& operator/=(const mint& rhs) { for (auto& x : f) x /= rhs; return *this; } mfps operator+(const mfps& rhs) const { return mfps{*this} += rhs; } mfps operator-(const mfps& rhs) const { return mfps{*this} -= rhs; } mfps operator*(const mfps& rhs) const { return mfps{*this} *= rhs; } mfps operator/(const mfps& rhs) const { return mfps{*this} /= rhs; } mfps operator+(const mint& rhs) const { return mfps{*this} += rhs; } mfps operator-(const mint& rhs) const { return mfps{*this} -= rhs; } mfps operator*(const mint& rhs) const { return mfps{*this} *= rhs; } mfps operator/(const mint& rhs) const { return mfps{*this} /= rhs; } mfps operator+() const { return mfps{*this}; } mfps operator-() const { return mfps{-f, base}; } friend bool operator==(const mfps& lhs, const mfps& rhs) { return lhs.f == rhs.f && lhs.base == rhs.base; } friend bool operator!=(const mfps& lhs, const mfps& rhs) { return lhs.f != rhs.f || lhs.base != rhs.base; } mfps diff() const { mfps g{*this}; for (int i = 0; i < (int)g.f.size(); i++) g.f[i] *= i; return g; } mfps integral() const { static Binomial<mint> binom; mfps g{*this}; for (int i = 1; i < (int)g.f.size(); i++) g.f[i] *= binom.inv(i); return g; } mfps inv() const { assert(f[0] != 0); if (base.empty()) return mfps{fps{f[0].inverse()}, base}; int n = f.size(), s = base.size(), W = 1; while (W < 2 * n) W *= 2; vector<int> chi(W); for (int i = 0; i < W; i++) { int x = i; for (int j = 0; j < s - 1; j++) chi[i] += (x /= base[j]); chi[i] %= s; } auto hadamard_prod = [&s](vector<fps>& F, vector<fps>& G, vector<fps>& H) { fps a(s); for (int k = 0; k < (int)F[0].size(); k++) { fill(begin(a), end(a), typename fps::value_type()); for (int i = 0; i < s; i++) for (int j = 0; j < s; j++) { a[i + j - (i + j >= s ? s : 0)] += F[i][k] * G[j][k]; } for (int i = 0; i < s; i++) H[i][k] = a[i]; } }; fps g(W); g[0] = f[0].inverse(); for (int d = 1; d < n; d *= 2) { vector<fps> F(s, fps(2 * d)), G(s, fps(2 * d)), H(s, fps(2 * d)); for (int j = 0; j < min((int)f.size(), 2 * d); j++) F[chi[j]][j] = f[j]; for (int j = 0; j < d; j++) G[chi[j]][j] = g[j]; for (auto& x : F) x.ntt(); for (auto& x : G) x.ntt(); hadamard_prod(F, G, H); for (auto& x : H) x.intt(); for (auto& x : F) fill(begin(x), end(x), typename fps::value_type()); for (int j = d; j < 2 * d; j++) F[chi[j]][j] = H[chi[j]][j]; for (auto& x : F) x.ntt(); hadamard_prod(F, G, H); for (auto& x : H) x.intt(); for (int j = d; j < 2 * d; j++) g[j] = -H[chi[j]][j]; } mfps res{*this}; res.f = fps{begin(g), begin(g) + n}; return res; } mfps log() const { assert(!f.empty() && f[0] == 1); return ((*this).diff() / *this).integral(); } mfps exp() const { assert(!f.empty() && f[0] == mint{0}); int n = f.size(); mfps g{fps{1}, base}; for (int d = 1; d < n; d *= 2) { int s = min(n, d * 2); g.f.resize(s, mint{0}); g *= mfps{fps{begin(f), begin(f) + s}, base} - g.log() + 1; } return g; } mfps pow(long long e) const { assert(!f.empty()); if (f[0] != 0) { mint f0inv = f[0].inverse(), coe = f[0].pow(e); return (((*this) * f0inv).log() * e).exp() * coe; } int n = f.size(); long long base_sum = 0; for (auto& b : base) base_sum += b - 1; if (e > base_sum) return mfps{fps(n), base}; mfps res{fps(n), base}, a{*this}; res.f[0] = 1; for (; e; a *= a, e >>= 1) { if (e & 1) res *= a; } return res; } }; /** * @brief 多変数形式的冪級数ライブラリ */
#line 2 "fps/multivariate-fps.hpp" #include <algorithm> #include <cassert> #include <vector> using namespace std; #line 2 "modulo/binomial.hpp" #line 4 "modulo/binomial.hpp" #include <type_traits> #line 6 "modulo/binomial.hpp" 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 "ntt/multivariate-multiplication.hpp" // 長さが等しい列同士の畳み込みしかしない template <typename fps> fps multivariate_multiplication(const fps& f, const fps& g, const vector<int>& base) { assert(f.size() == g.size()); int n = f.size(), s = base.size(), W = 1; if (s == 0) return fps{f[0] * g[0]}; while (W < 2 * n) W *= 2; vector<int> chi(n); for (int i = 0; i < n; i++) { int x = i; for (int j = 0; j < s - 1; j++) chi[i] += (x /= base[j]); chi[i] %= s; } vector<fps> F(s, fps(W)), G(s, fps(W)); for (int i = 0; i < n; i++) F[chi[i]][i] = f[i], G[chi[i]][i] = g[i]; for (auto& x : F) x.ntt(); for (auto& x : G) x.ntt(); fps a(s); for (int k = 0; k < W; k++) { fill(begin(a), end(a), typename fps::value_type()); for (int i = 0; i < s; i++) for (int j = 0; j < s; j++) { a[i + j - (i + j >= s ? s : 0)] += F[i][k] * G[j][k]; } for (int i = 0; i < s; i++) F[i][k] = a[i]; } for (auto& x : F) x.intt(); fps h(n); for (int i = 0; i < n; i++) h[i] = F[chi[i]][i]; return h; } /** * @brief Multivariate Multiplication * @docs docs/ntt/multivariate-multiplication.md */ #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 11 "fps/multivariate-fps.hpp" // FFT mod でないと使えない // f.size() != (base の総積) でも動くが, 二項演算は長さが等しい列同士しかダメ // 添え字アクセスは operator() を使う template <typename mint> struct MultivariateFormalPowerSeries { using mfps = MultivariateFormalPowerSeries<mint>; using fps = FormalPowerSeries<mint>; fps f; vector<int> base; MultivariateFormalPowerSeries() = default; MultivariateFormalPowerSeries(const vector<int>& _base) : base(_base) { int n = 1; for (auto& x : base) n *= x; f.resize(n); } MultivariateFormalPowerSeries(const fps& _f, const vector<int>& _base) : f(_f), base(_base) {} friend ostream& operator<<(ostream& os, const mfps& rhs) { os << "[ "; for (int i = 0; i < (int)rhs.f.size(); i++) { os << rhs.f[i] << (i + 1 == (int)rhs.f.size() ? "" : ", "); } return os << " ]"; } long long _id(int) { return 0; } template <typename Head, typename... Tail> long long _id(int i, Head&& head, Tail&&... tail) { assert(i < (int)base.size() && (int)head < base[i]); return head + _id(i + 1, forward<Tail>(tail)...) * base[i]; } template <typename... Args> long long id(Args&&... args) { return _id(0, forward<Args>(args)...); } template <typename... Args> mint& operator()(Args&&... args) { return f[id(forward<Args>(args)...)]; } mfps& operator+=(const mfps& rhs) { assert(base == rhs.base && f.size() == rhs.f.size()); for (int i = 0; i < (int)f.size(); i++) f[i] += rhs.f[i]; return *this; } mfps& operator-=(const mfps& rhs) { assert(base == rhs.base && f.size() == rhs.f.size()); for (int i = 0; i < (int)f.size(); i++) f[i] -= rhs.f[i]; return *this; } mfps& operator*=(const mfps& rhs) { assert(base == rhs.base && f.size() == rhs.f.size()); f = multivariate_multiplication(f, rhs.f, base); return *this; } mfps& operator/=(const mfps& rhs) { assert(base == rhs.base); assert(f.size() == rhs.f.size()); return (*this) *= rhs.inv(); } mfps& operator+=(const mint& rhs) { assert(!f.empty()); f[0] += rhs; return *this; } mfps& operator-=(const mint& rhs) { assert(!f.empty()); f[0] -= rhs; return *this; } mfps& operator*=(const mint& rhs) { for (auto& x : f) x *= rhs; return *this; } mfps& operator/=(const mint& rhs) { for (auto& x : f) x /= rhs; return *this; } mfps operator+(const mfps& rhs) const { return mfps{*this} += rhs; } mfps operator-(const mfps& rhs) const { return mfps{*this} -= rhs; } mfps operator*(const mfps& rhs) const { return mfps{*this} *= rhs; } mfps operator/(const mfps& rhs) const { return mfps{*this} /= rhs; } mfps operator+(const mint& rhs) const { return mfps{*this} += rhs; } mfps operator-(const mint& rhs) const { return mfps{*this} -= rhs; } mfps operator*(const mint& rhs) const { return mfps{*this} *= rhs; } mfps operator/(const mint& rhs) const { return mfps{*this} /= rhs; } mfps operator+() const { return mfps{*this}; } mfps operator-() const { return mfps{-f, base}; } friend bool operator==(const mfps& lhs, const mfps& rhs) { return lhs.f == rhs.f && lhs.base == rhs.base; } friend bool operator!=(const mfps& lhs, const mfps& rhs) { return lhs.f != rhs.f || lhs.base != rhs.base; } mfps diff() const { mfps g{*this}; for (int i = 0; i < (int)g.f.size(); i++) g.f[i] *= i; return g; } mfps integral() const { static Binomial<mint> binom; mfps g{*this}; for (int i = 1; i < (int)g.f.size(); i++) g.f[i] *= binom.inv(i); return g; } mfps inv() const { assert(f[0] != 0); if (base.empty()) return mfps{fps{f[0].inverse()}, base}; int n = f.size(), s = base.size(), W = 1; while (W < 2 * n) W *= 2; vector<int> chi(W); for (int i = 0; i < W; i++) { int x = i; for (int j = 0; j < s - 1; j++) chi[i] += (x /= base[j]); chi[i] %= s; } auto hadamard_prod = [&s](vector<fps>& F, vector<fps>& G, vector<fps>& H) { fps a(s); for (int k = 0; k < (int)F[0].size(); k++) { fill(begin(a), end(a), typename fps::value_type()); for (int i = 0; i < s; i++) for (int j = 0; j < s; j++) { a[i + j - (i + j >= s ? s : 0)] += F[i][k] * G[j][k]; } for (int i = 0; i < s; i++) H[i][k] = a[i]; } }; fps g(W); g[0] = f[0].inverse(); for (int d = 1; d < n; d *= 2) { vector<fps> F(s, fps(2 * d)), G(s, fps(2 * d)), H(s, fps(2 * d)); for (int j = 0; j < min((int)f.size(), 2 * d); j++) F[chi[j]][j] = f[j]; for (int j = 0; j < d; j++) G[chi[j]][j] = g[j]; for (auto& x : F) x.ntt(); for (auto& x : G) x.ntt(); hadamard_prod(F, G, H); for (auto& x : H) x.intt(); for (auto& x : F) fill(begin(x), end(x), typename fps::value_type()); for (int j = d; j < 2 * d; j++) F[chi[j]][j] = H[chi[j]][j]; for (auto& x : F) x.ntt(); hadamard_prod(F, G, H); for (auto& x : H) x.intt(); for (int j = d; j < 2 * d; j++) g[j] = -H[chi[j]][j]; } mfps res{*this}; res.f = fps{begin(g), begin(g) + n}; return res; } mfps log() const { assert(!f.empty() && f[0] == 1); return ((*this).diff() / *this).integral(); } mfps exp() const { assert(!f.empty() && f[0] == mint{0}); int n = f.size(); mfps g{fps{1}, base}; for (int d = 1; d < n; d *= 2) { int s = min(n, d * 2); g.f.resize(s, mint{0}); g *= mfps{fps{begin(f), begin(f) + s}, base} - g.log() + 1; } return g; } mfps pow(long long e) const { assert(!f.empty()); if (f[0] != 0) { mint f0inv = f[0].inverse(), coe = f[0].pow(e); return (((*this) * f0inv).log() * e).exp() * coe; } int n = f.size(); long long base_sum = 0; for (auto& b : base) base_sum += b - 1; if (e > base_sum) return mfps{fps(n), base}; mfps res{fps(n), base}, a{*this}; res.f[0] = 1; for (; e; a *= a, e >>= 1) { if (e & 1) res *= a; } return res; } }; /** * @brief 多変数形式的冪級数ライブラリ */