結果

問題 No.1781 LCM
コンテスト
ユーザー Nachia
提出日時 2024-08-28 22:47:00
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
WA  
実行時間 -
コード長 6,886 bytes
コンパイル時間 1,633 ms
コンパイル使用メモリ 101,664 KB
最終ジャッジ日時 2025-02-24 02:39:03
ジャッジサーバーID
（参考情報） 		judge5 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 2
other AC * 27 WA * 4
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <utility>

namespace nachia{

std::pair<std::vector<long long>, std::vector<long long>>
    CountingPrimesTable(long long maxval)
{
    using i64 = long long;
    i64 n = maxval;
    if(n <= 1) return { {0,0}, {0,0} };

    i64 sqrtn = 0; while(sqrtn * sqrtn <= n) sqrtn++; sqrtn--;
    
    std::vector<i64> prefix_sum_fairy(sqrtn+1, 0);
    for(i64 i=1; i<=sqrtn; i++) prefix_sum_fairy[i] = i-1;

    std::vector<i64> prefix_sum_devil(sqrtn+1, 0);
    for(i64 i=1; i<=sqrtn; i++) prefix_sum_devil[i] = n/i-1;

    for(i64 p=2; p<=sqrtn; p++){
        i64 prime_count_p = prefix_sum_fairy[p];
        i64 prime_count_p_minus1 = prefix_sum_fairy[p-1];
        if(prime_count_p == prime_count_p_minus1) continue;
        for(i64 devil_id = 1; devil_id <= sqrtn; devil_id++){
            if(devil_id * p <= sqrtn){
                prefix_sum_devil[devil_id] -= prefix_sum_devil[devil_id * p] - prime_count_p_minus1;
            }
            else{
                i64 tg_fairy = n / (devil_id * p);
                if(tg_fairy < p) break;
                prefix_sum_devil[devil_id] -= prefix_sum_fairy[tg_fairy] - prime_count_p_minus1;
            }
        }
        for(i64 fairy_id = sqrtn/p; fairy_id >= p; fairy_id--){
            i64 dc = prefix_sum_fairy[fairy_id] - prime_count_p_minus1;
            i64 max_tg = std::min(fairy_id * p + p - 1, sqrtn);
            for(i64 tg_fairy = fairy_id * p; tg_fairy <= max_tg; tg_fairy++) prefix_sum_fairy[tg_fairy] -= dc;
        }
    }

    return std::make_pair(std::move(prefix_sum_fairy), std::move(prefix_sum_devil));
}

long long CountingPrimes(long long maxval){
    return CountingPrimesTable(maxval).second[1];
}

} // namespace nachia

