萌新刚学 OI，求调数论模板题

#include <bits/stdc++.h>
#define mod 998244353
using namespace std;
namespace modulo_int {
	#define mint_int long long
	void exgcd(mint_int& x, mint_int& y, mint_int n, mint_int m) {
		if (m) {
			exgcd(y, x, m, n % m);
			y -= n / m * x;
		} else {
			x = 1;
			y = 0;
		}
	}
	mint_int mint_exgcd_a, mint_exgcd_b;
	#define mint_mod mod
	struct mint {
		mint_int n;
		mint() : n(0) {}
		mint(const mint_int &_n) { n = (_n % mint_mod + mint_mod) % mint_mod; }
		template<typename T> const mint operator+(const T &x) const { return n + mint(x).n; }
		template<typename T> const mint operator-(const T &x) const { return n - mint(x).n; }
		template<typename T> const mint operator*(const T &x) const { return n * mint(x).n; }
		template<typename T> const mint operator%(const T &x) const { return n % mint(x).n; }
		template<typename T> const mint operator|(const T &x) const { return n | mint(x).n; }
		template<typename T> const mint operator&(const T &x) const { return n & mint(x).n; }
		template<typename T> const mint operator^(const T &x) const { return n ^ mint(x).n; }
		template<typename T> const mint operator/(const T &x) const {
			exgcd(mint_exgcd_a, mint_exgcd_b, mint(x).n, mint_mod);
			return n * mint_exgcd_a;
		}
		template<typename T> mint &operator+=(const T x) { return *this = *this + (mint)x; }
		template<typename T> mint &operator-=(const T x) { return *this = *this - (mint)x; }
		template<typename T> mint &operator*=(const T x) { return *this = *this * (mint)x; }
		template<typename T> mint &operator/=(const T x) { return *this = *this / (mint)x; }
		template<typename T> mint &operator%=(const T x) { return *this = *this % (mint)x; }
		template<typename T> mint &operator|=(const T x) { return *this = *this | (mint)x; }
		template<typename T> mint &operator&=(const T x) { return *this = *this & (mint)x; }
		template<typename T> mint &operator^=(const T x) { return *this = *this ^ (mint)x; }
		mint &operator++() { return *this = n + 1; }
		mint &operator--() { return *this = n - 1; }
		mint operator++(signed) {
			n = ((n + 1 == mint_mod) ? 0 : (n + 1));
			return n - 1;
		}
		mint operator--(signed) {
			n = (n ? (n - 1) : (mint_mod - 1));
			return n + 1;
		}
		template<typename T> bool  operator<(const T &b) { return n <  mint(b).n; }
		template<typename T> bool operator<=(const T &b) { return n <= mint(b).n; }
		template<typename T> bool  operator>(const T &b) { return n >  mint(b).n; }
		template<typename T> bool operator>=(const T &b) { return n >= mint(b).n; }
		template<typename T> bool operator==(const T &b) { return n == mint(b).n; }
		template<typename T> bool operator!=(const T &b) { return n != mint(b).n; }
		// operator mint_int() const { return n; }
		mint operator-() const { return -n; }
		const mint &operator+() const { return *this; }
	};
	istream &operator>>(istream &A, mint &b) {
		mint_int c;
		A >> c;
		b = c;
		return A;
	}
	ostream &operator<<(ostream &A, const mint &b) {
		return A << b.n;
	}
	#undef mint_mod
	#undef mint_int
}
using modulo_int::mint;
ostream &cans(cout);
istream &operator>>(istream &a, __int128 &b) {
	string s;
	a >> s;
	b = 0;
	for (char i : s) {
		b = (b << 1) + (b << 3) + (i ^ '0');
	}
	return a;
}
void print(ostream &a, __int128 b) {
	if (b == 0) {
		return;
	}
	print(a, b / 10);
	a << int(b % 10);
}
ostream &operator<<(ostream &a, __int128 b) {
	print(a, b / 10);
	return a << int(b % 10);
}
#define int long long
#define lll __int128
#define MAXN 1000005
int prime[MAXN], minp[MAXN], phi[MAXN];
mint sum_phi_small[MAXN], sum_phiid_small[MAXN];
bool flag[MAXN];
unordered_map<int, mint> solved_phi, solved_phiid;;
mint solve_phi(int n) {
	if (n < MAXN) {
		return sum_phi_small[n];
	}
	if (solved_phi.count(n)) {
		return solved_phi[n];
	}
	mint ans = n * (n + 1) / 2;
	for (int l = 2, r = n / (n / 2); r <= n; l = r + 1, r = l <= n ? n / (n / l) : n + 1) {
		ans -= mint(r - l + 1) * solve_phi(n / l);
	}
	return solved_phi[n] = ans;
}
mint solve_phiid(int n) {
	if (n < MAXN) {
		return sum_phiid_small[n];
	}
	if (solved_phiid.count(n)) {
		return solved_phiid[n];
	}
	mint ans = mint(n) * (n + 1) * (2 * n + 1) / 6;
	for (int l = 2, r = n / (n / 2); r <= n; l = r + 1, r = l <= n ? n / (n / l) : n + 1) {
		ans -= mint(l + r) * mint(r - l + 1) / 2 * solve_phiid(n / l);
	}
	return solved_phiid[n] = ans;
}
signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.tie(0);
	flag[0] = flag[1] = minp[1] = phi[1] = 1;
	for (int i = 2; i < MAXN; i++) {
		if (!flag[i]) {
			prime[++prime[0]] = i;
			minp[i] = i;
		}
		if (i / minp[i] % minp[i] == 0) {
			phi[i] = phi[i / minp[i]] * minp[i];
		} else {
			phi[i] = phi[i / minp[i]] * (minp[i] - 1);
		}
		for (int j = 1; j <= prime[0] && prime[j] * (int)i < MAXN && prime[j] <= minp[i]; j++) {
			flag[i * prime[j]] = 1;
			minp[i * prime[j]] = prime[j];
		}
	}
	for (int i = 1; i < MAXN; i++) {
		sum_phi_small[i] = sum_phi_small[i - 1] + phi[i];
		sum_phiid_small[i] = sum_phiid_small[i - 1] + (mint)phi[i] * i;
	}
	lll n;
	cin >> n;
	int r = pow(n, 1.0 / 3.0) + 10;
	while ((lll)(r) * r * r > n) {
		r--;
	}
	set<int> r_prime_divisor, r_divisor;
	vector<pair<int, int> > p_factor;
	map<int, int> phi_divisor;
	int tr = r;
	for (int i = 2; i * i <= tr; i++) {
		if (tr % i == 0) {
			r_prime_divisor.insert(i);
			p_factor.push_back({i, 0});
			while (tr % i == 0) {
				p_factor.back().second++;
				tr /= i;
			}
		}
	}
	if (tr != 1) {
		r_prime_divisor.insert(tr);
		p_factor.push_back({tr, 1});
	}
	function<void(int, int, int)> dfs_get_divisor = [&](int x, int step, int now_phi)->void {
		if (step >= p_factor.size()) {
			r_divisor.insert(x);
			phi_divisor[x] = now_phi;
			return;
		}
		for (int i = 0; i <= p_factor[step].second; i++) {
			dfs_get_divisor(x, step + 1, now_phi);
			x *= p_factor[step].first;
			if (i == 0) {
				now_phi *= p_factor[step].first - 1;
			} else {
				now_phi *= p_factor[step].first;
			}
		}
		return;
	};
	dfs_get_divisor(1, 0, 1);
	mint answer = 0;
	for (pair<int, int> i : phi_divisor) {
		answer += (n / i.first - ((lll)r * r * r - 1) / i.first) % mod * i.second;
	}
	r--;
	for (int l = 1, R = 1; R <= r; l = R + 1, R = l <= r ? r / (r / l) : l + 1) {
		mint q = r / l;
		answer += q * (q + 1) * (q * 2 + 1) / 2 /* / 6 * 3 */ * (solve_phiid(R) - solve_phiid(l - 1));
		answer += (q * (q + 1) / 2 * 3 + q) * (solve_phi(R) - solve_phi(l - 1));
	}
	cans << answer << endl;
	return 0;
}