QOJ.ac

QOJ

IDProblemSubmitterResultTimeMemoryLanguageFile sizeSubmit timeJudge time
#903420#10095. Producthos.lyric🐇AC ✓2680ms170744kbC++1442.4kb2025-02-17 12:04:342025-02-17 12:04:35

Judging History

This is the latest submission verdict.

  • [2025-02-17 12:04:35]
  • Judged
  • Verdict: AC
  • Time: 2680ms
  • Memory: 170744kb
  • [2025-02-17 12:04:34]
  • Submitted

answer

#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <bitset>
#include <complex>
#include <deque>
#include <functional>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

using namespace std;

using Int = long long;

template <class T1, class T2> ostream &operator<<(ostream &os, const pair<T1, T2> &a) { return os << "(" << a.first << ", " << a.second << ")"; };
template <class T> ostream &operator<<(ostream &os, const vector<T> &as) { const int sz = as.size(); os << "["; for (int i = 0; i < sz; ++i) { if (i >= 256) { os << ", ..."; break; } if (i > 0) { os << ", "; } os << as[i]; } return os << "]"; }
template <class T> void pv(T a, T b) { for (T i = a; i != b; ++i) cerr << *i << " "; cerr << endl; }
template <class T> bool chmin(T &t, const T &f) { if (t > f) { t = f; return true; } return false; }
template <class T> bool chmax(T &t, const T &f) { if (t < f) { t = f; return true; } return false; }
#define COLOR(s) ("\x1b[" s "m")

////////////////////////////////////////////////////////////////////////////////
template <unsigned M_> struct ModInt {
  static constexpr unsigned M = M_;
  unsigned x;
  constexpr ModInt() : x(0U) {}
  constexpr ModInt(unsigned x_) : x(x_ % M) {}
  constexpr ModInt(unsigned long long x_) : x(x_ % M) {}
  constexpr ModInt(int x_) : x(((x_ %= static_cast<int>(M)) < 0) ? (x_ + static_cast<int>(M)) : x_) {}
  constexpr ModInt(long long x_) : x(((x_ %= static_cast<long long>(M)) < 0) ? (x_ + static_cast<long long>(M)) : x_) {}
  ModInt &operator+=(const ModInt &a) { x = ((x += a.x) >= M) ? (x - M) : x; return *this; }
  ModInt &operator-=(const ModInt &a) { x = ((x -= a.x) >= M) ? (x + M) : x; return *this; }
  ModInt &operator*=(const ModInt &a) { x = (static_cast<unsigned long long>(x) * a.x) % M; return *this; }
  ModInt &operator/=(const ModInt &a) { return (*this *= a.inv()); }
  ModInt pow(long long e) const {
    if (e < 0) return inv().pow(-e);
    ModInt a = *this, b = 1U; for (; e; e >>= 1) { if (e & 1) b *= a; a *= a; } return b;
  }
  ModInt inv() const {
    unsigned a = M, b = x; int y = 0, z = 1;
    for (; b; ) { const unsigned q = a / b; const unsigned c = a - q * b; a = b; b = c; const int w = y - static_cast<int>(q) * z; y = z; z = w; }
    assert(a == 1U); return ModInt(y);
  }
  ModInt operator+() const { return *this; }
  ModInt operator-() const { ModInt a; a.x = x ? (M - x) : 0U; return a; }
  ModInt operator+(const ModInt &a) const { return (ModInt(*this) += a); }
  ModInt operator-(const ModInt &a) const { return (ModInt(*this) -= a); }
  ModInt operator*(const ModInt &a) const { return (ModInt(*this) *= a); }
  ModInt operator/(const ModInt &a) const { return (ModInt(*this) /= a); }
  template <class T> friend ModInt operator+(T a, const ModInt &b) { return (ModInt(a) += b); }
  template <class T> friend ModInt operator-(T a, const ModInt &b) { return (ModInt(a) -= b); }
  template <class T> friend ModInt operator*(T a, const ModInt &b) { return (ModInt(a) *= b); }
  template <class T> friend ModInt operator/(T a, const ModInt &b) { return (ModInt(a) /= b); }
  explicit operator bool() const { return x; }
  bool operator==(const ModInt &a) const { return (x == a.x); }
  bool operator!=(const ModInt &a) const { return (x != a.x); }
  friend std::ostream &operator<<(std::ostream &os, const ModInt &a) { return os << a.x; }
};
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
constexpr unsigned MO = 998244353U;
constexpr unsigned MO2 = 2U * MO;
constexpr int FFT_MAX = 23;
using Mint = ModInt<MO>;
constexpr Mint FFT_ROOTS[FFT_MAX + 1] = {1U, 998244352U, 911660635U, 372528824U, 929031873U, 452798380U, 922799308U, 781712469U, 476477967U, 166035806U, 258648936U, 584193783U, 63912897U, 350007156U, 666702199U, 968855178U, 629671588U, 24514907U, 996173970U, 363395222U, 565042129U, 733596141U, 267099868U, 15311432U};
constexpr Mint INV_FFT_ROOTS[FFT_MAX + 1] = {1U, 998244352U, 86583718U, 509520358U, 337190230U, 87557064U, 609441965U, 135236158U, 304459705U, 685443576U, 381598368U, 335559352U, 129292727U, 358024708U, 814576206U, 708402881U, 283043518U, 3707709U, 121392023U, 704923114U, 950391366U, 428961804U, 382752275U, 469870224U};
constexpr Mint FFT_RATIOS[FFT_MAX] = {911660635U, 509520358U, 369330050U, 332049552U, 983190778U, 123842337U, 238493703U, 975955924U, 603855026U, 856644456U, 131300601U, 842657263U, 730768835U, 942482514U, 806263778U, 151565301U, 510815449U, 503497456U, 743006876U, 741047443U, 56250497U, 867605899U};
constexpr Mint INV_FFT_RATIOS[FFT_MAX] = {86583718U, 372528824U, 373294451U, 645684063U, 112220581U, 692852209U, 155456985U, 797128860U, 90816748U, 860285882U, 927414960U, 354738543U, 109331171U, 293255632U, 535113200U, 308540755U, 121186627U, 608385704U, 438932459U, 359477183U, 824071951U, 103369235U};

// as[rev(i)] <- \sum_j \zeta^(ij) as[j]
void fft(Mint *as, int n) {
  assert(!(n & (n - 1))); assert(1 <= n); assert(n <= 1 << FFT_MAX);
  int m = n;
  if (m >>= 1) {
    for (int i = 0; i < m; ++i) {
      const unsigned x = as[i + m].x;  // < MO
      as[i + m].x = as[i].x + MO - x;  // < 2 MO
      as[i].x += x;  // < 2 MO
    }
  }
  if (m >>= 1) {
    Mint prod = 1U;
    for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
      for (int i = i0; i < i0 + m; ++i) {
        const unsigned x = (prod * as[i + m]).x;  // < MO
        as[i + m].x = as[i].x + MO - x;  // < 3 MO
        as[i].x += x;  // < 3 MO
      }
      prod *= FFT_RATIOS[__builtin_ctz(++h)];
    }
  }
  for (; m; ) {
    if (m >>= 1) {
      Mint prod = 1U;
      for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
        for (int i = i0; i < i0 + m; ++i) {
          const unsigned x = (prod * as[i + m]).x;  // < MO
          as[i + m].x = as[i].x + MO - x;  // < 4 MO
          as[i].x += x;  // < 4 MO
        }
        prod *= FFT_RATIOS[__builtin_ctz(++h)];
      }
    }
    if (m >>= 1) {
      Mint prod = 1U;
      for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
        for (int i = i0; i < i0 + m; ++i) {
          const unsigned x = (prod * as[i + m]).x;  // < MO
          as[i].x = (as[i].x >= MO2) ? (as[i].x - MO2) : as[i].x;  // < 2 MO
          as[i + m].x = as[i].x + MO - x;  // < 3 MO
          as[i].x += x;  // < 3 MO
        }
        prod *= FFT_RATIOS[__builtin_ctz(++h)];
      }
    }
  }
  for (int i = 0; i < n; ++i) {
    as[i].x = (as[i].x >= MO2) ? (as[i].x - MO2) : as[i].x;  // < 2 MO
    as[i].x = (as[i].x >= MO) ? (as[i].x - MO) : as[i].x;  // < MO
  }
}

