Loj #6682. 梦中的数论（数论，min_25筛）

---

相当于求 $\displaystyle\sum_{i = 1}^nC(\sigma(i),2)$i=1∑n​C(σ(i),2)，拆开就是 $\displaystyle\frac{1}{2}(\sum_{i = 1}^n\sigma(i)^2 - \sum_{i = 1}^n\sigma(i))$21​(i=1∑n​σ(i)2−i=1∑n​σ(i))

前面一半可以用25筛筛出来，后面一半可以直接分块（也可以25筛筛出来）

---

代码：

#include<bits/stdc++.h>
using namespace std;
const int maxn = 2e6 + 10;
typedef long long ll;
const int mod = 998244353;
int num,tot;
bool ispri[maxn];
ll sum[maxn],pri[maxn],w[maxn],g[maxn],id1[maxn],id2[maxn];
ll n,sqr;
void sieve(int x) {
	num = 0; ispri[0] = ispri[1] = true;
	for (int i = 2; i <= x; i++) {
		if (!ispri[i]) {
			pri[++num] = i;
			sum[num] = (sum[num - 1] + 1) % mod;
		}
		for (int j = 1; j <= num && i * pri[j] <= x; j++) {
			ispri[i * pri[j]] = true;
			if (i % pri[j] == 0) break;
		}
	}
}
ll fpow(ll a,ll b) {
	ll r = 1;
	while (b) {
		if (b & 1) r = r * a % mod;
		b >>= 1;
		a = a * a % mod;
	}
	return r;
}
ll S1(ll x,int y) {
	if (pri[y] >= x) return 0;
	ll k = x <= sqr ? id1[x] : id2[n / x];
	ll ans = 2 * (g[k] - sum[y] + mod) % mod;
	for (int i = y + 1; i <= num && pri[i] * pri[i] <= x; i++) {
		for (ll e = 1, pe = pri[i]; pe <= x; e++, pe = pe * pri[i]) {
			ans += (e + 1) * (S1(x / pe,i) + (e != 1)) % mod;
			ans %= mod;
		}
	}
	return ans;
}
ll S2(ll x,int y) {
	if (pri[y] >= x) return 0;
	ll k = x <= sqr ? id1[x] : id2[n / x];
	ll ans = 4 * (g[k] - sum[y] + mod) % mod;
	for (int i = y + 1; i <= num && pri[i] * pri[i] <= x; i++) {
		for (ll e = 1, pe = pri[i]; pe <= x; e++, pe = pe * pri[i]) {
			ans += (e + 1) * (e + 1) % mod * (S2(x / pe,i) + (e != 1)) % mod;
			ans %= mod;
		}
	}
	return ans;
}
int main() {
	sieve(maxn - 10);
	scanf("%lld",&n);
	sqr = sqrt(n);
	for (ll i = 1,j; i <= n; i = j + 1) {
		j = n / (n / i);
		w[++tot] = n / i;
		ll p = n / i % mod;
		g[tot] = (p - 1) % mod;
		if (n / i <= sqr) id1[n / i] = tot;
		else id2[j] = tot;
	}
	for (int i = 1; i <= num; i++) {
		for (int j = 1; j <= tot && 1ll * pri[i] * pri[i] <= w[j]; j++) {
			ll x = w[j] / pri[i];
			x = x <= sqr ? id1[x] : id2[n / x];
			g[j] -= (g[x] - sum[i - 1] + mod) % mod;
			if (g[j] < 0) g[j] += mod;
		}
	}
	ll res = (S2(n,0) - S1(n,0) + mod) * fpow(2,mod - 2) % mod;
	printf("%lld\n",res);
	return 0;
}