#include <iostream>
#include <string>
#include <vector>
#include <array>
#include <queue>
#include <deque>
#include <algorithm>
#include <set>
#include <map>
#include <bitset>
#include <cmath>
#include <functional>
#include <cassert>
#include <iomanip>
#define vll vector<ll>
#define vvvl vector<vvl>
#define vvl vector<vector<ll>>
#define VV(a, b, c, d) vector<vector<d>>(a, vector<d>(b, c))
#define VVV(a, b, c, d) vector<vvl>(a, vvl(b, vll (c, d)));
#define re(c, b) for(ll c=0;c<b;c++)
#define all(obj) (obj).begin(), (obj).end()
typedef long long int ll;
typedef long double ld;
using namespace std;

#include <numeric>
#include <type_traits>

namespace internal{
  template <class T> using is_integral = typename std::is_integral<T>;
  template <class T>
  using is_signed_int =
  typename std::conditional<is_integral<T>::value && std::is_signed<T>::value, std::true_type, std::false_type>::type;
  template <class T>
  using is_unsigned_int =
  typename std::conditional<is_integral<T>::value && std::is_unsigned<T>::value,
  std::true_type, std::false_type>::type;
  template <class T>
  using to_unsigned = typename std::conditional<is_signed_int<T>::value,
  std::make_unsigned<T>, std::common_type<T>>::type;
  template <class T>
  using is_signed_int_t = std::enable_if_t<is_signed_int<T>::value>;
  template <class T>
  using is_unsigned_int_t = std::enable_if_t<is_unsigned_int<T>::value>;
  template <class T> using to_unsigned_t = typename to_unsigned<T>::type;
  //<internal_math>
  constexpr long long safe_mod(long long x, long long m) {
    x %= m;
    if (x < 0) x += m;
    return x;
  }
  struct barrett {
    unsigned int _m;
    unsigned long long im;
    barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {}
    unsigned int umod() const { return _m; }
    unsigned int mul(unsigned int a, unsigned int b) const {
      unsigned long long z = a;
      z *= b;
      unsigned long long x = (unsigned long long)(((unsigned __int128)(z)*im) >> 64);
      unsigned int v = (unsigned int)(z - x * _m);
      if (_m <= v) v += _m;
      return v;
    }
  };
  constexpr long long pow_mod_constexpr(long long x, long long n, int m) {
    if (m == 1) return 0;
    unsigned int _m = (unsigned int)(m);
    unsigned long long r = 1;
    unsigned long long y = safe_mod(x, m);
    while (n) {
      if (n & 1) r = (r * y) % _m;
      y = (y * y) % _m;
      n >>= 1;
    }
    return r;
  }

  constexpr bool is_prime_constexpr(int n) {
    if (n <= 1) return false;
    if (n == 2 || n == 7 || n == 61) return true;
    if (n % 2 == 0) return false;
    long long d = n - 1;
    while (d % 2 == 0) d /= 2;
    constexpr long long bases[3] = {2, 7, 61};
    for (long long a : bases) {
      long long t = d;
      long long y = pow_mod_constexpr(a, t, n);
      while (t != n - 1 && y != 1 && y != n - 1) {
        y = y * y % n;
        t <<= 1;
      }
      if (y != n - 1 && t % 2 == 0) {
        return false;
      }
    }
    return true;
  }
  template <int n> constexpr bool is_prime = is_prime_constexpr(n);

  constexpr std::pair<long long, long long> inv_gcd(long long a, long long b) {
    a = safe_mod(a, b);
    if (a == 0) return {b, 0};
    long long s = b, t = a;
    long long m0 = 0, m1 = 1;

    while (t) {
      long long u = s / t;
      s -= t * u;
      m0 -= m1 * u;
      auto tmp = s;
      s = t;
      t = tmp;
      tmp = m0;
      m0 = m1;
      m1 = tmp;
    }
    if (m0 < 0) m0 += b / s;
    return {s, m0};
  }
  constexpr int primitive_root_constexpr(int m) {
    if (m == 2) return 1;
    if (m == 167772161) return 3;
    if (m == 469762049) return 3;
    if (m == 754974721) return 11;
    if (m == 998244353) return 3;
    int divs[20] = {};
    divs[0] = 2;
    int cnt = 1;
    int x = (m - 1) / 2;
    while (x % 2 == 0) x /= 2;
    for (int i = 3; (long long)(i)*i <= x; i += 2) {
      if (x % i == 0) {
        divs[cnt++] = i;
        while (x % i == 0) {
          x /= i;
        }
      }
    }
    if (x > 1) divs[cnt++] = x;
    for (int g = 2;; g++) {
      bool ok = true;
      for (int i = 0; i < cnt; i++) {
        if (pow_mod_constexpr(g, (m - 1) / divs[i], m) == 1) {
          ok = false;
          break;
        }
      }
      if (ok) return g;
    }
  }
  template <int m> constexpr int primitive_root = primitive_root_constexpr(m);
  //<internal_bit>
  int ceil_pow2(int n) {
    int x = 0;
    while ((1U << x) < (unsigned int)(n)) x++;
    return x;
  }
  int bsf(unsigned int n) {
    return __builtin_ctz(n);
  }

