#include <bits/stdc++.h>
#ifndef ATCODER_INTERNAL_BITOP_HPP
#define ATCODER_INTERNAL_BITOP_HPP 1
#ifdef _MSC_VER
#include <intrin.h>
#endif
namespace atcoder
{
    namespace internal
    {
        // @param n `0 <= n`
        // @return minimum non-negative `x` s.t. `n <= 2**x`
        int ceil_pow2(int n)
        {
            int x = 0;
            while ((1U << x) < (unsigned int)(n))
                x++;
            return x;
        }
        // @param n `1 <= n`
        // @return minimum non-negative `x` s.t. `(n & (1 << x)) != 0`
        int bsf(unsigned int n)
        {
#ifdef _MSC_VER
            unsigned long index;
            _BitScanForward(&index, n);
            return index;
#else
            return __builtin_ctz(n);
#endif
        }
    } // namespace internal
} // namespace atcoder
#endif // ATCODER_INTERNAL_BITOP_HPP
#ifndef ATCODER_INTERNAL_MATH_HPP
#define ATCODER_INTERNAL_MATH_HPP 1
#include <utility>
namespace atcoder
{
    namespace internal
    {
        // @param m `1 <= m`
        // @return x mod m
        constexpr long long safe_mod(long long x, long long m)
        {
            x %= m;
            if (x < 0)
                x += m;
            return x;
        }
        // Fast moduler by barrett reduction
        // Reference: https://en.wikipedia.org/wiki/Barrett_reduction
        // NOTE: reconsider after Ice Lake
        struct barrett
        {
            unsigned int _m;
            unsigned long long im;
            // @param m `1 <= m`
            barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {}
            // @return m
            unsigned int umod() const { return _m; }
            // @param a `0 <= a < m`
            // @param b `0 <= b < m`
            // @return `a * b % m`
            unsigned int mul(unsigned int a, unsigned int b) const
            {
                // [1] m = 1
                // a = b = im = 0, so okay
                // [2] m >= 2
                // im = ceil(2^64 / m)
                // -> im * m = 2^64 + r (0 <= r < m)
                // let z = a*b = c*m + d (0 <= c, d < m)
                // a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
                // c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
                // ((ab * im) >> 64) == c or c + 1
                unsigned long long z = a;
                z *= b;
#ifdef _MSC_VER
                unsigned long long x;
                _umul128(z, im, &x);
#else
                unsigned long long x =
                    (unsigned long long)(((unsigned __int128)(z)*im) >> 64);
#endif
                unsigned int v = (unsigned int)(z - x * _m);
                if (_m <= v)
                    v += _m;
                return v;
            }
        };
        // @param n `0 <= n`
        // @param m `1 <= m`
        // @return `(x ** n) % m`
        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;
        }
        // Reference:
        // M. Forisek and J. Jancina,
        // Fast Primality Testing for Integers That Fit into a Machine Word
        // @param n `0 <= n`
        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;
            for (long long a : {2, 7, 61})
            {
                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);
        // @param b `1 <= b`
        // @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g
        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};
            // Contracts:
            // [1] s - m0 * a = 0 (mod b)
            // [2] t - m1 * a = 0 (mod b)
            // [3] s * |m1| + t * |m0| <= b
            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; // |m1 * u| <= |m1| * s <= b
                // [3]:
                // (s - t * u) * |m1| + t * |m0 - m1 * u|
                // <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u)
                // = s * |m1| + t * |m0| <= b
                auto tmp = s;
                s = t;
                t = tmp;
                tmp = m0;
                m0 = m1;
                m1 = tmp;
            }
            // by [3]: |m0| <= b/g
            // by g != b: |m0| < b/g
            if (m0 < 0)
                m0 += b / s;
            return {s, m0};
        }
        // Compile time primitive root
        // @param m must be prime
        // @return primitive root (and minimum in now)
        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);
    } // namespace internal
} // namespace atcoder
#endif // ATCODER_INTERNAL_MATH_HPP
#ifndef ATCODER_INTERNAL_QUEUE_HPP
#define ATCODER_INTERNAL_QUEUE_HPP 1
#include <vector>
namespace atcoder
{
    namespace internal
    {
        template <class T>
        struct simple_queue
        {
            std::vector<T> payload;
            int pos = 0;
            void reserve(int n) { payload.reserve(n); }
            int size() const { return int(payload.size()) - pos; }
            bool empty() const { return pos == int(payload.size()); }
            void push(const T &t) { payload.push_back(t); }
            T &front() { return payload[pos]; }
            void clear()
            {
                payload.clear();
                pos = 0;
            }
            void pop() { pos++; }
        };
    } // namespace internal
} // namespace atcoder
#endif // ATCODER_INTERNAL_QUEUE_HPP
#ifndef ATCODER_INTERNAL_SCC_HPP
#define ATCODER_INTERNAL_SCC_HPP 1
#include <algorithm>
#include <utility>
#include <vector>
namespace atcoder
{
    namespace internal
    {
        template <class E>
        struct csr
        {
            std::vector<int> start;
            std::vector<E> elist;
            csr(int n, const std::vector<std::pair<int, E>> &edges)
                : start(n + 1), elist(edges.size())
            {
                for (auto e : edges)
                {
                    start[e.first + 1]++;
                }
                for (int i = 1; i <= n; i++)
                {
                    start[i] += start[i - 1];
                }
                auto counter = start;
                for (auto e : edges)
                {
                    elist[counter[e.first]++] = e.second;
                }
            }
        };
        // Reference:
        // R. Tarjan,
        // Depth-First Search and Linear Graph Algorithms
        struct scc_graph
        {
        public:
            scc_graph(int n) : _n(n) {}
            int num_vertices() { return _n; }
            void add_edge(int from, int to) { edges.push_back({from, {to}}); }
            // @return pair of (# of scc, scc id)
            std::pair<int, std::vector<int>> scc_ids()
            {
                auto g = csr<edge>(_n, edges);
                int now_ord = 0, group_num = 0;
                std::vector<int> visited, low(_n), ord(_n, -1), ids(_n);
                visited.reserve(_n);
                auto dfs = [&](auto self, int v) -> void {
                    low[v] = ord[v] = now_ord++;
                    visited.push_back(v);
                    for (int i = g.start[v]; i < g.start[v + 1]; i++)
                    {
                        auto to = g.elist[i].to;
                        if (ord[to] == -1)
                        {
                            self(self, to);
                            low[v] = std::min(low[v], low[to]);
                        }
                        else
                        {
                            low[v] = std::min(low[v], ord[to]);
                        }
                    }
                    if (low[v] == ord[v])
                    {
                        while (true)
                        {
                            int u = visited.back();
                            visited.pop_back();
                            ord[u] = _n;
                            ids[u] = group_num;
                            if (u == v)
                                break;
                        }
                        group_num++;
                    }
                };
                for (int i = 0; i < _n; i++)
                {
                    if (ord[i] == -1)
                        dfs(dfs, i);
                }
                for (auto &x : ids)
                {
                    x = group_num - 1 - x;
                }
                return {group_num, ids};
            }
            std::vector<std::vector<int>> scc()
            {
                auto ids = scc_ids();
                int group_num = ids.first;
                std::vector<int> counts(group_num);
                for (auto x : ids.second)
                    counts[x]++;
                std::vector<std::vector<int>> groups(ids.first);
                for (int i = 0; i < group_num; i++)
                {
                    groups[i].reserve(counts[i]);
                }
                for (int i = 0; i < _n; i++)
                {
                    groups[ids.second[i]].push_back(i);
                }
                return groups;
            }

