LOJ 6682 梦中的数论

题目链接
\(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{k=1}^{n}[(j\mid i) \land ((j+k)\mid i)]\)
\(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{d=j+1}^{j+n}[(j\mid i) \land (d\mid i)]\)

把枚举后验证变为直接计算贡献，显然只有当\(j,d\)是\(i\)的约数时才有贡献。
从\(i\)的约数中选两个大小不同的约数的选法有\(\binom{\sigma-0 i}{2}\)种。
即求\(\sum\limits_{i=1}^n\frac{\sigma-0^2 i-\sigma-0 i}{2}\)

因为\(\sigma\)是积性函数，用数论分块和\(\text{min25}\) 分别求解即可。

#include <cstdio>
#include <cmathjax>

typedef long long ll;
const int N = 1e6 + 5,  MOD = 998244353;
const ll i2 = 499122177;
ll n, id1[N], id2[N], w[N], g[N], p[N], sqr;
int m, cnt;
bool vis[N];

ll add(ll a, ll b) {a %= MOD, b %= MOD; return (a + b >= MOD) ? a + b - MOD : a + b;}
ll mul(ll a, ll b) {a %= MOD, b %= MOD; return a * b % MOD;}
ll dec(ll a, ll b) {a %= MOD, b %= MOD; return ((a - b) % MOD + MOD) % MOD;}

void init(int m) {
    for (int i = 2; i <= m; i++) {
		if (!vis[i]) p[++cnt] = i;
		for (int j = 1; j <= cnt && i * p[j] <= m; j++) {
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) break;
		}
	}
}

ll solve() {
	ll res = 0;
	for (ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		res = add(res, mul(r - l + 1, n / l));
	}
	return res;
}
ll S(ll x, int y) {
	if (p[y] > x || x <= 1 ) return 0;
	int k = (x <= sqr) ? id1[x] : id2[n / x];
	ll res = dec(g[k], mul(4, y - 1));
	for (int i = y; i <= cnt && p[i] * p[i] <= x; i++) {
		ll pow1 = p[i], pow2 = p[i] * p[i];
		for (int e = 1; pow2 <= x; pow1 = pow2, pow2 *= p[i], e++) {
			res =  add(res, add(mul(mul(e + 1, e + 1), S(x / pow1, i + 1)), mul(e + 2, e + 2)));
		}
	}
	return res;
}

int main() {
    scanf("%lld", &n);
	sqr = int(sqrt(n)) + 1;
	init(sqr);
	for (ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		w[++m] = n / l;
		g[m] = dec(w[m], 1);
	    (w[m] <= sqr) ? id1[w[m]] = m : id2[r] = m;
	}
	for (int j = 1; j <= cnt; j++)
		for (int i = 1; i <= m && p[j] * p[j] <= w[i]; i++) {
			int k = (w[i] / p[j] <= sqr) ? id1[w[i] / p[j]] : id2[n / (w[i] / p[j])];
			g[i] = dec(g[i], dec(g[k], j - 1));
		}
	for (int i = 1; i <= m; i++) g[i] = mul(g[i], 4);
	printf("%lld\n", mul(dec(add(S(n, 1), 1), solve()), i2));
	return 0;
}