// as[i] <- (1/n) \sum_j \zeta^(-ij) as[rev(j)]
void invFft(Mint *as, int n) {
  assert(!(n & (n - 1))); assert(1 <= n); assert(n <= 1 << FFT_MAX);
  int m = 1;
  if (m < n >> 1) {
    Mint prod = 1U;
    for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
      for (int i = i0; i < i0 + m; ++i) {
        const unsigned long long y = as[i].x + MO - as[i + m].x;  // < 2 MO
        as[i].x += as[i + m].x;  // < 2 MO
        as[i + m].x = (prod.x * y) % MO;  // < MO
      }
      prod *= INV_FFT_RATIOS[__builtin_ctz(++h)];
    }
    m <<= 1;
  }
  for (; m < n >> 1; m <<= 1) {
    Mint prod = 1U;
    for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
      for (int i = i0; i < i0 + (m >> 1); ++i) {
        const unsigned long long y = as[i].x + MO2 - as[i + m].x;  // < 4 MO
        as[i].x += as[i + m].x;  // < 4 MO
        as[i].x = (as[i].x >= MO2) ? (as[i].x - MO2) : as[i].x;  // < 2 MO
        as[i + m].x = (prod.x * y) % MO;  // < MO
      }
      for (int i = i0 + (m >> 1); i < i0 + m; ++i) {
        const unsigned long long y = as[i].x + MO - as[i + m].x;  // < 2 MO
        as[i].x += as[i + m].x;  // < 2 MO
        as[i + m].x = (prod.x * y) % MO;  // < MO
      }
      prod *= INV_FFT_RATIOS[__builtin_ctz(++h)];
    }
  }
  if (m < n) {
    for (int i = 0; i < m; ++i) {
      const unsigned y = as[i].x + MO2 - as[i + m].x;  // < 4 MO
      as[i].x += as[i + m].x;  // < 4 MO
      as[i + m].x = y;  // < 4 MO
    }
  }
  const Mint invN = Mint(n).inv();
  for (int i = 0; i < n; ++i) {
    as[i] *= invN;
  }
}

void fft(vector<Mint> &as) {
  fft(as.data(), as.size());
}
void invFft(vector<Mint> &as) {
  invFft(as.data(), as.size());
}

vector<Mint> convolve(vector<Mint> as, vector<Mint> bs) {
  if (as.empty() || bs.empty()) return {};
  const int len = as.size() + bs.size() - 1;
  int n = 1;
  for (; n < len; n <<= 1) {}
  as.resize(n); fft(as);
  bs.resize(n); fft(bs);
  for (int i = 0; i < n; ++i) as[i] *= bs[i];
  invFft(as);
  as.resize(len);
  return as;
}
vector<Mint> square(vector<Mint> as) {
  if (as.empty()) return {};
  const int len = as.size() + as.size() - 1;
  int n = 1;
  for (; n < len; n <<= 1) {}
  as.resize(n); fft(as);
  for (int i = 0; i < n; ++i) as[i] *= as[i];
  invFft(as);
  as.resize(len);
  return as;
}
// m := |as|, n := |bs|
// cs[k] = \sum[i-j=k] as[i] bs[j]  (0 <= k <= m-n)
// transpose of ((multiply by bs): K^[0,m-n] -> K^[0,m-1])
vector<Mint> middle(vector<Mint> as, vector<Mint> bs) {
  const int m = as.size(), n = bs.size();
  assert(m >= n); assert(n >= 1);
  int len = 1;
  for (; len < m; len <<= 1) {}
  as.resize(len, 0);
  fft(as);
  std::reverse(bs.begin(), bs.end());
  bs.resize(len, 0);
  fft(bs);
  for (int i = 0; i < len; ++i) as[i] *= bs[i];
  invFft(as);
  as.resize(m);
  as.erase(as.begin(), as.begin() + (n - 1));
  return as;
}
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
// inv: log, exp, pow
// fac: shift
// invFac: shift
constexpr int LIM_INV = 1 << 20;  // @
Mint inv[LIM_INV], fac[LIM_INV], invFac[LIM_INV];
struct ModIntPreparator {
  ModIntPreparator() {
    inv[1] = 1;
    for (int i = 2; i < LIM_INV; ++i) inv[i] = -((Mint::M / i) * inv[Mint::M % i]);
    fac[0] = 1;
    for (int i = 1; i < LIM_INV; ++i) fac[i] = fac[i - 1] * i;
    invFac[0] = 1;
    for (int i = 1; i < LIM_INV; ++i) invFac[i] = invFac[i - 1] * inv[i];
  }
} preparator;

// polyWork0: *, inv, div, divAt, log, exp, pow, sqrt, shift
// polyWork1: inv, div, divAt, log, exp, pow, sqrt, shift
// polyWork2: divAt, exp, pow, sqrt
// polyWork3: exp, pow, sqrt
static constexpr int LIM_POLY = 1 << 20;  // @
static_assert(LIM_POLY <= 1 << FFT_MAX, "Poly: LIM_POLY <= 1 << FFT_MAX must hold.");
static Mint polyWork0[LIM_POLY], polyWork1[LIM_POLY], polyWork2[LIM_POLY], polyWork3[LIM_POLY];

struct Poly : public vector<Mint> {
  Poly() {}
  explicit Poly(int n) : vector<Mint>(n) {}
  Poly(const vector<Mint> &vec) : vector<Mint>(vec) {}
  Poly(std::initializer_list<Mint> il) : vector<Mint>(il) {}
  int size() const { return vector<Mint>::size(); }
  Mint at(long long k) const { return (0 <= k && k < size()) ? (*this)[k] : 0U; }
  int ord() const { for (int i = 0; i < size(); ++i) if ((*this)[i]) return i; return -1; }
  int deg() const { for (int i = size(); --i >= 0; ) if ((*this)[i]) return i; return -1; }
  Poly mod(int n) const { return Poly(vector<Mint>(data(), data() + min(n, size()))); }
  friend std::ostream &operator<<(std::ostream &os, const Poly &fs) {
    os << "[";
    for (int i = 0; i < fs.size(); ++i) { if (i > 0) os << ", "; os << fs[i]; }
    return os << "]";
  }

  Poly &operator+=(const Poly &fs) {
    if (size() < fs.size()) resize(fs.size());
    for (int i = 0; i < fs.size(); ++i) (*this)[i] += fs[i];
    return *this;
  }
  Poly &operator-=(const Poly &fs) {
    if (size() < fs.size()) resize(fs.size());
    for (int i = 0; i < fs.size(); ++i) (*this)[i] -= fs[i];
    return *this;
  }
  // 3 E(|t| + |f|)
  Poly &operator*=(const Poly &fs) {
    if (empty() || fs.empty()) return *this = {};
    const int nt = size(), nf = fs.size();
    int n = 1;
    for (; n < nt + nf - 1; n <<= 1) {}
    assert(n <= LIM_POLY);
    resize(n);
    fft(data(), n);  // 1 E(n)
    memcpy(polyWork0, fs.data(), nf * sizeof(Mint));
    memset(polyWork0 + nf, 0, (n - nf) * sizeof(Mint));
    fft(polyWork0, n);  // 1 E(n)
    for (int i = 0; i < n; ++i) (*this)[i] *= polyWork0[i];
    invFft(data(), n);  // 1 E(n)
    resize(nt + nf - 1);
    return *this;
  }
  // 13 E(deg(t) - deg(f) + 1)
  // rev(t) = rev(f) rev(q) + x^(deg(t)-deg(f)+1) rev(r)
  Poly &operator/=(const Poly &fs) {
    const int m = deg(), n = fs.deg();
    assert(n != -1);
    if (m < n) return *this = {};
    Poly tsRev(m - n + 1), fsRev(min(m - n, n) + 1);
    for (int i = 0; i <= m - n; ++i) tsRev[i] = (*this)[m - i];
    for (int i = 0, i0 = min(m - n, n); i <= i0; ++i) fsRev[i] = fs[n - i];
    const Poly qsRev = tsRev.div(fsRev, m - n + 1);  // 13 E(m - n + 1)
    resize(m - n + 1);
    for (int i = 0; i <= m - n; ++i) (*this)[i] = qsRev[m - n - i];
    return *this;
  }
  // 13 E(deg(t) - deg(f) + 1) + 3 E(|t|)
  Poly &operator%=(const Poly &fs) {
    const Poly qs = *this / fs;  // 13 E(deg(t) - deg(f) + 1)
    *this -= fs * qs;  // 3 E(|t|)
    resize(deg() + 1);
    return *this;
  }
  Poly &operator*=(const Mint &a) {
    for (int i = 0; i < size(); ++i) (*this)[i] *= a;
    return *this;
  }
  Poly &operator/=(const Mint &a) {
    const Mint b = a.inv();
    for (int i = 0; i < size(); ++i) (*this)[i] *= b;
    return *this;
  }
  Poly operator+() const { return *this; }
  Poly operator-() const {
    Poly fs(size());
    for (int i = 0; i < size(); ++i) fs[i] = -(*this)[i];
    return fs;
  }
  Poly operator+(const Poly &fs) const { return (Poly(*this) += fs); }
  Poly operator-(const Poly &fs) const { return (Poly(*this) -= fs); }
  Poly operator*(const Poly &fs) const { return (Poly(*this) *= fs); }
  Poly operator/(const Poly &fs) const { return (Poly(*this) /= fs); }
  Poly operator%(const Poly &fs) const { return (Poly(*this) %= fs); }
  Poly operator*(const Mint &a) const { return (Poly(*this) *= a); }
  Poly operator/(const Mint &a) const { return (Poly(*this) /= a); }
  friend Poly operator*(const Mint &a, const Poly &fs) { return fs * a; }

