結果
問題 | No.2645 Sum of Divisors? |
ユーザー |
|
提出日時 | 2024-02-19 22:30:37 |
言語 | C++14 (gcc 13.3.0 + boost 1.87.0) |
結果 |
AC
|
実行時間 | 1,903 ms / 2,000 ms |
コード長 | 4,019 bytes |
コンパイル時間 | 808 ms |
コンパイル使用メモリ | 88,588 KB |
実行使用メモリ | 27,668 KB |
最終ジャッジ日時 | 2024-09-29 02:23:01 |
合計ジャッジ時間 | 13,270 ms |
ジャッジサーバーID (参考情報) |
judge4 / judge3 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
sample | AC * 3 |
other | AC * 32 |
ソースコード
#include<iostream>#include<iomanip>#include<cmath>#include<vector>#include<numeric>#include<cassert>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<int> prime_enumerate(int N) {vector<bool> 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<int> 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=(int)1e6+1;double H[LIM];double calc(long long x){if(x<LIM)return H[x];return (log(x)+0.5772156649+0.5/x-1)*2;}// f(p, c) : f(p^c) 縺ョ蛟、繧定ソ斐☆template <typename T, T (*f)(long long, long long)>struct mf_prefix_sum {using i64 = long long;i64 M, sq, s;vector<int> p;int ps;vector<T> 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{};}}// 邏謨ー縺ョ prefix sum 縺ォ髢「縺吶k繝・・繝悶Νvector<T> prime_sum_table() {if (M == 0) return {};i64 hls = md(M, sq);if (hls != 1 && md(M, hls - 1) == sq) hls--;vector<T> h(s);for (int i = 1; i < hls; i++) {T x = md(M, i);h[i] = calc(x);}for (int i = 1; i <= sq; i++) {T x = i;h[s - i] = calc(x);}for (auto& x : p) {T xt = (T)1/(T)x;T pi = h[s - x + 1];i64 x2 = i64(x) * x;i64 imax = min<i64>(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, T pc, i64 prod, T cur) {ans += cur * (c+2)*(pc/(T)p[i]);//f(p[i], c + 1);i64 lim = md(M, prod);if (lim >= 1LL * p[i] * p[i]) dfs(i, c + 1, pc/(T)p[i], p[i] * prod, cur);//cur *= f(p[i], c);cur *= (c+1)*pc;ans += cur * (buf[idx(lim)] - buf[idx(p[i])]);int j = i + 1;// M < 2**42 -> p_j < 2**21 -> (p_j)^3 < 2**63for (; j < ps && 1LL * p[j] * p[j] * p[j] <= lim; j++) {dfs(j, 1, (T)1/(T)p[j], prod * p[j], cur);}for (; j < ps && 1LL * p[j] * p[j] <= lim; j++) {T sm = (T)3/((T)p[j]*p[j]);//f(p[j], 2);int id1 = idx(md(lim, p[j])), id2 = idx(p[j]);sm += (T)2/(T)p[j] * (buf[id1] - buf[id2]);ans += cur * sm;}}// fprime 遐エ螢顔噪T run(vector<T>& 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, (T)1/(T)p[i], 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<T>& _buf) { swap(buf, _buf); }};/*** @brief 荵玲ウ慕噪髢「謨ー縺ョprefix sum* @docs docs/multiplicative-function/sum-of-multiplicative-function.md*/double f(long long p,long long c){assert(false);return(double)(c+1)/pow((double)p,c);}int main(){for(int i=2;i<LIM;i++)H[i]=H[i-1]+(double)2/i;long long N;cin>>N;mf_prefix_sum<double,f>solve(N);vector<double>t=solve.prime_sum_table();cout<<fixed<<setprecision(16)<<solve.run(t)<<endl;}