#include <bits/stdc++.h>
using namespace std;
template <typename T1, typename T2>
istream &operator>>(istream &in, pair<T1, T2> &a) {
  in >> a.first >> a.second;
  return in;
}
template <typename T1, typename T2>
ostream &operator<<(ostream &out, pair<T1, T2> a) {
  out << a.first << " " << a.second;
  return out;
}
template <typename T>
pair<T, T> operator+=(pair<T, T> &a, pair<T, T> b) {
  a.first += b.first;
  a.second += b.second;
  return a;
}
template <typename T>
pair<T, T> operator-=(pair<T, T> &a, pair<T, T> b) {
  a.first -= b.first;
  a.second -= b.second;
  return a;
}
template <typename T>
pair<T, T> operator+(pair<T, T> a, pair<T, T> b) {
  a.first += b.first;
  a.second += b.second;
  return a;
}
template <typename T>
pair<T, T> operator-(pair<T, T> a, pair<T, T> b) {
  a.first -= b.first;
  a.second -= b.second;
  return a;
}
template <typename T, typename T1>
T amax(T &a, T1 b) {
  if (b > a) a = b;
  return a;
}
template <typename T, typename T1>
T amin(T &a, T1 b) {
  if (b < a) a = b;
  return a;
}
const int INF = 1e18;
const int N = 1e6 + 5;
int pw(int a, int b = 1000000007 - 2, int p = 1000000007) {
  if (b == 0) return 1;
  int t = pw(a, b / 2, p);
  if (b % 2)
    return (((1LL * a * t) % p) * t) % p;
  else
    return ((1LL * t * t) % p);
}
struct mod_int {
  int val;
  mod_int(long long v = 0) {
    if (v < 0) v = v % 1000000007 + 1000000007;
    if (v >= 1000000007) v %= 1000000007;
    val = v;
  }
  static int mod_inv(int a, int m = 1000000007) {
    int g = m, r = a, x = 0, y = 1;
    while (r != 0) {
      int q = g / r;
      g %= r;
      swap(g, r);
      x -= q * y;
      swap(x, y);
    }
    return x < 0 ? x + m : x;
  }
  explicit operator int() const { return val; }
  mod_int &operator+=(const mod_int &other) {
    val += other.val;
    if (val >= 1000000007) val -= 1000000007;
    return *this;
  }
  mod_int &operator-=(const mod_int &other) {
    val -= other.val;
    if (val < 0) val += 1000000007;
    return *this;
  }
  static unsigned fast_mod(uint64_t x, unsigned m = 1000000007) {
    return x % m;
  }
  mod_int &operator*=(const mod_int &other) {
    val = fast_mod((uint64_t)val * other.val);
    return *this;
  }
  mod_int &operator/=(const mod_int &other) { return *this *= other.inv(); }
  friend mod_int operator+(const mod_int &a, const mod_int &b) {
    return mod_int(a) += b;
  }
  friend mod_int operator-(const mod_int &a, const mod_int &b) {
    return mod_int(a) -= b;
  }
  friend mod_int operator*(const mod_int &a, const mod_int &b) {
    return mod_int(a) *= b;
  }
  friend mod_int operator/(const mod_int &a, const mod_int &b) {
    return mod_int(a) /= b;
  }
  mod_int &operator++() {
    val = val == 1000000007 - 1 ? 0 : val + 1;
    return *this;
  }
  mod_int &operator--() {
    val = val == 0 ? 1000000007 - 1 : val - 1;
    return *this;
  }
  mod_int operator++(int32_t) {
    mod_int before = *this;
    ++*this;
    return before;
  }
  mod_int operator--(int32_t) {
    mod_int before = *this;
    --*this;
    return before;
  }
  mod_int operator-() const { return val == 0 ? 0 : 1000000007 - val; }
  bool operator==(const mod_int &other) const { return val == other.val; }
  bool operator!=(const mod_int &other) const { return val != other.val; }
  mod_int inv() const { return mod_inv(val); }
  mod_int pow(long long p = 1000000007 - 2) const {
    assert(p >= 0);
    mod_int a = *this, result = 1;
    while (p > 0) {
      if (p & 1) result *= a;
      a *= a;
      p >>= 1;
    }
    return result;
  }
  friend ostream &operator<<(ostream &stream, const mod_int &m) {
    return stream << m.val;
  }
  friend istream &operator>>(istream &stream, mod_int &m) {
    return stream >> m.val;
  }
};
mod_int P[N];
mod_int ele[N];
void solve() {
  int n, k;
  cin >> n >> k;
  for (int i = 1; i < k + 4; i++) {
    mod_int j = i;
    P[i] = P[i - 1] + j.pow(k);
  }
  if (n <= k + 2) {
    cout << P[n];
    return;
  }
  ele[1] = 1;
  for (int i = 1; i < k + 2; i++) {
    ele[1] *= (1000000007 + n - i - 1);
    ele[1] /= mod_int(1000000007 - i);
  }
  mod_int ans = 0;
  ans += ele[1] * P[1];
  for (int i = 2; i < k + 3; i++) {
    ele[i] = ele[i - 1] * (1000000007 - k + i - 3);
    ele[i] /= mod_int(i - 1);
    ele[i] /= mod_int(1000000007 + n - i);
    ele[i] *= (1000000007 + n - i + 1);
    ans += ele[i] * P[i];
  }
  cout << ans;
}
signed main() {
  ios_base::sync_with_stdio(false);
  cin.tie(0);
  cout.tie(0);
  int t = 1;
  while (t--) solve();
  return 0;
}

