結果

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

ソースコード

diff # 

#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <map>
#include <cassert>

using namespace std;

vector<int> divisors (int N) {
    vector<int> res;
    for (int i = 1; i <= N; i++) {
        if (N < 1LL * i * i) break;
        if (N % i == 0) {
            res.push_back(i);
            if (1LL * i * i < N) res.push_back(N / i);
        }
    }

    sort(res.begin(), res.end());
    return res;
}

int main () {
    int N; cin >> N;
    // いつものお気持ちdpで立式は行ける。あとはあるxに対してE(gcd(x, v))の総和をとるパートをどうにかしたい。
    // gcd(x, v)の値で主客転倒して考えてみると、約数系除原理でgcd(x, v) = xからはじめて後ろ向きに求まる。
    // これであるxに対してO(sqrt(x) + (約数の個数)^2)で次のステップに進める。
    // さて、xの約数の約数もxの約数であるため、sqrt(x)で計算した分はつかいまわせる。また、「次のステップ」の最悪ケースは一つ小さな約数に行くことであるから、結局O((約数の個数)^3)で抑えることが出来て、全体でO(sqrt(x) + (約数の個数)^3)
    // 補足: xの約数dに対して、gcd(x, v) = dの必要条件はvがdの倍数であること。すなわち上限はfloor(x / d)個。ここからd < gcd(x, v)なる場合の数を除けばよい。

    auto div = divisors(N);

    map<int, double> memo;
    auto E = [&] (auto self, int x) -> double {
        if (memo.find(x) != memo.end()) return memo[x];
        if (x == 1) return 0;
        double res = 1;

        // gcd(x, v)の値で主客転倒
        int index = 0;
        while (div[index] < x) index++;
        map<int, int> mp;

        for (int i = index; 0 <= i; i--) {
            // ここのcontinueがないとダメ(N未満の数に対してはすべての場所が約数になるとは限らないので)
            if (x % div[i] != 0) continue;
            mp[div[i]] = x / div[i];

            // 約数系除原理
            for (int j = 1; i + j < div.size(); j++) {
                if (div[i + j] % div[i] == 0 && mp.find(div[i + j]) != mp.end()) {
                    mp[div[i]] -= mp[div[i + j]];
                }
            }
        }

        for (auto it : mp) {
            if (it.first == x) continue;
            res += (1 + self(self, it.first)) * it.second;
        }

        return memo[x] = res / (x - 1);
    };

    cout << setprecision(10);
    cout << E(E, N) << "\n";
}
0