結果

問題 No.502 階乗を計算するだけ
ユーザー miscalcmiscalc
提出日時 2023-12-02 04:01:20
言語 C++17
(gcc 12.3.0 + boost 1.83.0)
結果
AC  
実行時間 101 ms / 1,000 ms
コード長 50,107 bytes
コンパイル時間 4,450 ms
コンパイル使用メモリ 268,452 KB
実行使用メモリ 9,548 KB
最終ジャッジ日時 2023-12-02 04:01:32
合計ジャッジ時間 10,356 ms
ジャッジサーバーID
（参考情報） 		judge11 / judge15
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 68 ms
9,548 KB
testcase_01 AC 68 ms
9,548 KB
testcase_02 AC 68 ms
9,548 KB
testcase_03 AC 69 ms
9,548 KB
testcase_04 AC 68 ms
9,548 KB
testcase_05 AC 70 ms
9,548 KB
testcase_06 AC 69 ms
9,548 KB
testcase_07 AC 69 ms
9,548 KB
testcase_08 AC 68 ms
9,548 KB
testcase_09 AC 68 ms
9,548 KB
testcase_10 AC 68 ms
9,548 KB
testcase_11 AC 68 ms
9,548 KB
testcase_12 AC 68 ms
9,548 KB
testcase_13 AC 70 ms
9,548 KB
testcase_14 AC 68 ms
9,548 KB
testcase_15 AC 69 ms
9,548 KB
testcase_16 AC 68 ms
9,548 KB
testcase_17 AC 71 ms
9,548 KB
testcase_18 AC 67 ms
9,548 KB
testcase_19 AC 68 ms
9,548 KB
testcase_20 AC 68 ms
9,548 KB
testcase_21 AC 69 ms
9,548 KB
testcase_22 AC 68 ms
9,548 KB
testcase_23 AC 68 ms
9,548 KB
testcase_24 AC 68 ms
9,548 KB
testcase_25 AC 68 ms
9,548 KB
testcase_26 AC 68 ms
9,548 KB
testcase_27 AC 69 ms
9,548 KB
testcase_28 AC 69 ms
9,548 KB
testcase_29 AC 69 ms
9,548 KB
testcase_30 AC 68 ms
9,548 KB
testcase_31 AC 69 ms
9,548 KB
testcase_32 AC 68 ms
9,548 KB
testcase_33 AC 68 ms
9,548 KB
testcase_34 AC 68 ms
9,548 KB
testcase_35 AC 68 ms
9,548 KB
testcase_36 AC 101 ms
9,548 KB
testcase_37 AC 68 ms
9,548 KB
testcase_38 AC 67 ms
9,548 KB
testcase_39 AC 68 ms
9,548 KB
testcase_40 AC 68 ms
9,548 KB
testcase_41 AC 69 ms
9,548 KB
testcase_42 AC 68 ms
9,548 KB
testcase_43 AC 69 ms
9,548 KB
testcase_44 AC 72 ms
9,548 KB
testcase_45 AC 68 ms
9,548 KB
testcase_46 AC 69 ms
9,548 KB
testcase_47 AC 69 ms
9,548 KB
testcase_48 AC 68 ms
9,548 KB
testcase_49 AC 68 ms
9,548 KB
testcase_50 AC 69 ms
9,548 KB
testcase_51 AC 68 ms
9,548 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using ld = long double;
using ull = unsigned long long;
using pll = pair<ll, ll>;
using tlll = tuple<ll, ll, ll>;
constexpr ll INF = 1LL << 60;
template<class T> bool chmin(T& a, T b) {if (a > b) {a = b; return true;} return false;}
template<class T> bool chmax(T& a, T b) {if (a < b) {a = b; return true;} return false;}
ll safemod(ll A, ll M) {ll res = A % M; if (res < 0) res += M; return res;}
ll divfloor(ll A, ll B) {if (B < 0) A = -A, B = -B; return (A - safemod(A, B)) / B;}
ll divceil(ll A, ll B) {if (B < 0) A = -A, B = -B; return divfloor(A + B - 1, B);}
ll pow_ll(ll A, ll B) {if (A == 0 || A == 1) {return A;} if (A == -1) {return B & 1 ? -1 : 1;} ll res = 1; for (int i = 0; i < B; i++) {res *= A;} return res;}
ll mul_limited(ll A, ll B, ll M = INF) { return B == 0 ? 0 : A > M / B ? M : A * B; }
ll pow_limited(ll A, ll B, ll M = INF) { if (A == 0 || A == 1) {return A;} ll res = 1; for (int i = 0; i < B; i++) {if (res > M / A) return M; res *= A;} return res;}
template<class T> void unique(vector<T> &V) {V.erase(unique(V.begin(), V.end()), V.end());}
template<class T> void sortunique(vector<T> &V) {sort(V.begin(), V.end()); V.erase(unique(V.begin(), V.end()), V.end());}
#define FINALANS(A) do {cout << (A) << '\n'; exit(0);} while (false)
template<class T> void printvec(const vector<T> &V) {int _n = V.size(); for (int i = 0; i < _n; i++) cout << V[i] << (i == _n - 1 ? "" : " ");cout << '\n';}
template<class T> void printvect(const vector<T> &V) {for (auto v : V) cout << v << '\n';}
template<class T> void printvec2(const vector<vector<T>> &V) {for (auto &v : V) printvec(v);}
//*
#include <atcoder/modint>
#include <atcoder/math>
#include <atcoder/convolution>
#include <atcoder/internal_math>
using namespace atcoder;
//*/

// http://drken1215.hatenablog.com/entry/2018/06/08/210000
template <class T>
class binom
{
public:
  vector<T> fac, finv, inv;

  binom(int M)
  {
    fac.resize(M + 1);
    finv.resize(M + 1);
    inv.resize(M + 1);

    //*
    fac[0] = T(1);
    for (int i = 1; i <= M; i++)
      fac[i] = fac[i - 1] * T::raw(i);
    finv[M] = fac[M].inv();
    for (int i = M - 1; i >= 0; i--)
      finv[i] = finv[i + 1] * T::raw(i + 1);
    for (int i = 1; i <= M; i++)
      inv[i] = fac[i - 1] * finv[i];
    //*/

    /*
    fac[0] = T(1), finv[0] = T(1);
    fac[1] = T(1), finv[1] = T(1), inv[1] = T(1);
    for (int i = 2; i <= M; i++)
    {
      fac[i] = fac[i - 1] * i;
      inv[i] = -inv[T::mod() % i] * (T::mod() / i);
      finv[i] = finv[i - 1] * inv[i];
    }
    //*/
  }

  T P(int N, int K)
  {
    if (N < K) return 0;
    if (N < 0 || K < 0) return 0;
    return fac[N] * finv[N - K];
  }

  T C(int N, int K)
  {
    if (N < K) return 0;
    if (N < 0 || K < 0) return 0;
    return fac[N] * finv[K] * finv[N - K];
  }

  T H(int N, int K)
  {
    if (N == 0 && K == 0)
      return 1;
    return C(N + K - 1, K);
  }
};

// http://drken1215.hatenablog.com/entry/2018/06/08/210000
template <class T>
class binom_mut
{
private:
  vector<T> fac, finv, inv;
  void calc(int n)
  {
    int i = fac.size();
    if (n < i)
      return;
    fac.resize(n + 1), finv.resize(n + 1), inv.resize(n + 1);
    for (; i <= n; i++)
    {
      fac[i] = fac[i - 1] * i;
      inv[i] = -inv[T::mod() % i] * (T::mod() / i);
      finv[i] = finv[i - 1] * inv[i];
    }
  }
public:
  binom_mut() { fac = {1, 1}, finv = {1, 1}, inv = {0, 1}; }

  T get_fac(int n) { assert(n >= 0); calc(n); return fac[n]; }
  T get_finv(int n) { assert(n >= 0); calc(n); return finv[n]; }
  T get_inv(int n) { assert(n > 0); calc(n); return inv[n]; }

