模板

namespace flow {
int h[maxn], e[maxn << 1], ne[maxn << 1], idx = 1;
int cph[maxn], w[maxn << 1];
void clear () {
	memset (h, 0, sizeof(int) * (idx + 1)), idx = 1;
}
void add (int a, int b, int c) {
	idx++, e[idx] = b, ne[idx] = h[a], w[idx] = c, h[a] = idx;
	idx++, e[idx] = a, ne[idx] = h[b], w[idx] = 0, h[b] = idx;
}
int dist[maxn]; queue<int> q;
bool bfs (int N, int S, int T) {
	memset (dist, 0, sizeof(int) * (N + 1));
	dist[S] = 1; q.push (S);
	while (q.size()) {
		int u = q.front(); q.pop();
		for (int i = cph[u], v; i; i = ne[i])
			(!dist[(v = e[i])] && w[i] > 0 ? q.push (v), dist[v] = dist[u] + 1 : 0);
	}
	return dist[T];
}
int dfs (int T, int u, int now) {
	if (u == T || now == 0) return now;
	int cnt = 0;
	for (int& i = h[u], v, d; i; i = ne[i]) if (w[i] > 0 && dist[(v = e[i])] == dist[u] + 1 && (d = dfs (T, v, min (now - cnt, w[i])))) {
		cnt += d, w[i] -= d, w[i ^ 1] += d; if (cnt == now) return cnt;
	}
	return cnt;
}
const int inf = 2e9;
int calc (int N, int S, int T) {
	int cnt = 0; memcpy (cph, h, sizeof(int) * (N + 1));
	while (bfs (N, S, T)) memcpy (h, cph, sizeof(int) * (N + 1)), cnt += dfs (T, S, inf);
	return cnt;
}
}

namespace flow {
int h[maxn], e[maxn << 1], ne[maxn << 1], idx = 1;
int f[maxn << 1], w[maxn << 1];
void add (int a, int b, int c, int d) {
	idx++, e[idx] = b, ne[idx] = h[a], f[idx] = c, w[idx] = d, h[a] = idx;
	idx++, e[idx] = a, ne[idx] = h[b], f[idx] = 0, w[idx] = -d, h[b] = idx;
}
int dist[maxn], flow[maxn], pos[maxn]; bool st[maxn];
bool spfa (int N, int S, int T) {
	queue<int> q;
	memset (flow, 0, sizeof(int) * (N + 1));
	memset (dist, 0x3f, sizeof(int) * (N + 1));
	while (q.size()) q.pop();
	q.push (S); dist[S] = 0, flow[S] = 1e9;
	st[S] = 1;
	while (q.size()) {
		int u = q.front(); q.pop(); st[u] = 0;
		for (int i = h[u], v; i; i = ne[i]) {
			if (f[i] && dist[(v = e[i])] > dist[u] + w[i]) {
				dist[v] = dist[u] + w[i], pos[v] = i;
				flow[v] = min (flow[u], f[i]);
				if (!st[v]) q.push (v), st[v] = 1;
			}
		}
	}
	return flow[T] > 0;
}
pii EK (int N, int S, int T) {
	int cnt1 = 0, cnt2 = 0, tmp = 0;
	while (spfa(N, S, T)) {
		cnt1 += tmp = flow[T], cnt2 += tmp * dist[T];
		for (int i = T; i != S; i = e[pos[i] ^ 1])
			f[pos[i]] -= tmp, f[pos[i] ^ 1] += tmp;
	}
	return make_pair (cnt1, cnt2);
}
}

const double pi = acos(-1);
struct Complex {
  double x, y;
  Complex(double xx = 0, double yy = 0) : x(xx), y(yy) {}
};
Complex operator + (Complex a, Complex b) { return Complex(a.x + b.x, a.y + b.y); }
Complex operator - (Complex a, Complex b) { return Complex(a.x - b.x, a.y - b.y); }
Complex operator * (Complex a, Complex b) { 
  return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); 
}
void fft(vector<Complex>& A, int opt) {
  int n = A.size();
  vector<int> rev(n);
  for (int i = 0; i < n; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1) rev[i] |= n >> 1;
  }
  for (int i = 0; i < n; ++i) {
    if (i < rev[i]) swap(A[i], A[rev[i]]);
  }
  for (int mid = 1; mid < n; mid <<= 1) {
    Complex w1(cos(pi / mid), opt * sin(pi / mid));
    for (int j = 0; j < n; j += (mid << 1)) {
      Complex w(1, 0);
      for (int k = 0; k < mid; ++k, w = w * w1) {
        Complex x = A[j + k], y = w * A[j + k + mid];
        A[j + k] = x + y;
        A[j + k + mid] = x - y;
      }
    }
  }
  if (opt == -1) {
    for (int i = 0; i < n; ++i) {
      A[i].x /= n, A[i].y /= n;
    }
  }
}
vector<int> mul(vector<int> A, vector<int> B) {
  int sz = A.size() + B.size() - 1;
  int n = 1;
  while (n < A.size() + B.size()) n <<= 1;
  vector<Complex> tA(n), tB(n);
  for (int i = 0; i < A.size(); i++) tA[i].x = A[i];
  for (int i = 0; i < B.size(); i++) tB[i].x = B[i];
  fft(tA, 1); fft(tB, 1);
  for (int i = 0; i < n; i++) tA[i] = tA[i] * tB[i];
  fft(tA, -1);
  vector<int> C(sz);
  for (int i = 0; i < sz; i++)  C[i] = (int)(tA[i].x + 0.1);
  return C;
}

