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

• 問題: No.3607 Sum of Powers of GCDs
• ユーザー: pockyny
• 提出日時: 2026-08-02 13:37:41
• 言語: C++23 (gcc 15.2.0 + boost 1.90.0)
• 結果: AC
• 実行時間: 1,361 ms / 2,500 ms
• コード長: 4,514 bytes
• コンパイル時間: 2,108 ms
• コンパイル使用メモリ: 200,608 KB
• 実行使用メモリ: 59,392 KB
• 最終ジャッジ日時: 2026-08-02 13:37:59
• 合計ジャッジ時間: 9,587 ms
• ジャッジサーバーID （参考情報）: judge1_0 / judge2_0

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

ソースコード

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

#include <iostream>
#include <vector>
#include <cmath>
#include <atcoder/modint>

using namespace std;
using namespace atcoder;
using mint = modint998244353;
const int MX = 1000000;
bool isPrime[MX + 10];
mint moe[MX + 10],moe_sum[MX + 10];
// k乗和の累積和 sum[i][j] = 0 ~ i - 1 までのj乗和
mint sum[MX + 10][11];
mint pw(mint a,int x){
    mint ret = 1;
    while(x){
        if(x&1) ret *= a;
        a *= a; x /= 2;
    }
    return ret;
}
mint pre_calc_pow[MX + 10];

// floor(n/i) = qとなる l<=i<=rを列挙してあれこれするライブラリ
// 除算回数が4√nくらい
template<class F>
inline void quotient_ranges(int n, F&& f) {
    for (int l = 1; l <= n; ) {
        int q = n / l;
        int r = n / q;
        // l<=x<=r で、floor(n/l) = qとなる
        // lambda式とかで処理を渡す
        f(l, r, q);
        l = r + 1;
    }
}

// floor(n/i)の種類数は少ない
// i <= sqrt(n)に対して、iで全列挙をする
// floor(n/i) <= sqrt(n) かつ i > sqrt(n) に対して、区間の列挙ができる
// 除算回数が2√nくらい
template<class SmallIndexF, class SmallQuotientF>
inline void floor_quotient_split(
    int n,
    SmallIndexF&& small_index,
    SmallQuotientF&& small_quotient
) {
    const int s = sqrt(n);

    // i <= sqrt(n) を直接列挙
    for (int i = 1; i <= s; ++i) {
        // (floor(n/i),i) で計算
        small_index(n/i, i);
    }

    // i > sqrt(n) を商 q ごとに列挙
    const int qmax = n / (s + 1);
    int right = n;

    for (int q = 1; q <= qmax; ++q) {
        const int left = n / (q + 1);
        // max(left,s) < i <= right
        small_quotient(q, std::max(left, s), right);
        right = left;
    }
}

int main(){
    std::ios_base::sync_with_stdio(false);
    std::cin.tie(nullptr);
    {
        int i,j;
        for(i=2;i<=MX;i++){
            if(isPrime[i]) continue;
            for(j=2*i;j<=MX;j+=i) isPrime[j] = true;
        }
        for(i=2;i<=MX;i++) isPrime[i] = !isPrime[i];
        for(i=1;i<=MX;i++) moe[i] = 1;
        for(i=2;i<=MX;i++){
            if(!isPrime[i]) continue;
            for(j=1;j*i<=MX;j++){
                if(j%i==0) moe[j*i] = 0;
                moe[j*i] *= -1;
            } 
        }
        for(i=1;i<=MX + 1;i++){
            mint val = i - 1;
            for(j=1;j<=10;j++){
                sum[i][j] = sum[i - 1][j] + val;
                val *= (i - 1);
            }
        }
        for(i=1;i<=MX + 1;i++) moe_sum[i] = moe_sum[i - 1] + moe[i - 1];
    }
    int t; cin >> t;
    while(t){
        t--;
        int n,m,k; cin >> n >> m >> k;
        // ans = \sum_i\sum_{floor(m/x) = i}x^k(\sum_{A<=i}moe[A]*floor(i/A)^n)
        // g(i) := (\sum_{A<=i}moe[A]*floor(i/A)^n) として、floor(i/A)が同じ値ごとにまとめて計算する

        // 毎回pw(floor(i/A),n)を計算するとTLEするので、候補のO(√M)通りを先に計算 (floor(M/x)のn乗数を先に計算)
        // ただ、前計算するだけだと、まだTLEしたので、floor(m/x)を先に計算して、割り算の回数を減らす
        // {floor(m,x),(l,r]}で格納
        vector<pair<int,pair<int,int>>> v;
        floor_quotient_split(
            m,
            [&](int q,int i){
                v.push_back({q,{i - 1,i}});
                pre_calc_pow[q] = pw(q,n);
            },
            [&](int q,int l,int r){
                v.push_back({q,{l,r}});
                pre_calc_pow[q] = pw(q,n);    
            }
        );
        // quotient_ranges(m,[&](int l,int r,int q){
        //     pre_calc_pow[q] = pw(q,n);
        //     v.push_back({q,{l,r}});
        // });
        auto calcG = [&](int i){
            mint ret = 0;
            floor_quotient_split(
                i,
                [&](int q,int i){
                    ret += (moe_sum[i + 1] - moe_sum[i])*pre_calc_pow[q];
                },
                [&](int q,int l,int r){
                    ret += (moe_sum[r + 1] - moe_sum[l + 1])*pre_calc_pow[q];
                }
            );
            // quotient_ranges(i,[&](int l,int r,int q){
            //     ret += (moe_sum[r + 1] - moe_sum[l])*pre_calc_pow[q];
            // });
            return ret;
        };
        mint ans = 0;
        for(auto [i,lr]:v){
            auto [l,r] = lr;
            ans += (sum[r + 1][k] - sum[l + 1][k])*calcG(i);
        }
        cout << ans.val() << "\n";
    }
}