  T P(int N, int K)
  {
    if (N < K) return 0;
    if (N < 0 || K < 0) return 0;
    calc(N);
    return fac[N] * finv[N - K];
  }
  T C(int N, int K)
  {
    if (N < K) return 0;
    if (N < 0 || K < 0) return 0;
    calc(N);
    return fac[N] * finv[K] * finv[N - K];
  }
  T H(int N, int K)
  {
    if (N == 0 && K == 0)
      return 1;
    return C(N + K - 1, K);
  }
};

// https://qiita.com/taiyaki8926/items/f62f534d43ff006129f7
ll sqrt_mod(ll n, int p) // p は素数
{
  n %= p;
  if (n == 0)
    return 0;
  if (p == 2)
    return n;
  if (pow_mod(n, (p - 1) / 2, p) == p - 1) // 平方非剰余
    return -1;
  if (p % 4 == 3)
    return pow_mod(n, (p + 1) / 4, p);

  internal::barrett ba(p);
  int q = p - 1, s = 0;
  while (q % 2 == 0)
    q /= 2, s++;
  int z = 2;
  while (pow_mod(z, (p - 1) / 2, p) != p - 1)
    z++;
  int m = s;
  ll c = pow_mod(z, q, p);
  ll t = pow_mod(n, q, p);
  ll r = pow_mod(n, (q + 1) / 2, p);
  while (t != 1)
  {
    int m2 = 1;
    for (ll tmp = ba.mul(t, t); tmp != 1; tmp = ba.mul(tmp, tmp), m2++);
    ll b = pow_mod(c, 1 << (m - m2 - 1), p);
    m = m2, c = ba.mul(b, b), t = ba.mul(t, c), r = ba.mul(r, b);
  }
  return r;
}

template<const int MOD = 1000000007, class T>
vector<T> convolution_anymod(const vector<T> &A, const vector<T> &B)
{
  int N = A.size(), M = B.size();
  if (min(N, M) <= 300)
  {
    using mint = static_modint<MOD>;
    vector<mint> A2(N), B2(M);
    for (int i = 0; i < N; i++)
      A2[i] = A[i];
    for (int j = 0; j < M; j++)
      B2[j] = B[j];
    vector<mint> C2(N + M - 1, 0);
    for (int i = 0; i < N; i++)
      for (int j = 0; j < M; j++)
        C2[i + j] += A2[i] * B2[j];
    vector<T> C(N + M - 1);
    for (int i = 0; i < N + M - 1; i++)
      C[i] = C2[i].val();
    return C;
  }

  constexpr ll MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
  using mint2 = static_modint<MOD2>;
  using mint3 = static_modint<MOD3>;
  using mint4 = static_modint<MOD>;
  constexpr int i1_2 = internal::inv_gcd(MOD1, MOD2).second;
  constexpr int i12_3 = internal::inv_gcd(MOD1 * MOD2, MOD3).second;
  constexpr int m12_4 = MOD1 * MOD2 % MOD;

  auto C1 = convolution<MOD1>(A, B);
  auto C2 = convolution<MOD2>(A, B);
  auto C3 = convolution<MOD3>(A, B);

  vector<T> C(N + M - 1);
  for (ll i = 0; i < N + M - 1; i++)
  {
    int c1 = C1[i], c2 = C2[i], c3 = C3[i];
    int t1 = (mint2(c2 - c1) * mint2::raw(i1_2)).val();
    mint3 x2_m3 = mint3::raw(c1) + mint3::raw(t1) * mint3::raw(MOD1);
    mint4 x2_m = mint4::raw(c1) + mint4::raw(t1) * mint4::raw(MOD1);
    int t2 = ((mint3::raw(c3) - x2_m3) * mint3::raw(i12_3)).val();
    C[i] = (x2_m + mint4::raw(t2) * mint4::raw(m12_4)).val();
  }
  return C;
}
template<class T>
vector<T> convolution_anymod(const vector<T> &A, const vector<T> &B, const int MOD)
{
  int N = A.size(), M = B.size();
  if (min(N, M) <= 300)
  {
    using mint = dynamic_modint<100>;
    mint::set_mod(MOD);
    vector<mint> A2(N), B2(M);
    for (int i = 0; i < N; i++)
      A2[i] = A[i];
    for (int j = 0; j < M; j++)
      B2[j] = B[j];
    vector<mint> C2(N + M - 1, 0);
    for (int i = 0; i < N; i++)
      for (int j = 0; j < M; j++)
        C2[i + j] += A2[i] * B2[j];
    vector<T> C(N + M - 1);
    for (int i = 0; i < N + M - 1; i++)
      C[i] = C2[i].val();
    return C;
  }

  constexpr ll MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
  using mint2 = static_modint<MOD2>;
  using mint3 = static_modint<MOD3>;
  using mint4 = dynamic_modint<100>;
  mint4::set_mod(MOD);
  constexpr int i1_2 = internal::inv_gcd(MOD1, MOD2).second;
  constexpr int i12_3 = internal::inv_gcd(MOD1 * MOD2, MOD3).second;

  auto C1 = convolution<MOD1>(A, B);
  auto C2 = convolution<MOD2>(A, B);
  auto C3 = convolution<MOD3>(A, B);

  vector<T> C(N + M - 1);
  for (ll i = 0; i < N + M - 1; i++)
  {
    int c1 = C1[i], c2 = C2[i], c3 = C3[i];
    int t1 = (mint2(c2 - c1) * mint2::raw(i1_2)).val();
    mint3 x2_m3 = mint3::raw(c1) + mint3::raw(t1) * mint3::raw(MOD1);
    mint4 x2_m = mint4::raw(c1) + mint4::raw(t1) * mint4::raw(MOD1);
    int t2 = ((mint3::raw(c3) - x2_m3) * mint3::raw(i12_3)).val();
    C[i] = (x2_m + mint4::raw(t2) * mint4::raw(MOD1) * mint4::raw(MOD2)).val();
  }
  return C;
}
template<const int MOD>
vector<static_modint<MOD>> convolution_anymod(const vector<static_modint<MOD>> &A, const vector<static_modint<MOD>> &B)
{
  int N = A.size(), M = B.size();
  vector<int> A2(N), B2(M);
  for (int i = 0; i < N; i++)
    A2[i] = A[i].val();
  for (int i = 0; i < M; i++)
    B2[i] = B[i].val();
  vector<int> C2 = convolution_anymod<MOD>(A2, B2);
  vector<static_modint<MOD>> C(N + M - 1);
  for (int i = 0; i < N + M - 1; i++)
    C[i] = static_modint<MOD>::raw(C2[i]);
  return C;
}
template<const int id>
vector<dynamic_modint<id>> convolution_anymod(const vector<dynamic_modint<id>> &A, const vector<dynamic_modint<id>> &B)
{
  int N = A.size(), M = B.size();
  vector<int> A2(N), B2(M);
  for (int i = 0; i < N; i++)
    A2[i] = A[i].val();
  for (int i = 0; i < M; i++)
    B2[i] = B[i].val();
  vector<int> C2 = convolution_anymod(A2, B2, dynamic_modint<id>::mod());
  vector<dynamic_modint<id>> C(N + M - 1);
  for (int i = 0; i < N + M - 1; i++)
    C[i] = dynamic_modint<id>::raw(C2[i]);
  return C;
}

template<class T1>
struct LagrangeInterpolation
{
  int D;
  vector<T1> Y, fac, finv, prodl, prodr;
  
  template<class T2>
  LagrangeInterpolation(const vector<T2> &y)
  {
    D = (int)y.size() - 1;
    Y.resize(D + 1);
    for (int i = 0; i <= D; i++)
    {
      Y[i] = y[i];
    }

    fac.resize(D + 1), finv.resize(D + 1);
    fac[0] = 1;
    for (int i = 1; i <= D; i++)
      fac[i] = fac[i - 1] * i;
    finv[D] = fac[D].inv();
    for (int i = D - 1; i >= 0; i--)
      finv[i] = finv[i + 1] * (i + 1);

    prodl.resize(D + 2), prodr.resize(D + 2);
  }