        private:
            int _n;
            struct edge
            {
                int to;
            };
            std::vector<std::pair<int, edge>> edges;
        };
    } // namespace internal
} // namespace atcoder
#endif // ATCODER_INTERNAL_SCC_HPP
#ifndef ATCODER_INTERNAL_TYPE_TRAITS_HPP
#define ATCODER_INTERNAL_TYPE_TRAITS_HPP 1
#include <cassert>
#include <numeric>
#include <type_traits>
namespace atcoder
{
    namespace internal
    {
#ifndef _MSC_VER
        template <class T>
        using is_signed_int128 =
            typename std::conditional<std::is_same<T, __int128_t>::value ||
                                          std::is_same<T, __int128>::value,
                                      std::true_type,
                                      std::false_type>::type;
        template <class T>
        using is_unsigned_int128 =
            typename std::conditional<std::is_same<T, __uint128_t>::value ||
                                          std::is_same<T, unsigned __int128>::value,
                                      std::true_type,
                                      std::false_type>::type;
        template <class T>
        using make_unsigned_int128 =
            typename std::conditional<std::is_same<T, __int128_t>::value,
                                      __uint128_t,
                                      unsigned __int128>;
        template <class T>
        using is_integral = typename std::conditional<std::is_integral<T>::value ||
                                                          is_signed_int128<T>::value ||
                                                          is_unsigned_int128<T>::value,
                                                      std::true_type,
                                                      std::false_type>::type;
        template <class T>
        using is_signed_int = typename std::conditional<(is_integral<T>::value &&
                                                         std::is_signed<T>::value) ||
                                                            is_signed_int128<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) ||
                                          is_unsigned_int128<T>::value,
                                      std::true_type,
                                      std::false_type>::type;
        template <class T>
        using to_unsigned = typename std::conditional<
            is_signed_int128<T>::value,
            make_unsigned_int128<T>,
            typename std::conditional<std::is_signed<T>::value,
                                      std::make_unsigned<T>,
                                      std::common_type<T>>::type>::type;
#else
        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;
#endif
        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;
    } // namespace internal
} // namespace atcoder
#endif // ATCODER_INTERNAL_TYPE_TRAITS_HPP
#ifndef ATCODER_MODINT_HPP
#define ATCODER_MODINT_HPP 1
#include <cassert>
#include <numeric>
#include <type_traits>
#ifdef _MSC_VER
#include <intrin.h>
#endif
namespace atcoder
{
    namespace internal
    {
        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>;
    } // namespace internal
    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>;
    namespace internal
    {
        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>;
    } // namespace internal
} // namespace atcoder
#endif // ATCODER_MODINT_HPP
#pragma GCC optimize("Ofast")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#pragma GCC target("avx")
#include <bits/stdc++.h>
#define _GLIBCXX_DEBUG
using namespace std;
using namespace atcoder;
using ll = long long;
using vec = vector<ll>;
using vect = vector<double>;
using Graph = vector<vector<ll>>;
#define endl '\n'
#define loop(i, n) for (register int i = 0; i < n; i++)
#define Loop(i, m, n) for (int i = m; i < n; i++)
#define pool(i, n) for (int i = n; i >= 0; i--)
#define Pool(i, m, n) for (int i = n; i >= m; i--)
#define modd 1000000007ll
//#define modd 998244353ll
#define flagcount(bit) __builtin_popcount(bit)
#define flag(x) (1ll << x)
#define flagadd(bit, x) bit |= flag(x)
#define flagpop(bit, x) bit &= ~flag(x)
#define flagon(bit, i) bit &flag(i)
#define flagoff(bit, i) !(bit & (1ll << i))
#define all(v) v.begin(), v.end()
#define low2way(v, x) lower_bound(all(v), x)
#define high2way(v, x) upper_bound(all(v), x)
#define idx_lower(v, x) (distance(v.begin(), low2way(v, x)))  //配列vでx未満の要素数を返す
#define idx_upper(v, x) (distance(v.begin(), high2way(v, x))) //配列vでx以下の要素数を返す
#define idx_lower2(v, x) (v.size() - idx_lower(v, x))         //配列vでx以上の要素数を返す
#define idx_upper2(v, x) (v.size() - idx_upper(v, x))         //配列vでxより大きい要素の数を返す
#define putout(a) cout << a << '\n'
#define Sum(v) accumulate(all(v), 0ll)
ll ctoi(char c)
{
    if (c >= '0' && c <= '9')
    {
        return c - '0';
    }
    return -1;
}
template <typename T>
string make_string(T N)
{
    string ret;
    T now = N;
    while (now > 0)
    {
        T x = now % 10;
        ret += (char)('0' + x);
        now /= 10;
    }
    reverse(all(ret));
    return ret;
}
template <typename T>
T gcd(T a, T b)
{
    if (a % b == 0)
    {
        return (b);
    }
    else
    {
        return (gcd(b, a % b));
    }
}
template <typename T>
T lcm(T x, T y)
{
    T z = gcd(x, y);
    return x * y / z;
}
template <typename T>
bool primejudge(T n)
{
    if (n < 2)
        return false;
    else if (n == 2)
        return true;
    else if (n % 2 == 0)
        return false;
    double sqrtn = sqrt(n);
    for (T i = 3; i < sqrtn + 1; i++)
    {
        if (n % i == 0)
        {
            return false;
        }
        i++;
    }
    return true;
}
template <typename T>
bool chmax(T &a, const T &b)
{
    if (a < b)
    {
        a = b; // aをbで更新
        return true;
    }
    return false;
}
template <typename T>
bool chmin(T &a, const T &b)
{
    if (a > b)
    {
        a = b; // aをbで更新
        return true;
    }
    return false;
}
//場合によって使い分ける
//const ll dx[4]={1,0,-1,0};
//const ll dy[4]={0,1,0,-1};
const ll dx[8] = {1, 1, 0, -1, -1, -1, 0, 1};
const ll dy[8] = {0, 1, 1, 1, 0, -1, -1, -1};
//cout << fixed << setprecision(10);
//vector<vector<ll>> field(h, vector<ll>(w));
using mint = modint1000000007;
mint modpow(mint a, long long n)
{
    mint res = 1;
    while (n > 0)
    {
        if (n & 1)
            res = res * a;
        a *= a;
        n >>= 1;
    }
    return res;
}
mint modinv(mint a)
{
    return modpow(a, modd - 2);
}
vector<mint> fact(3000001);    //fact[i]=(i!)
vector<mint> factinv(3000001); //factinv[i]=(i!)^-1
void COMinit()
{
    mint now = 1;
    fact[0] = 1;
    factinv[0] = modinv(0);
    for (long long i = 1; i < 3000001; i++)
    {
        now *= i;
        fact[i] = now;
        factinv[i] = modinv(now);
    }
}
mint COM(long long n, long long r)
{
    if (n < r)
        return 0;
    if (n < 0 || r < 0)
        return 0;
    if (n == r)
        return 1;
    if (r == 0)
        return 1;
    mint ans = fact[n];
    ans *= factinv[r];
    ans *= factinv[n - r];
    return ans;
}
//1^K+2^K+...+N^Kを求める
mint Faulhaber(ll N, ll K)
{
    N++;
    vector<mint> B(K + 1, -1);
    B[0] = 1;
    for (ll i = 1; i <= K; i++)
    {
        B[i] *= modinv(i + 1);
        mint S = 0;
        for (ll k = 0; k <= i - 1; k++)
        {
            mint x = COM(1 + i, k);
            x *= B[k];
            S += x;
        }
        B[i] *= S;
    }
    vector<mint> pown(K + 2);
    pown[0] = 1;
    loop(i, K + 1) pown[i + 1] = pown[i] * N;
    mint ans = modinv(K + 1);
    mint SUM = 0;
    for (ll i = 0; i <= K; i++)
    {
        mint x = COM(K + 1, i);
        x *= B[i];
        x *= pown[K + 1 - i];
        SUM += x;
    }
    ans *= SUM;
    return ans;
}
int main()
{
    COMinit();
    ll N, K;
    cin >> N >> K;
    mint ANS = Faulhaber(N, K);
    putout(ANS.val());
    return 0;
}