// warm heart, wagging tail,and a smile just for you! // // ▒█████▒▒ // ██████████▒ // ▒████████████▒ // ██████████████████ // ████████████████████▒ // ▒██████████████████████▒ // ▒███████████████████████ // ▒████▒▒▒▒▒▒█████████████████▒ // ███▒▒▒▒▒▒██████████████████████▒▒▒ // ▒██▒▒███████████████████████▒▒▒▒▒██████ // ▒█████████████████████████▒▒▒▒▒▒█████████▒ // ▒█████████████████████▒▒▒▒▒▒██████████████ // ▒████ ████▒▒▒▒▒████ ████▒ // ▒█████▒ ████ ▒▒▒▒███████ ████ ██████▒ // ▒██▒▒▒▒▒ ██████ █████████ ██████ ██▒▒▒██▒ // █████████ ████████ █████████ ████████ ▒▒▒▒█████ // ▒█████████ ██████ ████████▒ ██████ █████████ // ▒██████████ ████ █████▒▒▒▒▒▒ ████ ██████████ // ████████████ ▒▒▒▒▒▒▒████████ ███████████▒ // ▒██████████▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒███████████████████████████████████▒ // ███▒▒▒▒▒▒▒▒▒▒▒▒█████████████████████████████████████████▒▒████████▒ // ▒▒▒▒▒▒▒▒▒██████████████ ███████▒▒▒▒███████████ // █████████████████████████ ███████▒▒▒██████████████▒ // █████████████████████████████ ███████▒▒▒██████████████████▒ // ██████████████████████████████████████████████████████████████████████ // ██████████████████████████████████████████████████████████████████▒ // ▒█████████████████▒▒▒▒▒▒▒██████████████████████████████████▒▒▒ // #include "bits/stdc++.h" using namespace std; #define MOD 1000000007 //#define MOD 998244353 const double EPS = 1e-9; #define INF (1LL<<60) #define D double #define fs first #define sc second #define int long long #define FOR(i,a,b) for(int i=(a);i<(b);++i) #define RFOR(i,a,b) for(int i = (b-1);i>=(a);--i) #define REP(i,n) FOR(i,0,(n)) #define RREP(i,n) RFOR(i,0,(n)) #define ITR(itr,mp) for(auto itr = (mp).begin(); itr != (mp).end(); ++itr) #define RITR(itr,mp) for(auto itr = (mp).rbegin(); itr != (mp).rend(); ++itr) #define range(i,a,b) ((a)<=(i) && (i)<(b)) #define debug(x) cout << #x << " = " << (x) << endl; #define SP << " " << typedef pair P; #include template class modint { using u64 = std::uint_fast64_t; u64 a; public: constexpr modint(const u64 x = 0) noexcept : a(x % Modulus) {} constexpr u64 val() const noexcept { return a; } constexpr modint operator+(const modint rhs) const noexcept { return modint(*this) += rhs; } constexpr modint operator-(const modint rhs) const noexcept { return modint(*this) -= rhs; } constexpr modint operator*(const modint rhs) const noexcept { return modint(*this) *= rhs; } constexpr modint operator/(const modint rhs) const noexcept { return modint(*this) /= rhs; } constexpr bool operator==(const modint rhs) const noexcept { return modint(*this).val() == rhs.val(); } modint &operator+=(const modint rhs) noexcept { a += rhs.a; if (a >= Modulus) { a -= Modulus; } return *this; } modint &operator-=(const modint rhs) noexcept { if (a < rhs.a) { a += Modulus; } a -= rhs.a; return *this; } modint &operator*=(const modint rhs) noexcept { a = a * rhs.a % Modulus; return *this; } modint &operator/=(modint rhs) noexcept { u64 exp = Modulus - 2; while (exp) { if (exp % 2) { *this *= rhs; } rhs *= rhs; exp /= 2; } return *this; } }; using mint = modint; typedef vector vec; typedef vector> mat; int m; vec matmul(vec &dp, mat &mt){ vec ret(m,0); REP(i,m) REP(j,m) ret[i] += mt[i][j]*dp[j]; return ret; } mat update(mat &mt){ mat ret(m,vec(m,0)); REP(i,m) REP(j,m) REP(k,m) ret[i][j] += mt[i][k]*mt[k][j]; return ret; } void matpow(vec &dp, mat &mt, int k){ m = dp.size(); while(k){ if(k%2) dp = matmul(dp,mt); mt = update(mt); k /= 2; } } vector divisor(const int n){ vector ret; for(int i=1;i*i<=n;i++){ if(n % i == 0){ ret.push_back(i); if(i*i!= n) ret.push_back(n/i); } } sort(ret.begin(),ret.end()); return ret; } signed main(){ ios::sync_with_stdio(false); cin.tie(0); int n,k; cin >> n >> k; int g = __gcd(n,k); vector list; for(int i=2;i<=n;i+=2){ bool flag = true; for(int j=3; j<=sqrt(i);j+=2){ if(i%j == 0){ flag = false; break; } } if(flag || i == 2){ int cnt = 0; while(g%i==0) g /= i, cnt++; if(cnt) g *= i,list.push_back(i); } if(i == 2) i--; } vector a = divisor(g); vec fact(n+1,1); REP(i,n) fact[i+1] = fact[i]*(i+1); mint ans = 0; FOR(i,1,a.size()){ int cnt = 0; REP(j,list.size()) if(a[i]%list[j]==0) cnt++; if(cnt%2) ans += fact[n/a[i]]/fact[n/a[i]-k/a[i]]/fact[k/a[i]]; else ans -= fact[n/a[i]]/fact[n/a[i]-k/a[i]]/fact[k/a[i]]; } cout << ans.val() << endl; return 0; }