  T1 eval(T1 x)
  {
    prodl[0] = 1;
    for (int i = 0; i <= D; i++)
    {
      prodl[i + 1] = prodl[i] * (x - i);
    }
    prodr[D + 1] = 1;
    for (int i = D; i >= 0; i--)
    {
      prodr[i] = prodr[i + 1] * (x - i);
    }

    T1 res = 0;
    for (int i = 0; i <= D; i++)
    {
      T1 tmp = Y[i] * prodl[i] * prodr[i + 1] * finv[i] * finv[D - i];
      if ((D - i) % 2 == 0)
        res += tmp;
      else
        res -= tmp;
    }
    return res;
  }
};

// https://opt-cp.com/fps-implementation/
// https://qiita.com/hotman78/items/f0e6d2265badd84d429a
// https://opt-cp.com/fps-fast-algorithms/
// https://maspypy.com/%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%83%BB%E5%BD%A2%E5%BC%8F%E7%9A%84%E3%81%B9%E3%81%8D%E7%B4%9A%E6%95%B0-%E9%AB%98%E9%80%9F%E3%81%AB%E8%A8%88%E7%AE%97%E3%81%A7%E3%81%8D%E3%82%8B%E3%82%82%E3%81%AE
template<class T, bool is_ntt_friendly>
struct FormalPowerSeries : vector<T>
{
private:
  static vector<T> fac, finv, invmint;
  void calc(int n)
  {
    while ((int)fac.size() <= n)
    {
      int i = fac.size();
      fac.emplace_back(fac[i - 1] * i);
      invmint.emplace_back(-invmint[T::mod() % i] * (T::mod() / i));
      finv.emplace_back(finv[i - 1] * invmint[i]);
    }
  }

public:
  T get_fac(int n) { calc(n); return fac[n]; }
  T get_finv(int n) { calc(n); return finv[n]; }
  T get_invmint(int n) { calc(n); return invmint[n]; }

  using vector<T>::vector;
  using vector<T>::operator=;
  using F = FormalPowerSeries;
  using S = vector<pair<ll, T>>;

  FormalPowerSeries(const S &f, int n = -1)
  {
    if (n == -1)
      n = f.back().first + 1;
    (*this).assign(n, T(0));
    for (auto [d, a] : f)
      (*this)[d] += a;
  }

  F operator-() const
  {
    F res(*this);
    for (auto &a : res)
      a = -a;
    return res;
  }

  F operator*=(const T &k)
  {
    for (auto &a : *this)
      a *= k;
    return *this;
  }
  F operator*(const T &k) const { return F(*this) *= k; }
  friend F operator*(const T k, const F &f) { return f * k; }
  F operator/=(const T &k)
  {
    *this *= k.inv();
    return *this;
  }
  F operator/(const T &k) const { return F(*this) /= k; }

  F &operator+=(const F &g)
  {
    int n = (*this).size(), m = g.size();
    (*this).resize(max(n, m), T(0));
    for (int i = 0; i < m; i++)
      (*this)[i] += g[i];
    return *this;
  }
  F operator+(const F &g) const { return F(*this) += g; }
  F &operator-=(const F &g)
  {
    int n = (*this).size(), m = g.size();
    (*this).resize(max(n, m), T(0));
    for (int i = 0; i < m; i++)
      (*this)[i] -= g[i];
    return *this;
  }
  F operator-(const F &g) const { return F(*this) -= g; }

  F &operator<<=(const ll d)
  {
    int n = (*this).size();
    (*this).insert((*this).begin(), min(ll(n), d), T(0));
    (*this).resize(n);
    return *this;
  }
  F operator<<(const ll d) const { return F(*this) <<= d; }
  F &operator>>=(const ll d)
  {
    int n = (*this).size();
    (*this).erase((*this).begin(), (*this).begin() + min(ll(n), d));
    (*this).resize(n, T(0));
    return *this;
  }
  F operator>>(const ll d) const { return F(*this) >>= d; }

  F &operator*=(const S &g)
  {
    int n = (*this).size();
    auto [d, c] = g.front();
    if (d != 0)
      c = 0;
    for (int i = n - 1; i >= 0; i--)
    {
      (*this)[i] *= c;
      for (auto &[j, b] : g)
      {
        if (j == 0)
          continue;
        if (j > i)
          break;
        (*this)[i] += (*this)[i - j] * b;
      }
    }
    return *this;
  }
  F operator*(const S &g) const { return F(*this) *= g; }
  F &operator/=(const S &g)
  {
    int n = (*this).size();
    auto [d, c] = g.front();
    assert(d == 0 && c != T(0));
    T inv_c = c.inv();
    for (int i = 0; i < n; i++)
    {
      for (auto &[j, b] : g)
      {
        if (j == 0)
          continue;
        if (j > i)
          break;
        (*this)[i] -= (*this)[i - j] * b;
      }
      (*this)[i] *= inv_c;
    }
    return *this;
  }
  F operator/(const S &g) const { return F(*this) /= g; }

  // (1 + cx^d) を掛ける
  F multiply(const int d, const T c)
  {
    int n = (*this).size();
    if (c == T(1))
    {
      for (int i = n - 1 - d; i >= 0; i--)
        (*this)[i + d] += (*this)[i];
    }
    else if (c == T(-1))
    {
      for (int i = n - 1 - d; i >= 0; i--)
        (*this)[i + d] -= (*this)[i];
    }
    else
    {
      for (int i = n - 1 - d; i >= 0; i--)
        (*this)[i + d] += (*this)[i] * c;
    }
    return *this;
  }
  F multiplication(const int d, const T c) const { return multiply(F(*this)); }
  // (1 + cx^d) で割る
  F divide(const int d, const T c)
  {
    int n = (*this).size();
    if (c == T(1))
    {
      for (int i = 0; i < n - d; i++)
        (*this)[i + d] -= (*this)[i];
    }
    else if (c == T(-1))
    {
      for (int i = 0; i < n - d; i++)
        (*this)[i + d] += (*this)[i];
    }
    else
    {
      for (int i = 0; i < n - d; i++)
        (*this)[i + d] -= (*this)[i] * c;
    }
    return *this;
  }
  F division(const int d, const T c) const { return divide(F(*this)); }

  template<const int MOD>
  F convolution2(const vector<static_modint<MOD>> &A, const vector<static_modint<MOD>> &B, const int d = -1) const
  {
    F res;
    if (is_ntt_friendly)
      res = convolution(A, B);
    else
      res = convolution_anymod(A, B);
    if (d != -1 && (int)res.size() > d)
      res.resize(d);
    return res;
  }
  template<const int id>
  F convolution2(const vector<dynamic_modint<id>> &A, const vector<dynamic_modint<id>> &B, const int d = -1) const
  {
    F res;
    res = convolution_anymod(A, B);
    if (d != -1 && (int)res.size() > d)
      res.resize(d);
    return res;
  }

  F &operator*=(const F &g)
  {
    int n = (*this).size();
    if (n == 0)
      return *this;
    *this = convolution2(*this, g, n);
    return *this;
  }
  F operator*(const F &g) const { return F(*this) *= g; }

  template <const int MOD>
  void butterfly2(FormalPowerSeries<static_modint<MOD>, true> &A) const { internal::butterfly(A); }
  template <const int MOD>
  void butterfly2(FormalPowerSeries<static_modint<MOD>, false> &A) const { assert(false); }
  template <const int id>
  void butterfly2(FormalPowerSeries<dynamic_modint<id>, false> &A) const { assert(false); }
  template <const int MOD>
  void butterfly_inv2(FormalPowerSeries<static_modint<MOD>, true> &A) const { internal::butterfly_inv(A); }
  template <const int MOD>
  void butterfly_inv2(FormalPowerSeries<static_modint<MOD>, false> &A) const { assert(false); }
  template <const int id>
  void butterfly_inv2(FormalPowerSeries<dynamic_modint<id>, false> &A) const { assert(false); }

  // mod (x^n - 1) をとったものを返す
  F circular_mod(int n) const
  {
    F res(n, T(0));
    for (int i = 0; i < (int)(*this).size(); i++)
      res[i % n] += (*this)[i];
    return res;
  }

