#line 2 "matrix/matrix-tree.hpp"
#line 2 "matrix/matrix.hpp"
#line 2 "matrix/inverse-matrix.hpp"
#line 2 "matrix/gauss-elimination.hpp"
#include <utility>
#include <vector>
using namespace std;
// {rank, det(非正方行列の場合は未定義)} を返す
// 型が double や Rational でも動くはず？(未検証)
//
// pivot 候補 : [0, pivot_end)
template <typename T>
std::pair<int, T> GaussElimination(vector<vector<T>> &a, int pivot_end = -1,
bool diagonalize = false) {
int H = a.size(), W = a[0].size(), rank = 0;
if (pivot_end == -1) pivot_end = W;
T det = 1;
for (int j = 0; j < pivot_end; j++) {
int idx = -1;
for (int i = rank; i < H; i++) {
if (a[i][j] != T(0)) {
idx = i;
break;
}
}
if (idx == -1) {
det = 0;
continue;
}
if (rank != idx) det = -det, swap(a[rank], a[idx]);
det *= a[rank][j];
if (diagonalize && a[rank][j] != T(1)) {
T coeff = T(1) / a[rank][j];
for (int k = j; k < W; k++) a[rank][k] *= coeff;
}
int is = diagonalize ? 0 : rank + 1;
for (int i = is; i < H; i++) {
if (i == rank) continue;
if (a[i][j] != T(0)) {
T coeff = a[i][j] / a[rank][j];
for (int k = j; k < W; k++) a[i][k] -= a[rank][k] * coeff;
}
}
rank++;
}
return make_pair(rank, det);
}
#line 4 "matrix/inverse-matrix.hpp"
template <typename mint>
vector<vector<mint>> inverse_matrix(const vector<vector<mint>>& a) {
int N = a.size();
assert(N > 0);
assert(N == (int)a[0].size());
vector<vector<mint>> m(N, vector<mint>(2 * N));
for (int i = 0; i < N; i++) {
copy(begin(a[i]), end(a[i]), begin(m[i]));
m[i][N + i] = 1;
}
auto [rank, det] = GaussElimination(m, N, true);
if (rank != N) return {};
vector<vector<mint>> b(N);
for (int i = 0; i < N; i++) {
copy(begin(m[i]) + N, end(m[i]), back_inserter(b[i]));
}
return b;
}
#line 4 "matrix/matrix.hpp"
template <class T>
struct Matrix {
vector<vector<T> > A;
Matrix() = default;
Matrix(int n, int m) : A(n, vector<T>(m, T())) {}
Matrix(int n) : A(n, vector<T>(n, T())){};
int H() const { return A.size(); }
int W() const { return A[0].size(); }
int size() const { return A.size(); }
inline const vector<T> &operator[](int k) const { return A[k]; }
inline vector<T> &operator[](int k) { return A[k]; }
static Matrix I(int n) {
Matrix mat(n);
for (int i = 0; i < n; i++) mat[i][i] = 1;
return (mat);
}
Matrix &operator+=(const Matrix &B) {
int n = H(), m = W();
assert(n == B.H() && m == B.W());
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++) (*this)[i][j] += B[i][j];
return (*this);
}
Matrix &operator-=(const Matrix &B) {
int n = H(), m = W();
assert(n == B.H() && m == B.W());
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++) (*this)[i][j] -= B[i][j];
return (*this);
}
Matrix &operator*=(const Matrix &B) {
int n = H(), m = B.W(), p = W();
assert(p == B.H());
vector<vector<T> > C(n, vector<T>(m, T{}));
for (int i = 0; i < n; i++)
for (int k = 0; k < p; k++)
for (int j = 0; j < m; j++) C[i][j] += (*this)[i][k] * B[k][j];
A.swap(C);
return (*this);
}
Matrix &operator^=(long long k) {
Matrix B = Matrix::I(H());
while (k > 0) {
if (k & 1) B *= *this;
*this *= *this;
k >>= 1LL;
}
A.swap(B.A);
return (*this);
}
Matrix operator+(const Matrix &B) const { return (Matrix(*this) += B); }
Matrix operator-(const Matrix &B) const { return (Matrix(*this) -= B); }
Matrix operator*(const Matrix &B) const { return (Matrix(*this) *= B); }
Matrix operator^(const long long k) const { return (Matrix(*this) ^= k); }
bool operator==(const Matrix &B) const {
assert(H() == B.H() && W() == B.W());
for (int i = 0; i < H(); i++)
for (int j = 0; j < W(); j++)
if (A[i][j] != B[i][j]) return false;
return true;
}
bool operator!=(const Matrix &B) const {
assert(H() == B.H() && W() == B.W());
for (int i = 0; i < H(); i++)
for (int j = 0; j < W(); j++)
if (A[i][j] != B[i][j]) return true;
return false;
}
Matrix inverse() const {
assert(H() == W());
Matrix B(H());
B.A = inverse_matrix(A);
return B;
}
friend ostream &operator<<(ostream &os, const Matrix &p) {
int n = p.H(), m = p.W();
for (int i = 0; i < n; i++) {
os << (i ? " " : "") << "[";
for (int j = 0; j < m; j++) {
os << p[i][j] << (j + 1 == m ? "]\n" : ",");
}
}
return (os);
}
T determinant() const {
Matrix B(*this);
assert(H() == W());
T ret = 1;
for (int i = 0; i < H(); i++) {
int idx = -1;
for (int j = i; j < W(); j++) {
if (B[j][i] != 0) {
idx = j;
break;
}
}
if (idx == -1) return 0;
if (i != idx) {
ret *= T(-1);
swap(B[i], B[idx]);
}
ret *= B[i][i];
T inv = T(1) / B[i][i];
for (int j = 0; j < W(); j++) {
B[i][j] *= inv;
}
for (int j = i + 1; j < H(); j++) {
T a = B[j][i];
if (a == 0) continue;
for (int k = i; k < W(); k++) {
B[j][k] -= B[i][k] * a;
}
}
}
return ret;
}
};
/**
* @brief 行列ライブラリ
*/
#line 2 "matrix/polynomial-matrix-determinant.hpp"
#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 2 "fps/polynomial-interpolation.hpp"
#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 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 8 "matrix/polynomial-matrix-determinant.hpp"
template <typename mint>
FormalPowerSeries<mint> PolynomialMatrixDeterminant(
const Matrix<FormalPowerSeries<mint>> &m) {
int N = m.size();
int deg = 0;
for (int i = 0; i < N; ++i) deg += max<int>(1, m[i][i].size()) - 1;
vector<mint> xs(deg + 1);
vector<mint> ys(deg + 1);
Matrix<mint> M(N);
for (int x = 0; x <= deg; x++) {
xs[x] = x;
for (int i = 0; i < N; ++i)
for (int j = 0; j < N; ++j) M[i][j] = m[i][j].eval(x);
ys[x] = M.determinant();
}
return PolynomialInterpolation<mint>(xs, ys);
}
/**
* @brief 多項式行列の行列式
* @docs docs/matrix/polynomial-matrix-determinant.md
*/
#line 7 "matrix/matrix-tree.hpp"
template <typename T>
struct MatrixTree {
int n;
Matrix<T> m;
MatrixTree(int _n) : n(_n), m(_n) { assert(n > 0); }
void add(int i, int j, const T& x) {
if (i < n) m[i][i] += x;
if (j < n) m[j][j] += x;
if (i < n and j < n) {
m[i][j] -= x;
m[j][i] -= x;
}
}
Matrix<T> get() const { return m; }
template <typename U, typename = void>
struct has_value_type : false_type {};
template <typename U>
struct has_value_type<
U, typename conditional<false, typename U::value_type, void>::type>
: true_type {};
template <typename U = T,
enable_if_t<has_value_type<U>::value, nullptr_t> = nullptr>
T calc() {
return PolynomialMatrixDeterminant(m);
}
template <typename U = T,
enable_if_t<!has_value_type<U>::value, nullptr_t> = nullptr>
T calc() {
return m.determinant();
}
};
/**
* @brief 行列木定理(ラプラシアン行列)
*/
行列木定理(ラプラシアン行列)