  // 10 E(n)
  // f <- f - (t f - 1) f
  Poly inv(int n) const {
    assert(!empty()); assert((*this)[0]); assert(1 <= n);
    assert(n == 1 || 1 << (32 - __builtin_clz(n - 1)) <= LIM_POLY);
    Poly fs(n);
    fs[0] = (*this)[0].inv();
    for (int m = 1; m < n; m <<= 1) {
      memcpy(polyWork0, data(), min(m << 1, size()) * sizeof(Mint));
      memset(polyWork0 + min(m << 1, size()), 0, ((m << 1) - min(m << 1, size())) * sizeof(Mint));
      fft(polyWork0, m << 1);  // 2 E(n)
      memcpy(polyWork1, fs.data(), min(m << 1, n) * sizeof(Mint));
      memset(polyWork1 + min(m << 1, n), 0, ((m << 1) - min(m << 1, n)) * sizeof(Mint));
      fft(polyWork1, m << 1);  // 2 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
      invFft(polyWork0, m << 1); // 2 E(n)
      memset(polyWork0, 0, m * sizeof(Mint));
      fft(polyWork0, m << 1); // 2 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
      invFft(polyWork0, m << 1); // 2 E(n)
      for (int i = m, i0 = min(m << 1, n); i < i0; ++i) fs[i] = -polyWork0[i];
    }
    return fs;
  }
  // 9 E(n)
  // Need (4 m)-th roots of unity to lift from (mod x^m) to (mod x^(2m)).
  // f <- f - (t f - 1) f
  // (t f^2) mod ((x^(2m) - 1) (x^m - 1^(1/4)))
  /*
  Poly inv(int n) const {
    assert(!empty()); assert((*this)[0]); assert(1 <= n);
    assert(n == 1 || 3 << (31 - __builtin_clz(n - 1)) <= LIM_POLY);
    assert(n <= 1 << (FFT_MAX - 1));
    Poly fs(n);
    fs[0] = (*this)[0].inv();
    for (int h = 2, m = 1; m < n; ++h, m <<= 1) {
      const Mint a = FFT_ROOTS[h], b = INV_FFT_ROOTS[h];
      memcpy(polyWork0, data(), min(m << 1, size()) * sizeof(Mint));
      memset(polyWork0 + min(m << 1, size()), 0, ((m << 1) - min(m << 1, size())) * sizeof(Mint));
      {
        Mint aa = 1;
        for (int i = 0; i < m; ++i) { polyWork0[(m << 1) + i] = aa * polyWork0[i]; aa *= a; }
        for (int i = 0; i < m; ++i) { polyWork0[(m << 1) + i] += aa * polyWork0[m + i]; aa *= a; }
      }
      fft(polyWork0, m << 1);  // 2 E(n)
      fft(polyWork0 + (m << 1), m);  // 1 E(n)
      memcpy(polyWork1, fs.data(), min(m << 1, n) * sizeof(Mint));
      memset(polyWork1 + min(m << 1, n), 0, ((m << 1) - min(m << 1, n)) * sizeof(Mint));
      {
        Mint aa = 1;
        for (int i = 0; i < m; ++i) { polyWork1[(m << 1) + i] = aa * polyWork1[i]; aa *= a; }
        for (int i = 0; i < m; ++i) { polyWork1[(m << 1) + i] += aa * polyWork1[m + i]; aa *= a; }
      }
      fft(polyWork1, m << 1);  // 2 E(n)
      fft(polyWork1 + (m << 1), m);  // 1 E(n)
      for (int i = 0; i < (m << 1) + m; ++i) polyWork0[i] *= polyWork1[i] * polyWork1[i];
      invFft(polyWork0, m << 1);  // 2 E(n)
      invFft(polyWork0 + (m << 1), m);  // 1 E(n)
      // 2 f0 + (-f2), (-f1) + (-f3), 1^(1/4) (-f1) - (-f2) - 1^(1/4) (-f3)
      {
        Mint bb = 1;
        for (int i = 0, i0 = min(m, n - m); i < i0; ++i) {
          unsigned x = polyWork0[i].x + (bb * polyWork0[(m << 1) + i]).x + MO2 - (fs[i].x << 1);  // < 4 MO
          fs[m + i] = Mint(static_cast<unsigned long long>(FFT_ROOTS[2].x) * x) - polyWork0[m + i];
          fs[m + i].x = ((fs[m + i].x & 1) ? (fs[m + i].x + MO) : fs[m + i].x) >> 1;
          bb *= b;
        }
      }
    }
    return fs;
  }
  */
  // 13 E(n)
  // g = (1 / f) mod x^m
  // h <- h - (f h - t) g
  Poly div(const Poly &fs, int n) const {
    assert(!fs.empty()); assert(fs[0]); assert(1 <= n);
    if (n == 1) return {at(0) / fs[0]};
    // m < n <= 2 m
    const int m = 1 << (31 - __builtin_clz(n - 1));
    assert(m << 1 <= LIM_POLY);
    Poly gs = fs.inv(m);  // 5 E(n)
    gs.resize(m << 1);
    fft(gs.data(), m << 1);  // 1 E(n)
    if (size()) memcpy(polyWork0, data(), min(m, size()) * sizeof(Mint));
    memset(polyWork0 + min(m, size()), 0, ((m << 1) - min(m, size())) * sizeof(Mint));
    fft(polyWork0, m << 1);  // 1 E(n)
    for (int i = 0; i < m << 1; ++i) polyWork0[i] *= gs[i];
    invFft(polyWork0, m << 1);  // 1 E(n)
    Poly hs(n);
    memcpy(hs.data(), polyWork0, m * sizeof(Mint));
    memset(polyWork0 + m, 0, m * sizeof(Mint));
    fft(polyWork0, m << 1);  // 1 E(n)
    memcpy(polyWork1, fs.data(), min(m << 1, fs.size()) * sizeof(Mint));
    memset(polyWork1 + min(m << 1, fs.size()), 0, ((m << 1) - min(m << 1, fs.size())) * sizeof(Mint));
    fft(polyWork1, m << 1);  // 1 E(n)
    for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
    invFft(polyWork0, m << 1);  // 1 E(n)
    memset(polyWork0, 0, m * sizeof(Mint));
    for (int i = m, i0 = min(m << 1, size()); i < i0; ++i) polyWork0[i] -= (*this)[i];
    fft(polyWork0, m << 1);  // 1 E(n)
    for (int i = 0; i < m << 1; ++i) polyWork0[i] *= gs[i];
    invFft(polyWork0, m << 1);  // 1 E(n)
    for (int i = m; i < n; ++i) hs[i] = -polyWork0[i];
    return hs;
  }
  // (4 (floor(log_2 k) - ceil(log_2 |f|)) + 16) E(|f|)  for  |t| < |f|
  // [x^k] (t(x) / f(x)) = [x^k] ((t(x) f(-x)) / (f(x) f(-x))
  // polyWork0: half of (2 m)-th roots of unity, inversed, bit-reversed
  Mint divAt(const Poly &fs, long long k) const {
    assert(k >= 0);
    if (size() >= fs.size()) {
      const Poly qs = *this / fs;  // 13 E(deg(t) - deg(f) + 1)
      Poly rs = *this - fs * qs;  // 3 E(|t|)
      rs.resize(rs.deg() + 1);
      return qs.at(k) + rs.divAt(fs, k);
    }
    int h = 0, m = 1;
    for (; m < fs.size(); ++h, m <<= 1) {}
    if (k < m) {
      const Poly gs = fs.inv(k + 1);  // 10 E(|f|)
      Mint sum;
      for (int i = 0, i0 = min<int>(k + 1, size()); i < i0; ++i) sum += (*this)[i] * gs[k - i];
      return sum;
    }
    assert(m << 1 <= LIM_POLY);
    polyWork0[0] = Mint(2U).inv();
    for (int hh = 0; hh < h; ++hh) for (int i = 0; i < 1 << hh; ++i) polyWork0[1 << hh | i] = polyWork0[i] * INV_FFT_ROOTS[hh + 2];
    const Mint a = FFT_ROOTS[h + 1];
    if (size()) memcpy(polyWork2, data(), size() * sizeof(Mint));
    memset(polyWork2 + size(), 0, ((m << 1) - size()) * sizeof(Mint));
    fft(polyWork2, m << 1);  // 2 E(|f|)
    memcpy(polyWork1, fs.data(), fs.size() * sizeof(Mint));
    memset(polyWork1 + fs.size(), 0, ((m << 1) - fs.size()) * sizeof(Mint));
    fft(polyWork1, m << 1);  // 2 E(|f|)
    for (; ; ) {
      if (k & 1) {
        for (int i = 0; i < m; ++i) polyWork2[i] = polyWork0[i] * (polyWork2[i << 1 | 0] * polyWork1[i << 1 | 1] - polyWork2[i << 1 | 1] * polyWork1[i << 1 | 0]);
      } else {
        for (int i = 0; i < m; ++i) {
          polyWork2[i] = polyWork2[i << 1 | 0] * polyWork1[i << 1 | 1] + polyWork2[i << 1 | 1] * polyWork1[i << 1 | 0];
          polyWork2[i].x = ((polyWork2[i].x & 1) ? (polyWork2[i].x + MO) : polyWork2[i].x) >> 1;
        }
      }
      for (int i = 0; i < m; ++i) polyWork1[i] = polyWork1[i << 1 | 0] * polyWork1[i << 1 | 1];
      if ((k >>= 1) < m) {
        invFft(polyWork2, m);  // 1 E(|f|)
        invFft(polyWork1, m);  // 1 E(|f|)
        // Poly::inv does not use polyWork2
        const Poly gs = Poly(vector<Mint>(polyWork1, polyWork1 + k + 1)).inv(k + 1);  // 10 E(|f|)
        Mint sum;
        for (int i = 0; i <= k; ++i) sum += polyWork2[i] * gs[k - i];
        return sum;
      }
      memcpy(polyWork2 + m, polyWork2, m * sizeof(Mint));
      invFft(polyWork2 + m, m);  // (floor(log_2 k) - ceil(log_2 |f|)) E(|f|)
      memcpy(polyWork1 + m, polyWork1, m * sizeof(Mint));
      invFft(polyWork1 + m, m);  // (floor(log_2 k) - ceil(log_2 |f|)) E(|f|)
      Mint aa = 1;
      for (int i = m; i < m << 1; ++i) { polyWork2[i] *= aa; polyWork1[i] *= aa; aa *= a; }
      fft(polyWork2 + m, m);  // (floor(log_2 k) - ceil(log_2 |f|)) E(|f|)
      fft(polyWork1 + m, m);  // (floor(log_2 k) - ceil(log_2 |f|)) E(|f|)
    }
  }
  // 13 E(n)
  // D log(t) = (D t) / t
  Poly log(int n) const {
    assert(!empty()); assert((*this)[0].x == 1U); assert(n <= LIM_INV);
    Poly fs = mod(n);
    for (int i = 0; i < fs.size(); ++i) fs[i] *= i;
    fs = fs.div(*this, n);
    for (int i = 1; i < n; ++i) fs[i] *= ::inv[i];
    return fs;
  }
  // (16 + 1/2) E(n)
  // f = exp(t) mod x^m  ==>  (D f) / f == D t  (mod x^m)
  // g = (1 / exp(t)) mod x^m
  // f <- f - (log f - t) / (1 / f)
  //   =  f - (I ((D f) / f) - t) f
  //   == f - (I ((D f) / f + (f g - 1) ((D f) / f - D (t mod x^m))) - t) f  (mod x^(2m))
  //   =  f - (I (g (D f - f D (t mod x^m)) + D (t mod x^m)) - t) f
  // g <- g - (f g - 1) g
  // polyWork1: DFT(f, 2 m), polyWork2: g, polyWork3: DFT(g, 2 m)
  Poly exp(int n) const {
    assert(!empty()); assert(!(*this)[0]); assert(1 <= n);
    assert(n == 1 || 1 << (32 - __builtin_clz(n - 1)) <= min(LIM_INV, LIM_POLY));
    if (n == 1) return {1U};
    if (n == 2) return {1U, at(1)};
    Poly fs(n);
    fs[0].x = polyWork1[0].x = polyWork1[1].x = polyWork2[0].x = 1U;
    int m;
    for (m = 1; m << 1 < n; m <<= 1) {
      for (int i = 0, i0 = min(m, size()); i < i0; ++i) polyWork0[i] = i * (*this)[i];
      memset(polyWork0 + min(m, size()), 0, (m - min(m, size())) * sizeof(Mint));
      fft(polyWork0, m);  // (1/2) E(n)
      for (int i = 0; i < m; ++i) polyWork0[i] *= polyWork1[i];
      invFft(polyWork0, m);  // (1/2) E(n)
      for (int i = 0; i < m; ++i) polyWork0[i] -= i * fs[i];
      memset(polyWork0 + m, 0, m * sizeof(Mint));
      fft(polyWork0, m << 1);  // 1 E(n)
      memcpy(polyWork3, polyWork2, m * sizeof(Mint));
      memset(polyWork3 + m, 0, m * sizeof(Mint));
      fft(polyWork3, m << 1);  // 1 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork3[i];
      invFft(polyWork0, m << 1);  // 1 E(n)
      for (int i = 0; i < m; ++i) polyWork0[i] *= ::inv[m + i];
      for (int i = 0, i0 = min(m, size() - m); i < i0; ++i) polyWork0[i] += (*this)[m + i];
      memset(polyWork0 + m, 0, m * sizeof(Mint));
      fft(polyWork0, m << 1);  // 1 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
      invFft(polyWork0, m << 1);  // 1 E(n)
      memcpy(fs.data() + m, polyWork0, m * sizeof(Mint));
      memcpy(polyWork1, fs.data(), (m << 1) * sizeof(Mint));
      memset(polyWork1 + (m << 1), 0, (m << 1) * sizeof(Mint));
      fft(polyWork1, m << 2);  // 2 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork0[i] = polyWork1[i] * polyWork3[i];
      invFft(polyWork0, m << 1);  // 1 E(n)
      memset(polyWork0, 0, m * sizeof(Mint));
      fft(polyWork0, m << 1);  // 1 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork3[i];
      invFft(polyWork0, m << 1);  // 1 E(n)
      for (int i = m; i < m << 1; ++i) polyWork2[i] = -polyWork0[i];
    }
    for (int i = 0, i0 = min(m, size()); i < i0; ++i) polyWork0[i] = i * (*this)[i];
    memset(polyWork0 + min(m, size()), 0, (m - min(m, size())) * sizeof(Mint));
    fft(polyWork0, m);  // (1/2) E(n)
    for (int i = 0; i < m; ++i) polyWork0[i] *= polyWork1[i];
    invFft(polyWork0, m);  // (1/2) E(n)
    for (int i = 0; i < m; ++i) polyWork0[i] -= i * fs[i];
    memcpy(polyWork0 + m, polyWork0 + (m >> 1), (m >> 1) * sizeof(Mint));
    memset(polyWork0 + (m >> 1), 0, (m >> 1) * sizeof(Mint));
    memset(polyWork0 + m + (m >> 1), 0, (m >> 1) * sizeof(Mint));
    fft(polyWork0, m);  // (1/2) E(n)
    fft(polyWork0 + m, m);  // (1/2) E(n)
    memcpy(polyWork3 + m, polyWork2 + (m >> 1), (m >> 1) * sizeof(Mint));
    memset(polyWork3 + m + (m >> 1), 0, (m >> 1) * sizeof(Mint));
    fft(polyWork3 + m, m);  // (1/2) E(n)
    for (int i = 0; i < m; ++i) polyWork0[m + i] = polyWork0[i] * polyWork3[m + i] + polyWork0[m + i] * polyWork3[i];
    for (int i = 0; i < m; ++i) polyWork0[i] *= polyWork3[i];
    invFft(polyWork0, m);  // (1/2) E(n)
    invFft(polyWork0 + m, m);  // (1/2) E(n)
    for (int i = 0; i < m >> 1; ++i) polyWork0[(m >> 1) + i] += polyWork0[m + i];
    for (int i = 0; i < m; ++i) polyWork0[i] *= ::inv[m + i];
    for (int i = 0, i0 = min(m, size() - m); i < i0; ++i) polyWork0[i] += (*this)[m + i];
    memset(polyWork0 + m, 0, m * sizeof(Mint));
    fft(polyWork0, m << 1);  // 1 E(n)
    for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
    invFft(polyWork0, m << 1);  // 1 E(n)
    memcpy(fs.data() + m, polyWork0, (n - m) * sizeof(Mint));
    return fs;
  }
  // (29 + 1/2) E(n)
  // g <- g - (log g - a log t) g
  Poly pow1(Mint a, int n) const {
    assert(!empty()); assert((*this)[0].x == 1U); assert(1 <= n);
    return (a * log(n)).exp(n);  // 13 E(n) + (16 + 1/2) E(n)
  }
  // (29 + 1/2) E(n - a ord(t))
  Poly pow(long long a, int n) const {
    assert(a >= 0); assert(1 <= n);
    if (a == 0) { Poly gs(n); gs[0].x = 1U; return gs; }
    const int o = ord();
    if (o == -1 || o > (n - 1) / a) return Poly(n);
    const Mint b = (*this)[o].inv(), c = (*this)[o].pow(a);
    const int ntt = min<int>(n - a * o, size() - o);
    Poly tts(ntt);
    for (int i = 0; i < ntt; ++i) tts[i] = b * (*this)[o + i];
    tts = tts.pow1(a, n - a * o);  // (29 + 1/2) E(n - a ord(t))
    Poly gs(n);
    for (int i = 0; i < n - a * o; ++i) gs[a * o + i] = c * tts[i];
    return gs;
  }
  // (10 + 1/2) E(n)
  // f = t^(1/2) mod x^m,  g = 1 / t^(1/2) mod x^m
  // f <- f - (f^2 - h) g / 2
  // g <- g - (f g - 1) g
  // polyWork1: DFT(f, m), polyWork2: g, polyWork3: DFT(g, 2 m)
  Poly sqrt(int n) const {
    assert(!empty()); assert((*this)[0].x == 1U); assert(1 <= n);
    assert(n == 1 || 1 << (32 - __builtin_clz(n - 1)) <= LIM_POLY);
    if (n == 1) return {1U};
    if (n == 2) return {1U, at(1) / 2};
    Poly fs(n);
    fs[0].x = polyWork1[0].x = polyWork2[0].x = 1U;
    int m;
    for (m = 1; m << 1 < n; m <<= 1) {
      for (int i = 0; i < m; ++i) polyWork1[i] *= polyWork1[i];
      invFft(polyWork1, m);  // (1/2) E(n)
      for (int i = 0, i0 = min(m, size()); i < i0; ++i) polyWork1[i] -= (*this)[i];
      for (int i = 0, i0 = min(m, size() - m); i < i0; ++i) polyWork1[i] -= (*this)[m + i];
      memset(polyWork1 + m, 0, m * sizeof(Mint));
      fft(polyWork1, m << 1);  // 1 E(n)
      memcpy(polyWork3, polyWork2, m * sizeof(Mint));
      memset(polyWork3 + m, 0, m * sizeof(Mint));
      fft(polyWork3, m << 1);  // 1 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork1[i] *= polyWork3[i];
      invFft(polyWork1, m << 1);  // 1 E(n)
      for (int i = 0; i < m; ++i) { polyWork1[i] = -polyWork1[i]; fs[m + i].x = ((polyWork1[i].x & 1) ? (polyWork1[i].x + MO) : polyWork1[i].x) >> 1; }
      memcpy(polyWork1, fs.data(), (m << 1) * sizeof(Mint));
      fft(polyWork1, m << 1);  // 1 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork0[i] = polyWork1[i] * polyWork3[i];
      invFft(polyWork0, m << 1);  // 1 E(n)
      memset(polyWork0, 0, m * sizeof(Mint));
      fft(polyWork0, m << 1);  // 1 E(n)
      for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork3[i];
      invFft(polyWork0, m << 1);  // 1 E(n)
      for (int i = m; i < m << 1; ++i) polyWork2[i] = -polyWork0[i];
    }
    for (int i = 0; i < m; ++i) polyWork1[i] *= polyWork1[i];
    invFft(polyWork1, m);  // (1/2) E(n)
    for (int i = 0, i0 = min(m, size()); i < i0; ++i) polyWork1[i] -= (*this)[i];
    for (int i = 0, i0 = min(m, size() - m); i < i0; ++i) polyWork1[i] -= (*this)[m + i];
    memcpy(polyWork1 + m, polyWork1 + (m >> 1), (m >> 1) * sizeof(Mint));
    memset(polyWork1 + (m >> 1), 0, (m >> 1) * sizeof(Mint));
    memset(polyWork1 + m + (m >> 1), 0, (m >> 1) * sizeof(Mint));
    fft(polyWork1, m);  // (1/2) E(n)
    fft(polyWork1 + m, m);  // (1/2) E(n)
    memcpy(polyWork3 + m, polyWork2 + (m >> 1), (m >> 1) * sizeof(Mint));
    memset(polyWork3 + m + (m >> 1), 0, (m >> 1) * sizeof(Mint));
    fft(polyWork3 + m, m);  // (1/2) E(n)
    // for (int i = 0; i < m << 1; ++i) polyWork1[i] *= polyWork3[i];
    for (int i = 0; i < m; ++i) polyWork1[m + i] = polyWork1[i] * polyWork3[m + i] + polyWork1[m + i] * polyWork3[i];
    for (int i = 0; i < m; ++i) polyWork1[i] *= polyWork3[i];
    invFft(polyWork1, m);  // (1/2) E(n)
    invFft(polyWork1 + m, m);  // (1/2) E(n)
    for (int i = 0; i < m >> 1; ++i) polyWork1[(m >> 1) + i] += polyWork1[m + i];
    for (int i = 0; i < n - m; ++i) { polyWork1[i] = -polyWork1[i]; fs[m + i].x = ((polyWork1[i].x & 1) ? (polyWork1[i].x + MO) : polyWork1[i].x) >> 1; }
    return fs;
  }
  // (10 + 1/2) E(n)
  // modSqrt must return a quadratic residue if exists, or anything otherwise.
  // Return {} if *this does not have a square root.
  template <class F> Poly sqrt(int n, F modSqrt) const {
    assert(1 <= n);
    const int o = ord();
    if (o == -1) return Poly(n);
    if (o & 1) return {};
    const Mint c = modSqrt((*this)[o]);
    if (c * c != (*this)[o]) return {};
    if (o >> 1 >= n) return Poly(n);
    const Mint b = (*this)[o].inv();
    const int ntt = min(n - (o >> 1), size() - o);
    Poly tts(ntt);
    for (int i = 0; i < ntt; ++i) tts[i] = b * (*this)[o + i];
    tts = tts.sqrt(n - (o >> 1));  // (10 + 1/2) E(n)
    Poly gs(n);
    for (int i = 0; i < n - (o >> 1); ++i) gs[(o >> 1) + i] = c * tts[i];
    return gs;
  }
  // 6 E(|t|)
  // x -> x + a
  Poly shift(const Mint &a) const {
    if (empty()) return {};
    const int n = size();
    int m = 1;
    for (; m < n; m <<= 1) {}
    for (int i = 0; i < n; ++i) polyWork0[i] = fac[i] * (*this)[i];
    memset(polyWork0 + n, 0, ((m << 1) - n) * sizeof(Mint));
    fft(polyWork0, m << 1);  // 2 E(|t|)
    {
      Mint aa = 1;
      for (int i = 0; i < n; ++i) { polyWork1[n - 1 - i] = invFac[i] * aa; aa *= a; }
    }
    memset(polyWork1 + n, 0, ((m << 1) - n) * sizeof(Mint));
    fft(polyWork1, m << 1);  // 2 E(|t|)
    for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
    invFft(polyWork0, m << 1);  // 2 E(|t|)
    Poly fs(n);
    for (int i = 0; i < n; ++i) fs[i] = invFac[i] * polyWork0[n - 1 + i];
    return fs;
  }
};

