/** * date : 2023-04-07 23:53:12 */ #define NDEBUG using namespace std; // intrinstic #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include // utility namespace Nyaan { using ll = long long; using i64 = long long; using u64 = unsigned long long; using i128 = __int128_t; using u128 = __uint128_t; template using V = vector; template using VV = vector>; using vi = vector; using vl = vector; using vd = V; using vs = V; using vvi = vector>; using vvl = vector>; template struct P : pair { template P(Args... args) : pair(args...) {} using pair::first; using pair::second; P &operator+=(const P &r) { first += r.first; second += r.second; return *this; } P &operator-=(const P &r) { first -= r.first; second -= r.second; return *this; } P &operator*=(const P &r) { first *= r.first; second *= r.second; return *this; } template P &operator*=(const S &r) { first *= r, second *= r; return *this; } P operator+(const P &r) const { return P(*this) += r; } P operator-(const P &r) const { return P(*this) -= r; } P operator*(const P &r) const { return P(*this) *= r; } template P operator*(const S &r) const { return P(*this) *= r; } P operator-() const { return P{-first, -second}; } }; using pl = P; using pi = P; using vp = V; constexpr int inf = 1001001001; constexpr long long infLL = 4004004004004004004LL; template int sz(const T &t) { return t.size(); } template inline bool amin(T &x, U y) { return (y < x) ? (x = y, true) : false; } template inline bool amax(T &x, U y) { return (x < y) ? (x = y, true) : false; } template inline T Max(const vector &v) { return *max_element(begin(v), end(v)); } template inline T Min(const vector &v) { return *min_element(begin(v), end(v)); } template inline long long Sum(const vector &v) { return accumulate(begin(v), end(v), 0LL); } template int lb(const vector &v, const T &a) { return lower_bound(begin(v), end(v), a) - begin(v); } template int ub(const vector &v, const T &a) { return upper_bound(begin(v), end(v), a) - begin(v); } constexpr long long TEN(int n) { long long ret = 1, x = 10; for (; n; x *= x, n >>= 1) ret *= (n & 1 ? x : 1); return ret; } template pair mkp(const T &t, const U &u) { return make_pair(t, u); } template vector mkrui(const vector &v, bool rev = false) { vector ret(v.size() + 1); if (rev) { for (int i = int(v.size()) - 1; i >= 0; i--) ret[i] = v[i] + ret[i + 1]; } else { for (int i = 0; i < int(v.size()); i++) ret[i + 1] = ret[i] + v[i]; } return ret; }; template vector mkuni(const vector &v) { vector ret(v); sort(ret.begin(), ret.end()); ret.erase(unique(ret.begin(), ret.end()), ret.end()); return ret; } template vector mkord(int N,F f) { vector ord(N); iota(begin(ord), end(ord), 0); sort(begin(ord), end(ord), f); return ord; } template vector mkinv(vector &v) { int max_val = *max_element(begin(v), end(v)); vector inv(max_val + 1, -1); for (int i = 0; i < (int)v.size(); i++) inv[v[i]] = i; return inv; } vector mkiota(int n) { vector ret(n); iota(begin(ret), end(ret), 0); return ret; } template T mkrev(const T &v) { T w{v}; reverse(begin(w), end(w)); return w; } template bool nxp(vector &v) { return next_permutation(begin(v), end(v)); } template using minpq = priority_queue, greater>; } // namespace Nyaan // bit operation namespace Nyaan { __attribute__((target("popcnt"))) inline int popcnt(const u64 &a) { return _mm_popcnt_u64(a); } inline int lsb(const u64 &a) { return a ? __builtin_ctzll(a) : 64; } inline int ctz(const u64 &a) { return a ? __builtin_ctzll(a) : 64; } inline int msb(const u64 &a) { return a ? 63 - __builtin_clzll(a) : -1; } template inline int gbit(const T &a, int i) { return (a >> i) & 1; } template inline void sbit(T &a, int i, bool b) { if (gbit(a, i) != b) a ^= T(1) << i; } constexpr long long PW(int n) { return 1LL << n; } constexpr long long MSK(int n) { return (1LL << n) - 1; } } // namespace Nyaan // inout namespace Nyaan { template ostream &operator<<(ostream &os, const pair &p) { os << p.first << " " << p.second; return os; } template istream &operator>>(istream &is, pair &p) { is >> p.first >> p.second; return is; } template ostream &operator<<(ostream &os, const vector &v) { int s = (int)v.size(); for (int i = 0; i < s; i++) os << (i ? " " : "") << v[i]; return os; } template istream &operator>>(istream &is, vector &v) { for (auto &x : v) is >> x; return is; } istream &operator>>(istream &is, __int128_t &x) { string S; is >> S; x = 0; int flag = 0; for (auto &c : S) { if (c == '-') { flag = true; continue; } x *= 10; x += c - '0'; } if (flag) x = -x; return is; } istream &operator>>(istream &is, __uint128_t &x) { string S; is >> S; x = 0; for (auto &c : S) { x *= 10; x += c - '0'; } return is; } ostream &operator<<(ostream &os, __int128_t x) { if (x == 0) return os << 0; if (x < 0) os << '-', x = -x; string S; while (x) S.push_back('0' + x % 10), x /= 10; reverse(begin(S), end(S)); return os << S; } ostream &operator<<(ostream &os, __uint128_t x) { if (x == 0) return os << 0; string S; while (x) S.push_back('0' + x % 10), x /= 10; reverse(begin(S), end(S)); return os << S; } void in() {} template void in(T &t, U &...u) { cin >> t; in(u...); } void out() { cout << "\n"; } template void out(const T &t, const U &...u) { cout << t; if (sizeof...(u)) cout << sep; out(u...); } struct IoSetupNya { IoSetupNya() { cin.tie(nullptr); ios::sync_with_stdio(false); cout << fixed << setprecision(15); cerr << fixed << setprecision(7); } } iosetupnya; } // namespace Nyaan // debug #ifdef NyaanDebug #define trc(...) (void(0)) #else #define trc(...) (void(0)) #endif #ifdef NyaanLocal #define trc2(...) (void(0)) #else #define trc2(...) (void(0)) #endif // macro #define each(x, v) for (auto&& x : v) #define each2(x, y, v) for (auto&& [x, y] : v) #define all(v) (v).begin(), (v).end() #define rep(i, N) for (long long i = 0; i < (long long)(N); i++) #define repr(i, N) for (long long i = (long long)(N)-1; i >= 0; i--) #define rep1(i, N) for (long long i = 1; i <= (long long)(N); i++) #define repr1(i, N) for (long long i = (N); (long long)(i) > 0; i--) #define reg(i, a, b) for (long long i = (a); i < (b); i++) #define regr(i, a, b) for (long long i = (b)-1; i >= (a); i--) #define fi first #define se second #define ini(...) \ int __VA_ARGS__; \ in(__VA_ARGS__) #define inl(...) \ long long __VA_ARGS__; \ in(__VA_ARGS__) #define ins(...) \ string __VA_ARGS__; \ in(__VA_ARGS__) #define in2(s, t) \ for (int i = 0; i < (int)s.size(); i++) { \ in(s[i], t[i]); \ } #define in3(s, t, u) \ for (int i = 0; i < (int)s.size(); i++) { \ in(s[i], t[i], u[i]); \ } #define in4(s, t, u, v) \ for (int i = 0; i < (int)s.size(); i++) { \ in(s[i], t[i], u[i], v[i]); \ } #define die(...) \ do { \ Nyaan::out(__VA_ARGS__); \ return; \ } while (0) namespace Nyaan { void solve(); } int main() { Nyaan::solve(); } // #ifdef _MSC_VER #include #endif namespace atcoder { namespace internal { // @param m `1 <= m` // @return x mod m constexpr long long safe_mod(long long x, long long m) { x %= m; if (x < 0) x += m; return x; } // Fast modular multiplication by barrett reduction // Reference: https://en.wikipedia.org/wiki/Barrett_reduction // NOTE: reconsider after Ice Lake struct barrett { unsigned int _m; unsigned long long im; // @param m `1 <= m < 2^31` barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {} // @return m unsigned int umod() const { return _m; } // @param a `0 <= a < m` // @param b `0 <= b < m` // @return `a * b % m` unsigned int mul(unsigned int a, unsigned int b) const { // [1] m = 1 // a = b = im = 0, so okay // [2] m >= 2 // im = ceil(2^64 / m) // -> im * m = 2^64 + r (0 <= r < m) // let z = a*b = c*m + d (0 <= c, d < m) // a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im // c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2 // ((ab * im) >> 64) == c or c + 1 unsigned long long z = a; z *= b; #ifdef _MSC_VER unsigned long long x; _umul128(z, im, &x); #else unsigned long long x = (unsigned long long)(((unsigned __int128)(z)*im) >> 64); #endif unsigned int v = (unsigned int)(z - x * _m); if (_m <= v) v += _m; return v; } }; // @param n `0 <= n` // @param m `1 <= m` // @return `(x ** n) % m` constexpr long long pow_mod_constexpr(long long x, long long n, int m) { if (m == 1) return 0; unsigned int _m = (unsigned int)(m); unsigned long long r = 1; unsigned long long y = safe_mod(x, m); while (n) { if (n & 1) r = (r * y) % _m; y = (y * y) % _m; n >>= 1; } return r; } // Reference: // M. Forisek and J. Jancina, // Fast Primality Testing for Integers That Fit into a Machine Word // @param n `0 <= n` constexpr bool is_prime_constexpr(int n) { if (n <= 1) return false; if (n == 2 || n == 7 || n == 61) return true; if (n % 2 == 0) return false; long long d = n - 1; while (d % 2 == 0) d /= 2; constexpr long long bases[3] = {2, 7, 61}; for (long long a : bases) { long long t = d; long long y = pow_mod_constexpr(a, t, n); while (t != n - 1 && y != 1 && y != n - 1) { y = y * y % n; t <<= 1; } if (y != n - 1 && t % 2 == 0) { return false; } } return true; } template constexpr bool is_prime = is_prime_constexpr(n); // @param b `1 <= b` // @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g constexpr std::pair inv_gcd(long long a, long long b) { a = safe_mod(a, b); if (a == 0) return {b, 0}; // Contracts: // [1] s - m0 * a = 0 (mod b) // [2] t - m1 * a = 0 (mod b) // [3] s * |m1| + t * |m0| <= b long long s = b, t = a; long long m0 = 0, m1 = 1; while (t) { long long u = s / t; s -= t * u; m0 -= m1 * u; // |m1 * u| <= |m1| * s <= b // [3]: // (s - t * u) * |m1| + t * |m0 - m1 * u| // <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u) // = s * |m1| + t * |m0| <= b auto tmp = s; s = t; t = tmp; tmp = m0; m0 = m1; m1 = tmp; } // by [3]: |m0| <= b/g // by g != b: |m0| < b/g if (m0 < 0) m0 += b / s; return {s, m0}; } // Compile time primitive root // @param m must be prime // @return primitive root (and minimum in now) constexpr int primitive_root_constexpr(int m) { if (m == 2) return 1; if (m == 167772161) return 3; if (m == 469762049) return 3; if (m == 754974721) return 11; if (m == 998244353) return 3; int divs[20] = {}; divs[0] = 2; int cnt = 1; int x = (m - 1) / 2; while (x % 2 == 0) x /= 2; for (int i = 3; (long long)(i)*i <= x; i += 2) { if (x % i == 0) { divs[cnt++] = i; while (x % i == 0) { x /= i; } } } if (x > 1) { divs[cnt++] = x; } for (int g = 2;; g++) { bool ok = true; for (int i = 0; i < cnt; i++) { if (pow_mod_constexpr(g, (m - 1) / divs[i], m) == 1) { ok = false; break; } } if (ok) return g; } } template constexpr int primitive_root = primitive_root_constexpr(m); } // namespace internal } // namespace atcoder namespace atcoder { long long pow_mod(long long x, long long n, int m) { assert(0 <= n && 1 <= m); if (m == 1) return 0; internal::barrett bt((unsigned int)(m)); unsigned int r = 1, y = (unsigned int)(internal::safe_mod(x, m)); while (n) { if (n & 1) r = bt.mul(r, y); y = bt.mul(y, y); n >>= 1; } return r; } long long inv_mod(long long x, long long m) { assert(1 <= m); auto z = internal::inv_gcd(x, m); assert(z.first == 1); return z.second; } // (rem, mod) std::pair crt(const std::vector& r, const std::vector& m) { assert(r.size() == m.size()); int n = int(r.size()); // Contracts: 0 <= r0 < m0 long long r0 = 0, m0 = 1; for (int i = 0; i < n; i++) { assert(1 <= m[i]); long long r1 = internal::safe_mod(r[i], m[i]), m1 = m[i]; if (m0 < m1) { std::swap(r0, r1); std::swap(m0, m1); } if (m0 % m1 == 0) { if (r0 % m1 != r1) return {0, 0}; continue; } // assume: m0 > m1, lcm(m0, m1) >= 2 * max(m0, m1) // (r0, m0), (r1, m1) -> (r2, m2 = lcm(m0, m1)); // r2 % m0 = r0 // r2 % m1 = r1 // -> (r0 + x*m0) % m1 = r1 // -> x*u0*g % (u1*g) = (r1 - r0) (u0*g = m0, u1*g = m1) // -> x = (r1 - r0) / g * inv(u0) (mod u1) // im = inv(u0) (mod u1) (0 <= im < u1) long long g, im; std::tie(g, im) = internal::inv_gcd(m0, m1); long long u1 = (m1 / g); // |r1 - r0| < (m0 + m1) <= lcm(m0, m1) if ((r1 - r0) % g) return {0, 0}; // u1 * u1 <= m1 * m1 / g / g <= m0 * m1 / g = lcm(m0, m1) long long x = (r1 - r0) / g % u1 * im % u1; // |r0| + |m0 * x| // < m0 + m0 * (u1 - 1) // = m0 + m0 * m1 / g - m0 // = lcm(m0, m1) r0 += x * m0; m0 *= u1; // -> lcm(m0, m1) if (r0 < 0) r0 += m0; } return {r0, m0}; } long long floor_sum(long long n, long long m, long long a, long long b) { long long ans = 0; if (a >= m) { ans += (n - 1) * n * (a / m) / 2; a %= m; } if (b >= m) { ans += n * (b / m); b %= m; } long long y_max = (a * n + b) / m, x_max = (y_max * m - b); if (y_max == 0) return ans; ans += (n - (x_max + a - 1) / a) * y_max; ans += floor_sum(y_max, a, m, (a - x_max % a) % a); return ans; } } // namespace atcoder // { (q, l, r) : forall x in (l,r], floor(N/x) = q } // を引数に取る関数f(q, l, r)を渡す。範囲が左に半開なのに注意 template void enumerate_quotient(T N, const F& f) { T sq = sqrt(N), upper = N, quo = 0; while (upper > sq) { T thres = N / (++quo + 1); f(quo, thres, upper); upper = thres; } while (upper > 0) { f(N / upper, upper - 1, upper); upper--; } } /** * @brief 商の列挙 */ struct Rational { using R = Rational; using i128 = __int128_t; using i64 = long long; using u64 = unsigned long long; long long x, y; Rational() : x(0), y(1) {} Rational(long long _x, long long _y = 1) : x(_x), y(_y) { assert(y != 0); if (_y != 1) { long long g = gcd(x, y); if (g != 0) x /= g, y /= g; if (y < 0) x = -x, y = -y; } } u64 gcd(i64 A, i64 B) { u64 a = A >= 0 ? A : -A; u64 b = B >= 0 ? B : -B; if (a == 0 || b == 0) return a + b; int n = __builtin_ctzll(a); int m = __builtin_ctzll(b); a >>= n; b >>= m; while (a != b) { int d = __builtin_ctzll(a - b); bool f = a > b; u64 c = f ? a : b; b = f ? b : a; a = (c - b) >> d; } return a << min(n, m); } friend R operator+(const R& l, const R& r) { return R(l.x * r.y + l.y * r.x, l.y * r.y); } friend R operator-(const R& l, const R& r) { return R(l.x * r.y - l.y * r.x, l.y * r.y); } friend R operator*(const R& l, const R& r) { return R(l.x * r.x, l.y * r.y); } friend R operator/(const R& l, const R& r) { assert(r.x != 0); return R(l.x * r.y, l.y * r.x); } R& operator+=(const R& r) { return (*this) = (*this) + r; } R& operator-=(const R& r) { return (*this) = (*this) - r; } R& operator*=(const R& r) { return (*this) = (*this) * r; } R& operator/=(const R& r) { return (*this) = (*this) / r; } R operator-() const { R r; r.x = -x, r.y = y; return r; } R inverse() const { assert(x != 0); R r; r.x = y, r.y = x; if (x < 0) r.x = -r.x, r.y = -r.y; return r; } R pow(long long p) const { R res(1), base(*this); while (p) { if (p & 1) res *= base; base *= base; p >>= 1; } return res; } friend bool operator==(const R& l, const R& r) { return l.x == r.x && l.y == r.y; }; friend bool operator!=(const R& l, const R& r) { return l.x != r.x || l.y != r.y; }; friend bool operator<(const R& l, const R& r) { return i128(l.x) * r.y < i128(l.y) * r.x; }; friend bool operator<=(const R& l, const R& r) { return l < r || l == r; } friend bool operator>(const R& l, const R& r) { return i128(l.x) * r.y > i128(l.y) * r.x; }; friend bool operator>=(const R& l, const R& r) { return l > r || l == r; } friend ostream& operator<<(ostream& os, const R& r) { os << r.x; if (r.x != 0 && r.y != 1) os << "/" << r.y; return os; } long long toMint(long long mod) { assert(mod != 0); i64 a = y, b = mod, u = 1, v = 0, t; while (b > 0) { t = a / b; swap(a -= t * b, b); swap(u -= t * v, v); } return i128((u % mod + mod) % mod) * x % mod; } }; template struct Binomial { vector fc; Binomial(int = 0) { fc.emplace_back(1); } void extend() { int n = fc.size(); R nxt = fc.back() * n; fc.push_back(nxt); } R fac(int n) { while ((int)fc.size() <= n) extend(); return fc[n]; } R finv(int n) { return fac(n).inverse(); } R inv(int n) { return R{1, max(n, 1)}; } R C(int n, int r) { if (n < 0 or r < 0 or n < r) return R{0}; return fac(n) * finv(n - r) * finv(r); } R operator()(int n, int r) { return C(n, r); } template R multinomial(const vector& r) { static_assert(is_integral::value == true); int n = 0; for (auto& x : r) { if (x < 0) return R{0}; n += x; } R res = fac(n); for (auto& x : r) res *= finv(x); return res; } template R operator()(const vector& r) { return multinomial(r); } }; // Prime Sieve {2, 3, 5, 7, 11, 13, 17, ...} vector prime_enumerate(int N) { vector sieve(N / 3 + 1, 1); for (int p = 5, d = 4, i = 1, sqn = sqrt(N); p <= sqn; p += d = 6 - d, i++) { if (!sieve[i]) continue; for (int q = p * p / 3, r = d * p / 3 + (d * p % 3 == 2), s = 2 * p, qe = sieve.size(); q < qe; q += r = s - r) sieve[q] = 0; } vector ret{2, 3}; for (int p = 5, d = 4, i = 1; p <= N; p += d = 6 - d, i++) if (sieve[i]) ret.push_back(p); while (!ret.empty() && ret.back() > N) ret.pop_back(); return ret; } struct divisor_transform { template static void zeta_transform(vector &a) { int N = a.size() - 1; auto sieve = prime_enumerate(N); for (auto &p : sieve) for (int k = 1; k * p <= N; ++k) a[k * p] += a[k]; } template static void mobius_transform(T &a) { int N = a.size() - 1; auto sieve = prime_enumerate(N); for (auto &p : sieve) for (int k = N / p; k > 0; --k) a[k * p] -= a[k]; } template static void zeta_transform(map &a) { for (auto p = rbegin(a); p != rend(a); p++) for (auto &x : a) { if (p->first == x.first) break; if (p->first % x.first == 0) p->second += x.second; } } template static void mobius_transform(map &a) { for (auto &x : a) { for (auto p = rbegin(a); p != rend(a); p++) { if (x.first == p->first) break; if (p->first % x.first == 0) p->second -= x.second; } } } }; struct multiple_transform { template static void zeta_transform(vector &a) { int N = a.size() - 1; auto sieve = prime_enumerate(N); for (auto &p : sieve) for (int k = N / p; k > 0; --k) a[k] += a[k * p]; } template static void mobius_transform(vector &a) { int N = a.size() - 1; auto sieve = prime_enumerate(N); for (auto &p : sieve) for (int k = 1; k * p <= N; ++k) a[k] -= a[k * p]; } template static void zeta_transform(map &a) { for (auto &x : a) for (auto p = rbegin(a); p->first != x.first; p++) if (p->first % x.first == 0) x.second += p->second; } template static void mobius_transform(map &a) { for (auto p1 = rbegin(a); p1 != rend(a); p1++) for (auto p2 = rbegin(a); p2 != p1; p2++) if (p2->first % p1->first == 0) p1->second -= p2->second; } }; /** * @brief 倍数変換・約数変換 * @docs docs/multiplicative-function/divisor-multiple-transform.md */ // f(p, c) : f(p^c) の値を返す template struct mf_prefix_sum { using i64 = long long; i64 M, sq, s; vector p; int ps; vector buf; T ans; mf_prefix_sum(i64 m) : M(m) { assert(m < (1LL << 42)); sq = sqrt(M); while (sq * sq > M) sq--; while ((sq + 1) * (sq + 1) <= M) sq++; if (M != 0) { i64 hls = md(M, sq); if (hls != 1 && md(M, hls - 1) == sq) hls--; s = hls + sq; p = prime_enumerate(sq); ps = p.size(); ans = T{}; } } // 素数の個数関数に関するテーブル vector pi_table() { if (M == 0) return {}; i64 hls = md(M, sq); if (hls != 1 && md(M, hls - 1) == sq) hls--; vector hl(hls); for (int i = 1; i < hls; i++) hl[i] = md(M, i) - 1; vector hs(sq + 1); iota(begin(hs), end(hs), -1); int pi = 0; for (auto& x : p) { i64 x2 = i64(x) * x; i64 imax = min(hls, md(M, x2) + 1); for (i64 i = 1, ix = x; i < imax; ++i, ix += x) { hl[i] -= (ix < hls ? hl[ix] : hs[md(M, ix)]) - pi; } for (int n = sq; n >= x2; n--) hs[n] -= hs[md(n, x)] - pi; pi++; } vector res; res.reserve(2 * sq + 10); for (auto& x : hl) res.push_back(x); for (int i = hs.size(); --i;) res.push_back(hs[i]); assert((int)res.size() == s); return res; } // 素数の prefix sum に関するテーブル vector prime_sum_table() { if (M == 0) return {}; i64 hls = md(M, sq); if (hls != 1 && md(M, hls - 1) == sq) hls--; vector h(s); T inv2 = T{2}.inverse(); for (int i = 1; i < hls; i++) { T x = md(M, i); h[i] = x * (x + 1) * inv2 - 1; } for (int i = 1; i <= sq; i++) { T x = i; h[s - i] = x * (x + 1) / 2 - 1; } for (auto& x : p) { T xt = x; T pi = h[s - x + 1]; i64 x2 = i64(x) * x; i64 imax = min(hls, md(M, x2) + 1); i64 ix = x; for (i64 i = 1; i < imax; ++i, ix += x) { h[i] -= ((ix < hls ? h[ix] : h[s - md(M, ix)]) - pi) * xt; } for (int n = sq; n >= x2; n--) { h[s - n] -= (h[s - md(n, x)] - pi) * xt; } } assert((int)h.size() == s); return h; } void dfs(int i, int c, i64 prod, T cur) { ans += cur * f(p[i], c + 1); i64 lim = md(M, prod); if (lim >= 1LL * p[i] * p[i]) dfs(i, c + 1, p[i] * prod, cur); cur *= f(p[i], c); ans += cur * (buf[idx(lim)] - buf[idx(p[i])]); int j = i + 1; // M < 2**42 -> p_j < 2**21 -> (p_j)^3 < 2**63 for (; j < ps && 1LL * p[j] * p[j] * p[j] <= lim; j++) { dfs(j, 1, prod * p[j], cur); } for (; j < ps && 1LL * p[j] * p[j] <= lim; j++) { T sm = f(p[j], 2); int id1 = idx(md(lim, p[j])), id2 = idx(p[j]); sm += f(p[j], 1) * (buf[id1] - buf[id2]); ans += cur * sm; } } // fprime 破壊的 T run(vector& fprime) { if (M == 0) return {}; set_buf(fprime); assert((int)buf.size() == s); ans = buf[idx(M)] + 1; for (int i = 0; i < ps; i++) dfs(i, 1, p[i], 1); return ans; } i64 md(i64 n, i64 d) { return double(n) / d; } i64 idx(i64 n) { return n <= sq ? s - n : md(M, n); } void set_buf(vector& _buf) { swap(buf, _buf); } }; /** * @brief 乗法的関数のprefix sum * @docs docs/multiplicative-function/sum-of-multiplicative-function.md */ // f(n, p, c) : n = pow(p, c), f is multiplicative function template struct enamurate_multiplicative_function { enamurate_multiplicative_function(int _n) : ps(prime_enumerate(_n)), a(_n + 1, T()), n(_n), p(ps.size()) {} vector run() { a[1] = 1; dfs(-1, 1, 1); return a; } private: vector ps; vector a; int n, p; void dfs(int i, long long x, T y) { a[x] = y; if (y == T()) return; for (int j = i + 1; j < p; j++) { long long nx = x * ps[j]; long long dx = ps[j]; if (nx > n) break; for (int c = 1; nx <= n; nx *= ps[j], dx *= ps[j], ++c) { dfs(j, nx, y * f(dx, ps[j], c)); } } } }; /** * @brief 乗法的関数の列挙 */ namespace multiplicative_function { template T moebius(int, int, int c) { return c == 0 ? 1 : c == 1 ? -1 : 0; } template T sigma0(int, int, int c) { return c + 1; } template T sigma1(int n, int p, int) { return (n - 1) / (p - 1) + n; } template T totient(int n, int p, int) { return n - n / p; } } // namespace multiplicative_function template static constexpr vector mobius_function(int n) { enamurate_multiplicative_function> em( n); return em.run(); } template static constexpr vector sigma0(int n) { enamurate_multiplicative_function> em( n); return em.run(); } template static constexpr vector sigma1(int n) { enamurate_multiplicative_function> em( n); return em.run(); } template static constexpr vector totient(int n) { enamurate_multiplicative_function> em( n); return em.run(); } /** * @brief 有名な乗法的関数 * @docs docs/multiplicative-function/mf-famous-series.md */ template T sum_of_totient(long long N) { if (N < 2) return N; using i64 = long long; auto f = [](i64 v, i64 p, i64) -> i64 { return v / p * (p - 1); }; vector ns{0}, p; for (i64 i = N; i > 0; i = N / (N / i + 1)) ns.push_back(i); i64 s = ns.size(), sq = sqrt(N); auto idx = [&](i64 n) { return n <= sq ? s - n : N / n; }; vector h0(s), h1(s), buf(s); for (int i = 0; i < s; i++) { T x = ns[i]; h0[i] = x - 1; h1[i] = x * (x + 1) / 2 - 1; } for (i64 x = 2; x <= sq; ++x) { if (h0[s - x] == h0[s - x + 1]) continue; p.push_back(x); i64 x2 = x * x; for (i64 i = 1, n = ns[i]; i < s && n >= x2; n = ns[++i]) { int id = (i * x <= sq ? i * x : s - n / x); h0[i] -= h0[id] - h0[s - x + 1]; h1[i] -= (h1[id] - h1[s - x + 1]) * x; } } for (int i = 0; i < s; i++) buf[i] = h1[i] - h0[i]; T ans = buf[idx(N)] + 1; auto dfs = [&](auto rec, int i, int c, i64 v, i64 lim, T cur) -> void { ans += cur * f(p[i] * v, p[i], c + 1); if (lim >= p[i] * p[i]) rec(rec, i, c + 1, p[i] * v, lim / p[i], cur); cur *= f(v, p[i], c); ans += cur * (buf[idx(lim)] - buf[idx(p[i])]); for (int j = i + 1; j < (int)p.size() && p[j] * p[j] <= lim; j++) { rec(rec, j, 1, p[j], lim / p[j], cur); } }; for (int i = 0; i < (int)p.size(); i++) dfs(dfs, i, 1, p[i], N / p[i], 1); return ans; } /** * @brief トーシェント関数の和 */ using namespace Nyaan; V mo; void precalc(ll) { if (mo.empty()) { mo = mobius_function(TEN(8) + 200); rep(i, sz(mo) - 1) mo[i + 1] += mo[i]; } } // k 番目に小さい pl calc(ll N, ll K) { precalc(N); auto cnt = [&](Rational f) -> ll { ll s = 0; enumerate_quotient(N, [&](ll q, ll l, ll r) { ll x = 0; x += atcoder::floor_sum(r + 1, f.y, f.x, 0); x -= atcoder::floor_sum(l + 1, f.y, f.x, 0); s += x * mo[q]; }); /* each(i, mop) { if (i > N) break; s += i * ll(f.x) / int(f.y); } */ trc2(f, s); return s; }; Rational L{0, 1}; Rational M{1, 2}; Rational R{1, 1}; while (true) { // trc2(L.x, L.y, M.x, M.y, R.x, R.y); ll c = cnt(M); if (c == K) { break; } if (c < K) { for (ll i = 1;; i *= 2) { Rational f{L.x + R.x * i, L.y + R.y * i}; if (max(f.x, f.y) > N) break; if (cnt(f) == K) return {f.x, f.y}; if (cnt(f) < K) { L = f; } else { break; } } } else { for (ll i = 1;; i *= 2) { Rational f{L.x * i + R.x, L.y * i + R.y}; if (max(f.x, f.y) > N) break; if (cnt(f) == K) return {f.x, f.y}; if (cnt(f) >= K) { R = f; } else { break; } } } M = Rational{L.x + R.x, L.y + R.y}; } return {M.x, M.y}; } void q() { inl(N, K); ll s = sum_of_totient(N) - 1; trc(s); ll p = -1, q = -1; if (K <= s) { tie(p, q) = calc(N, K); } else if (K == s + 1) { p = q = 1; } else if (K <= s * 2 + 1) { tie(q, p) = calc(N, 2 * s + 1 - (K - 1)); } else { // do nothing } if (p == -1) { out(-1); } else { cout << p << "/" << q << "\n"; } } void Nyaan::solve() { int t = 1; // in(t); while (t--) q(); }