#1152701 (C++23) No.3505 Sum of Prod of Root

• 問題: No.3505 Sum of Prod of Root
• ユーザー: Naru820
• 提出日時: 2026-03-17 15:53:08
• 言語: C++23 (gcc 15.2.0 + boost 1.89.0)
• 結果: AC
• 実行時間: 128 ms / 2,000 ms
• コード長: 3,532 bytes
• コンパイル時間: 2,628 ms
• コンパイル使用メモリ: 341,568 KB
• 実行使用メモリ: 23,552 KB
• 最終ジャッジ日時: 2026-04-17 19:40:37
• 合計ジャッジ時間: 4,458 ms
• ジャッジサーバーID （参考情報）: judge1_1 / judge2_0

ファイルパターン	結果
sample	AC * 1
other	AC * 13

ソースコード

//GPT 5.4-pro
#include <bits/stdc++.h>
using namespace std;

using u64 = unsigned long long;
using u128 = __uint128_t;

static constexpr int MOD = 998244353;
static constexpr int INV2  = 499122177; // 2^{-1}
static constexpr int INV6  = 166374059; // 6^{-1}
static constexpr int INV30 = 432572553; // 30^{-1}

int addmod(int a, int b) {
    int s = a + b;
    if (s >= MOD) s -= MOD;
    return s;
}
int submod(int a, int b) {
    int s = a - b;
    if (s < 0) s += MOD;
    return s;
}
int mulmod(long long a, long long b) {
    return int(a * b % MOD);
}

// a^k が limit を超えたら limit+1 を返す
u64 ipow_limit(u64 a, int k, u64 limit) {
    u128 v = 1;
    for (int i = 0; i < k; ++i) {
        v *= a;
        if (v > limit) return limit + 1;
    }
    return (u64)v;
}

u64 isqrt_u64(u64 x) {
    long double d = sqrt((long double)x);
    u64 r = (u64)d;
    while ((u128)(r + 1) * (r + 1) <= x) ++r;
    while ((u128)r * r > x) --r;
    return r;
}

// G(n) = sum_{i=1}^n i * floor(sqrt(i)) mod MOD
int G(u64 n) {
    if (n == 0) return 0;

    u64 m = isqrt_u64(n);
    u64 u = m - 1; // 完全ブロックは t=1..u

    int uu = int(u % MOD);

    // sum t^2
    int s2 = mulmod(mulmod(uu, uu + 1LL), mulmod(2LL * uu + 1, INV6));

    // sum t^3
    int half = mulmod(mulmod(uu, uu + 1LL), INV2);
    int s3 = mulmod(half, half);

    // sum t^4
    int poly = ((3LL * uu % MOD) * uu + 3LL * uu - 1) % MOD;
    if (poly < 0) poly += MOD;
    int s4 = mulmod(
        mulmod(mulmod(uu, uu + 1LL), 2LL * uu + 1),
        mulmod(poly, INV30)
    );

    int full = (2LL * s4 + 3LL * s3 + s2) % MOD;

    // 末尾の不完全ブロック: m * sum_{i=m^2}^n i
    u64 a = m * m;
    int aa = int(a % MOD);
    int nn = int(n % MOD);
    int cnt = int((n - a + 1) % MOD);
    int tail_sum = mulmod(mulmod((aa + 1LL * nn) % MOD, cnt), INV2);

    return addmod(full, mulmod(int(m % MOD), tail_sum));
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    u64 N;
    cin >> N;

    // base は最大でも floor((1e18)^(1/3)) = 1e6
    const int MAXB = 1000000 + 5;

    // inv[x] = x^{-1} mod MOD
    vector<int> inv(MAXB + 1);
    inv[1] = 1;
    for (int i = 2; i <= MAXB; ++i) {
        inv[i] = MOD - 1LL * (MOD / i) * inv[MOD % i] % MOD;
    }

    // events: (x=a^k, multiply by a/(a-1))
    vector<pair<u64, int>> events;
    events.reserve(1100000);

    for (int k = 3; k <= 60; ++k) {
        if (ipow_limit(2, k, N) > N) break; // これ以上の k は全部 1
        for (u64 a = 2;; ++a) {
            u64 x = ipow_limit(a, k, N);
            if (x > N) break;
            events.push_back({x, mulmod(a % MOD, inv[a - 1])});
        }
    }

    sort(events.begin(), events.end(),
         [](const auto& L, const auto& R) { return L.first < R.first; });

    int ans = 0;
    int curH = 1;   // 現在の H(i)
    u64 prev = 1;

    size_t i = 0;
    while (i < events.size()) {
        u64 x = events[i].first;

        // [prev, x-1] は H が一定
        if (prev <= x - 1) {
            int seg = submod(G(x - 1), G(prev - 1));
            ans = addmod(ans, mulmod(curH, seg));
        }

        // x で起こるすべての更新をまとめて反映
        while (i < events.size() && events[i].first == x) {
            curH = mulmod(curH, events[i].second);
            ++i;
        }
        prev = x;
    }

    if (prev <= N) {
        int seg = submod(G(N), G(prev - 1));
        ans = addmod(ans, mulmod(curH, seg));
    }

    cout << ans << '\n';
    return 0;
}