Mint linearRecurrenceAt(const vector<Mint> &as, const vector<Mint> &cs, long long k) {
  assert(!cs.empty()); assert(cs[0]);
  const int d = cs.size() - 1;
  assert(as.size() >= static_cast<size_t>(d));
  return (Poly(vector<Mint>(as.begin(), as.begin() + d)) * cs).mod(d).divAt(cs, k);
}

struct SubproductTree {
  int logN, n, nn;
  vector<Mint> xs;
  // [DFT_4((X-xs[0])(X-xs[1])(X-xs[2])(X-xs[3]))] [(X-xs[0])(X-xs[1])(X-xs[2])(X-xs[3])mod X^4]
  // [         DFT_4((X-xs[0])(X-xs[1]))         ] [         DFT_4((X-xs[2])(X-xs[3]))         ]
  // [   DFT_2(X-xs[0])   ] [   DFT_2(X-xs[1])   ] [   DFT_2(X-xs[2])   ] [   DFT_2(X-xs[3])   ]
  vector<Mint> buf;
  vector<Mint *> gss;
  // (1 - xs[0] X) ... (1 - xs[nn-1] X)
  Poly all;
  // (ceil(log_2 n) + O(1)) E(n)
  SubproductTree(const vector<Mint> &xs_) {
    n = xs_.size();
    for (logN = 0, nn = 1; nn < n; ++logN, nn <<= 1) {}
    xs.assign(nn, 0U);
    memcpy(xs.data(), xs_.data(), n * sizeof(Mint));
    buf.assign((logN + 1) * (nn << 1), 0U);
    gss.assign(nn << 1, nullptr);
    for (int h = 0; h <= logN; ++h) for (int u = 1 << h; u < 1 << (h + 1); ++u) {
      gss[u] = buf.data() + (h * (nn << 1) + ((u - (1 << h)) << (logN - h + 1)));
    }
    for (int i = 0; i < nn; ++i) {
      gss[nn + i][0] = -xs[i] + 1;
      gss[nn + i][1] = -xs[i] - 1;
    }
    if (nn == 1) gss[1][1] += 2;
    for (int h = logN; --h >= 0; ) {
      const int m = 1 << (logN - h);
      for (int u = 1 << (h + 1); --u >= 1 << h; ) {
        for (int i = 0; i < m; ++i) gss[u][i] = gss[u << 1][i] * gss[u << 1 | 1][i];
        memcpy(gss[u] + m, gss[u], m * sizeof(Mint));
        invFft(gss[u] + m, m);  // ((1/2) ceil(log_2 n) + O(1)) E(n)
        if (h > 0) {
          gss[u][m] -= 2;
          const Mint a = FFT_ROOTS[logN - h + 1];
          Mint aa = 1;
          for (int i = m; i < m << 1; ++i) { gss[u][i] *= aa; aa *= a; };
          fft(gss[u] + m, m);  // ((1/2) ceil(log_2 n) + O(1)) E(n)
        }
      }
    }
    all.resize(nn + 1);
    all[0] = 1;
    for (int i = 1; i < nn; ++i) all[i] = gss[1][nn + nn - i];
    all[nn] = gss[1][nn] - 1;
  }
  // ((3/2) ceil(log_2 n) + O(1)) E(n) + 10 E(|f|) + 3 E(|f| + 2^(ceil(log_2 n)))
  vector<Mint> multiEval(const Poly &fs) const {
    vector<Mint> work0(nn), work1(nn), work2(nn);
    {
      const int m = max(fs.size(), 1);
      auto invAll = all.inv(m);  // 10 E(|f|)
      std::reverse(invAll.begin(), invAll.end());
      int mm;
      for (mm = 1; mm < m - 1 + nn; mm <<= 1) {}
      invAll.resize(mm, 0U);
      fft(invAll);  // E(|f| + 2^(ceil(log_2 n)))
      vector<Mint> ffs(mm, 0U);
      if (fs.size()) memcpy(ffs.data(), fs.data(), fs.size() * sizeof(Mint));
      fft(ffs);  // E(|f| + 2^(ceil(log_2 n)))
      for (int i = 0; i < mm; ++i) ffs[i] *= invAll[i];
      invFft(ffs);  // E(|f| + 2^(ceil(log_2 n)))
      memcpy(((logN & 1) ? work1 : work0).data(), ffs.data() + m - 1, nn * sizeof(Mint));
    }
    for (int h = 0; h < logN; ++h) {
      const int m = 1 << (logN - h);
      for (int u = 1 << h; u < 1 << (h + 1); ++u) {
        Mint *hs = (((logN - h) & 1) ? work1 : work0).data() + ((u - (1 << h)) << (logN - h));
        Mint *hs0 = (((logN - h) & 1) ? work0 : work1).data() + ((u - (1 << h)) << (logN - h));
        Mint *hs1 = hs0 + (m >> 1);
        fft(hs, m);  // ((1/2) ceil(log_2 n) + O(1)) E(n)
        for (int i = 0; i < m; ++i) work2[i] = gss[u << 1 | 1][i] * hs[i];
        invFft(work2.data(), m);  // ((1/2) ceil(log_2 n) + O(1)) E(n)
        memcpy(hs0, work2.data() + (m >> 1), (m >> 1) * sizeof(Mint));
        for (int i = 0; i < m; ++i) work2[i] = gss[u << 1][i] * hs[i];
        invFft(work2.data(), m);  // ((1/2) ceil(log_2 n) + O(1)) E(n)
        memcpy(hs1, work2.data() + (m >> 1), (m >> 1) * sizeof(Mint));
      }
    }
    work0.resize(n);
    return work0;
  }
  // ((5/2) ceil(log_2 n) + O(1)) E(n)
  Poly interpolate(const vector<Mint> &ys) const {
    assert(static_cast<int>(ys.size()) == n);
    Poly gs(n);
    for (int i = 0; i < n; ++i) gs[i] = (i + 1) * all[n - (i + 1)];
    const vector<Mint> denoms = multiEval(gs);  // ((3/2) ceil(log_2 n) + O(1)) E(n)
    vector<Mint> work(nn << 1, 0U);
    for (int i = 0; i < n; ++i) {
      // xs[0], ..., xs[n - 1] are not distinct
      assert(denoms[i]);
      work[i << 1] = work[i << 1 | 1] = ys[i] / denoms[i];
    }
    for (int h = logN; --h >= 0; ) {
      const int m = 1 << (logN - h);
      for (int u = 1 << (h + 1); --u >= 1 << h; ) {
        Mint *hs = work.data() + ((u - (1 << h)) << (logN - h + 1));
        for (int i = 0; i < m; ++i) hs[i] = gss[u << 1 | 1][i] * hs[i] + gss[u << 1][i] * hs[m + i];
        if (h > 0) {
          memcpy(hs + m, hs, m * sizeof(Mint));
          invFft(hs + m, m);  // ((1/2) ceil(log_2 n) + O(1)) E(n)
          const Mint a = FFT_ROOTS[logN - h + 1];
          Mint aa = 1;
          for (int i = m; i < m << 1; ++i) { hs[i] *= aa; aa *= a; };
          fft(hs + m, m);  // ((1/2) ceil(log_2 n) + O(1)) E(n)
        }
      }
    }
    invFft(work.data(), nn);  // E(n)
    return Poly(vector<Mint>(work.data() + nn - n, work.data() + nn));
  }
};

