#1174778 (C++23) No.3607 Sum of Powers of GCDs

• 問題: No.3607 Sum of Powers of GCDs
• ユーザー: pockyny
• 提出日時: 2026-08-02 02:11:11
• 言語: C++23 (gcc 15.2.0 + boost 1.90.0)
• 結果: AC
• 実行時間: 1,657 ms / 2,500 ms
• コード長: 6,220 bytes
• コンパイル時間: 2,013 ms
• コンパイル使用メモリ: 197,448 KB
• 実行使用メモリ: 62,412 KB
• 最終ジャッジ日時: 2026-08-02 02:11:24
• 合計ジャッジ時間: 11,235 ms
• ジャッジサーバーID （参考情報）: judge1_1 / judge2_0

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

ソースコード

// 3/4乗方針を定数倍高速化で通すために、AI様にガン頼り

#pragma GCC target("avx2")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")

#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>

#include <atcoder/modint>

using namespace std;
using namespace atcoder;

using mint = modint998244353;

constexpr int MX = 1'000'000;
constexpr int KMAX = 10;

/*
 * mu[x] = Möbius 関数
 *
 * mu_prefix[i] = sum_{x=1}^{i-1} mu[x]
 */
int mu[MX + 2];
int mu_prefix[MX + 2];
int least_prime[MX + 2];

/*
 * power_sum[k][i] = sum_{x=0}^{i-1} x^k
 */
mint power_sum[KMAX + 1][MX + 2];

/*
 * quotient_power[q] = q^N
 *
 * 各テストケースについて、
 * q = floor(M / x) として現れるものだけ設定する。
 */
mint quotient_power[MX + 2];

inline mint mod_pow(mint a, int exponent) {
    mint result = 1;

    while (exponent > 0) {
        if (exponent & 1) {
            result *= a;
        }

        a *= a;
        exponent >>= 1;
    }

    return result;
}

inline int integer_sqrt(int x) {
    int result = static_cast<int>(sqrt(static_cast<double>(x)));

    while (1LL * (result + 1) * (result + 1) <= x) {
        ++result;
    }

    while (1LL * result * result > x) {
        --result;
    }

    return result;
}

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

    /*
     * 線形篩で Möbius 関数を前計算する。
     */
    {
        vector<int> primes;
        primes.reserve(80'000);

        mu[1] = 1;

        for (int i = 2; i <= MX; ++i) {
            if (least_prime[i] == 0) {
                least_prime[i] = i;
                primes.push_back(i);
                mu[i] = -1;
            }

            for (int p : primes) {
                int64_t next = static_cast<int64_t>(i) * p;

                if (next > MX) {
                    break;
                }

                least_prime[next] = p;

                if (i % p == 0) {
                    mu[next] = 0;
                    break;
                }

                mu[next] = -mu[i];
            }
        }

        for (int i = 1; i <= MX + 1; ++i) {
            mu_prefix[i] = mu_prefix[i - 1] + mu[i - 1];
        }
    }

    /*
     * power_sum[k][i] = 0^k + 1^k + ... + (i-1)^k
     *
     * i^1, i^2, ..., i^10 を順に掛けて作る。
     */
    {
        for (int x = 1; x <= MX; ++x) {
            mint current_power = 1;
            mint value = x;

            for (int k = 1; k <= KMAX; ++k) {
                current_power *= value;
                power_sum[k][x + 1] =
                    power_sum[k][x] + current_power;
            }
        }
    }

    int T;
    cin >> T;

    while (T--) {
        int N, M, K;
        cin >> N >> M >> K;

        /*
         * q = floor(M / x) として現れる q について、
         * q^N を前計算する。
         *
         * 商の異なる値は O(sqrt(M)) 個。
         */
        {
            int q = 1;

            while (true) {
                int left_boundary = M / (q + 1);

                quotient_power[q] = mod_pow(q, N);

                if (left_boundary == 0) {
                    break;
                }

                q = M / left_boundary;
            }
        }

        /*
         * G(x) =
         * sum_{d=1}^x mu[d] * floor(x/d)^N
         *
         * s = floor(sqrt(x)) として、
         *
         * 1. d <= s は d を直接列挙
         * 2. d > s は q = floor(x/d) を直接列挙
         *
         * と分割する。
         */
        auto calc_G = [&](int x) -> mint {
            mint result = 0;

            const int s = integer_sqrt(x);

            /*
             * d <= sqrt(x) の部分。
             *
             * 商 floor(x/d) は、必ず floor(M/z) として
             * 現れる値なので quotient_power に存在する。
             */
            for (int d = 1; d <= s; ++d) {
                const int coefficient = mu[d];

                if (coefficient != 0) {
                    result +=
                        quotient_power[x / d] * coefficient;
                }
            }

            /*
             * d > sqrt(x) の部分。
             *
             * q = floor(x/d) は
             * 1 <= q <= floor(x/(s+1))。
             *
             * floor(x/d) = q となる d の範囲は
             *
             * floor(x/(q+1)) < d <= floor(x/q)。
             *
             * 右端 floor(x/q) は、一つ前のループの左端と
             * 等しいので持ち回る。
             */
            const int quotient_max = x / (s + 1);
            int right = x;

            for (int q = 1; q <= quotient_max; ++q) {
                const int left = x / (q + 1);

                /*
                 * d > s のみを扱う。
                 *
                 * 区間は (max(left,s), right]。
                 */
                const int effective_left =
                    left > s ? left : s;

                const int coefficient =
                    mu_prefix[right + 1] -
                    mu_prefix[effective_left + 1];

                if (coefficient != 0) {
                    result += quotient_power[q] * coefficient;
                }

                right = left;
            }

            return result;
        };

        mint answer = 0;

        /*
         * floor(M/x) が等しい x の区間ごとにまとめる。
         *
         * 現在の区間は
         *
         * left_boundary < x <= right_boundary
         *
         * で、商は q。
         *
         * right_boundary は前の left_boundary を
         * 持ち回せる。
         */
        int q = 1;
        int right_boundary = M;

        while (true) {
            const int left_boundary = M / (q + 1);

            const mint x_power_sum =
                power_sum[K][right_boundary + 1] -
                power_sum[K][left_boundary + 1];

            answer += x_power_sum * calc_G(q);

            if (left_boundary == 0) {
                break;
            }

            right_boundary = left_boundary;
            q = M / left_boundary;
        }

        cout << answer.val() << '\n';
    }
}