Verify/Sum_of_Totient_Function.test.cpp

View this file on GitHub
Last update: 2023-02-11 22:28:05+08:00
Problem: https://judge.yosupo.jp/problem/sum_of_totient_function

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/sum_of_totient_function"

#include "../ModInt/Modint32.hpp"
#include "../Number_Theory/Mf_Sieve.hpp"
#include "../Template/Template.hpp"

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

    constexpr int mod = 998244353;
    using Z = mint<mod>;

    i64 n;
    std::cin >> n;

    int TH = (int)pow(n, 2.0 / 3.0);

    auto f = [&](i64 n, i64 p, i64 c) -> Z { return Z(n - n / p); };
    auto phi = mf_sieve<Z>(TH, f);
    std::vector<Z> dphi(phi.size());
    for (int i = 1; i <= TH; i++) dphi[i] = dphi[i - 1] + phi[i];

    const Z inv2 = Z(2).inverse();
    std::unordered_map<i64, Z> mp;
    std::function<Z(i64)> dfs = [&](i64 x) -> Z {
        if (mp.count(x)) return mp[x];
        if (x <= TH) return dphi[x];
        Z ans = Z(x) * Z(x + 1) * inv2;
        for (i64 l = 2, r; l <= x; l = r + 1) {
            r = x / (x / l);
            ans -= dfs(x / l) * Z(r - l + 1);
        }
        return mp[x] = ans;
    };

    std::cout << dfs(n) << '\n';

    return 0;
}

#line 1 "Verify/Sum_of_Totient_Function.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/sum_of_totient_function"

#line 2 "ModInt/Modint32.hpp"

#line 2 "Template/Template.hpp"

#include <bits/stdc++.h>

using i64 = std::int64_t;
#line 4 "ModInt/Modint32.hpp"

template <int mod>
struct mint {
    int x;
    mint() : x(0) {}
    mint(int64_t y) : x(y >= 0 ? y % mod : (mod - (-y) % mod) % mod) {}
    mint &operator+=(const mint &p) {
        if ((x += p.x) >= mod) x -= mod;
        return *this;
    }
    mint &operator-=(const mint &p) {
        if ((x += mod - p.x) >= mod) x -= mod;
        return *this;
    }
    mint &operator*=(const mint &p) {
        x = (int)(1LL * x * p.x % mod);
        return *this;
    }
    mint &operator/=(const mint &p) {
        *this *= p.inverse();
        return *this;
    }
    mint operator-() const { return mint(-x); }
    mint operator+(const mint &p) const { return mint(*this) += p; }
    mint operator-(const mint &p) const { return mint(*this) -= p; }
    mint operator*(const mint &p) const { return mint(*this) *= p; }
    mint operator/(const mint &p) const { return mint(*this) /= p; }
    bool operator==(const mint &p) const { return x == p.x; }
    bool operator!=(const mint &p) const { return x != p.x; }
    mint inverse() const {
        int a = x, b = mod, u = 1, v = 0, t;
        while (b > 0) {
            t = a / b;
            std::swap(a -= t * b, b);
            std::swap(u -= t * v, v);
        }
        return mint(u);
    }
    friend std::ostream &operator<<(std::ostream &os, const mint &p) {
        return os << p.x;
    }
    friend std::istream &operator>>(std::istream &is, mint &a) {
        int64_t t;
        is >> t;
        a = mint<mod>(t);
        return (is);
    }
    int get() const { return x; }
    static constexpr int get_mod() { return mod; }
};
#line 2 "Number_Theory/Mf_Sieve.hpp"

#line 2 "Number_Theory/Prime_Sieve.hpp"

#line 4 "Number_Theory/Prime_Sieve.hpp"

std::vector<int> prime_sieve(int N) {
    std::vector<bool> sieve(N / 3 + 1, 1);
    for (int p = 5, d = 4, i = 1, sqn = sqrt(N); p <= sqn;
         p += d = 6 - d, i++) {
        if (!sieve[i]) continue;
        for (int q = p * p / 3, r = d * p / 3 + (d * p % 3 == 2), s = 2 * p,
                 qe = sieve.size();
             q < qe; q += r = s - r)
            sieve[q] = 0;
    }
    std::vector<int> ret{2, 3};
    for (int p = 5, d = 4, i = 1; p <= N; p += d = 6 - d, i++)
        if (sieve[i]) ret.push_back(p);
    while (!ret.empty() && ret.back() > N) ret.pop_back();
    return ret;
}
#line 4 "Number_Theory/Mf_Sieve.hpp"

template <class T>
std::vector<T> mf_sieve(int n, std::function<T(i64, i64, i64)> f) {
    /*
        ##pragma f is a multiplicative-function
        f(n,p,c) <-> n=p^c
    */
    std::vector<T> ans(n + 1, T());
    std::vector<int> ps = prime_sieve(n);
    int p(ps.size());

    std::function<void(int, i64, T)> dfs = [&](int i, i64 x, T y) {
        ans[x] = y;
        if (y == T()) return;
        for (int j = i + 1; j < p; j++) {
            i64 nx = x * ps[j];
            i64 dx = ps[j];
            if (nx > n) break;
            for (int c = 1; nx <= n; nx *= ps[j], dx *= ps[j], ++c)
                dfs(j, nx, y * f(dx, ps[j], c));
        }
    };
    ans[1] = 1;
    dfs(-1, 1, 1);
    return ans;
};
#line 6 "Verify/Sum_of_Totient_Function.test.cpp"

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

    constexpr int mod = 998244353;
    using Z = mint<mod>;

    i64 n;
    std::cin >> n;

    int TH = (int)pow(n, 2.0 / 3.0);

    auto f = [&](i64 n, i64 p, i64 c) -> Z { return Z(n - n / p); };
    auto phi = mf_sieve<Z>(TH, f);
    std::vector<Z> dphi(phi.size());
    for (int i = 1; i <= TH; i++) dphi[i] = dphi[i - 1] + phi[i];

    const Z inv2 = Z(2).inverse();
    std::unordered_map<i64, Z> mp;
    std::function<Z(i64)> dfs = [&](i64 x) -> Z {
        if (mp.count(x)) return mp[x];
        if (x <= TH) return dphi[x];
        Z ans = Z(x) * Z(x + 1) * inv2;
        for (i64 l = 2, r; l <= x; l = r + 1) {
            r = x / (x / l);
            ans -= dfs(x / l) * Z(r - l + 1);
        }
        return mp[x] = ans;
    };

    std::cout << dfs(n) << '\n';

    return 0;
}