// q: rev([0, m]) * [0, n], [t^0] q(t, x) = 1 omitted
// ret: [0, m-1] * [0, n]
vector<Mint> comRec(int m, int n, const vector<Mint> &as, const vector<Mint> &qss) {
  if (!n) { auto ret = as; ret.resize(m, 0); return ret; }
  // reuse DFT(q(t, -x)); (2n+2) instead of (2n+1)
  int len;
  for (len = 2; len < (2*m) * (2*n+2); len <<= 1) {}
  vector<Mint> qs(len, 0);
  for (int i = 0; i < m; ++i) for (int j = 0; j <= n; ++j) qs[i * (2*n+2) + j] = qss[i * (n+1) + j];
  fft(qs);
  vector<Mint> work(len >> 1, 0);
  for (int k = 0; k < len >> 1; ++k) { work[k] = qs[k << 1] * qs[k << 1 | 1]; swap(qs[k << 1], qs[k << 1 | 1]); }
  invFft(work);
  vector<Mint> qqss((2*m) * (n/2+1), 0);
  for (int i = 0; i < 2*m-1; ++i) for (int j = 0; j <= n/2; ++j) qqss[i * (n/2+1) + j] = work[i * (n+1) + j];
  for (int i = 0; i < m; ++i) for (int j = 0; j <= n/2; ++j) qqss[(m+i) * (n/2+1) + j] += qss[i * (n+1) + 2*j] + qss[i * (n+1) + 2*j];
  const auto res = comRec(2*m, n/2, as, qqss);
  vector<Mint> ps(len, 0);
  for (int i = 0; i < 2*m; ++i) for (int j = 0; j <= n/2; ++j) ps[i * (2*n+2) + (2*n+1) - (2*j+(n&1))] = res[i * (n/2+1) + j];
  fft(ps);
  for (int k = 0; k < len; ++k) ps[k] *= qs[k];
  invFft(ps);
  vector<Mint> ret(m * (n+1));
  for (int i = 0; i < m; ++i) for (int j = 0; j <= n; ++j) ret[i * (n+1) + j] = ps[(m+i) * (2*n+2) + (2*n+1) - j];
  for (int i = 0; i < m; ++i) for (int j = 0; j <= n/2; ++j) ret[i * (n+1) + (2*j+(n&1))] += res[i * (n/2+1) + j];
  return ret;
}
// a(b(x)) mod x^n
// transpose and rev: p(x) -> [x^(n-1)] p(x) b(x)^i for each 0 <= i < n
// [x^(n-1)] p(x) / (1 - t b(x))
vector<Mint> com(const vector<Mint> &as, const vector<Mint> &bs, int n) {
  assert(bs.size() == 0 || !bs[0]);
  if (n == 0) return {};
  vector<Mint> qss(n, 0);
  for (int j = 0; j < min<int>(bs.size(), n); ++j) qss[j] = -bs[j];
  auto cs = comRec(1, n - 1, as, qss);
  reverse(cs.begin(), cs.end());
  return cs;
}
// [x^(n-1)] a(x) b(x)^i for each 0 <= i < n
// [x^(n-1)] a(x) / (1 - t b(x))
vector<Mint> powProj(const vector<Mint> &as, const vector<Mint> &bs, int n) {
  assert(bs.size() == 0 || !bs[0]);
  assert(n >= 1);
  // p(t, x): [0, m-1] * [0, n]
  // q(t, x): rev([0, m]) * [0, n],  [t^0] q(t, x) = 1 omitted
  vector<Mint> pss(n, 0), qss(n, 0);
  for (int j = 0; j < min<int>(as.size(), n); ++j) pss[j] = as[j];
  for (int j = 0; j < min<int>(bs.size(), n); ++j) qss[j] = -bs[j];
  const int n0 = n--;
  for (int m = 1; n; m *= 2, n /= 2) {
    // reuse DFT(q(t, -x)); (2n+2) instead of (2n+1)
    int len;
    for (len = 2; len < (m*2) * (n*2+2); len <<= 1) {}
    vector<Mint> qs(len, 0);
    for (int i = 0; i < m; ++i) for (int j = 0; j <= n; ++j) qs[i * (n*2+2) + j] = qss[i * (n+1) + j];
    fft(qs);
    vector<Mint> work(len >> 1, 0);
    for (int k = 0; k < len >> 1; ++k) { work[k] = qs[k << 1] * qs[k << 1 | 1]; swap(qs[k << 1], qs[k << 1 | 1]); }
    invFft(work);
    vector<Mint> qqss((m*2) * (n/2+1), 0);
    for (int i = 0; i <= m*2-2; ++i) for (int j = 0; j <= n/2; ++j) qqss[i * (n/2+1) + j] = work[i * (n+1) + j];
    for (int i = 0; i < m; ++i) for (int j = 0; j <= n/2; ++j) qqss[(m+i) * (n/2+1) + j] += qss[i * (n+1) + j*2] + qss[i * (n+1) + j*2];
    vector<Mint> ps(len, 0);
    for (int i = 0; i < m; ++i) for (int j = 0; j <= n; ++j) ps[i * (n*2+2) + j] = pss[i * (n+1) + j];
    fft(ps);
    for (int k = 0; k < len; ++k) ps[k] *= qs[k];
    invFft(ps);
    vector<Mint> ppss((m*2) * (n/2+1), 0);
    for (int i = 0; i <= m*2-2; ++i) for (int j = 0; j <= n/2; ++j) ppss[i * (n/2+1) + j] = ps[i * (n*2+2) + (j*2+(n&1))];
    for (int i = 0; i < m; ++i) for (int j = 0; j <= n/2; ++j) ppss[(m+i) * (n/2+1) + j] += pss[i * (n+1) + (j*2+(n&1))];
    pss.swap(ppss);
    qss.swap(qqss);
  }
  std::reverse(pss.begin(), pss.end());
  pss.resize(n0, 0);
  return pss;
}
// a^<-1>(x) mod x^n
// (n-1) [x^(n-1)] a(x)^i = i [x^(n-1-i)] (x/a^<-1>(x))^(n-1)
Poly comInv(const Poly &as, int n) {
  assert(as.size() >= 2); assert(!as[0]); assert(as[1]);
  assert(n >= 0);
  if (n <= 1) return Poly(n);
  // reduce to [x^1] a(x) = 1 in order to take (n-1)-th root
  const Mint t = as[1].inv();
  const auto res = powProj({1}, t * as, n);
  Poly ret(n - 1);
  for (int i = 1; i < n; ++i) ret[n - 1 - i] = inv[i] * (n - 1) * res[i];
  ret = ret.pow1(-inv[n - 1], n - 1);
  ret.insert(ret.begin(), 0);
  { Mint tt = 1; for (int i = 0; i < n; ++i) { ret[i] *= tt; tt *= t; } }
  return ret;
}
////////////////////////////////////////////////////////////////////////////////

