結果

問題 No.2829 GCD Divination
コンテスト
ユーザー Harui
提出日時 2024-08-03 11:01:59
言語 C++23
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 44 ms / 2,000 ms
コード長 3,716 bytes
コンパイル時間 1,559 ms
コンパイル使用メモリ 113,260 KB
実行使用メモリ 6,944 KB
最終ジャッジ日時 2024-08-03 11:02:02
合計ジャッジ時間 3,435 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 35
権限があれば一括ダウンロードができます

ソースコード

diff # 

#line 1 "playground_A.cpp"
#include <iostream>
#include <iomanip>
#include <cassert>
#include <map>
using namespace std;
using ll = long long;
template<class T> inline bool chmax(T& a, const T& b) {if (a<b) {a=b; return true;} return false;}
template<class T> inline bool chmin(T& a, const T& b) {if (b<a) {a=b; return true;} return false;}
const int INTINF = 1000001000;
const int INTMAX = 2147483647;
const ll LLMAX = 9223372036854775807;
const ll LLINF = 1000000000000000000;

#line 1 "/home/samejima/CompetitiveProgramming/library/convolution/multiple-zeta-moebius-transform.hpp"



#include <vector>

namespace multiple {

  // 倍数についてのゼータ変換。 g_n = \Sigma_{n|m} f_m なる g を求める。
  // n|mというのは、m%n==0という意味。
  // O(N log N) (調和級数)
  // うまくやるとO(Nlog(log(N)))にできることがよく知られているが、難しいしlogは定数なので妥協する。
  template <typename T, T (*op)(T, T)>
  std::vector<T> zeta_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;
    for (int i = 1; i <= N; i++) {
      for (int j = 2 * i; j <= N; j += i) {
        g[i] = op(g[i], f[j]);
      }
    }

    return g;
  }

  // 倍数についてのメビウス変換
  // f_n = \Sigma_{n|m} g_m なる g を求める。
  // O(N log N) (調和級数)
  template <typename T, T(*invop)(T, T)>
  std::vector<T> moebius_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;
    for (int i = N; i >= 1; i--) {
      for (int j = 2 * i; j <= N; j += i) {
        g[i] = invop(g[i], g[j]);
      }
    }
    return g;
  }


  template <typename I, typename T, T(*op)(T, T)>
  std::map<I, T> zeta_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (std::pair<I, T> pit : ret) {
      for (auto p2itr = ret.rbegin(); (*p2itr).first != pit.first; p2itr++) {
        if ((*p2itr).first % pit.first == 0) {
          ret[pit.first] = op(ret[pit.first], (*p2itr).second);
        }
      }
    }

    return ret;
  }


  template <typename I, typename T, T (*invop)(T, T)>
  std::map<I, T> moebius_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (auto p1itr = ret.rbegin(); p1itr != ret.rend(); p1itr++) {
      for (auto p2itr = ret.rbegin(); p2itr != p1itr; p2itr++) {
        if ((*p2itr).first % (*p1itr).first == 0) {
          (*p1itr).second = invop((*p1itr).second, (*p2itr).second);
        }
      }
    }

    return ret;
  }

} // namespace multiple


#line 15 "playground_A.cpp"

ll llpow(ll a, ll power) {
  ll ret =1;
  ll base =a;
  while (power > 0) {
    if (power & 1) {
      ret *= base;
    }

    base *= base;
    power >>= 1;
  }
  return ret;
}

ll op(ll a, ll b) {return a+b;}
ll invop(ll a, ll b) {return a-b;}

int main() {
  ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
  ll N; cin >> N;
  map<ll,double> memo;

  memo[1] = 0.0;

  auto solve = [&memo] (auto self, ll n) -> double {
    if (memo.contains(n)) return memo.at(n);

    double divsum = 0;

    map<ll,ll> g;
    for (int d=1; d*d <= n; d++) {
      if (n % d == 0) {
        g[d] = n/d;
      
        if (d != n/d) {
          g[n/d] = d;
        }
      }
    }
    map<ll,ll>f = multiple::moebius_transform<ll,ll,invop>(g);

    for (int d=2; d*d <= n; d++) {
      if (n % d == 0) {
        divsum += self(self, d) * f[d];
      
        if (d != n/d) {
          divsum += self(self, n/d) * f[n/d];
        }
      }
    }

    double ret = double(n) / double(n-1) + divsum / double(n-1);
    memo[n] = ret;
    return ret;
  };
  cout << fixed << setprecision(10);
  cout << solve(solve,N) << endl;
  return 0;
}
0