  //<modint>
  struct modint_base {};
  struct static_modint_base : modint_base {};
  template <class T> using is_modint = std::is_base_of<modint_base, T>;
  template <class T> using is_modint_t = std::enable_if_t<is_modint<T>::value>;
  template <int m, std::enable_if_t<(1 <= m)>* = nullptr>
  struct static_modint : internal::static_modint_base {
    using mint = static_modint;
  public:
    static constexpr int mod() { return m; }
    static mint raw(int v) {
      mint x;
      x._v = v;
      return x;
    }
    static_modint() : _v(0) {}
    template <class T, internal::is_signed_int_t<T>* = nullptr>
    static_modint(T v) {
      long long x = (long long)(v % (long long)(umod()));
      if (x < 0) x += umod();
      _v = (unsigned int)(x);
    }
    template <class T, internal::is_unsigned_int_t<T>* = nullptr>
    static_modint(T v) {
      _v = (unsigned int)(v % umod());
    }
    static_modint(bool v) { _v = ((unsigned int)(v) % umod()); }

    unsigned int val() const { return _v; }

    mint& operator++() {
      _v++;
      if (_v == umod()) _v = 0;
      return *this;
    }
    mint& operator--() {
      if (_v == 0) _v = umod();
      _v--;
      return *this;
    }
    mint operator++(int) {
      mint result = *this;
      ++*this;
      return result;
    }
    mint operator--(int) {
      mint result = *this;
      --*this;
      return result;
    }