Mint binom(Int n, Int k) {
  if (n < 0) {
    if (k >= 0) {
      return ((k & 1) ? -1 : +1) * binom(-n + k - 1, k);
    } else if (n - k >= 0) {
      return (((n - k) & 1) ? -1 : +1) * binom(-k - 1, n - k);
    } else {
      return 0;
    }
  } else {
    if (0 <= k && k <= n) {
      assert(n < LIM_INV);
      return fac[n] * invFac[k] * invFac[n - k];
    } else {
      return 0;
    }
  }
}

////////////////////////////////////////////////////////////////////////////////


int main() {
  int N, K;
  Int M;
  for (; ~scanf("%d%lld%d", &N, &M, &K); ) {
    Poly ys(K + 1);
    for (int k = 0; k <= K; ++k) ys[k] = binom(K, k) * (k&1?-1:+1);
    ys[0] -= 1;
    for (int k = 1; k <= K; ++k) ys[k] *= inv[k];
    ys = ys.exp(N);
    ys.insert(ys.begin(), 0);
    const Poly xs = comInv(ys, N + 1);
    Poly gs = xs.div(Poly{1} - xs, N + 1);
    for (int i = 0; i <= N; ++i) gs[i] *= Mint(i).pow(M);
    Poly ans = com(gs, ys, N + 1);
    for (int i = 1; i <= N; ++i) printf("%u\n", ans[i].x);
  }
  return 0;
}

