[Luogu3768]简单的数学题

Portal

\[求 \sum_{i = 1}^{n}\sum_{j = 1}^{n}(i,j)ij \]

---

好像直接利用\(\varphi\)很好做诶:

\[\sum_{i = 1}^{n}\sum_{j = 1}^{n}(i,j)ij\\ = \sum_{i = 1}^{n} i\sum_{j = 1}^{n}j \sum_{d | (i,j)} \varphi(d)\\ = \sum_{i = 1}^{n} i\sum_{j = 1}^{n}j \sum_{d | i,d | j} \varphi(d)\\ = \sum_{d = 1}^{n} d ^ 2\varphi(d) \sum_{i = 1}^{\lfloor \frac{n}{d}\rfloor}i\sum_{j = 1}^{\lfloor \frac{n}{d}\rfloor}j\\ \]

令$Sum(n) = \sum_{i = 1}^{n} i $

那么有：

\[原式 = \sum_{d = 1}^{n}d^2\varphi(d) Sum^2(\lfloor\frac{n}{d}\rfloor) \]

后面一部分可以整除分块。考虑前面一部分：

构造：\(g(n) = n ^ 2\)

\[(f * g)(n) = \sum_{d | n} d^2\varphi(d) * \frac{n ^ 2}{d ^ 2} \\ = n ^ 2 \sum_{d | n} \varphi(d) = n ^ 3 \]

有：\(\sum{n ^ 3} = (\sum{n})^2\)

然后杜教筛筛出前面的和，然后整除分块即可。

---

这里有莫比乌斯反演推式子的办法,可以参考一下

#include <bits/stdc++.h>
#include <bits/extc++.h>
using namespace std;
#define rep(i, a, b) for(LL i = (a), i##_end_ = (b); i <= i##_end_; ++i)
#define drep(i, a, b) for(LL i = (a), i##_end_ = (b); i >= i##_end_; --i)
#define clar(a, b) memset((a), (b), sizeof(a))
#define debug(...) fprintf(stderr, __VA_ARGS__)
typedef long long LL;
typedef long double LD;
LL read() {
    char ch = getchar();
    LL x = 0, flag = 1;
    for (;!isdigit(ch); ch = getchar()) if (ch == '-') flag *= -1;
    for (;isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
    return x * flag;
}

__gnu_pbds :: gp_hash_table <LL, LL> Sum;
const int Maxn = 10000009;
LL isnprime[Maxn], prime[Maxn], tot, phi[Maxn];
LL prefix[Maxn], Mod, n;

LL fpm(LL a, LL tims) {
	LL r = 1;
	for (; tims; tims >>= 1, a = a * a % Mod) 
		if (tims & 1) r = r * a % Mod;
	return r;
}

LL inv(LL a) { return fpm(a, Mod - 2); }

void init() {
	Mod = read();
	phi[1] = 1;
	rep (i, 2, Maxn - 1) {
		if (!isnprime[i]) 
			prime[++tot] = i, phi[i] = i - 1;

		for (int j = 1, k; j <= tot && (k = i * prime[j]) < Maxn; ++j) {
			isnprime[k] = 1;
			if (i % prime[j] == 0) {
				phi[k] = phi[i] * prime[j];
				break;
			} else phi[k] = phi[prime[j]] * phi[i];
		}
	}

	rep (i, 1, Maxn - 1) prefix[i] = (prefix[i - 1] + (phi[i] * i % Mod * i % Mod)) % Mod;
}

int inv2, inv6;
inline LL singleSum(LL val) { return val % Mod * (val % Mod + 1) % Mod * inv2 % Mod; }
inline LL doubleSum(LL val) { return (val % Mod * (val % Mod + 1) % Mod * (val % Mod * 2 + 1) % Mod) * inv6 % Mod; }
inline LL tripleSum(LL val) { return singleSum(val) * singleSum(val) % Mod; }

LL calcSum(LL val) {
	if (val < Maxn) return prefix[val];
	if (Sum[val]) return Sum[val];

	LL ans = tripleSum(val);
	for (LL l = 2, r; l <= val; l = r + 1) {
		r = val / (val / l);
		LL res = ((doubleSum(r) - doubleSum(l - 1)) % Mod + Mod) % Mod;
		(ans += (Mod - res * calcSum(val / l) % Mod) % Mod) %= Mod;
	}

	return Sum[val] = (ans + Mod) % Mod;
}

void solve() {
	inv2 = inv(2), inv6 = inv(6);
	n = read();

	LL ans = 0; 
	for (LL l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		(ans += tripleSum(n / l) * ((calcSum(r) - calcSum(l - 1) + Mod) % Mod) % Mod) %= Mod;
	}

	cout << ans << endl;
}

int main() {
	freopen("LG3768.in", "r", stdin);
	freopen("LG3768.out", "w", stdout);

	init();
	solve();

#ifdef Qrsikno
    debug("\nRunning time: %.3lf(s)\n", clock() * 1.0 / CLOCKS_PER_SEC);
#endif
    return 0;
}