  F inv(int d = -1) const
  {
    int n = (*this).size();
    assert(!(*this).empty() && (*this).at(0) != T(0));
    if (d == -1)
      d = n;
    //assert(d > 0);
    F f, g2;
    F g{(*this).front().inv()};
    while ((int)g.size() < d)
    {
      if (is_ntt_friendly)
      {
        int m = g.size();
        f = F{(*this).begin(), (*this).begin() + min(n, 2 * m)};
        g2 = F(g);
        f.resize(2 * m, T(0)), butterfly2(f);
        g2.resize(2 * m, T(0)), butterfly2(g2);
        for (int i = 0; i < 2 * m; i++)
          f[i] *= g2[i];
        butterfly_inv2(f);
        f.erase(f.begin(), f.begin() + m);
        f.resize(2 * m, T(0)), butterfly2(f);
        for (int i = 0; i < 2 * m; i++)
          f[i] *= g2[i];
        butterfly_inv2(f);
        T iz = T(2 * m).inv();
        iz *= -iz;
        for (int i = 0; i < m; i++)
          f[i] *= iz;
        g.insert(g.end(), f.begin(), f.begin() + m);
      }
      else
      {
        g.resize(2 * g.size(), T(0));
        g *= F{T(2)} - g * (*this);
      }
    }
    return {g.begin(), g.begin() + d};
  }
  F &operator/=(const F &g)
  {
    *this *= g.inv((*this).size());
    return *this;
  }
  F operator/(const F &g) const { return F(*this) *= g.inv((*this).size()); }

  F differentiate()
  {
    *this >>= 1;
    for (int i = 0; i < int((*this).size()) - 1; i++)
      (*this)[i] *= i + 1;
    return *this;
  }
  F differential() const { return F(*this).differentiate(); }
  F integrate()
  {
    int n = (*this).size();
    vector<T> minv(n);
    minv[1] = T(1);
    *this <<= 1;
    for (int i = 2; i < n; i++)
    {
      minv[i] = -minv[T::mod() % i] * (T::mod() / i);
      (*this)[i] *= minv[i];
    }
    return *this;
  }
  F integral() const { return F(*this).integrate(); }

  F log() const
  {
    assert((*this).front() == T(1));
    return ((*this).differential() / (*this)).integral();
  }
  F exp() const // https://arxiv.org/pdf/1301.5804.pdf
  {
    int n = (*this).size();
    assert(n != 0 && (*this).front() == T(0));

    //*
    if (is_ntt_friendly)
    {
      F f{T(1)}, g{T(1)};
      F dh = (*this).differential();
      F f2, g2, f3, q, s, h, u;
      g2 = {T(0)};
      while ((int)f.size() < n)
      {
        int m = f.size();
        T im = T(m).inv(), i2m = T(2 * m).inv();
        f2 = F(f);
        f2.resize(2 * m), butterfly2(f2);

        // a
        F f3(f);
        butterfly2(f3);
        for (int i = 0; i < m; i++)
          f3[i] *= g2[i];
        butterfly_inv2(f3);
        f3.erase(f3.begin(), f3.begin() + m / 2);
        f3.resize(m, T(0)), butterfly2(f3);
        for (int i = 0; i < m; i++)
          f3[i] *= g2[i];
        butterfly_inv2(f3);
        for (int i = 0; i < m / 2; i++)
          f3[i] *= -im * im;
        g.insert(g.end(), f3.begin(), f3.begin() + m / 2);

        g2 = F(g), g2.resize(2 * m), butterfly2(g2);

        // b, c
        q = F(dh);
        q.resize(2 * m);
        for (int i = m - 1; i < 2 * m; i++)
          q[i] = T(0);
        butterfly2(q);
        for (int i = 0; i < 2 * m; i++)
          q[i] *= f2[i];
        butterfly_inv2(q);
        q = q.circular_mod(m);
        for (int i = 0; i < m; i++)
          q[i] *= i2m;

        // d, e
        q.resize(m + 1);
        s = ((f.differential() - q) << 1).circular_mod(m);
        s.resize(2 * m);
        butterfly2(s);
        for (int i = 0; i < 2 * m; i++)
          s[i] *= g2[i];
        butterfly_inv2(s);
        for (int i = 0; i < m; i++)
          s[i] *= i2m;
        s.resize(m);

        // f, g
        h = (*this);
        h.resize(2 * m), s.resize(2 * m);
        u = (h - (s << (m - 1)).integral()) >> m;
        butterfly2(u);
        for (int i = 0; i < 2 * m; i++)
          u[i] *= f2[i];
        butterfly_inv2(u);
        for (int i = 0; i < m; i++)
          u[i] *= i2m;
        u.resize(m);

        // h
        f.insert(f.end(), u.begin(), u.end());
      }
      return {f.begin(), f.begin() + n};
    }
    else
    //*/
    {
      F f{T(1)}, g{T(1)};
      while ((int)f.size() < n)
      {
        int m = f.size();
        g = convolution2(g, F{T(2)} - f * g, m);
        F q = (*this).differential();
        q.resize(m - 1);
        F r = f.convolution2(f, q).circular_mod(m);
        r.resize(m + 1);
        F s = ((f.differential() - r) << 1).circular_mod(m);
        F t = g * s;
        F h = (*this);
        h.resize(2 * m), t.resize(2 * m);
        F u = (h - (t << (m - 1)).integral()) >> m;
        F v = f * u;
        f.insert(f.end(), v.begin(), v.end());
      }
      return {f.begin(), f.begin() + n};
      /*
      F f{T(1)};
      while ((int)f.size() < n)
      {
        int m = f.size();
        f.resize(min(n, 2 * m), T(0));
        f *= (*this) + F{T(1)} - f.log();
      }
      return f;
      //*/
    }
  }

  F pow(const ll k) const
  {
    if (k == 0)
    {
      F res((*this).size(), T(0));
      res[0] = T(1);
      return res;
    }
    int n = (*this).size(), d;
    for (d = 0; d < n; d++)
    {
      if ((*this)[d] != T(0))
        break;
    }
    if (d == n)
      return F(n, 0);
    F res = F(*this) >> d;
    T c = res[0];
    res /= c;
    res = (res.log() * T(k)).exp();
    res *= c.pow(k), res <<= (d != 0 && k > n ? n : d * k);
    return res;
  }

  F powmod(ll k, const F &g) const
  {
    F res(2 * g.size(), 0);
    res.front() = 1;
    F tmp = (*this) % g;
    tmp.resize(g.size());
    while (k > 0)
    {
      if (k & 1)
      {
        res *= tmp;
        res %= g;
        res.resize(2 * g.size());
      }
      tmp = tmp.convolution2(tmp, tmp);
      tmp %= g;
      tmp.resize(g.size());
      k >>= 1;
    }
    return res;
  }

  // f(x)^k mod (x^n - 1)
  F powmod_circular(ll k, ll n) const
  {
    F res(n, 0);
    res.front() = 1;
    F tmp = (*this).circular_mod(n);
    while (k > 0)
    {
      if (k & 1)
        res = res.convolution2(res, tmp).circular_mod(n);
      tmp = tmp.convolution2(tmp, tmp).circular_mod(n);
      k >>= 1;
    }
    return res;
  }

  // 素数 mod を要求
  // 存在しないなら空配列を返す
  F sqrt() const
  {
    int n = (*this).size(), d;
    for (d = 0; d < n; d += 2)
    {
      if ((*this)[d] != 0)
        break;
      if (d + 1 < n && (*this)[d + 1] != 0)
        return F(0);
    }
    if (d >= n)
      return F(n, 0);

    T a = (*this)[d];
    int p = T::mod();
    int r = sqrt_mod(a.val(), p);
    if (r == -1)
      return F(0);

    T inv_2 = T(2).inv();
    F f = F(*this) >> d, res = F{r};
    while (res.size() < f.size())
    {
      res.resize(min(f.size(), 2 * res.size()), T(0));
      res = (res + res.inv() * f) * inv_2;
    }
    res <<= d / 2;
    return res;
  }