这程序好像有点Bug,我给组数据试试?

Details

Tip: Click on the bar to expand more detailed information

Test #1:

score: 100
Accepted
time: 13ms
memory: 22200kb

input:

2 2 2

output:

1
10

result:

ok 2 number(s): "1 10"

Test #2:

score: 0
Accepted
time: 13ms
memory: 22820kb

input:

100 100 100

output:

1
288724001
477404412
847240327
388844119
60081063
382511336
656778732
715208757
458347567
173839136
330412968
871594390
54163986
520130283
894861257
697751100
191809807
596519881
264752553
380475772
774753694
205003415
258884583
507657187
224789401
441690786
733910671
722165795
522099197
529537234
...

result:

ok 100 numbers

Test #3:

score: 0
Accepted
time: 111ms
memory: 30496kb

input:

10000 10000 10000

output:

1
841616425
879018678
828464494
824353479
934437594
684261651
283427416
176711404
593462157
341689268
550748568
97309592
457838656
136079766
94711704
154856521
701974251
861457551
899123401
767934812
783681548
901411309
833556675
895864303
387820295
442453423
769529393
88723206
733802073
167561588
8...

result:

ok 10000 numbers

Test #4:

score: 0
Accepted
time: 2646ms
memory: 170604kb

input:

250000 250000 250000