    mint& operator+=(const mint& rhs) {
      _v += rhs._v;
      if (_v >= umod()) _v -= umod();
      return *this;
    }
    mint& operator-=(const mint& rhs) {
      _v -= rhs._v;
      if (_v >= umod()) _v += umod();
      return *this;
    }
    mint& operator*=(const mint& rhs) {
      unsigned long long z = _v;
      z *= rhs._v;
      _v = (unsigned int)(z % umod());
      return *this;
    }
    mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); }

    mint operator+() const { return *this; }
    mint operator-() const { return mint() - *this; }

    mint pow(long long n) const {
      assert(0 <= n);
      mint x = *this, r = 1;
      while (n) {
        if (n & 1) r *= x;
        x *= x;
        n >>= 1;
      }
      return r;
    }
    mint inv() const {
      if (prime) {
        assert(_v);
        return pow(umod() - 2);
      } else {
        auto eg = internal::inv_gcd(_v, m);
        assert(eg.first == 1);
        return eg.second;
      }
    }

    friend mint operator+(const mint& lhs, const mint& rhs) {
      return mint(lhs) += rhs;
    }
    friend mint operator-(const mint& lhs, const mint& rhs) {
      return mint(lhs) -= rhs;
    }
    friend mint operator*(const mint& lhs, const mint& rhs) {
      return mint(lhs) *= rhs;
    }
    friend mint operator/(const mint& lhs, const mint& rhs) {
      return mint(lhs) /= rhs;
    }
    friend bool operator==(const mint& lhs, const mint& rhs) {
      return lhs._v == rhs._v;
    }
    friend bool operator!=(const mint& lhs, const mint& rhs) {
      return lhs._v != rhs._v;
    }

  private:
    unsigned int _v;
    static constexpr unsigned int umod() { return m; }
    static constexpr bool prime = internal::is_prime<m>;
  };

  template <int id> struct dynamic_modint : internal::modint_base {
    using mint = dynamic_modint;
  public:
    static int mod() { return (int)(bt.umod()); }
    static void set_mod(int m) {
      assert(1 <= m);
      bt = internal::barrett(m);
    }
    static mint raw(int v) {
      mint x;
      x._v = v;
      return x;
    }
    dynamic_modint() : _v(0) {}
    template <class T, internal::is_signed_int_t<T>* = nullptr>
    dynamic_modint(T v) {
      long long x = (long long)(v % (long long)(mod()));
      if (x < 0) x += mod();
      _v = (unsigned int)(x);
    }
    template <class T, internal::is_unsigned_int_t<T>* = nullptr>
    dynamic_modint(T v) {
      _v = (unsigned int)(v % mod());
    }
    dynamic_modint(bool v) { _v = ((unsigned int)(v) % mod()); }

    unsigned int val() const { return _v; }

    mint& operator++() {
      _v++;
      if (_v == umod()) _v = 0;
      return *this;
    }
    mint& operator--() {
      if (_v == 0) _v = umod();
      _v--;
      return *this;
    }
    mint operator++(int) {
      mint result = *this;
      ++*this;
      return result;
    }
    mint operator--(int) {
      mint result = *this;
      --*this;
      return result;
    }

    mint& operator+=(const mint& rhs) {
      _v += rhs._v;
      if (_v >= umod()) _v -= umod();
      return *this;
    }
    mint& operator-=(const mint& rhs) {
      _v += mod() - rhs._v;
      if (_v >= umod()) _v -= umod();
      return *this;
    }
    mint& operator*=(const mint& rhs) {
      _v = bt.mul(_v, rhs._v);
      return *this;
    }
    mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); }

    mint operator+() const { return *this; }
    mint operator-() const { return mint() - *this; }

    mint pow(long long n) const {
      assert(0 <= n);
      mint x = *this, r = 1;
      while (n) {
        if (n & 1) r *= x;
        x *= x;
        n >>= 1;
      }
      return r;
    }
    mint inv() const {
      auto eg = internal::inv_gcd(_v, mod());
      assert(eg.first == 1);
      return eg.second;
    }

    friend mint operator+(const mint& lhs, const mint& rhs) {
      return mint(lhs) += rhs;
    }
    friend mint operator-(const mint& lhs, const mint& rhs) {
      return mint(lhs) -= rhs;
    }
    friend mint operator*(const mint& lhs, const mint& rhs) {
      return mint(lhs) *= rhs;
    }
    friend mint operator/(const mint& lhs, const mint& rhs) {
      return mint(lhs) /= rhs;
    }
    friend bool operator==(const mint& lhs, const mint& rhs) {
      return lhs._v == rhs._v;
    }
    friend bool operator!=(const mint& lhs, const mint& rhs) {
      return lhs._v != rhs._v;
    }

  private:
    unsigned int _v;
    static internal::barrett bt;
    static unsigned int umod() { return bt.umod(); }
  };
  template <int id> internal::barrett dynamic_modint<id>::bt = 998244353;
  //using modint998244353 = static_modint<998244353>;
  //using modint1000000007 = static_modint<1000000007>;
  //using modint = dynamic_modint<-1>;
  template <class T>
  using is_static_modint = std::is_base_of<internal::static_modint_base, T>;

  template <class T>
  using is_static_modint_t = std::enable_if_t<is_static_modint<T>::value>;

  template <class> struct is_dynamic_modint : public std::false_type {};
  template <int id>
  struct is_dynamic_modint<dynamic_modint<id>> : public std::true_type {};

  template <class T>
  using is_dynamic_modint_t = std::enable_if_t<is_dynamic_modint<T>::value>;

  //<conbovution>
  template <class mint, internal::is_static_modint_t<mint>* = nullptr>
  void butterfly(std::vector<mint>& a) {
    static constexpr int g = internal::primitive_root<mint::mod()>;
    int n = int(a.size());
    int h = internal::ceil_pow2(n);
    static bool first = true;
    static mint sum_e[30];  // sum_e[i] = ies[0] * ... * ies[i - 1] * es[i]
    if (first) {
      first = false;
      mint es[30], ies[30];  // es[i]^(2^(2+i)) == 1
      int cnt2 = bsf(mint::mod() - 1);
      mint e = mint(g).pow((mint::mod() - 1) >> cnt2), ie = e.inv();
      for (int i = cnt2; i >= 2; i--) {
        // e^(2^i) == 1
        es[i - 2] = e;
        ies[i - 2] = ie;
        e *= e;
        ie *= ie;
      }
      mint now = 1;
      for (int i = 0; i <= cnt2 - 2; i++) {
        sum_e[i] = es[i] * now;
        now *= ies[i];
      }
    }
    for (int ph = 1; ph <= h; ph++) {
      int w = 1 << (ph - 1), p = 1 << (h - ph);
      mint now = 1;
      for (int s = 0; s < w; s++) {
        int offset = s << (h - ph + 1);
        for (int i = 0; i < p; i++) {
          auto l = a[i + offset];
          auto r = a[i + offset + p] * now;
          a[i + offset] = l + r;
          a[i + offset + p] = l - r;
        }
        now *= sum_e[internal::bsf(~(unsigned int)(s))];
      }
    }
  }

  template <class mint, internal::is_static_modint_t<mint>* = nullptr>
  void butterfly_inv(std::vector<mint>& a) {
    static constexpr int g = internal::primitive_root<mint::mod()>;
    int n = int(a.size());
    int h = internal::ceil_pow2(n);

    static bool first = true;
    static mint sum_ie[30];  // sum_ie[i] = es[0] * ... * es[i - 1] * ies[i]
    if (first) {
      first = false;
      mint es[30], ies[30];  // es[i]^(2^(2+i)) == 1
      int cnt2 = bsf(mint::mod() - 1);
      mint e = mint(g).pow((mint::mod() - 1) >> cnt2), ie = e.inv();
      for (int i = cnt2; i >= 2; i--) {
        // e^(2^i) == 1
        es[i - 2] = e;
        ies[i - 2] = ie;
        e *= e;
        ie *= ie;
      }
      mint now = 1;
      for (int i = 0; i <= cnt2 - 2; i++) {
        sum_ie[i] = ies[i] * now;
        now *= es[i];
      }
    }

    for (int ph = h; ph >= 1; ph--) {
      int w = 1 << (ph - 1), p = 1 << (h - ph);
      mint inow = 1;
      for (int s = 0; s < w; s++) {
        int offset = s << (h - ph + 1);
        for (int i = 0; i < p; i++) {
          auto l = a[i + offset];
          auto r = a[i + offset + p];
          a[i + offset] = l + r;
          a[i + offset + p] =
            (unsigned long long)(mint::mod() + l.val() - r.val()) *
            inow.val();
        }
        inow *= sum_ie[internal::bsf(~(unsigned int)(s))];
      }
    }
  }

  template <class mint, internal::is_static_modint_t<mint>* = nullptr>
  std::vector<mint> _convolution(std::vector<mint> a, std::vector<mint> b) {
    int n = int(a.size()), m = int(b.size());
    if (!n || !m) return {};
    if (std::min(n, m) <= 60) {
      if (n < m) {
        std::swap(n, m);
        std::swap(a, b);
      }
      std::vector<mint> ans(n + m - 1);
      for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
          ans[i + j] += a[i] * b[j];
        }
      }
      return ans;
    }
    int z = 1 << internal::ceil_pow2(n + m - 1);
    a.resize(z);
    internal::butterfly(a);
    b.resize(z);
    internal::butterfly(b);
    for (int i = 0; i < z; i++) {
      a[i] *= b[i];
    }
    internal::butterfly_inv(a);
    a.resize(n + m - 1);
    mint iz = mint(z).inv();
    for (int i = 0; i < n + m - 1; i++) a[i] *= iz;
    return a;
  }

  template <unsigned int mod = 998244353, class T,
  std::enable_if_t<internal::is_integral<T>::value>* = nullptr>
  std::vector<T> _convolution(const std::vector<T>& a, const std::vector<T>& b) {
    int n = int(a.size()), m = int(b.size());
    if (!n || !m) return {};

    using mint = internal::static_modint<mod>;
    std::vector<mint> a2(n), b2(m);
    for (int i = 0; i < n; i++) {
      a2[i] = mint(a[i]);
    }
    for (int i = 0; i < m; i++) {
      b2[i] = mint(b[i]);
    }
    auto c2 = _convolution(move(a2), move(b2));
    std::vector<T> c(n + m - 1);
    for (int i = 0; i < n + m - 1; i++) {
      c[i] = c2[i].val();
    }
    return c;
  }

  std::vector<long long> convolution_ll(const std::vector<long long>& a,
                      const std::vector<long long>& b) {
    int n = int(a.size()), m = int(b.size());
    if (!n || !m) return {};
    static constexpr unsigned long long MOD1 = 754974721;  // 2^24
    static constexpr unsigned long long MOD2 = 167772161;  // 2^25
    static constexpr unsigned long long MOD3 = 469762049;  // 2^26
    static constexpr unsigned long long M2M3 = MOD2 * MOD3;
    static constexpr unsigned long long M1M3 = MOD1 * MOD3;
    static constexpr unsigned long long M1M2 = MOD1 * MOD2;
    static constexpr unsigned long long M1M2M3 = MOD1 * MOD2 * MOD3;
    static constexpr unsigned long long i1 =
      internal::inv_gcd(MOD2 * MOD3, MOD1).second;
    static constexpr unsigned long long i2 =
      internal::inv_gcd(MOD1 * MOD3, MOD2).second;
    static constexpr unsigned long long i3 =
      internal::inv_gcd(MOD1 * MOD2, MOD3).second;
    auto c1 = internal::_convolution<MOD1>(a, b);
    auto c2 = internal::_convolution<MOD2>(a, b);
    auto c3 = internal::_convolution<MOD3>(a, b);
    std::vector<long long> c(n + m - 1);
    for (int i = 0; i < n + m - 1; i++) {
      unsigned long long x = 0;
      x += (c1[i] * i1) % MOD1 * M2M3;
      x += (c2[i] * i2) % MOD2 * M1M3;
      x += (c3[i] * i3) % MOD3 * M1M2;
      long long diff =
        c1[i] - internal::safe_mod((long long)(x), (long long)(MOD1));
      if (diff < 0) diff += MOD1;
      static constexpr unsigned long long offset[5] = {
        0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3};
      x -= offset[diff % 5];
      c[i] = x;
    }
    return c;
  }
}

