#include using namespace std; using i64 = int64_t; #define rep(i, x, y) for (i64 i = i64(x), i##_max_for_repmacro = i64(y); i < i##_max_for_repmacro; ++i) #define debug(x) #x << "=" << (x) #ifdef DEBUG #define _GLIBCXX_DEBUG #define print(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl #else #define print(x) #endif template class fp { public: i64 x; fp() : x(0) {} fp(i64 x_) : x((x_ % p + p) % p) {} fp operator+() const { return fp(x); } fp operator-() const { return fp(-x); } fp& operator+=(const fp& y) { x += y.x; if (x >= p) x -= p; return *this; } fp& operator-=(const fp& y) { return *this += -y; } fp& operator*=(const fp& y) { x = x * y.x % p; return *this; } fp& operator/=(const fp& y) { return *this *= fp(inverse(y.x)); } fp operator+(const fp& y) const { return fp(x) += y; } fp operator-(const fp& y) const { return fp(x) -= y; } fp operator*(const fp& y) const { return fp(x) *= y; } fp operator/(const fp& y) const { return fp(x) /= y; } bool operator==(const fp& y) const { return x == y.x; } bool operator!=(const fp& y) const { return !(*this == y); } i64 extgcd(i64 a, i64 b, i64& x, i64& y) { i64 d = a; if (b != 0) { d = extgcd(b, a % b, y, x); y -= (a / b) * x; } else { x = 1; y = 0; } return d; } i64 inverse(i64 a) { i64 x, y; extgcd(a, p, x, y); return (x % p + p) % p; } }; template i64 abs(const fp

& x) { return x.x; } template istream& operator>>(istream& is, fp

& x) { is >> x.x; return is; } template ostream& operator<<(ostream& os, const fp

& x) { os << x.x; return os; } template ostream& operator<<(ostream& os, const pair& p) { os << "(" << p.first << ", " << p.second << ")"; return os; } template ostream& operator<<(ostream& os, const vector& vec) { os << "["; for (const auto& v : vec) { os << v << ","; } os << "]"; return os; } template bool chmin(T& a, const T& b) { if (a > b) { a = b; return true; } return false; } template bool chmax(T& a, const T& b) { if (a < b) { a = b; return true; } return false; } template void fill(A (&ary)[size], const T& val) { fill((T*)ary, (T*)(ary + size), val); } constexpr int inf = 1.01e9; constexpr i64 inf64 = 4.01e18; constexpr long double eps = 1e-9; // double(64bit浮動小数)のn分探索のループ回数の上限(2分探索なら50でも十分かもしれない). long double(80ビットの x87 浮動小数点型?)だと, 2分探索であってもこれだと足りないケースがある気がするので, もうちょっと余裕を持たせた方が良さそう. constexpr i64 max_loop = 100; template class matrix { public: int m, n; vector> a; matrix(int m_, int n_) : m(m_), n(n_), a(m_, vector(n_)) {} matrix(const vector>& v) : m(v.size()), n((v.size() == 0 or v[0].size() == 0) ? 0 : v[0].size()), a(v) {} //matrix(const matrix& other)=default; //~matrix()=default; vector& operator[](int i) { assert(0 <= i); assert(i < m); return a[i]; } matrix operator+() const { return *this; } matrix operator-() const { return K(-1) * *this; } K at(int i, int j) const { return a[i][j]; } matrix& operator=(const matrix& other) = default; matrix& operator+=(const matrix& other) { assert(m == other.m); assert(n == other.n); for (int i = 0; i < m; ++i) for (int j = 0; j < n; ++j) (*this)[i][j] += other[i][j]; return *this; } matrix& operator-=(const matrix& other) { return (*this) += -other; } matrix& operator*=(const matrix& other) { assert(n == other.m); vector> b(m, vector(other.n)); for (int i = 0; i < m; ++i) for (int j = 0; j < other.n; ++j) for (int k = 0; k < n; ++k) b[i][j] += (*this)[i][k] * other.at(k, j); a = b; n = other.n; return *this; } matrix& operator/=(const K& k) { assert(k != K(0)); for (int i = 0; i < m; ++i) for (int j = 0; j < n; ++j) (*this)[i][j] /= k; return *this; } matrix& operator%=(const K& k) { assert(k != K(0)); for (int i = 0; k < m; ++i) for (int j = 0; j < n; ++j) (*this)[i][j] = ((*this)[i][j] + k) % k; return *this; } matrix operator+(const matrix& other) const { return matrix(*this) += other; } matrix operator-(const matrix& other) const { return matrix(*this) -= other; } matrix operator*(const matrix& other) const { return matrix(*this) *= other; } matrix operator/(const K& k) const { return matrix(*this) /= k; } matrix operator%(const K& k) const { return matrix(*this) %= k; } bool operator==(const matrix& other) const { if (m != other.m or n != other.n) return false; for (int i = 0; i < m; ++i) for (int j = 0; j < n; ++j) if ((*this)[i][j] != other[i][j]) return false; return true; } bool operator!=(const matrix& other) const { return !((*this) == other); } static matrix E(int m) { assert(0 <= m); matrix E_(m, m); for (int i = 0; i < m; ++i) E_[i][i] = K(1); return E_; } static matrix e(int m, int i) { assert(0 <= i); assert(i < m); matrix e_; e_[i][0] = K(1); return e_; } }; // Kの単位元がK(1)で取得できる場合こっちを使用する template matrix rep_pow(matrix x, int n /* = xのサイズ */, int64_t y) { matrix res(matrix::E(n)); while (y > 0) { if (y & 1) res *= x; x *= x; y >>= 1; } return res; } void solve() { constexpr i64 mod = 1'000'000'007; i64 M,K; cin >> M >> K; matrix> A(M,M); rep(i,0,M){ rep(j,0,M){ A[i*j%M][i]+=1; A[(i+j)%M][i]+=1; } } auto Ak=rep_pow(A,M,K); matrix> v(M,1); v[0][0]=1; auto w=Ak*v; cout << w[0][0] << endl; } int main() { std::cin.tie(0); std::ios::sync_with_stdio(false); cout.setf(ios::fixed); cout.precision(16); solve(); return 0; }