  F div_poly(const F &g) const
  {
    F f2 = F(*this), g2 = F(g);
    while (!f2.empty() && f2.back() == T(0))
      f2.pop_back();
    while (!g2.empty() && g2.back() == T(0))
      g2.pop_back();
    
    int n = f2.size(), m = g2.size();
    int k = n - m + 1;
    if (k <= 0)
      return F{};

    reverse(f2.begin(), f2.end());
    reverse(g2.begin(), g2.end());
    f2.resize(k, T(0)), g2.resize(k, T(0));
    F q = f2 / g2;
    reverse(q.begin(), q.end());
    while (!q.empty() && q.back() == T(0))
      q.pop_back();
    return q;
  }
  pair<F, F> divmod(const F &g) const
  {
    int m = g.size();
    assert(m != 0);
    F q = (*this).div_poly(g);

    F f3 = F(*this), g3 = F(g), q3 = F(q);
    f3.resize(m - 1, T(0)), g3.resize(m - 1, T(0)), q3.resize(m - 1, T(0));
    F r = f3 - q3 * g3;
    while (!r.empty() && r.back() == T(0))
      r.pop_back();

    return make_pair(q, r);
  }
  F operator%(const F &g) const { return (*this).divmod(g).second; }
  F &operator%=(const F &g) { return (*this) = (*this) % g; }

  F div_poly(const S &g) const
  {
    F f2 = F(*this);
    while (!f2.empty() && f2.back() == T(0))
      f2.pop_back();
    assert(!g.empty());

    int n = f2.size(), m = g.back().first + 1;
    int k = n - m + 1;
    if (k <= 0)
      return F{};

    reverse(f2.begin(), f2.end());
    S g2(g.size());
    for (int i = 0; i < (int)g.size(); i++)
      g2[(int)g.size() - 1 - i] = make_pair(m - 1 - g[i].first, g[i].second);
    f2.resize(k, T(0));
    F q = f2 / g2;
    reverse(q.begin(), q.end());
    while (!q.empty() && q.back() == T(0))
      q.pop_back();
    return q;
  }
  pair<F, F> divmod(const S &g) const
  {
    assert(!g.empty());
    int m = g.back().first + 1;
    F q = (*this).div_poly(g);

    F f3 = F(*this), q3 = F(q);
    f3.resize(m - 1, T(0)), q3.resize(m - 1, T(0));
    F r = f3 - q3 * g;
    while (!r.empty() && r.back() == T(0))
      r.pop_back();

    return make_pair(q, r);
  }
  F operator%(const S &g) const { return (*this).divmod(g).second; }
  F &operator%=(const S &g) { return (*this) = (*this) % g; }

  T eval(const T &x) const
  {
    T res(0);
    for (int i = (int)(*this).size() - 1; i >= 0; i--)
    {
      res *= x;
      res += (*this)[i];
    }
    return res;
  }

  // 各係数 a_n を n! で割る
  F to_egf()
  {
    for (int i = 0; i < (int)(*this).size(); i++)
      (*this)[i] *= get_finv(i);
    return (*this);
  }
  // 各係数 a_n を n! で割ったものを返す
  F get_egf() const { return F(*this).to_egf(); }
  // 各係数 a_n に n! をかける
  F to_ogf()
  {
    for (int i = 0; i < (int)(*this).size(); i++)
      (*this)[i] *= get_fac(i);
    return (*this);
  }
  // 各係数 a_n に n! をかけたものを返す
  F get_ogf() const { return F(*this).to_ogf(); }

  F taylor_shift(const T &c) const
  {
    int n = (*this).size();

    F f = F(*this).get_ogf();
    reverse(f.begin(), f.end());
    F g = F(n);
    g[0] = 1;
    for (int i = 1; i < n; i++)
      g[i] = c * g[i - 1];
    g.to_egf();
    F h = f * g;
    reverse(h.begin(), h.end());
    return h.to_egf();
  }

  vector<T> eval_multipoint(const vector<T> &xs) const
  {
    int m0 = xs.size(), m = 1;
    while (m < m0)
      m <<= 1;
    vector<F> node(2 * m, F{1});
    for (int i = 0; i < m0; i++)
      node[m + i] = {-xs[i], T(1)};
    for (int i = m - 1; i > 0; i--)
      node[i] = convolution2(node[i << 1], node[(i << 1) | 1]);

    node[1] = (*this).divmod(node[1]).second;
    for (int i = 2; i < m + m0; i++)
      node[i] = node[i >> 1].divmod(node[i]).second;

    vector<T> res(m0);
    for (int i = 0; i < m0; i++)
      res[i] = node[m + i].empty() ? T(0) : node[m + i][0];
    return res;
  }

  // i = 0..m-1 に対する f(ar^i)
  // https://noshi91.github.io/algorithm-encyclopedia/chirp-z-transform
  vector<T> eval_multipoint_geo(int m, T a, T r) const
  {
    if (r == 0)
    {
      vector<T> res(m, (*this).eval(0));
      res.front() = (*this).eval(a);
      return res;
    }

    auto calc_pw = [&](T x, int k) -> vector<T>
    {
      vector<T> res(k);
      res.front() = 1;
      for (int i = 1; i < k; i++)
        res[i] = res[i - 1] * x;
      return res;
    };
    auto get_pw_tri = [&](const vector<T> &pw, int i) -> T
    {
      if (i == 0)
        return 1;
      return i % 2 == 0 ? pw[i - 1].pow(i / 2) : pw[i].pow((i - 1) / 2); 
    };
    int n = (*this).size();
    T invr = r.inv();
    vector<T> pwa = calc_pw(a, n), pwr = calc_pw(r, n + m), pwir = calc_pw(invr, max(n, m));

    vector<T> s(n), t(n + m);
    for (int i = 0; i < n; i++)
      s[n - 1 - i] = (*this)[i] * pwa[i] * get_pw_tri(pwir, i);
    for (int i = 0; i < n + m; i++)
      t[i] = get_pw_tri(pwr, i);
    vector<T> u = convolution2(s, t, n + m - 1);
    u.erase(u.begin(), u.begin() + n - 1);
    for (int i = 0; i < m; i++)
      u[i] *= get_pw_tri(pwir, i);
    return u;
  }
};

// (次数, 係数) を昇順に並べたもの
template <class T, bool is_ntt_friendly>
struct SparseFormalPowerSeries : vector<pair<ll, T>>
{
  using vector<pair<ll, T>>::vector;
  using vector<pair<ll, T>>::operator=;
  using F = FormalPowerSeries<T, is_ntt_friendly>;
  using S = SparseFormalPowerSeries;

  F to_fps(int n) const
  {
    F res(n, T(0));
    for (auto [d, a] : (*this))
      res[d] += a;
    return res;
  }

  SparseFormalPowerSeries(const F &f)
  {
    (*this).clear();
    for (int i = 0; i < (int)f.size(); i++)
    {
      if (f[i] != T(0))
        (*this).emplace_back(make_pair(i, f[i]));
    }
  }

  S operator-() const
  {
    S res(*this);
    for (auto &[d, a] : res)
      a = -a;
    return res;
  }

  S operator*=(const T &k)
  {
    for (auto &[d, a] : (*this))
      a *= k;
    return (*this);
  }
  S operator/=(const T &k)
  {
    (*this) *= k.inv();
    return (*this);
  }
  S operator*(const T &k) const { return S(*this) *= k; }
  S operator/(const T &k) const { return S(*this) /= k; }
  friend S operator*(const T k, const S &f) { return f * k; }

  S operator+(const S &g) const
  {
    S res;
    int n = (*this).size(), m = g.size(), i = 0, j = 0;
    while (i < n || j < m)
    {
      pair<ll, T> tmp;
      if (j == m || (i != n && (*this)[i].first <= g[j].first))
        tmp = (*this)[i++];
      else
        tmp = g[j++];

      if (!res.empty() && res.back().first == tmp.first)
        res.back().second += tmp.second;
      else
        res.emplace_back(tmp);
    }
    return res;
  }
  S operator-(const S &g) const
  {
    S res;
    int n = (*this).size(), m = g.size(), i = 0, j = 0;
    while (i < n || j < m)
    {
      pair<ll, T> tmp;
      if (j == m || (i != n && (*this)[i].first <= g[j].first))
        tmp = (*this)[i++];
      else
      {
        tmp = g[j++];
        tmp.second = -tmp.second;
      }

      if (!res.empty() && res.back().first == tmp.first)
        res.back().second += tmp.second;
      else
        res.emplace_back(tmp);
    }
    return res;
  }
  S operator*(const S &g) const
  {
    S res;
    for (auto [d, a] : (*this))
      for (auto [e, b] : g)
        res.emplace_back(make_pair(d + e, a * b));
    sort(res.begin(), res.end(), [&](pair<ll, T> p1, pair<ll, T> p2)
         { return p1.first < p2.first; });
    S res2;
    for (auto da : res)
    {
      auto [d, a] = da;
      if (res2.empty() || res2.back().first != d)
        res2.emplace_back(da);
      else
        res2.back().second += a;
    }
    return res2;
  }
  S operator+=(const S &g) { return (*this) = (*this) + g; }
  S operator-=(const S &g) { return (*this) = (*this) - g; }
  S operator*=(const S &g) { return (*this) = (*this) * g; }

