#include #include #include #include #include #include #include #include #include #include #include static const int MOD = 1000000007; using ll = long long; using u32 = uint32_t; using namespace std; template constexpr T INF = ::numeric_limits::max()/32*15+208; template struct modint{ ll val; modint(): val(0){} template explicit modint(T t){val = t%M; if(val < 0) val += M;} modint pow(ll k){ modint res(1), x(val); while(k){ if(k&1) res *= x; x *= x; k >>= 1; } return res; } template modint& operator=(T a){ val = a%M; if(val < 0) val += M; return *this; } modint inv() {return pow(M-2);} modint& operator+=(modint a){ val += a.val; if(val >= M) val -= M; return *this;} modint& operator-=(modint a){ val += M-a.val; if(val >= M) val -= M; return *this;} modint& operator*=(modint a){ val = 1LL*val*a.val%M; return *this;} modint& operator/=(modint a){ return (*this) *= a.inv();} modint operator+(modint a) const {return modint(val) +=a;} modint operator-(modint a) const {return modint(val) -=a;} modint operator*(modint a) const {return modint(val) *=a;} modint operator/(modint a) const {return modint(val) /=a;} modint operator-(){ return modint(-val);} bool operator==(const modint a) const {return val == a.val;} bool operator!=(const modint a) const {return val != a.val;} bool operator<(const modint a) const {return val < a.val;} }; template struct matrix { vector> A; matrix() = default; matrix(size_t n, size_t m) : A(n, vector(m)) {} explicit matrix(size_t n) : A(n, vector (n)) {}; size_t height() const { return (A.size()); } size_t width() const { return (A[0].size()); } const vector &operator [] (int k) const { return A.at(k); } vector &operator[] (int k) { return A.at(k); } static matrix I(size_t n){ matrix mat(n); for (int i = 0; i < n; ++i) mat[i][i] = 1; return mat; } matrix &operator+= (const matrix &B){ size_t h = height(), w = width(); for (int i = 0; i < h; ++i) { for (int j = 0; j < w; ++j) { (*this)[i][j] += B[i][j]; } } } matrix &operator-= (const matrix &B){ size_t h = height(), w = width(); for (int i = 0; i < h; ++i) { for (int j = 0; j < w; ++j) { (*this)[i][j] -= B[i][j]; } } } matrix &operator*=(const matrix &B) { size_t n = height(), m = B.width(), p = width(); matrix C (n, m); for (int i = 0; i < n; ++i) { for (int j = 0; j < m; ++j) { for (int k = 0; k < p; ++k) { C[i][j] = (C[i][j] + (*this)[i][k] * B[k][j]); } } } A.swap(C.A); return (*this); } template matrix &operator%= (const U &m){ for (int i = 0; i < height(); ++i) { for (int j = 0; j < width(); ++j) { (*this)[i][j] %= m; } } } matrix pow(ll n) const { matrix a = (*this), res = I(height()); while(n > 0){ if(n & 1) res *= a; a *= a; n >>= 1; } return res; } matrix operator+(const matrix &A) const {return matrix(*this) += A;} matrix operator-(const matrix &A) const {return matrix(*this) -= A;} matrix operator*(const matrix &A) const {return matrix(*this) *= A;} template matrix operator%(const U &m) const {return matrix(*this) %= m;} }; using mint = modint<>; int main() { int n, k; cin >> n >> k; int k3 = k*k*k; matrix M(k3, k3); auto f = [&k](int a, int b, int c){ return a+(b+c*k)*k; }; for (int a = 0; a < k; ++a) { for (int b = 0; b < k; ++b) { for (int c = 0; c < k; ++c) { M[f(a, b, c)][f((a+1)%k, b, c)] += mint(1); M[f(a, b, c)][f(a, (a+b)%k, c)] += mint(1); M[f(a, b, c)][f(a, b, (b+c)%k)] += mint(1); } } } M = M.pow(n); mint ans(0); for (int a = 0; a < k; ++a) { for (int b = 0; b < k; ++b) { ans += M[0][f(a, b, 0)]; } } cout << ans.val << "\n"; return 0; }