简单的数学题

题面

Solution1:

[egin{aligned} &sum_{i=1}^nsum_{i=1}^nijgcd(i,j) \ =&sum_{d=1}^ddsum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{n}{d}}ijd^2[ gcd(i, j)=1 ]\ =&sum_{d=1}^nd^3sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{n}{d}}ijsum_{k|(i,j)}mu(k)\ =&sum_{d=1}^{n}d^3sum_{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^3sum_{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})^2sum_{d|T}d^3(frac{T}{d})^2mu(frac{T}{d}) (ps:T=dk)\ =&sum_{T=1}^{n}(frac{(1+T)T}{2})^2T^2sum_{d|T}dmu(frac{T}{d})\ =&sum_{T=1}^{n}(frac{(1+T)T}{2})^2T^2varphi(T) (ps:Id* mu=varphi) end{aligned} ]

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

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

(g(n)=n^2),

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

[egin{aligned} &g(1)S(n)\ &=sum_{i=1}^n(f*g)(i) - sum_{i=2}^ng(i)S(frac{n}{i})\ &=sum_{i=1}^nsum_{d|i}d^2varphi(d)(frac{i}{d})^2-sum_{i=2}^ng(i)S(frac{n}{i})\ &=sum_{i=1}^ni^2sum_{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) 构造成这样)

Solution2:

[egin{aligned} &sum_{i=1}^nsum_{j=1}^nijgcd(i,j) \ &=sum_{i=1}^nsum_{j=1}^nijsum_{d|(i,j)}varphi(d)\ &=sum_{d=1}^nvarphi(d)d^2sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{n}{d}}ij\ &=sum_{d=1}^nvarphi(d)d^2(sum_{i=1}^{frac{n}{d}}i)^2\ &=sum_{d=1}^nvarphi(d)d^2sum_{i=1}^{frac{n}{d}}i^3 end{aligned} ]

同样杜教筛前面的。

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;
}