51nod 1220 约数之和
$\sigma_{0}(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[x\bot y]$
$\sigma_{1}(ij)=\sum\limits_{x|i}\sum\limits_{y|j}\dfrac{xj}{y}[x\bot y]$
所以 $$\begin{aligned}&\sum_{i=1}^n\sum_{j=1}^n \sigma_{1}(ij)\\=& \sum_{i=1}^n\sum_{j=1}^n \sum_{x|i}\sum_{y|j}\dfrac{xj}{y}[x\bot y]\\=& \sum_{d=1}^n\mu(d)\sum_{i=1}^n\sum_{j=1}^n\sum_{x|i}\sum_{y|j}\dfrac{xj}{y}[d|(x,y)]\\=&\sum_{d=1}^n\mu(d)\sum_{d|x}\sum_{d|y} \frac{x}{y}\sum_{x|i}\sum_{y|j}j \\=& \sum_{d=1}^n \mu(d)\sum_{d|x}\sum_{d|y}\frac{x}{y}s_0(\lfloor\frac{n}{x} \rfloor)ys_1(\lfloor\frac{n}{y} \rfloor)\\=& \sum_{i=1}^n \mu(d)d\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}is_0(\lfloor \frac{n}{id} \rfloor) \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}s_1(\lfloor \frac{n}{jd} \rfloor) \end{aligned}$$
$\sum \limits_{i=1}^n i\lfloor \dfrac{n}{i} \rfloor=\sum\limits_{i=1}^n \sigma_1(i)$
$\sum \limits_{i=1}^n s_1(\lfloor \dfrac{n}{i} \rfloor)=\sum\limits_{i=1}^n \sum\limits_{j=1}^{\lfloor \frac{n}{i} \rfloor}j=\sum_{i=1}^n i\lfloor \dfrac{n}{i} \rfloor$
二者等价
式子变成 $$\sum_{i=1}^n \mu(d)d \left(\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sigma_1(i)\right)^2$$
对于 $f(n)=n\mu(n)$,有 $f * I = \epsilon$,杜教筛之后的式子就是 $S(n)=1-\sum\limits_{i=2}^n iS(\lfloor \dfrac{n}{i} \rfloor)$
对于 $\sigma_1$ 的前缀和,小范围预处理,大的用 $\sum\limits_{i=1}^n \sigma_1(i)=\sum \limits_{i=1}^n i\lfloor \dfrac{n}{i} \rfloor$ 整除分块解决。
#include <bits/stdc++.h> const int N = 1e6 + 7; const int MOD = 1000000007; int mu[N], prime[N], prin, d[N], t[N]; bool vis[N]; inline void M(int &ans) { if (ans >= MOD) ans -= MOD; if (ans < 0) ans += MOD; } int qp(int a, int b = MOD - 2) { int ans = 1; while (b) { if (b & 1) ans = 1LL * ans * a % MOD; b >>= 1; a = 1LL * a * a % MOD; } return ans; } const int inv6 = qp(6); void init(int n) { mu[1] = 1; for (int i = 2; i <= n; i++) { if (!vis[i]) { prime[++prin] = i; mu[i] = -1; } for (int j = 1; j <= prin && i * prime[j] < N; j++) { vis[i * prime[j]] = 1; if (i % prime[j] == 0) { mu[i * prime[j]] = 0; break; } mu[i * prime[j]] = -mu[i]; } } for (int i = 1; i <= n; i++) { for (int j = i; j <= n; j += i) M(d[j] += i); } for (int i = 1; i <= n; i++) M(mu[i] = (1LL * mu[i] * i + mu[i - 1]) % MOD), M(d[i] += d[i - 1]); } std::unordered_map<int, int> muu, dd; inline int sum1(int n) { return 1LL * n * (n + 1) / 2 % MOD; } inline int sum1(int i, int j) { int ans; M(ans = sum1(j) - sum1(i - 1)); return ans; } inline int sum2(int n) { return 1LL * n * (n + 1) % MOD * (2 * n + 1) % MOD * inv6 % MOD; } inline int sum2(int i, int j) { int ans = sum2(j) - sum2(i - 1); M(ans); return ans; } inline int Mu(int n) { if (n < N) return mu[n]; if (muu.count(n)) return muu[n]; int ans = 1; for (int i = 2, j; i <= n; i = j + 1) { j = n / (n / i); M(ans -= 1LL * sum1(i, j) * Mu(n / i) % MOD); } return muu[n] = ans; } inline int D(int n) { if (n < N) return d[n]; if (dd.count(n)) return dd[n]; int ans = 0; for (int i = 1, j; i <= n; i = j + 1) { j = n / (n / i); M(ans += 1LL * sum1(i, j) * (n / i) % MOD); } return dd[n] = ans; } int solve(int n) { int ans = 0; for (int i = 1, j; i <= n; i = j + 1) { j = n / (n / i); M(ans += 1LL * (Mu(j) - Mu(i - 1) + MOD) % MOD * D(n / i) % MOD * D(n / i) % MOD); } return ans; } int main() { init(N - 1); int n; scanf("%d", &n); printf("%d\n", solve(n)); return 0; }