ll mpow(ll a, ll b, ll MOD = -1){
  ll ret = 1, num = a;
  while(b>0){
    if(b%2) ret = (ret*num)%MOD;
    num = (num*num)%MOD;
    b /= 2;
  }
  return ret;
}
vector<ll> int32mod_conv(vector<ll> a, vector<ll> b, ll MOD, int MAXSIZE=-1){
  if(MAXSIZE!=-1){
    if(a.size()>MAXSIZE) a.resize(MAXSIZE);
    if(b.size()>MAXSIZE) b.resize(MAXSIZE);
  }
  if(MOD==998244353) return internal::_convolution<998244353, ll>(a, b);
  vector<ll> A = internal::_convolution<167772161>(a, b);
  vector<ll> B = internal::_convolution<469762049>(a, b);
  vector<ll> C = internal::_convolution<1224736769>(a, b);
  ll N = A.size();
  vector<ll> ret(N);
  ll x = 167772161, y = 469762049, z = 1224736769;
  ll ix = mpow(x, y-2, y);
  ll ixy = mpow((x*y)%z, z-2, z);
  for(int i=0;i<N;i++){
    ll v = ((B[i] - A[i])*ix)%y;
    if(v<0) v += y;
    ll xxv = A[i]+x*v;
    v = ((C[i] - (xxv%z))*ixy)%z;
    if(v<0) v += z;
    ret[i] = ((xxv%MOD) + ((x*y)%MOD)*v)%MOD;
  }
  if(MAXSIZE!=-1&&(int)ret.size()>MAXSIZE) ret.resize(MAXSIZE);
  return ret;
}