  S operator<<=(ll k)
  {
    for (auto &[d, a] : (*this))
      d += k;
    return (*this);
  }
  S operator<<(ll k) const { return (*this) <<= k; }
  S operator>>(ll k) const
  {
    S res;
    for (auto [d, a] : (*this))
    {
      d -= k;
      if (d >= 0)
        res.emplace_back(make_pair(d, a));
    }
    return res;
  }
  S operator>>=(ll k) { return (*this) = (*this) >> k; }

  F inv(int n) const
  {
    F f(n, T(0));
    f.front() = T(1);
    return f / (*this);
  }

  S differentiate()
  {
    for (auto &[d, a] : (*this))
      a *= d--;
    if (!(*this).empty() && (*this).front().first == -1)
      (*this).erase((*this).begin());
    return (*this);
  }
  S differential() const { return S(*this).differentiate(); }
  S integrate()
  {
    for (auto &[d, a] : (*this))
      a /= T(++d);
    return (*this);
  }
  S integral() const { return S(*this).integrate(); }

  F log(int n) const
  {
    F f = (*this).to_fps(n);
    return (f.differential() / (*this)).integral();
  }

  // 微分方程式 a(x)F'(x) + b(x)F(x) = 0, [x^0]F(x) = 1 を満たす F を n 項まで求める
  // [x^0]a(x) = 1 である必要がある
  F diffeq(const S &a, const S &b, int n) const
  {
    assert(a.front().first == 0 && a.front().second == 1);
    vector<T> minv(n);
    minv[1] = T(1);
    for (int i = 2; i < n; i++)
      minv[i] = -minv[T::mod() % i] * (T::mod() / i);
    F f(n, T(0));
    f[0] = T(1);
    for (int k = 0; k < n - 1; k++)
    {
      for (auto [i, ai] : a)
      {
        if (0 <= k - i + 1 && k - i + 1 < k + 1)
          f[k + 1] -= ai * (k - i + 1) * f[k - i + 1];
      }
      for (auto [j, bj] : b)
      {
        if (0 <= k - j && k - j < k + 1)
          f[k + 1] -= bj * f[k - j];
      }
      f[k + 1] *= minv[k + 1];
    }
    return f;
  }

  F exp(int n) const
  {
    return diffeq(S{{0, 1}}, -((*this).differential()), n);
  }
  
  // m >= 0 のときは O(nk)  (k: sparse の non-zero の個数)
  // m < 0 のときは O((n + d_0 m)k)
  F pow(ll m, int n) const
  {
    S f(*this);
    if (f.empty())
    {
      F res(n, T(0));
      if (m == 0)
        res.front() = T(1);
      return res;
    }
    auto [d0, a0] = f.front();
    T a0_inv = a0.inv();
    for (auto &[d, a] : f)
      d -= d0, a *= a0_inv;
    if (m >= 0)
    {
      F g = diffeq(f, -m * f.differential(), n);
      return (g * a0.pow(m)) << mul_limited(d0, m);
    }
    else
    {
      F g = diffeq(f, -m * f.differential(), n + (d0 * (-m)));
      F h = (g * a0_inv.pow(-m)) >> (d0 * (-m));
      h.resize(n);
      return h;
    }
  }

  // 素数 mod を要求
  // 存在しないなら空配列を返す
  F sqrt(int n) const
  {
    S f(*this);
    if (f.empty())
      return F(n, T(0));
    auto [d0, a0] = f.front();
    if (d0 % 2 != 0)
      return F(0);
    if (d0 >= n)
      return F(n, T(0));
    int p = T::mod();
    int r = sqrt_mod(a0.val(), p);
    if (r == -1)
      return F(0);

    T a0_inv = a0.inv();
    T inv_2 = T(2).inv();
    for (auto &[d, a] : f)
      d -= d0, a *= a0_inv;
    F g = diffeq(f, -inv_2 * f.differential(), n);
    return ((g * r) << (d0 / 2));
  }
};
template <class T, bool is_ntt_friendly>
vector<T> FormalPowerSeries<T, is_ntt_friendly>::fac{1, 1};
template <class T, bool is_ntt_friendly>
vector<T> FormalPowerSeries<T, is_ntt_friendly>::finv{1, 1};
template<class T, bool is_ntt_friendly>
vector<T> FormalPowerSeries<T, is_ntt_friendly>::invmint{0, 1};

template<class T, bool is_ntt_friendly>
struct RationalFormalPowerSeries
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;
  using R = RationalFormalPowerSeries;

  F num, den;

  R operator-() const
  {
    R res(*this);
    res.num = -res.num;
    return res;
  }

  R operator*=(const T &k)
  {
    (*this).num *= k;
    return *this;
  }
  R operator*(const T &k) const { return R(*this) *= k; }
  friend R operator*(const T k, const R &r) { return r * k; }
  R operator/=(const T &k)
  {
    (*this).den *= k;
    return k;
  }
  R operator/(const T &k) const { return R(*this) /= k; }

  R &operator+=(const R &r)
  {
    // ここうまくやると FFT の回数が節約できる気がする
    // うまくやらないと次数に偏りがある場合にかえって遅くなったりしそうで面倒
    F f, g;
    f = f.convolution2((*this).num, r.den);
    g = g.convolution2((*this).den, r.num);
    (*this).num = f + g;
    (*this).den = (*this).den.convolution2((*this).den, r.den);
    return *this;
  }
  R operator+(const R &r) const { return R(*this) += r; }
  R &operator-=(const R &r)
  {
    F f, g;
    f = f.convolution2((*this).num, r.den);
    g = g.convolution2((*this).den, r.num);
    (*this).num = f - g;
    (*this).den = (*this).den.convolution2((*this).den, r.den);
    return *this;
  }
  R operator-(const R &r) const { return R(*this) -= r; }
  
  R operator*=(const R &r)
  {
    (*this).num = (*this).num.convolution2((*this).num, r.num);
    (*this).den = (*this).den.convolution2((*this).den, r.den);
    return *this;
  }
  R operator*(const R &r) const { return R(*this) *= r; }
  R operator/=(const R &r)
  {
    (*this).num = (*this).num.convolution2((*this).num, r.den);
    (*this).den = (*this).den.convolution2((*this).den, r.num);
    return *this;
  }
  R operator/(const R &r) const { return R(*this) /= r; }

  R inv()
  {
    R res(*this);
    swap(res.num, res.den);
    return res;
  }
};

template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> convolution_many(const vector<FormalPowerSeries<T, is_ntt_friendly>> &fs, int d = -1)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  if (fs.empty())
  {
    if (d == -1)
      d = 0;
    F res(d + 1, T(0));
    res.front() = T(1);
    return res;
  }
  deque<F> deq;
  for (auto f : fs)
    deq.push_back(f);
  while ((int)deq.size() > 1)
  {
    F f = deq.front();
    deq.pop_front();
    F g = deq.front();
    deq.pop_front();
    f = f.convolution2(f, g, d);
    deq.push_back(f);
  }
  if (d != -1)
    deq.front().resize(d);
  return deq.front();
}

