#include #include #include #include #include #include using namespace std; #line 2 "multiplicative-function/sum-of-multiplicative-function.hpp" #line 2 "prime/prime-enumerate.hpp" // 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; } #line 4 "multiplicative-function/sum-of-multiplicative-function.hpp" const int LIM=1e7; long double H[LIM]; long double calc(long long x) { if(x 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{}; } } // 邏謨ー縺ョ蛟区焚髢「謨ー縺ォ髢「縺吶k繝・・繝悶Ν 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 縺ォ髢「縺吶k繝・・繝悶Ν 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); for (int i = 1; i < hls; i++) { T x = md(M, i); h[i] = calc(x); log(x)+0.5772156649+0.5*x-1.0/12/x/x-1; } for (int i = 1; i <= sq; i++) { T x = i; h[s - i] = calc(x); } 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 */ long double f(long long p,long long c){ return(long double)(c+1)/pow((long double)p,c); } int main() { for(int i=2;i>N; mf_prefix_sumsolve(N); vectort=solve.prime_sum_table(); cout<