//using modint998244353 = internal::static_modint<998244353>;
//using modint1000000007 = internal::static_modint<1000000007>;
//using modint = internal::dynamic_modint<-1>;

template<ll MOD, int MAXSIZE=-1>
struct StaticModFPS: vector<ll>{
  using vector<ll>::vector;
  using fps = StaticModFPS<MOD, MAXSIZE>;
  StaticModFPS(vector<ll> v){
    int n = v.size();
    this->resize(n);
    for(int i=0;i<n;i++){
      (*this)[i] = v[i] % MOD;
      if((*this)[i] < 0) (*this)[i] += MOD;
    }
  }
  fps operator *= (const fps &vr) {
    *this = int32mod_conv(*this, vr, MOD);
    return *this;
  }
  fps operator /= (fps &vr){
    return (*this) *= vr.inv();
  }
  fps operator += (const fps &vr){
    int n = this->size();
    int m = vr.size();
    if(n < m) this->resize(m);
    for(int i=0;i<m;i++) {
      (*this)[i] += vr[i];
      if((*this)[i] >= MOD) (*this)[i] -= MOD;
    }
    return *this;
  }
  fps operator -= (const fps &vr){
    int n = this->size();
    int m = vr.size();
    if(n < m) this->resize(m);
    for(int i=0;i<m;i++) {
      (*this)[i] -= vr[i];
      if((*this)[i] < 0) (*this)[i] += MOD;
    }
    return *this;
  }
  fps operator += (const ll &vr){
    int n = this->size();
    ll r = vr % MOD;
    if(r < 0) r += MOD;
    for(int i=0;i<n;i++){
      (*this)[i] += r;
      if((*this)[i] >= MOD) (*this)[i] -= MOD;
    }
    return *this;
  }
  fps operator -= (const ll &vr){
    int n = this->size();
    ll r = vr % MOD;
    if(r<0) r += MOD;
    for(int i=0;i<n;i++){
      (*this[i]) -= r;
      if((*this[i]) < 0) (*this)[i] += MOD;
    }
    return *this;
  }
  fps operator *= (const ll &vr){
    int n = this->size();
    ll r = vr % MOD;
    if(r<0) r += MOD;
    for(int i=0;i<n;i++){
      (*this)[i] = ((*this)[i] * r)%MOD;
    }
    return *this;
  }
  fps operator /= (const ll &vr){
    ll r = vr % MOD;
    if(r < 0) r += MOD;
    assert(r!=0);
    r = mpow(r, MOD-2, MOD);
    int n = (int)this->size();
    for(int i=0;i<n;i++) (*this)[i] = ((*this)[i] * r)%MOD;
    return *this;
  }
  fps operator + (const fps& vr){return fps(*this) += vr;}
  fps operator - (const fps& vr){return fps(*this) -= vr;}
  fps operator * (const fps& vr){return fps(*this) *= vr;}
  fps operator / (const fps& vr){return fps(*this) /= vr;}
  fps operator + (const ll& vr){return fps(*this) += vr;}
  fps operator - (const ll& vr){return fps(*this) -= vr;}
  fps operator * (const ll& vr){return fps(*this) *= vr;}
  fps operator / (const ll& vr){return fps(*this) /= vr;}

