結果

問題 No.2645 Sum of Divisors?
ユーザー kotatsugamekotatsugame
提出日時 2024-02-19 22:18:48
言語 C++14
(gcc 12.3.0 + boost 1.83.0)
結果
TLE  
実行時間 -
コード長 4,864 bytes
コンパイル時間 1,149 ms
コンパイル使用メモリ 89,556 KB
実行使用メモリ 198,264 KB
最終ジャッジ日時 2024-09-29 02:13:47
合計ジャッジ時間 7,150 ms
ジャッジサーバーID
（参考情報） 		judge3 / judge1
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 63 ms
167,032 KB
testcase_01 AC 61 ms
159,920 KB
testcase_02 TLE -
testcase_03 AC 63 ms
159,972 KB
testcase_04 AC 63 ms
159,824 KB
testcase_05 AC 63 ms
160,136 KB
testcase_06 TLE -
testcase_07 -- -
testcase_08 -- -
testcase_09 -- -
testcase_10 -- -
testcase_11 -- -
testcase_12 -- -
testcase_13 -- -
testcase_14 -- -
testcase_15 -- -
testcase_16 -- -
testcase_17 -- -
testcase_18 -- -
testcase_19 -- -
testcase_20 -- -
testcase_21 -- -
testcase_22 -- -
testcase_23 -- -
testcase_24 -- -
testcase_25 -- -
testcase_26 -- -
testcase_27 -- -
testcase_28 -- -
testcase_29 -- -
testcase_30 -- -
testcase_31 -- -
testcase_32 -- -
testcase_33 -- -
testcase_34 -- -
権限があれば一括ダウンロードができます

ソースコード

diff # 

#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=1e7;
long double H[LIM];

long double calc(long long x)
{
	if(x<LIM)return H[x];
	return (log(x)+0.5772156649+0.5/x-1.0/12/x/x+1.0/120/x/x/x/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{};
    }
  }

  // 邏謨ｰ縺ｮ蛟区焚髢｢謨ｰ縺ｫ髢｢縺吶ｋ繝・・繝悶Ν
  vector<T> pi_table() {
    if (M == 0) return {};
    i64 hls = md(M, sq);
    if (hls != 1 && md(M, hls - 1) == sq) hls--;

    vector<i64> hl(hls);
    for (int i = 1; i < hls; i++) hl[i] = md(M, i) - 1;

    vector<int> hs(sq + 1);
    iota(begin(hs), end(hs), -1);

    int pi = 0;
    for (auto& x : p) {
      i64 x2 = i64(x) * x;
      i64 imax = min<i64>(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<T> 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<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);
	  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<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, 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<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, 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
 */
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<LIM;i++)H[i]=H[i-1]+(long double)2/i;
	long long N;
	cin>>N;
	mf_prefix_sum<long double,f>solve(N);
	vector<long double>t=solve.prime_sum_table();
	cout<<fixed<<setprecision(16)<<solve.run(t)<<endl;
}
0