#include using namespace std; #define FOR(i,a,b) for(int i=(a);i<(b);++i) #define REP(i,n) FOR(i,0,n) #define ALL(v) begin(v),end(v) template inline bool chmax(A & a, const B & b) { if (a < b) { a = b; return true; } return false; } template inline bool chmin(A & a, const B & b) { if (a > b) { a = b; return true; } return false; } using ll = long long; using pii = pair; constexpr ll INF = 1ll<<30; constexpr ll longINF = 1ll<<60; constexpr ll MOD = 1000000007; constexpr bool debug = false; //---------------------------------// template struct ModInt { static_assert(M > 0); public: using value_type = int; using calc_type = std::int_fast64_t; private: value_type val_; public: constexpr ModInt(calc_type val = 0) : val_(val < 0 ? (val % M + M) % M : val % M) {} constexpr const value_type & val() const noexcept { return val_; } constexpr static decltype(M) mod() noexcept { return M; } explicit constexpr operator bool() const noexcept { return val_; } constexpr bool operator !() const noexcept { return !static_cast(*this); } constexpr ModInt operator +() const noexcept { return ModInt(*this); } constexpr ModInt operator -() const noexcept { return ModInt(-val_); } constexpr ModInt operator ++(int) noexcept { ModInt res = *this; ++*this; return res; } constexpr ModInt operator --(int) noexcept { ModInt res = *this; --*this; return res; } constexpr ModInt & operator ++() noexcept { val_ = val_ + 1 == M ? 0 : val_ + 1; return *this; } constexpr ModInt & operator --() noexcept { val_ = val_ == 0 ? M - 1 : val_ - 1; return *this; } constexpr ModInt & operator +=(const ModInt & rhs) noexcept { val_ += val_ < M - rhs.val_ ? rhs.val_ : rhs.val_ - M; return *this; } constexpr ModInt & operator -=(const ModInt & rhs) noexcept { val_ += val_ >= rhs.val_ ? -rhs.val_ : M - rhs.val_; return *this; } constexpr ModInt & operator *=(const ModInt & rhs) noexcept { val_ = static_cast(val_) * rhs.val_ % M; return *this; } constexpr ModInt & operator /=(const ModInt & rhs) noexcept { return *this *= rhs.inv(); } friend constexpr ModInt operator +(const ModInt & lhs, const ModInt & rhs) noexcept { return ModInt(lhs) += rhs; } friend constexpr ModInt operator -(const ModInt & lhs, const ModInt & rhs) noexcept { return ModInt(lhs) -= rhs; } friend constexpr ModInt operator *(const ModInt & lhs, const ModInt & rhs) noexcept { return ModInt(lhs) *= rhs; } friend constexpr ModInt operator /(const ModInt & lhs, const ModInt & rhs) noexcept { return ModInt(lhs) /= rhs; } friend constexpr bool operator ==(const ModInt & lhs, const ModInt & rhs) noexcept { return lhs.val_ == rhs.val_; } friend constexpr bool operator !=(const ModInt & lhs, const ModInt & rhs) noexcept { return !(lhs == rhs); } friend std::ostream & operator <<(std::ostream & os, const ModInt & rhs) { return os << rhs.val_; } friend std::istream & operator >>(std::istream & is, ModInt & rhs) { calc_type x; is >> x; rhs = ModInt(x); return is; } constexpr ModInt pow(calc_type n) const noexcept { ModInt res = 1, x = val_; if (n < 0) { x = x.inv(); n = -n; } while (n) { if (n & 1) res *= x; x *= x; n >>= 1; } return res; } constexpr ModInt inv() const noexcept { value_type a = val_, a1 = 1, a2 = 0, b = M, b1 = 0, b2 = 1; while (b > 0) { value_type q = a / b, r = a % b; value_type nb1 = a1 - q * b1, nb2 = a2 - q * b2; a = b; b = r; a1 = b1; b1 = nb1; a2 = b2; b2 = nb2; } assert(a == 1); return a1; } }; using mint = ModInt; template struct Matrix { public: using value_type = T; using size_type = std::size_t; Matrix() {} Matrix(size_type h, size_type w, const value_type & x = 0) : h(h), w(w), val(h, std::vector(w, x)) { assert(h > 0 && w > 0); } Matrix(std::vector> val) : h(val.size()), w(val.size() ? val[0].size() : 0), val(val) { assert(h > 0 && w > 0); for (size_type i = 1; i < h; ++i) assert(val[i].size() == w); } Matrix(std::initializer_list> init) : val(init.begin(), init.end()) { h = val.size(); w = val.size() ? val[0].size() : 0; assert(h > 0 && w > 0); for (size_type i = 1; i < h; ++i) assert(val[i].size() == w); } std::vector & operator [](size_type i) noexcept { return val[i]; } const std::vector & operator [](size_type i) const noexcept { return val[i]; }; value_type & operator ()(size_type i, size_type j) noexcept { return val[i][j]; }; const value_type & operator ()(size_type i, size_type j) const noexcept { return val[i][j]; } value_type & at(size_type i, size_type j) { assert(i < h && j < w); return val[i][j]; } const value_type & at(size_type i, size_type j) const { assert(i < h & j < w); return val[i][j]; } bool empty() const { return !(h || w); } std::pair type() const { return std::make_pair(h, w); } bool match_type(const Matrix & rhs) const noexcept { return h == rhs.h && w == rhs.w; } bool is_square() const { return h == w; } const std::vector> & get() const noexcept { return val; } bool operator ==(const Matrix & rhs) const noexcept { return match_type(rhs) && val == rhs.val; } bool operator !=(const Matrix & rhs) const noexcept { return !(*this == rhs); } Matrix operator +() const { return Matrix(*this); } Matrix operator -() const { return Matrix(h, w) - Matrix(*this); } Matrix operator +(const Matrix & rhs) const { return Matrix(*this) += rhs; } Matrix operator -(const Matrix & rhs) const { return Matrix(*this) -= rhs; } Matrix operator *(const Matrix & rhs) const { assert(w == rhs.h); Matrix res(h, rhs.w); for (size_type i = 0; i < h; ++i) for (size_type j = 0; j < rhs.w; ++j) for (size_type k = 0; k < w; ++k) res.val[i][j] += val[i][k] * rhs.val[k][j]; return res; } Matrix operator /(const Matrix & rhs) const { return Matrix(*this) /= rhs; } friend Matrix operator *(const value_type & lhs, const Matrix & rhs) { Matrix res(rhs.val); for (size_type i = 0; i < res.h; ++i) for (size_type j = 0; j < res.w; ++j) res.val[i][j] = lhs * res.val[i][j]; return res; } Matrix operator *(const value_type & rhs) const { Matrix res(val); for (size_type i = 0; i < h; ++i) for (size_type j = 0; j < w; ++j) res.val[i][j] *= rhs; return res; } Matrix operator /(const value_type & rhs) const { Matrix res(val); for (size_type i = 0; i < h; ++i) for (size_type j = 0; j < w; ++j) res.val[i][j] /= rhs; return res; } Matrix & operator +=(const Matrix & rhs) { assert(match_type(rhs)); for (size_type i = 0; i < h; ++i) for (size_type j = 0; j < w; ++j) val[i][j] += rhs.val[i][j]; return *this; } Matrix & operator -=(const Matrix & rhs) { assert(match_type(rhs)); for (size_type i = 0; i < h; ++i) for(size_type j = 0; j < w; ++j) val[i][j] -= rhs.val[i][j]; return *this; } Matrix & operator *=(const Matrix & rhs) { *this = *this * rhs; return *this; } Matrix & operator /=(const Matrix & rhs) { *this *= rhs.inverse(); return *this; } Matrix pow(long long n) const { Matrix res = identity(h), x(*this); if (n < 0) { x = x.inverse(); n = -n; } while (n) { if (n & 1) res *= x; x *= x; n >>= 1; } return res; } Matrix trans() const { Matrix res(w, h); for (size_type i = 0; i < h; ++i) for (size_type j = 0; j < w; ++j) res.val[j][i] = val[i][j]; return res; } Matrix inverse() const { assert(is_square()); Matrix aug_mat = this->hstack(identity(h)); if (aug_mat.gauss_jordan().first != h) return Matrix(); return aug_mat.submat(0, w, h, 2 * w); } Matrix vstack(const Matrix & A) const { assert(w == A.w); Matrix res(h + A.h, w); std::copy(val.begin(), val.end(), res.val.begin()); std::copy(A.val.begin(), A.val.end(), res.val.begin() + h); return res; } Matrix hstack(const Matrix & A) const { assert(h == A.h); Matrix res(h, w + A.w); for (int i = 0; i < h; ++i) { std::copy(val[i].begin(), val[i].end(), res.val[i].begin()); std::copy(A.val[i].begin(), A.val[i].end(), res.val[i].begin() + w); } return res; } Matrix submat(size_type i1, size_type j1, size_type i2, size_type j2) const { assert(i1 < i2 && j1 < j2 && i2 <= h && j2 <= w); Matrix res(i2 - i1, j2 - j1); for (size_type i = 0; i < i2 - i1; ++i) std::copy(val[i + i1].begin() + j1, val[i + i1].begin() + j2, res.val[i].begin()); return res; } static Matrix identity(size_type n) { Matrix res(n, n); for (size_type i = 0; i < n; ++i) res(i, i) = 1; return res; } std::pair gauss_jordan(size_type colnum = -1) { if (colnum == -1) colnum = w; size_type rank = 0; value_type det {}; bool done = false, rflag = false; for (size_type k = 0; k < colnum; ++k) { size_type pivot = -1; for (size_type i = rank; i < h; ++i) { if (val[i][k] != 0) { pivot = i; break; } } if (pivot == -1) continue; if (pivot != rank) { rflag ^= 1; std::swap(val[rank], val[pivot]); } if (!done) { det = val[rank][k]; done = true; } else det *= val[rank][k]; value_type div = static_cast(1) / val[rank][k]; for (size_type j = k; j < w; ++j) val[rank][j] *= div; for (size_type i = 0; i < h; ++i) if (i != rank) { for (size_type j = k + 1; j < w; ++j) val[i][j] -= val[rank][j] * val[i][k]; val[i][k] = 0; } ++rank; } if (!is_square() || rank != h) det = 0; else if (rflag) det *= -1; return {rank, det}; } friend std::ostream & operator <<(std::ostream & os, const Matrix & rhs) { os << "type = (" << rhs.h << "," << rhs.w << ") [\n"; for (size_type i = 0; i < rhs.h; ++i) for (size_type j = 0; j < rhs.w; ++j) os << (j == 0 ? " " : "") << rhs.val[i][j] << (j + 1 == rhs.w ? '\n' : ' '); return os << "]"; } private: size_type h, w; std::vector> val; }; int main() { int P, K; cin >> P >> K; Matrix dp(2, 1); // [0] := N = 0, [1] := N != 0 dp[0][0] = 1; Matrix A(2, 2); A[0][0] = P + 1; A[0][1] = 2; A[1][0] = P - 1; A[1][1] = 2 * (P - 1); mint ans = (A.pow(K) * dp)[0][0]; cout << ans << endl; }