  void debug(int printsize = 20){
    int n = min(20, (int)this->size());
    for(int i=0;i<n;i++){
      if(i==n-1) printf("%lld\n", (*this)[i]);
      else printf("%lld ", (*this)[i]);
    }
  }
  void print(){
    int n = (int)this->size();
    for(int i=0;i<n;i++){
      if(i==n-1) printf("%lld\n", (*this)[i]);
      else printf("%lld ", (*this)[i]);
    }
  }

  fps rev(int deg=-1) const {
   fps ret(*this);
   if(deg != -1) ret.resize(deg, 0);
   reverse(ret.begin(), ret.end());
   return ret;
 }

  fps prefix(int deg){
    int n = min((int)this->size(), deg);
    return fps(this->begin(), this->begin() + n);
  }
  // https://en.wikipedia.org/wiki/Formal_power_series#Multiplicative_inverse
  // invertible in R[X] if and only if its constant coefficient a_{0} is invertible in R
  fps inv(int deg=-1){
    assert((*this)[0]);
    int n = this->size();
    if(deg==-1) deg = n;
    fps ret({mpow((*this)[0], MOD-2, MOD)});
    for(int i=1;i<deg;i<<=1){
      ret = ((ret + ret) - (ret * ret * prefix(i << 1))).prefix(i << 1);
    }
    return ret.prefix(deg);
  }
  fps modulo(const fps &vr){
    int n = this->size();
    int m = vr.size();
    if(n<m) return (*this);
    n = n - m + 1;
    //n-m次多項式を返す
    fps r = ((rev().prefix(n) * vr.rev().inv(n)).prefix(n).rev(n)) * vr;
    return (*this - r).prefix(m-1);
  }
};
// n: 多項式の次数
// m: 評価する点の数
// n >> mの場合, n, m共に小さい場合は愚直に計算したほうが早い
// n < mの場合はn<-mとなるとして、
// n * m < 10^6 なら愚直
// n * m < 2*10^8　かつ n / m > 2000
template<ll MOD>
vector<ll> simpleMultipointEvaluation(const StaticModFPS<MOD> &F, const vector<ll> &v){
  int n = F.size();
  int m = v.size();
  vector<ll> ret(m, 0);
  for(int i=0;i<m;i++){
    ll x = 1;
    for(int j=0;j<n;j++){
      ret[i] = (ret[i] + (x * F[j])%MOD)%MOD;
      if(ret[i] < 0) ret[i] += MOD;
      x = (x * v[i])%MOD;
    }
  }
  return ret;
}
template<ll MOD>
vector<ll> MultipointEvaluation(StaticModFPS<MOD> F, const vector<ll> &v){
  typedef StaticModFPS<MOD> fps;
  ll n = F.size();
  ll m = v.size();
  int add_coef = -1;
  if(MOD==998244353){
    if(max(n, m)*m < 2000000LL || (max(n, m)*m<100000000LL&&n/m>2500LL)) {
      return simpleMultipointEvaluation<MOD>(F, v);
    }
  }else{
    if(max(n, m)*m < 15000000LL || (max(n, m)*m<100000000LL&&n/m>1000LL)) {
      return simpleMultipointEvaluation<MOD>(F, v);
    }
  }
  if(n < m) {
    add_coef = m-1;//m-1次に1を足す
    F.resize(m);
    F[m-1] = 1;
    n = m;
  }
  int N = 1;
  while(N < m) N *= 2;
  vector<fps> tree(2*N-1, fps{1});
  for(int i=0;i<m;i++) tree[N-1+i] = fps{(MOD-v[i]%MOD)%MOD, 1};
  for(int i=N-2;i>=0;i--){
    tree[i] = tree[i*2+1] * tree[i*2+2];
  }
  tree[0] = F.modulo(tree[0]);
  for(int i=1;i<2*N-1;i++){
    int par = (i-1)/2;
    tree[i] = tree[par].modulo(tree[i]);
  }
  vector<ll> ret(m);
  for(int i=0;i<m;i++) {
    ret[i] = tree[N-1+i][0];
    if(add_coef!=-1) ret[i] = (ret[i] - mpow(v[i], add_coef, MOD) + MOD)%MOD;
  }
  return ret;
}