namespace ploy {
const int G = 3, invG = 332748118;
int qmi (int a, int b) {
	int res = 1;
	while (b) {
		if (b & 1) res = 1ll * res * a % mod;
		a = 1ll * a * a % mod; b >>= 1; 
	}
	return res;
}
void ntt(vector<int>& A, int opt) {
  int n = A.size();
  vector<int> rev(n);
  for (int i = 0; i < n; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1) rev[i] |= n >> 1;
  }
  for (int i = 0; i < n; ++i) {
    if (i < rev[i]) swap(A[i], A[rev[i]]);
  }
  for (int mid = 1; mid < n; mid <<= 1) {
  	int w1 = qmi (opt == 1 ? G : invG, (mod - 1) / (mid * 2));
    for (int j = 0; j < n; j += (mid << 1)) {
    	int w = 1;
      for (int k = 0; k < mid; ++k, w = 1ll * w * w1 % mod) {
        int x = A[j + k], y = 1ll * w * A[j + k + mid] % mod;
        A[j + k] = (x + y) % mod;
        A[j + k + mid] = (x - y + mod) % mod;
      }
    }
  }
  if (opt == -1) {
  	int x = qmi (n, mod - 2);
    for (int i = 0; i < n; ++i) {
      A[i] = 1ll * A[i] * x % mod;
    }
  }
}
vector<int> mul(vector<int> A, vector<int> B) {
  int sz = A.size() + B.size() - 1;
  int n = 1;
  while (n < (int)A.size() + (int)B.size()) n <<= 1;
  vector<int> tA(n), tB(n);
  for (int i = 0; i < (int)A.size(); i++) tA[i] = A[i];
  for (int i = 0; i < (int)B.size(); i++) tB[i] = B[i];
  ntt(tA, 1); ntt(tB, 1);
  for (int i = 0; i < n; i++) tA[i] = 1ll * tA[i] * tB[i] % mod;
  ntt(tA, -1);
  return vector<int>(tA.begin(), tA.begin() + sz);
}
vector<int> inv(vector<int>& A, int n) {
  if (n == 1) {
    return {qmi(A[0], mod - 2)};
  }
  vector<int> B = inv(A, (n + 1) / 2);
  int sz = 1;
  while (sz < 2 * n) sz <<= 1;
  vector<int> tA(sz), tB(sz);
  for (int i = 0; i < min((int)A.size(), n); i++) tA[i] = A[i];
  for (int i = 0; i < (int)B.size(); i++) tB[i] = B[i];
  ntt(tA, 1); ntt(tB, 1);
  for (int i = 0; i < sz; i++) {
    tA[i] = 1ll * tB[i] * (2 - 1ll * tA[i] * tB[i] % mod + mod) % mod;
  }
  ntt(tA, -1);
  return vector<int>(tA.begin(), tA.begin() + n);
}
vector<int> inv(vector<int>& A) {
  return inv(A, A.size());
}
const int inv2 = 499122177;
int sqrt_mod(int a) {
  if (a == 0) return 0;
  if (qmi(a, (mod - 1) / 2) != 1) return -1;
  if (mod % 4 == 3) return qmi(a, (mod + 1) / 4);
  int s = mod - 1;
  int e = 0;
  while (s % 2 == 0) s /= 2, e++;
  int n = 2;
  while (qmi(n, (mod - 1) / 2) != mod - 1) n++;
  int x = qmi(a, (s + 1) / 2), b = qmi(a, s), g = qmi(n, s), r = e;
  while (true) {
    int t = b;
    int m = 0;
    for (; m < r; m++) {
      if (t == 1) break;
      t = 1ll * t * t % mod;
    }
    if (m == 0) return x;
    int gs = qmi(g, 1 << (r - m - 1));
    g = 1ll * gs * gs % mod, x = 1ll * x * gs % mod;
    b = 1ll * b * g % mod, r = m;
  }
}
vector<int> sqrt(vector<int>& A, int n) {
	if (n == 1) {
    int root = sqrt_mod(A[0]);
    if (root == -1) return{};
    return {min(root, mod - root)};
  }
  vector<int> B = sqrt(A, (n + 1) / 2);
  if (B.empty()) return {};
  B.resize(n);
  vector<int> B_inv = inv(B);
  int sz = 1;
  while (sz < 2 * n) sz <<= 1;
  vector<int> tA(sz), tB(sz), tB_inv(sz);
  for (int i = 0; i < min((int)A.size(), n); i++) tA[i] = A[i];
  for (int i = 0; i < (int)B.size(); i++) tB[i] = B[i];
  for (int i = 0; i < (int)B_inv.size(); i++) tB_inv[i] = B_inv[i];
  ntt(tA, 1); ntt(tB, 1); ntt(tB_inv, 1);
  for (int i = 0; i < sz; i++) {
    tA[i] = 1ll * (1ll * tA[i] * tB_inv[i] % mod + tB[i]) % mod * inv2 % mod;
  }
  ntt(tA, -1);
  return vector<int>(tA.begin(), tA.begin() + n);
}
vector<int> sqrt(vector<int>& A) {
  return sqrt(A, A.size());
}
vector<int> derivative(vector<int>& A) {
  int n = A.size();
  if (n == 1) return {0};
  vector<int> res(n - 1);
  for (int i = 1; i < n; i++) {
    res[i - 1] = 1ll * A[i] * i % mod;
  }
  return res;
}
vector<int> init(vector<int>& A) {
  int n = A.size();
  vector<int> res(n + 1);
  for (int i = 0; i < n; i++) {
    res[i + 1] = 1ll * A[i] * qmi(i + 1, mod - 2) % mod;
  }
  return res;
}
vector<int> ln(vector<int>& A) {
  if (A.empty() || A[0] != 1) return {};
  vector<int> A_deriv = derivative(A);
  vector<int> A_inv = inv(A);
  vector<int> p = mul(A_deriv, A_inv);
  p = init(p);
  if (p.size() > A.size()) {
    p.resize(A.size());
  } else if (p.size() < A.size()) {
    while (p.size() < A.size()) p.push_back(0);
  }
  return p;
}
vector<int> exp(vector<int>& A, int n) {
  if (n == 1) return {1};
  int m = (n + 1) / 2;
  vector<int> B = exp(A, m);
  vector<int> ln_B = B;
  ln_B.resize(n);
  ln_B = ln(ln_B);
  vector<int> d(n);
  for (int i = 0; i < n; ++i) {
    int a_val = (i < (int)A.size()) ? A[i] : 0;
    int ln_val = (i < (int)ln_B.size()) ? ln_B[i] : 0;
    d[i] = (a_val - ln_val + mod) % mod;
  }
  d[0] = (d[0] + 1) % mod;
  vector<int> p = mul(B, d);
  p.resize(n);
  return p;
}
vector<int> exp(vector<int>& A) {
  if (A.empty() || A[0] != 0)  return {};
  return exp(A, A.size());
}
vector<int> pow(vector<int> A, int k) {
  if (A.empty()) return {};
  int d = 0;
  while (d < (int)A.size() && A[d] == 0) ++d;
  if (d == (int)A.size()) {
    if (k == 0) return {1};
    else return vector<int>(A.size(), 0);
  }
  if (d * k >= (int)A.size()) return vector<int>(A.size(), 0);
  int coef = A[d], inv_coef = qmi(coef, mod - 2);
  vector<int> B(A.size() - d);
  for (int i = 0; i < (int)B.size(); ++i) B[i] = 1ll * A[i + d] * inv_coef % mod;
  vector<int> ln_B = ln(B);
  for (int i = 0; i < (int)ln_B.size(); ++i) ln_B[i] = 1ll * ln_B[i] * k % mod;
  vector<int> exp_ln = exp(ln_B);
  int pow_coef = qmi(coef, k);
  vector<int> res(A.size(), 0);
  for (int i = 0; i + d * k < (int)res.size() && i < (int)exp_ln.size(); ++i)
  	res[i + d * k] = 1ll * exp_ln[i] * pow_coef % mod;
  return res;
}
vector<int> pow(vector<int> A, string k_str) {
  if (A.empty()) return {};
  int d = 0;
  while (d < (int)A.size() && A[d] == 0) ++d;
  if (d == (int)A.size()) {
    if (k_str == "0") return {1};
    else return vector<int>(A.size(), 0);
  }
  i64 k_mod = 0, k_mod_phi = 0;
  for (char c : k_str) {
    k_mod = (k_mod * 10 + (c - '0')) % mod;
    k_mod_phi = (k_mod_phi * 10 + (c - '0')) % (mod - 1);
  }
  if ((i64)d * k_mod >= (i64)A.size()) return vector<int>(A.size(), 0);
  int coef = A[d];
  int inv_coef = qmi(coef, mod - 2);
  vector<int> B(A.size() - d);
  for (int i = 0; i < (int)B.size(); ++i) B[i] = 1ll * A[i + d] * inv_coef % mod;
  vector<int> ln_B = ln(B);
  for (int i = 0; i < (int)ln_B.size(); ++i) ln_B[i] = 1ll * ln_B[i] * k_mod % mod;
  vector<int> exp_ln = exp(ln_B);
  int pow_coef = qmi(coef, k_mod_phi);
  vector<int> res(A.size(), 0);
  for (int i = 0; i + d * k_mod < (int)res.size() && i < (int)exp_ln.size(); ++i)
    res[i + d * k_mod] = 1ll * exp_ln[i] * pow_coef % mod;
  return res;
}
} // using namespace ploy;