### Prompt

Please provide a cpp coded solution to the problem described below:

There are well-known formulas: <image>, <image>, <image>. Also mathematicians found similar formulas for higher degrees.

Find the value of the sum <image> modulo 109 + 7 (so you should find the remainder after dividing the answer by the value 109 + 7).

Input

The only line contains two integers n, k (1 ≤ n ≤ 109, 0 ≤ k ≤ 106).

Output

Print the only integer a — the remainder after dividing the value of the sum by the value 109 + 7.

Examples

Input

4 1


Output

10


Input

4 2


Output

30


Input

4 3


Output

100


Input

4 0


Output

4

### Response

```cpp
#include <bits/stdc++.h>
using namespace std;
template <typename T1, typename T2>
istream &operator>>(istream &in, pair<T1, T2> &a) {
  in >> a.first >> a.second;
  return in;
}
template <typename T1, typename T2>
ostream &operator<<(ostream &out, pair<T1, T2> a) {
  out << a.first << " " << a.second;
  return out;
}
template <typename T>
pair<T, T> operator+=(pair<T, T> &a, pair<T, T> b) {
  a.first += b.first;
  a.second += b.second;
  return a;
}
template <typename T>
pair<T, T> operator-=(pair<T, T> &a, pair<T, T> b) {
  a.first -= b.first;
  a.second -= b.second;
  return a;
}
template <typename T>
pair<T, T> operator+(pair<T, T> a, pair<T, T> b) {
  a.first += b.first;
  a.second += b.second;
  return a;
}
template <typename T>
pair<T, T> operator-(pair<T, T> a, pair<T, T> b) {
  a.first -= b.first;
  a.second -= b.second;
  return a;
}
template <typename T, typename T1>
T amax(T &a, T1 b) {
  if (b > a) a = b;
  return a;
}
template <typename T, typename T1>
T amin(T &a, T1 b) {
  if (b < a) a = b;
  return a;
}
const int INF = 1e18;
const int N = 1e6 + 5;
int pw(int a, int b = 1000000007 - 2, int p = 1000000007) {
  if (b == 0) return 1;
  int t = pw(a, b / 2, p);
  if (b % 2)
    return (((1LL * a * t) % p) * t) % p;
  else
    return ((1LL * t * t) % p);
}
struct mod_int {
  int val;
  mod_int(long long v = 0) {
    if (v < 0) v = v % 1000000007 + 1000000007;
    if (v >= 1000000007) v %= 1000000007;
    val = v;
  }
  static int mod_inv(int a, int m = 1000000007) {
    int g = m, r = a, x = 0, y = 1;
    while (r != 0) {
      int q = g / r;
      g %= r;
      swap(g, r);
      x -= q * y;
      swap(x, y);
    }
    return x < 0 ? x + m : x;
  }
  explicit operator int() const { return val; }
  mod_int &operator+=(const mod_int &other) {
    val += other.val;
    if (val >= 1000000007) val -= 1000000007;
    return *this;
  }
  mod_int &operator-=(const mod_int &other) {
    val -= other.val;
    if (val < 0) val += 1000000007;
    return *this;
  }
  static unsigned fast_mod(uint64_t x, unsigned m = 1000000007) {
    return x % m;
  }
  mod_int &operator*=(const mod_int &other) {
    val = fast_mod((uint64_t)val * other.val);
    return *this;
  }
  mod_int &operator/=(const mod_int &other) { return *this *= other.inv(); }
  friend mod_int operator+(const mod_int &a, const mod_int &b) {
    return mod_int(a) += b;
  }
  friend mod_int operator-(const mod_int &a, const mod_int &b) {
    return mod_int(a) -= b;
  }
  friend mod_int operator*(const mod_int &a, const mod_int &b) {
    return mod_int(a) *= b;
  }
  friend mod_int operator/(const mod_int &a, const mod_int &b) {
    return mod_int(a) /= b;
  }
  mod_int &operator++() {
    val = val == 1000000007 - 1 ? 0 : val + 1;
    return *this;
  }
  mod_int &operator--() {
    val = val == 0 ? 1000000007 - 1 : val - 1;
    return *this;
  }
  mod_int operator++(int32_t) {
    mod_int before = *this;
    ++*this;
    return before;
  }
  mod_int operator--(int32_t) {
    mod_int before = *this;
    --*this;
    return before;
  }
  mod_int operator-() const { return val == 0 ? 0 : 1000000007 - val; }
  bool operator==(const mod_int &other) const { return val == other.val; }
  bool operator!=(const mod_int &other) const { return val != other.val; }
  mod_int inv() const { return mod_inv(val); }
  mod_int pow(long long p = 1000000007 - 2) const {
    assert(p >= 0);
    mod_int a = *this, result = 1;
    while (p > 0) {
      if (p & 1) result *= a;
      a *= a;
      p >>= 1;
    }
    return result;
  }
  friend ostream &operator<<(ostream &stream, const mod_int &m) {
    return stream << m.val;
  }
  friend istream &operator>>(istream &stream, mod_int &m) {
    return stream >> m.val;
  }
};
mod_int P[N];
mod_int ele[N];
void solve() {
  int n, k;
  cin >> n >> k;
  for (int i = 1; i < k + 4; i++) {
    mod_int j = i;
    P[i] = P[i - 1] + j.pow(k);
  }
  if (n <= k + 2) {
    cout << P[n];
    return;
  }
  ele[1] = 1;
  for (int i = 1; i < k + 2; i++) {
    ele[1] *= (1000000007 + n - i - 1);
    ele[1] /= mod_int(1000000007 - i);
  }
  mod_int ans = 0;
  ans += ele[1] * P[1];
  for (int i = 2; i < k + 3; i++) {
    ele[i] = ele[i - 1] * (1000000007 - k + i - 3);
    ele[i] /= mod_int(i - 1);
    ele[i] /= mod_int(1000000007 + n - i);
    ele[i] *= (1000000007 + n - i + 1);
    ans += ele[i] * P[i];
  }
  cout << ans;
}
signed main() {
  ios_base::sync_with_stdio(false);
  cin.tie(0);
  cout.tie(0);
  int t = 1;
  while (t--) solve();
  return 0;
}

```