template <class T, bool is_ntt_friendly>
RationalFormalPowerSeries<T, is_ntt_friendly> rational_sum(const vector<RationalFormalPowerSeries<T, is_ntt_friendly>> &rs, int d = -1)
{
  using R = RationalFormalPowerSeries<T, is_ntt_friendly>;

  if (rs.empty())
    return R{{1}, {1}};

  deque<R> deq;
  for (auto &r : rs)
    deq.emplace_back(r);
  while ((int)deq.size() > 1)
  {
    R r1 = deq.front();
    deq.pop_front();
    R r2 = deq.front();
    deq.pop_front();
    R r3 = r1 + r2;
    if (d != -1)
    {
      if ((int)r3.num.size() > d)
        r3.num.resize(d);
      if ((int)r3.den.size() > d)
        r3.den.resize(d);
    }
    deq.emplace_back(r3);
  }
  if (d != -1)
    deq.front().num.resize(d), deq.front().den.resize(d);
  return deq.front();
}

template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> interpolation(const vector<T> &xs, const vector<T> &ys)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;
  using R = RationalFormalPowerSeries<T, is_ntt_friendly>;

  int n = xs.size();
  assert(n == ys.size());
  vector<F> fs(n);
  for (int i = 0; i < n; i++)
    fs[i] = F{-xs[i], T(1)};
  F g = convolution_many(fs);
  F h = g.differential();
  vector<T> a = h.eval_multipoint(xs);

  vector<R> rs(n);
  for (int i = 0; i < n; i++)
    rs[i] = R{F{ys[i] / a[i]}, fs[i]};
  R q = rational_sum(rs, n);
  return q.num;
}

// f(i) = ys[i] で定まる多項式 f(x) について f(c), …, f(c + M - 1) を求める
template<class T, bool is_ntt_friendly>
vector<T> sample_points_shift(const vector<T> &ys, int M, T c)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;
  F f;

  int N = ys.size();
  vector<T> a;
  {
    vector<T> p(N), q(N);
    for (int i = 0; i < N; i++)
    {
      p[i] = ys[i] * f.get_finv(i);
      q[i] = i % 2 == 0 ? f.get_finv(i) : -f.get_finv(i);
    }
    a = f.convolution2(p, q);
    a.resize(N);
  }
  vector<T> b;
  {
    vector<T> p(N), q(N);
    T tmp = 1;
    for (int i = 0; i < N; i++)
    {
      p[i] = a[i] * f.get_fac(i);
      q[i] = tmp * f.get_finv(i);
      tmp *= c - i;
    }
    reverse(q.begin(), q.end());
    b = f.convolution2(p, q);
    b.erase(b.begin(), b.begin() + N - 1);
    for (int i = 0; i < N; i++)
      b[i] *= f.get_finv(i);
  }
  vector<T> res;
  {
    vector<T> p(M);
    for (int i = 0; i < M; i++)
      p[i] = f.get_finv(i);
    res = f.convolution2(b, p);
    res.resize(M);
    for (int i = 0; i < M; i++)
      res[i] *= f.get_fac(i);
  }
  return res;
}

// https://suisen-kyopro.hatenablog.com/entry/2023/11/22/201600
// 前計算 O(K 2^K + (P/2^K) log K), クエリ O(2^K)
template<class T, bool is_ntt_friendly>
struct FactorialFast
{
private:
  const int P, K;
  vector<T> Y, Z, fac;

public:
  FactorialFast(const int K = 9) : P(T::mod()), K(K)
  {
    Y = {1};
    for (int i = 0; i < K; i++)
    {
      Z = sample_points_shift<T, is_ntt_friendly>(Y, (1 << (i + 2)) - (1 << i), 1 << i);
      Z.insert(Z.begin(), Y.begin(), Y.end());
      Y.resize(1 << (i + 1));
      for (int j = 0; j < (1 << (i + 1)); j++)
        Y[j] = Z[2 * j] * Z[2 * j + 1] * T::raw((1 << i) * (2 * j + 1));
    }
    if ((1 << K) <= P / (1 << K))
    {
      Z = sample_points_shift<T, is_ntt_friendly>(Y, P / (1 << K), 1 << K);
      Y.insert(Y.end(), Z.begin(), Z.end());
    }
    fac.resize(P / (1 << K) + 1);
    fac.at(0) = 1;
    for (int i = 0; i < P / (1 << K); i++)
      fac[i + 1] = fac[i] * Y[i] * T::raw((1 + i) * (1 << K));
  }
  T query(ll n)
  {
    if (n >= T::mod())
      return 0;
    T res = fac.at(n / (1 << K));
    for (int j = n / (1 << K) * (1 << K) + 1; j <= n; j++)
      res *= T::raw(j);
    return res;
  }
};

// f(x)/prod[i](1-a[i]x) = sum[i] c[i]/(1-a[i]x) なる c を求める
template<class T, bool is_ntt_friendly>
vector<T> partial_fraction_decomposition(const FormalPowerSeries<T, is_ntt_friendly> &f, const vector<T> &as)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  int N = as.size();
  assert((int)f.size() <= N);
  vector<T> bs(N);
  for (int i = 0; i < N; i++)
    bs[i] = as[i].inv();
  vector<F> gs(N);
  for (int i = 0; i < N; i++)
    gs[i] = F{T(1), -as[i]};
  F g = convolution_many(gs);
  F dg = g.differential();
  vector<T> ys = f.eval_multipoint(bs), zs = dg.eval_multipoint(bs);
  vector<T> cs(N);
  for (int i = 0; i < N; i++)
    cs[i] = -as[i] * ys[i] / zs[i];
  return cs;
}

// n = 0, 1, …, N-1 に対して n^k を列挙
template<class T>
vector<T> enum_pow(int N, int k)
{
  vector<int> minfactor(N, -1);
  for (int i = 2; i < N; i++)
  {
    if (minfactor[i] != -1)
      continue;
    for (int j = 2 * i; j < N; j += i)
      minfactor[j] = i;
  }
  vector<T> power(N);
  for (int i = 0; i < N; i++)
  {
    if (minfactor[i] == -1)
      power[i] = T(i).pow(k);
    else
      power[i] = power[minfactor[i]] * power[i / minfactor[i]];
  }
  return power;
}

// sum_[i = 0..infty] r^i i^d
template <class T, bool is_ntt_friendly>
T sum_of_exp_times_poly_limit(T r, int d)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  vector<T> pws = enum_pow<T>(d + 2, d);
  vector<T> pwr(d + 2, 1);
  for (int i = 0; i < d + 1; i++)
    pwr[i + 1] = pwr[i] * r;
  F f(d + 2), g(d + 2);
  f.front() = pws.front() * pwr.front();
  for (int i = 0; i < d + 1; i++)
    f[i + 1] = f[i] + pws[i + 1] * pwr[i + 1];
  for (int i = 0; i <= d + 1; i++)
    g[i] = (i % 2 == 0 ? 1 : -1) * pwr[i] * g.get_fac(d + 1) * g.get_finv(i) * g.get_finv(d + 1 - i);
  T c = 0;
  for (int i = 0; i <= d + 1; i++)
    c += f[i] * g[d + 1 - i];
  c /= accumulate(g.begin(), g.end(), T(0));
  return c;
}

// sum_[i = 0..n-1] r^i i^d
template <class T, bool is_ntt_friendly>
T sum_of_exp_times_poly(T r, int d, ll n)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  if (n == 0)
    return 0;
  if (r == 0)
    return d == 0 ? 1 : 0;

  vector<T> pws = enum_pow<T>(d + 2, d);
  vector<T> pwr(d + 2, 1);
  for (int i = 0; i < d + 1; i++)
    pwr[i + 1] = pwr[i] * r;
  F f(d + 2), g(d + 2);
  f.front() = pws.front() * pwr.front();
  for (int i = 0; i < d + 1; i++)
    f[i + 1] = f[i] + pws[i + 1] * pwr[i + 1];
  if (r == 1)
  {
    LagrangeInterpolation<T> lag(f);
    return lag.eval(n - 1);
  }
  for (int i = 0; i <= d + 1; i++)
    g[i] = (i % 2 == 0 ? 1 : -1) * pwr[i] * g.get_fac(d + 1) * g.get_finv(i) * g.get_finv(d + 1 - i);
  T c = 0;
  for (int i = 0; i <= d + 1; i++)
    c += f[i] * g[d + 1 - i];
  c /= accumulate(g.begin(), g.end(), T(0));
  F h = f - F(d + 2, c);
  {
    T rinv = r.inv();
    T pwrinv = 1;
    for (int i = 0; i <= d + 1; i++)
    {
      h[i] *= pwrinv;
      pwrinv *= rinv;
    }
  }
  LagrangeInterpolation<T> lag(h);
  return c + r.pow(n - 1) * lag.eval(n - 1);
}