template<ll MOD>
ll adjacentInterpolation(vector<ll> x, const vector<ll> &y, ll p){
  int n = x.size();
  p = ((p%MOD) + MOD)%MOD;
  for(ll i=0;i<n;i++){
    x[i] = ((x[i]%MOD) + MOD)%MOD;
    if(x[i]!=(x[0]+i)%MOD || x[i] == p) return -1;
  }
  ll fac[n+1], inv[n+1], finv[n+1];
  fac[0] = finv[0] = fac[1] = finv[1] = inv[1] = 1;
  for(int i=2;i<=n;i++){
    fac[i] = (fac[i-1]*i)%MOD;
    inv[i] = ((-(MOD/i)*inv[MOD%i])%MOD+MOD)%MOD;
    finv[i] = (finv[i-1]*inv[i])%MOD;
  }
  ll M = 1, ret = 0;
  for(int i=0;i<n;i++){
    ll diff = (p - x[i] + MOD)%MOD;
    M = (M * diff)%MOD;
  }
  for(int i=0;i<n;i++){
    ll iQ = (finv[i] * finv[n-1-i])%MOD;
    if((n-i-1)&1) iQ = (MOD - iQ)%MOD;
    ll iC = ((y[i]%MOD) * iQ)%MOD;
    if(iC<0) iC += MOD;
    ll diff = (p - x[i] + MOD)%MOD;
    if(diff<=n) ret = (ret + iC * inv[diff])%MOD;
    else ret = (ret + iC * mpow(diff, MOD-2, MOD))%MOD;
  }
  ret = (ret * M)%MOD;
  return ret;
}

template<ll MOD>
StaticModFPS<MOD> InterpolationbyDivideandConquer(
  const vector<ll> &fx, StaticModFPS<MOD> F, const vector<StaticModFPS<MOD>> &tree, int k, int l, int r){
  if(r-l==1) return StaticModFPS<MOD>{(fx[l]*mpow(F[0], MOD-2, MOD))%MOD};
  int mid = (l+r)/2;
  if(tree[k*2+2].size()==0) return InterpolationbyDivideandConquer<MOD>(fx, F, tree, k*2+1, l, mid);
  StaticModFPS<MOD> left = InterpolationbyDivideandConquer<MOD>(fx, F.modulo(tree[k*2+1]), tree, k*2+1, l, mid);
  StaticModFPS<MOD> right = InterpolationbyDivideandConquer<MOD>(fx, F.modulo(tree[k*2+2]), tree, k*2+2, mid, r);
  return left * tree[k*2+2] + right * tree[k*2+1];
}
template<ll MOD>
StaticModFPS<MOD> PolynomialInterpolation(const vector<ll> &xi, const vector<ll> &fx){
  typedef StaticModFPS<MOD> fps;
  int n = xi.size();
  int N = 1;
  while(N<n) N*=2;
  vector<fps> tree(2*N-1, fps{});
  for(int i=0;i<n;i++) tree[N-1+i] = fps{(MOD - xi[i]%MOD)%MOD, 1};
  for(int i=N-2;i>=0;i--){
    if(tree[i*2+2].size()==0) tree[i] = tree[i*2+1];
    else tree[i] = tree[i*2+1] * tree[i*2+2];
  }
  for(ll i=0;i<n;i++) tree[0][i] = (tree[0][i+1] * (i+1))%MOD;
  tree[0].pop_back();
  return InterpolationbyDivideandConquer<MOD>(fx, tree[0], tree, 0, 0, N).prefix(n);
}