output:

1
977538799
713071020
642656262
843437895
953939551
334515765
125548294
916390140
646942240
375064651
982842993
112721417
258413725
555579789
386912580
203353646
620154542
815908071
212216000
122585098
548938540
485418623
786547521
773855095
408984962
705041768
722485169
484323895
627828268
51253248...

result:

ok 250000 numbers

Test #5:

score: 0
Accepted
time: 2649ms
memory: 164704kb

input:

250000 1000000000000000000 250000

output:

1
674474200
356838072
423544639
620250749
472812708
467429466
425849422
544855209
138797438
719778389
306535642
986339502
436790422
266311982
415658715
242144610
253918786
742897815
823455763
914993586
706117605
639931316
875046605
162988809
135552036
369784489
175683457
731389739
194710792
60429203...

result:

ok 250000 numbers

Test #6:

score: 0
Accepted
time: 2656ms
memory: 170744kb

input:

250000 1000000000000000000 249999

output:

1
432274433
635272012
372335736
314816479
31062876
429528258
978201943
878737643
119637359
783592035
279481481
402908546
162484816
518316813
824200271
80678587
119112724
893737985
179077628
931294746
546897588
934868287
2843164
784278632
213754991
531691200
645062305
574638099
278555801
517596452
44...

result:

ok 250000 numbers

Test #7:

score: 0
Accepted
time: 2673ms
memory: 169424kb

input:

249999 1000000000000000000 250000

output:

1
674474200
356838072
423544639
620250749
472812708
467429466
425849422
544855209
138797438
719778389
306535642
986339502
436790422
266311982
415658715
242144610
253918786
742897815
823455763
914993586
706117605
639931316
875046605
162988809
135552036
369784489
175683457
731389739
194710792
60429203...

result:

ok 249999 numbers

Test #8:

score: 0
Accepted
time: 2680ms
memory: 163300kb

input:

241987 981121343321878999 238890

output:

1
951533200
944966565
737435680
26700008
201315166
13857231
855199761
496080089
183602747
50367126
250392634
806744764
31731796
265606981
763963523
659670280
428966498
489144773
350912463
114843760
585520736
200445151
69015747
953474979
151333103
506137448
299204044
745196763
576503791
417285710
778...

result:

ok 241987 numbers

Test #9:

score: 0
Accepted
time: 2664ms
memory: 164288kb

input:

250000 900000000000000000 250000

output:

1
88596811
299212258
411700355
485247006
865280274
41471638
3875210
134894519
390851180
888938548
801662167
740423714
508619978
663288694
854166871
283417700
434201304
784008787
565207810
375107087
244797356
825989537
45237315
684036726
791404940
943671358
108944323
598512816
598163457
842839538
681...

result:

ok 250000 numbers

Test #10:

score: 0
Accepted
time: 2656ms
memory: 163724kb

input:

250000 800000000000000000 250000

output:

1
33286348
879551824
577033758
910057381
782116970
545969399
166006368
610131824
581490833
20583070
814938489
690703255
812359682
380313911
835948097
278305226
311436254
989981734
46261489
951280591
433869480
714897746
55813990
328309210
343245000
317168015
887963472
507790568
200021753
292353163
88...

result:

ok 250000 numbers

Test #11:

score: 0
Accepted
time: 2654ms
memory: 164256kb

input:

250000 850000000000000000 250000

output:

1
850463981
919130890
594062815
853844081
340915288
69602994
423763703
744995209
210850487
993120441
964969301
141123493
187865084
532455528
172265983
950171432
220979336
880371218
18648835
39236248
290415692
318667152
951653052
966470245
819239424
548317720
215382408
569671462
404511634
275906691
5...

result:

ok 250000 numbers

Test #12:

score: 0
Accepted
time: 2630ms
memory: 164172kb

input:

250000 930000000000000000 233888

output:

1
913879354
443430322
736466644
397769128
959111419
252430611
887770075
242934908
439535113
251827104
557340629
557604248
312331037
335304552
420836877
851068915
578274735
950030714
260730488
788581921
259988973
179306065
850277354
644761072
472098191
124456383
5233243
76862373
533101079
521456435
7...

result:

ok 250000 numbers

Test #13:

score: 0
Accepted
time: 2624ms
memory: 154620kb

input:

140000 1000000000000000000 250000

output:

1
674474200
356838072
423544639
620250749
472812708
467429466
425849422
544855209
138797438
719778389
306535642
986339502
436790422
266311982
415658715
242144610
253918786
742897815
823455763
914993586
706117605
639931316
875046605
162988809
135552036
369784489
175683457
731389739
194710792
60429203...

result:

ok 140000 numbers

Test #14:

score: 0
Accepted
time: 12ms
memory: 23244kb

input:

101 100 102

output:

1
57234729
327110970
837424947
308436726
90521176
189415601
11967831
818517526
694746636
654362038
277172308
675497915
569788167
27780728
536402869
650201480
662086702
482770872
860152713
120258945
620375259
392523019
974084341
95322104
782869418
92301258
565905704
144429014
56540330
639211147
52457...

result:

ok 101 numbers

Test #15:

score: 0
Accepted
time: 12ms
memory: 23344kb

input:

123 456 789

output:

1
606065649
219325320
752236190
746458375
973373078
830581510
600287883
440906891
305388579
534354385
714656541
89412019
784882591
865510770
984400276
674074764
388411354
205689582
377495363
827083283
724922254
703920660
102797755
201740794
885679171
538258076
505743240
434068122
828364937
196472614...

result:

ok 123 numbers

Test #16:

score: 0
Accepted
time: 12ms
memory: 21132kb

input:

1 10 100

output:

1

result:

ok 1 number(s): "1"

Test #17:

score: 0
Accepted
time: 11ms
memory: 24520kb

input:

100 1 10

output:

1
12
78
364
1365
4368
12376
31824
75582
167960
352716
705432
1352078
2496144
4457400
7726160
13037895
21474180
34597290
54627300
84672315
129024480
193536720
286097760
417225900
600805296
854992152
205077935
677811691
315312734
164728909
287583964
760782584
681628661
168152380
363606607
446979616
63...

result:

ok 100 numbers

Test #18:

score: 0
Accepted
time: 14ms
memory: 22752kb

input:

100 10 1

output:

1
2047
261625
10391745
210766920
738437852
751995961
367882293
626598267
990684424
32946479
746153195
309367626
577393442
149727732
683395486
756615148
203162153
948422841
561114284
798562552
982203275
695555427
600675931
268185566
532424771
926160003
509859984
295191035
522289937
629317861
71197675...

result:

ok 100 numbers

Test #19:

score: 0
Accepted
time: 2662ms
memory: 163912kb

input:

250000 12121212121 250000

output:

1
245722720
107796473
237280018
692737912
416919748
631746435
711104024
415511334
702734699
976979524
725791457
801016883
47722471
631420909
462427749
359269088
874578178
480887167
708391129
967439571
630126451
648395143
40124209
850129997
313151332
958549730
593454741
967467417
567017114
814340975
...

result:

ok 250000 numbers

Extra Test:

score: 0
Extra Test Passed