icpc_library
This documentation is automatically generated by online-judge-tools/verification-helper
test/math/Sum_of_Totient_Function.test.cpp
- View this file on GitHub
- Last update: 2023-11-14 21:03:22+09:00
- Problem: https://judge.yosupo.jp/problem/sum_of_totient_function
Depends on
- Dirichlet Series
(math/dirichlet_series.hpp)
- テンプレート
(template/template.hpp)
- modint
(utility/modint.hpp)
Code
#define PROBLEM "https://judge.yosupo.jp/problem/sum_of_totient_function"
#include "../../math/dirichlet_series.hpp"
#include "../../utility/modint.hpp"
#include "../../template/template.hpp"
namespace lib {
using mint = modint998244353;
void main_() {
ll n;
std::cin >> n;
using DirichletSeries = DirichletSeries<mint, 0>;
DirichletSeries::set_size(n);
mint ans = (DirichletSeries::zeta1() / DirichletSeries::zeta()).get_sum();
std::cout << ans.val() << '\n';
}
} // namespace ebi
int main() {
int t = 1;
// std::cin >> t;
while (t--) {
lib::main_();
}
return 0;
}#line 1 "test/math/Sum_of_Totient_Function.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/sum_of_totient_function"
#line 2 "math/dirichlet_series.hpp"
#line 2 "template/template.hpp"
#include <bits/stdc++.h>
#define rep(i, s, n) for (int i = (int)(s); i < (int)(n); i++)
#define rrep(i, s, n) for (int i = (int)(n)-1; i >= (int)(s); i--)
#define all(v) v.begin(), v.end()
using ll = long long;
using ld = long double;
using ull = unsigned long long;
template <typename T> bool chmin(T &a, const T &b) {
if (a <= b) return false;
a = b;
return true;
}
template <typename T> bool chmax(T &a, const T &b) {
if (a >= b) return false;
a = b;
return true;
}
namespace lib {
using namespace std;
} // namespace lib
// using namespace lib;
#line 4 "math/dirichlet_series.hpp"
namespace lib {
template <class T>
std::vector<T> dirichlet_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
assert(a.size() == b.size());
int n = a.size() - 1;
std::vector<T> c(n + 1, 0);
for (int i = 1; i <= n; i++) {
for (int j = 1; i * j <= n; j++) {
c[i * j] += a[i] * b[j];
}
}
return c;
}
template <class T, int id> struct DirichletSeries {
private:
using Self = DirichletSeries<T, id>;
void set(std::function<T(ll)> f, std::function<T(ll)> F) {
for (int i = 1; i <= K; i++) {
a[i] = f(i);
}
for (int i = 1; i <= L; i++) {
A[i] = F(N / i);
}
}
public:
DirichletSeries() : a(K + 1), A(L + 1) {}
DirichletSeries(std::function<T(ll)> f, std::function<T(ll)> F)
: a(K + 1), A(L + 1) {
set(f, F);
}
Self operator+(const Self &rhs) const noexcept {
return Self(*this) += rhs;
}
Self operator-(const Self &rhs) const noexcept {
return Self(*this) -= rhs;
}
Self operator*(const Self &rhs) const noexcept {
return Self(*this) *= rhs;
}
Self operator/(const Self &rhs) const noexcept {
return Self(*this) /= rhs;
}
Self operator+=(const Self &rhs) noexcept {
for (int i = 1; i <= K; i++) {
a[i] += rhs.a[i];
}
for (int i = 1; i <= L; i++) {
A[i] += rhs.A[i];
}
return *this;
}
Self operator-=(const Self &rhs) noexcept {
for (int i = 1; i <= K; i++) {
a[i] -= rhs.a[i];
}
for (int i = 1; i <= L; i++) {
A[i] -= rhs.A[i];
}
return *this;
}
Self operator*=(const Self &rhs) noexcept {
Self ret;
ret.a = dirichlet_convolution(a, rhs.a);
std::vector<T> sum_a = a, sum_b = rhs.a;
for (int i = 1; i < K; i++) {
sum_a[i + 1] += sum_a[i];
sum_b[i + 1] += sum_b[i];
}
auto get_A = [&](ll x) -> T {
if (x <= K) {
return sum_a[x];
} else {
return A[N / x];
}
};
auto get_B = [&](ll x) -> T {
if (x <= K) {
return sum_b[x];
} else {
return rhs.A[N / x];
}
};
for (ll l = L, m = 1; l >= 1; l--) {
ll n = N / l;
while (m * m <= n) m++;
m--;
for (int i = 1; i <= m; i++) {
ret.A[l] +=
a[i] * get_B(n / i) + (get_A(n / i) - get_A(m)) * rhs.a[i];
}
}
return ret;
}
// c = a / b
Self operator/=(const Self &b) noexcept {
Self c = *this;
T inv_a = b.a[1].inv();
for (int i = 1; i <= K; i++) {
c.a[i] *= inv_a;
for (int j = 2; i * j <= K; j++) {
c.a[i * j] -= c.a[i] * b.a[j];
}
}
std::vector<T> sum_b = b.a, sum_c = c.a;
for (int i = 1; i < K; ++i) {
sum_b[i + 1] += sum_b[i];
sum_c[i + 1] += sum_c[i];
}
auto get_B = [&](ll x) -> T {
if (x <= K) {
return sum_b[x];
} else {
return b.A[N / x];
}
};
auto get_C = [&](ll x) -> T {
if (x <= K) {
return sum_c[x];
} else {
return c.A[N / x];
}
};
for (ll l = L, m = 1; l >= 1; l--) {
ll n = N / l;
while (m * m <= n) m++;
m--;
for (int i = 2; i <= m; i++) {
c.A[l] -= b.a[i] * get_C(n / i);
}
for (int i = 1; i <= m; i++) {
c.A[l] -= c.a[i] * (get_B(n / i) - get_B(m));
}
c.A[l] *= inv_a;
}
return *this = c;
}
Self pow(ll n) const {
Self res;
res.a[1] = 1;
std::fill(res.A.begin(), res.A.end(), 1);
Self x = *this;
while (n > 0) {
if (n & 1) res = res * x;
x = x * x;
n >>= 1;
}
return res;
}
T get_sum() {
return A[1];
}
static Self zeta() {
Self ret;
std::fill(ret.a.begin(), ret.a.end(), 1);
for (int i = 1; i <= L; i++) {
ret.A[i] = N / i;
}
return ret;
}
static Self zeta1() {
Self ret;
std::iota(ret.a.begin(), ret.a.end(), 0);
T inv2 = T(2).inv();
for (int i = 1; i <= L; i++) {
ll n = N / i;
ret.A[i] = T(n) * T(n + 1) * inv2;
}
return ret;
}
static Self zeta2() {
Self ret;
for (int i = 1; i <= K; i++) {
ret.a[i] = i * i;
}
T inv6 = T(6).inv();
for (int i = 1; i <= L; i++) {
ll n = N / i;
ret.A[i] = T(n) * T(n + 1) * T(2 * n + 1) * inv6;
}
}
static void set_size(ll n) {
N = n;
if (N <= 10) {
K = N;
L = 1;
} else if (N <= 5000) {
K = 1;
while (K * K < N) K++;
L = (N + K - 1) / K;
} else {
L = 1;
while (L * L * L / 50 < N) L++;
K = (N + L - 1) / L;
}
}
static void set_size_multiplicative(ll n) {
N = n;
L = 1;
while (L * L * L < N) L++;
K = L * L;
}
private:
static ll N, K, L;
static std::vector<std::pair<int, int>> prime_pow_table;
std::vector<T> a, A;
};
template <class T, int id> ll DirichletSeries<T, id>::N = 1000000;
template <class T, int id> ll DirichletSeries<T, id>::K = 10000;
template <class T, int id> ll DirichletSeries<T, id>::L = 100;
template <class T, int id>
std::vector<std::pair<int, int>> DirichletSeries<T, id>::prime_pow_table = {};
} // namespace lib
#line 2 "utility/modint.hpp"
#line 4 "utility/modint.hpp"
namespace lib {
template <ll m> struct modint {
using mint = modint;
ll a;
modint(ll x = 0) : a((x % m + m) % m) {}
static constexpr ll mod() {
return m;
}
ll val() const {
return a;
}
ll& val() {
return a;
}
mint pow(ll n) const {
mint res = 1;
mint x = a;
while (n) {
if (n & 1) res *= x;
x *= x;
n >>= 1;
}
return res;
}
mint inv() const {
return pow(m - 2);
}
mint& operator+=(const mint rhs) {
a += rhs.a;
if (a >= m) a -= m;
return *this;
}
mint& operator-=(const mint rhs) {
if (a < rhs.a) a += m;
a -= rhs.a;
return *this;
}
mint& operator*=(const mint rhs) {
a = a * rhs.a % m;
return *this;
}
mint& operator/=(mint rhs) {
*this *= rhs.inv();
return *this;
}
friend mint operator+(const mint& lhs, const mint& rhs) {
return mint(lhs) += rhs;
}
friend mint operator-(const mint& lhs, const mint& rhs) {
return mint(lhs) -= rhs;
}
friend mint operator*(const mint& lhs, const mint& rhs) {
return mint(lhs) *= rhs;
}
friend mint operator/(const mint& lhs, const mint& rhs) {
return mint(lhs) /= rhs;
}
friend bool operator==(const modint &lhs, const modint &rhs) {
return lhs.a == rhs.a;
}
friend bool operator!=(const modint &lhs, const modint &rhs) {
return !(lhs == rhs);
}
mint operator+() const {
return *this;
}
mint operator-() const {
return mint() - *this;
}
};
using modint998244353 = modint<998244353>;
using modint1000000007 = modint<1'000'000'007>;
} // namespace lib
#line 6 "test/math/Sum_of_Totient_Function.test.cpp"
namespace lib {
using mint = modint998244353;
void main_() {
ll n;
std::cin >> n;
using DirichletSeries = DirichletSeries<mint, 0>;
DirichletSeries::set_size(n);
mint ans = (DirichletSeries::zeta1() / DirichletSeries::zeta()).get_sum();
std::cout << ans.val() << '\n';
}
} // namespace ebi
int main() {
int t = 1;
// std::cin >> t;
while (t--) {
lib::main_();
}
return 0;
}