namespace nachia {

template<class Elem, class F1, class F2>
std::pair<std::vector<Elem>, std::vector<Elem>>
MultiplicativePrefixSum(long long N, F1 f1, F2 f2){
    using i64 = long long;
    i64 Nr2 = 0; while((Nr2+1) * (Nr2+1) <= N) Nr2++;
    i64 Nr3 = 0; while((Nr3+1) * (Nr3+1) * (Nr3+1) <= N) Nr3++;
    i64 Nr6 = 0; while((Nr6+1) * (Nr6+1) <= Nr3) Nr6++;
    i64 N2r3 = std::max(Nr2, N / (N / Nr6 / Nr2 + 1));
    i64 Dsz = N / (Nr2+1);
    i64 Asz = (Nr2+1) + Dsz;
    std::vector<Elem> A(Asz + 1, Elem(0));
    auto Aat = [Asz,Nr2,N](i64 x) -> i64 { return x <= Nr2 ? x : (Asz-N/x); };

    std::vector<int> sieve(Nr2+1, 1);
    sieve[0] = sieve[1] = 0;
    std::vector<i64> primes;
    for(i64 p=2; p<=Nr2; p++) if(sieve[p]){
        primes.push_back(p);
        for(i64 q=p*p; q<=Nr2; q+=p) sieve[q] = 0;
    }
    i64 Nr2Pi = i64(primes.size());
    i64 Nr3Pi = 0;
    while(Nr3Pi < Nr2Pi && primes[Nr3Pi] <= Nr3) Nr3Pi++;

    // 1
    A[1] += Elem(1);

    // p
    for(i64 p : primes) if(Nr3 < p) A[p] = f1(p, 1, p);
    for(i64 i=1; i<=Dsz; i++) A[Asz-i] = f2(N/i);
    for(i64 i=1; i<Dsz; i++) A[Asz-i] -= A[Asz-(i+1)];
    A[Asz-Dsz] -= f2(Nr2);

    // pq, p^2
    for(i64 px=Nr3Pi; px<Nr2Pi; px++){
        i64 p = primes[px];
        A[Aat(p*p)] += f1(p, 2, p*p); // p^2
        Elem fp = f1(p, 1, p);
        Elem ppre = f2(p);
        // pq (p<q)
        for(i64 i=N/(p*(p+1)); i>=1; i--){
            Elem pnx = f2(N/(i*p));
            A[Asz-i] += (pnx - ppre) * fp;
            ppre = pnx;
        }
    }

    for(i64 i=1; i+1<Asz; i++) A[i+1] += A[i];

    auto getfpeBuf = [&](i64 p) -> std::pair<std::vector<i64>, std::vector<Elem>> {
        std::vector<Elem> fpebuf;
        std::vector<i64> pe;
        for(i64 pp=p, e=1; pp<=N; pp*=p, e++){
            fpebuf.push_back(f1(p, e, pp));
            pe.push_back(pp);
        }
        pe.push_back(N+1);
        return std::make_pair(std::move(pe), std::move(fpebuf));
    };
    std::vector<std::pair<std::vector<i64>, std::vector<Elem>>> fpebufall(Nr2Pi);
    for(i64 px=0; px<Nr2Pi; px++) fpebufall[px] = getfpeBuf(primes[px]);

    i64 px = Nr3Pi - 1;

    // setup fenwick tree
    i64 fwtSz = Aat(N2r3) + 1;
    auto fwAdd = [&](i64 i, Elem v) -> void {
        while(i < fwtSz){ A[i] += v; i += i & -i; }
    };
    auto fwSum = [&](i64 r) -> Elem {
        Elem x = Elem(0);
        while(r > 0){ x += A[r]; r -= r & -r; }
        return x;
    };
    for(i64 i=fwtSz-2; i>=1; i--) A[i+1] -= A[i];
    for(i64 i=1; i<fwtSz; i++) if(i+(i&-i) < fwtSz) A[i+(i&-i)] += A[i];

    auto getAsum = [&](i64 v) -> Elem { return v < fwtSz ? fwSum(v) : A[v]; };

    for( ; px>=0 && primes[px]>Nr6; px--){
        auto& [xpe, xfpe] = fpebufall[px];
        // N2r3 <= n
        for(i64 x=1; Asz-x>=fwtSz; x++){
            i64 n = N / x;
            for(i64 e=0; xpe[e]<=n; e++){
                A[Asz-x] += xfpe[e] * getAsum(Aat(n/xpe[e]));
            }
        }
        // search v <= N2r3 && primes[px] == lpf(v)
        std::vector<std::pair<i64, Elem>> bfs;
        for(i64 e=0; xpe[e] <= N2r3; e++){
            fwAdd(Aat(xpe[e]), xfpe[e]);
            bfs.push_back({ xpe[e], xfpe[e] });
        }
        for(i64 pxv=px+1; pxv<Nr2Pi; pxv++){
            auto& [vpe, vfpe] = fpebufall[pxv];
            i64 p = vpe[0];
            i64 bfsz = bfs.size();
            for(i64 i=0; i<bfsz; i++){
                auto [v, fv] = bfs[i];
                for(i64 e=0; e<i64(vpe.size()); e++) if(v * vpe[e] <= N2r3){
                    i64 nx = v * vpe[e];
                    Elem nxval = fv * vfpe[e];
                    fwAdd(Aat(nx), nxval);
                    if(nx*p <= N2r3) bfs.push_back({ nx, nxval });
                }
            }
        }
    }

    for(i64 i=fwtSz-1; i>=1; i--) if(i+(i&-i) < fwtSz) A[i+(i&-i)] -= A[i];
    for(i64 i=1; i<fwtSz-1; i++) A[i+1] += A[i];

    for( ; px>=0; px--){
        auto [pe, fpebuf] = std::move(fpebufall[px]);
        for(i64 i=1; i<=Dsz; i++){
            for(i64 t=0; i*pe[t]<=N; t++){
                A[Asz-i] += A[Aat(N/(i*pe[t]))] * fpebuf[t];
            }
        }
        for(i64 i=Nr2; i>=pe[0]; i--){
            for(i64 t=0; pe[t]<=i; t++){
                A[i] += A[i/pe[t]] * fpebuf[t];
            }
        }
    }

    return std::make_pair(
        std::vector<Elem>(A.begin(), A.begin() + Nr2),
        std::vector<Elem>(A.rbegin(), A.rbegin() + (Dsz + 1))
    );
}

} // namespace nachia
#include <atcoder/modint>
using Modint = atcoder::static_modint<998244353>;
using namespace std;

int main(){
    ios::sync_with_stdio(false); cin.tie(nullptr);
    using i64 = long long;
    i64 N; cin >> N;
    i64 M; cin >> M;
    
    vector<Modint> X(100);
    for(i64 i=1; i<100; i++) X[i] = Modint(i+1).pow(N);

    auto [pcf, pcd] = nachia::CountingPrimesTable(M);
    i64 s = i64(pcf.size() - 1);

    auto ans = nachia::MultiplicativePrefixSum<Modint>(
        M,
        [&](i64 p, i64 e, i64 pe){ return X[e]; },
        [&](i64 x){ return Modint((x <= s ? pcf[x] : pcd[M/x]) * X[1]); });

    cout << ans.second[1].val() << '\n';
    return 0;
}
0