// O(N^2logN) (deg(f)=deg(g)=deg(ans)=Nとして)
template<ll MOD>
StaticModFPS<MOD> simpleComposition(const StaticModFPS<MOD> &f, const StaticModFPS<MOD> &g, int deg){
  int n = f.size();
  StaticModFPS<MOD> c(deg, 0), gpower{1};
  for(int i=0;i<min(deg, n);i++){
    int d = min(deg, (int)gpower.size());
    for(int j=0;j<d;j++){
      c[j] = (c[j] + f[i] * gpower[j])%MOD;
    }
    gpower *= g;
    if(gpower.size() > deg) gpower.resize(deg);
  }
  return c;
}
// reference: http://www.eecs.harvard.edu/~htk/publication/1978-jacm-brent-kung.pdf
// N-1次多項式 f(x)にM-1次多項式 g(x)を合成した結果を求める
// deg(f) = deg(g) = deg(ans) = Nとして、
// f(x)をk:=ceil(√N+1)ブロックに平方分割すると　f(x) = f_0(x) + f_1(x) x^k ... と約k項になる
// f((g(x))) = f_0(g(x)) + f_1(g(x))g(x)^k + ...

// 1. f_i(g(x))はk項とN項の合成なので、g(x)^i (0<=i<=k)を前処理しておくとO(Nk)
// 2. g(x)^kiを求める、かけるのは共に一回あたりO(NlogN)
// (1, 2)をkブロック分行うのでO(k × (Nk + NlogN)) = O(N^2 + N^1.5 logN)

template<ll MOD>
StaticModFPS<MOD> Composition(const StaticModFPS<MOD> &f, const StaticModFPS<MOD> &g, int deg){
  typedef StaticModFPS<MOD> fps;
  int n = f.size();
  int k = (int)sqrt(n);
  if(k*k<n) k++;
  int d = n / k;
  if(k*d<n) d++;
  vector<fps> gpower(k+1, {1});
  for(int i=1;i<=k;i++){
    gpower[i] = gpower[i-1] * g;
    if(gpower[i].size()>deg) gpower[i].resize(deg);
  }
  vector<fps> fi(k, fps(deg, 0));
  for(int i=0;i<k;i++){
    for(int j=0;j<d;j++){
      int idx = i*d+j;
      if(idx>=n) break;
      int sz = gpower[j].size();
      for(int t=0;t<sz;t++){
        fi[i][t] = (fi[i][t] + gpower[j][t] * f[idx])%MOD;
      }
    }
  }
  fps ret(deg, 0), gd = {1};
  for(int i=0;i<k;i++){
    fi[i] *= gd;
    int sz = min(deg, (int)fi[i].size());
    for(int j=0;j<sz;j++) ret[j] = (ret[j] + fi[i][j])%MOD;
    gd *= gpower[d];
    if(gd.size() > deg) gd.resize(deg);
  }
  return ret;
}

template<ll MOD>
ll PowerSum(ll n, ll k){
  vector<ll> x(k+2), y(k+2);
  for(int i=0;i<k+2;i++){
    x[i] = i, y[i] = (i==0?0:y[i-1]) + mpow(i, k, MOD);
    if(y[i]>=MOD) y[i] -= MOD;
    if(i==n) return y[i];
  }
  return adjacentInterpolation<MOD>(x, y, n);
}

int main(){
  ll n, k;scanf("%lld %lld", &n, &k);
  std::cout << PowerSum<1000000007>(n, k) << '\n';
  /*
  ll r, d, n;scanf("%lld %lld %lld", &r, &d, &n);
  vector<ll> x(d+2), y(d+2);
  ll P = 998244353, tmp = 0;
  ll rs = 1;
  for(int i=0;i<d+2;i++){
    tmp += (rs * mpow(i, d, P))%P;
    if(tmp >= P) tmp -= P;
    x[i] = i, y[i] = tmp;
    rs = (rs * r)%P;
  }
  if(n==0) {
    std::cout << 0 << '\n';
    return 0;
  }
  if(n<d+3) {
    std::cout << y[n-1] << '\n';
    return 0;
  }
  std::cout << adjacentInterpolation<998244353>(x, y, n-1) << '\n';
  */
}