题面

Solution1:

\[ \begin{aligned} &\sum_{i=1}^n\sum_{i=1}^nijgcd(i,j) \\ =&\sum_{d=1}^dd\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ijd^2[\ gcd(i, j)=1\ ]\\ =&\sum_{d=1}^nd^3\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij\sum_{k|(i,j)}\mu(k)\\ =&\sum_{d=1}^{n}d^3\sum_{k=1}^{\frac{n}{d}}\mu(k)\sum_{i=1}^{\frac{n}{dk}}\sum_{i=1}^{\frac{n}{dk}}ijk^2\\ =&\sum_{d=1}^{n}d^3\sum_{k=1}^{\frac{n}{d}}\mu(k)k^2(1+2+3+...+\frac{n}{dk})^2\\ =&\sum_{T=1}^{n}(\frac{(1+T)T}{2})^2\sum_{d|T}d^3(\frac{T}{d})^2\mu(\frac{T}{d})\ \ \ \ \ \ \ (ps:T=dk)\\ =&\sum_{T=1}^{n}(\frac{(1+T)T}{2})^2T^2\sum_{d|T}d\mu(\frac{T}{d})\\ =&\sum_{T=1}^{n}(\frac{(1+T)T}{2})^2T^2\varphi(T)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (ps:Id* \mu=\varphi) \end{aligned} \]

对于前面的整数分块，后面的考虑用杜教筛快速筛出来。

设 \(f(n)=n^2\varphi(n)\),

\(g(n)=n^2\),

\(S(n)=\sum_{i=1}^nf(i)\)

\[ \begin{aligned} &g(1)S(n)\\ &=\sum_{i=1}^n(f*g)(i) - \sum_{i=2}^ng(i)S(\frac{n}{i})\\ &=\sum_{i=1}^n\sum_{d|i}d^2\varphi(d)(\frac{i}{d})^2-\sum_{i=2}^ng(i)S(\frac{n}{i})\\ &=\sum_{i=1}^ni^2\sum_{d|i}\varphi(d)-\sum_{i=2}^ng(i)S(\frac{n}{i})\ \ \ \ \ \ \ \ \ \ \ \ (ps:1*\varphi=Id)\\ &=\sum_{i=1}^ni^3-\sum_{i=1}^ni^2S(\frac{n}{i}) \end{aligned} \]

然后直接杜教筛.(请仔细思考为什么将 \(g\) 构造成这样)

Source:

#include <iostream>
#include <cstdio>
#include <map>
#include <cstring>

using namespace std;

#define LL long long
#define int long long

const int maxN = 1e6 + 1;

LL Mod, inv2, inv4, n, inv6;
LL qpow(LL a, LL b, LL p) {
    LL res = 1;
    while (b) { if (b & 1) res = res * a % p; a = a * a % p; b >>= 1; }
    return res ;
}
bool vis[maxN];
int prime[348514], tot;
LL f[maxN];
void Init() {
    f[1] = 1;
    for (int i = 2; i <= maxN; ++ i) {
        if (!vis[i]) {
            prime[++tot] = i;
            f[i] = i - 1;
        }
        for (int j = 1; j <= tot && prime[j] * i <= maxN; ++ j) {
            vis[prime[j] * i] = 1;
            if (i % prime[j] == 0) {
                f[i * prime[j]] = f[i] * prime[j];
                break;
            } else 
                f[i * prime[j]] = f[i] * f[prime[j]];
        }
    }
    for (int i = 1; i <= maxN; ++ i) f[i] = (f[i - 1] + f[i] * i % Mod * i % Mod) % Mod;
}

LL Sum1(LL x) { x %= Mod; return (1 + x) * x % Mod * inv2 % Mod; }
LL Sum2(LL x) { x %= Mod; return (1 + x) * x % Mod * (x + x + 1) % Mod * inv6 % Mod; } 
LL Sum3(LL x) { x %= Mod; return x * x % Mod * (x + 1) % Mod * (x + 1) % Mod * inv4 % Mod; }

map<LL, LL> F;

LL Sf(LL x) {
    if (x <= maxN) return f[x];
    if (F[x]) return F[x];
    LL res = 0;
    for (int l = 2, r; l <= x; l = r + 1) {
        r = x / (x / l);
        res = (res + (Sum2(r) - Sum2(l - 1) + Mod) % Mod * Sf(x / l)) % Mod;
    }
    return F[x] = (Sum3(x) - res + Mod) % Mod;
}

signed main() {
#ifndef ONLINE_JUDGE
    freopen("P3768.in", "r", stdin);
    freopen("P3768.out", "w", stdout);
#endif
    cin >> Mod >> n;
    Init();
    inv2 = qpow(2, Mod - 2, Mod);
    inv4 = qpow(4, Mod - 2, Mod);
    inv6 = qpow(6, Mod - 2, Mod);
    LL ans = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        LL a = Sum1(n / l);
        a = a * a % Mod;
        ans = (ans + a * (Sf(r) - Sf(l - 1) + Mod) % Mod) % Mod;
    }
    cout << ans << endl;
}