// prod[d in D](1 + cx^d) を M 次の項まで求める
template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> multiply_many(const int &M, const T &c, const vector<int> &D)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  vector<int> cnt(M + 1, 0);
  for (auto d : D)
  {
    if (d < 0 || M < d)
      continue;
    cnt[d]++;
  }

  F f(M + 1, 0);
  for (int k = 1; k <= M; k++)
  {
    T pw = 1;
    for (int i = 1; k * i <= M; i++)
    {
      pw *= c;
      if (i & 1)
        f[k * i] += T::raw(cnt[k]) * pw * f.get_invmint(i);
      else
        f[k * i] -= T::raw(cnt[k]) * pw * f.get_invmint(i);
    }
  }
  return f.exp();
}

// 多重集合 S の要素から何個か選んで総和を 0, 1, …, M にする方法の数
template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> subset_sum(const int &M, const vector<int> &S)
{
  return multiply_many<T, is_ntt_friendly>(M, T(1), S);
}

// 集合 S の各要素が無限個ある集合 T から何個か選んで総和を 0, 1, …, M にする方法の数
template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> partition(const int &M, const vector<int> &S)
{
  return multiply_many<T, is_ntt_friendly>(M, T(-1), S).inv();
}

template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> stirling1(const int &N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;
  using S = SparseFormalPowerSeries<T, is_ntt_friendly>;

  if (N == 0)
    return {1};
  if (N == 1)
    return {0, 1};
  if (N & 1)
  {
    F f = stirling1<T, is_ntt_friendly>(N - 1);
    f.resize(N + 1, T(0));
    return f * S{{0, 1 - N}, {1, 1}};
  }
  else
  {
    F f = stirling1<T, is_ntt_friendly>(N / 2);
    f.resize(N + 1, T(0));
    F g = f.taylor_shift(-(N / 2));
    return f * g;
  }
}

template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> stirling2(const int &N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  vector<T> power = enum_pow<T>(N + 1, N);
  F A(N + 1), B(N + 1);
  for (int i = 0; i <= N; i++)
  {
    A[i] = power[i] * A.get_finv(i);
    B[i] = (i & 1) ? -A.get_finv(i) : A.get_finv(i);
  }
  return A * B;
}

template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> bernoulli_number(const int &N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  F f(N + 1, T(0));
  for (int i = 0; i <= N; i++)
    f[i] = f.get_finv(i + 1);
  return f.inv().to_ogf();
}

// [x^N] P(x)/Q(x) を求める（P の次数は Q の次数より小さい）
template<class T, bool is_ntt_friendly>
T bostan_mori(const FormalPowerSeries<T, is_ntt_friendly> &P, const FormalPowerSeries<T, is_ntt_friendly> &Q, ll N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  int d = (int)Q.size() - 1;
  assert((int)P.size() <= d);
  if (is_ntt_friendly)
  {
    int z = 1;
    while (z < 2 * d + 1)
      z <<= 1;
    T iz = T(z).inv();
    F U = F(P), V = F(Q);
    U.resize(z), V.resize(z);
    while (N > 0)
    {
      U.butterfly2(U), V.butterfly2(V);
      for (int i = 0; i < z; i += 2)
      {
        T x = V[i + 1], y = V[i];
        U[i] *= x, V[i] *= x;
        U[i + 1] *= y, V[i + 1] *= y;
      }
      U.butterfly_inv2(U), V.butterfly_inv2(V);
      for (int i = 0; i < (z >> 1); i++)
      {
        U[i] = U[2 * i + (N & 1)] * iz;
        V[i] = V[2 * i] * iz;
      }
      for (int i = (z >> 1); i < z; i++)
        U[i] = 0, V[i] = 0;
      N >>= 1;
    }
    return U.front() / V.front();
  }
  else
  {
    F U = F(P), V = F(Q);
    U.resize(d), V.resize(d + 1);
    while (N > 0)
    {
      F U2 = F(U), V2 = F(V), V3 = F(V);
      for (int i = 1; i <= d; i += 2)
        V3[i] = -V3[i];
      U2 *= V3, V2 *= V3;
      for (int i = 0; i <= d; i++)
      {
        U[i] = U2[2 * i + (N & 1)];
        V[i] = V2[2 * i];
      }
      N >>= 1;
    }
    return U.front() / V.front();
  }
}
// a_n = sum[i = 1..d] c_i a_{n-i}（n ≥ d）を満たすとき、a_N を求める（A は 0-indexed で C は 1-indexed）
template<class T, bool is_ntt_friendly>
T linear_recurrence(const vector<T> &A, const vector<T> &C, ll N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  int d = C.size();
  assert((int)A.size() >= d);

  F Ga(d), Q(d + 1);
  Q[0] = 1;
  for (int i = 0; i < d; i++)
    Ga[i] = A[i], Q[i + 1] = -C[i];
  F P = Ga * Q;
  return bostan_mori(P, Q, N);
}

// (P の次数) < (Q の次数) とする
// P/Q = R + x^N (P'/Q) を満たす P' (R は N 次未満、P' は d 次未満)
// [x^{N+n}](P/Q) = [x^n](P'/Q)　線形漸化的数列のシフト
// 高速化の余地あり
template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> fiduccia(const FormalPowerSeries<T, is_ntt_friendly> &P, const FormalPowerSeries<T, is_ntt_friendly> &Q, ll N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  assert(P.size() < Q.size());

  F xinv = -(Q >> 1);
  if (Q[0] != 1)
    xinv /= Q[0];
  return xinv.powmod(N, Q) * P % Q;
}
// a_n = sum[i = 1..d] c_i a_{n-i}（n ≥ d）を満たすとき、a_N, …, a_{N+d-1} を求める（A は 0-indexed で C は 1-indexed）
template<class T, bool is_ntt_friendly>
vector<T> linear_recurrence_many(const vector<T> &A, const vector<T> &C, ll N)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;

  int d = C.size();
  assert((int)A.size() >= d);

  F Ga(d), Q(d + 1);
  Q[0] = 1;
  for (int i = 0; i < d; i++)
    Ga[i] = A[i], Q[i + 1] = -C[i];
  F P = Ga * Q;
  F P2 = fiduccia(P, Q, N);
  P2.resize(d);
  F Gb = P2 / Q + (P.div_poly(Q) >> N);
  return Gb;
}

// Σ[i = 0..M-1] a_i exp(b_i x) を N 項まで求める
template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> exp_sum(int N, const vector<T> &a, const vector<T> &b)
{
  using F = FormalPowerSeries<T, is_ntt_friendly>;
  using R = RationalFormalPowerSeries<T, is_ntt_friendly>;
  assert(a.size() == b.size());
  int M = a.size();
  vector<R> gs(M);
  for (int i = 0; i < M; i++)
    gs[i] = R{F{a[i]}, F{1, -b[i]}};
  R g = rational_sum(gs, N);
  return (g.num / g.den).to_egf();
}

// f(exp(kx)) を N 項まで求める 
template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> eval_exp(FormalPowerSeries<T, is_ntt_friendly> &f, T k, int N = -1)
{
  if (N == -1)
    N = (int)f.size();
  vector<T> b(f.size());
  for (int i = 0; i < (int)f.size(); i++)
    b[i] = k * i;
  return exp_sum<T, is_ntt_friendly>(N, f, b);
}

/*
using mint = modint998244353;
const bool ntt = true;
//*/
//*
using mint = modint1000000007;
const bool ntt = false;
//*/
/*
using mint = modint;
const bool ntt = false;
//*/

using fps = FormalPowerSeries<mint, ntt>;
using sfps = SparseFormalPowerSeries<mint, ntt>;
using rfps = RationalFormalPowerSeries<mint, ntt>;

int main()
{
  ll N;
  cin >> N;
  FactorialFast<mint, ntt> fac(13);
  cout << fac.query(N).val() << endl;
}
0