[洛谷 P3768] 简单的数学题 (莫比乌斯反演 狄利克雷卷积 杜教筛)

题意：

求

$$\sum _{i=1}^n\sum _{j=1}^{n}ij\gcd (i,j)$$

对 $p$ 取模，$n\le 10^{10},5\times 10^8\le p\le 1.1\times 10^9$

分析：

$$\sum _{i=1}^n\sum _{j=1}^{n}ij\gcd (i,j)$$

枚举 $\gcd (i,j)$

$$\sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^{n}ij[\gcd (i,j)=d]$$

利用 $\gcd$ 的性质

$$\sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^{n}ij[\gcd (\frac{i}{d},\frac{j}{d})=1]$$

让式子除 $d^2$ 再乘 $d^2$

$$\sum _{d=1}^n{d}^3\sum _{i=1}^n\sum _{j=1}^n\frac{i}{d}\cdot \frac{j}{d}[\gcd (\frac{i}{d},\frac{j}{d})=1]$$

换一下上界

$$\sum _{d=1}^n{d}^3\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }ij[\gcd (i,j)=1]$$

莫比乌斯反演

$$\sum _{d=1}^n{d}^3\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }ij\sum _{k\mid \gcd (i,j)}\mu (k)$$

交换求和次序

$$\sum _{d=1}^n{d}^3\sum _{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu (k)\sum _{i=1}^{\lfloor \frac{n}{dk}\rfloor }ik\sum _{j=1}^{\lfloor \frac{n}{dk}\rfloor }jk$$

化简后半部分

$$\sum _{d=1}^n{d}^3\sum _{k=1}^{\lfloor \frac{n}{d}\rfloor }{k}^2\mu (k)(\sum _{i=1}^{\lfloor \frac{n}{dk}\rfloor }i{)}^2$$

设 $T=dk$

$$\sum _{T=1}^n{T}^2\sum _{d\mid T}d\mu (\frac{T}{d})(\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }i{)}^2$$

因为 $\mu \ast Id=\phi$，替换 $\sum _{d\mid T}d\mu (\frac{T}{d})$ 得

$$\sum _{T=1}^n{T}^2\phi (T)(\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }i{)}^2$$

根据 $(\frac{n^2+n}{2}{)}^2=1^3+2^3+\cdots +n^3$

$$\sum _{T=1}^n{T}^2\phi (T)\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }{i}^3$$

设 $f(x)=x^2\phi (x)$ 套用杜教筛，$S(n)$ 为 $\sum _{i=1}^{n}f(i)$

$$g(1)S(n)=\sum _{i=1}^{n}h(i)-\sum _{i=1}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor )$$

那么 $h=f\ast g$，也就是 $h(n)=\sum _{d\mid n}f(d)g(\frac{n}{d})$，带入 $f(n)$ 得

$$h(n)=\sum _{d\mid n}{d}^2\phi (d)g(\frac{n}{d})$$

考虑把 $d^2$ 消去，所以设 $g(x)=x^2$，故

$$h(n)=n^2\sum _{d\mid n}\phi (d)$$

根据 $\phi \ast I=Id$

$$h(n)=n^3$$

那么 $f(x)$ 的前缀和就是

$$S(n)=(\frac{n^2+n}{2}{)}^2-\sum _{i=2}^n{i}^{2}S(\lfloor \frac{n}{i}\rfloor )$$

代码：

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 5e6 + 5;
int mod, n, euler[N], primes[N], cnt, sum[N], inv, res;
bool st[N];
unordered_map<int, int> mp;
int qmi(int a, int b) {
    int res = 1;
    while (b) {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}
void get_eulers(int n) {
    euler[1] = 1;
    for (int i = 2; i <= n; i ++) {
        if (!st[i]) {
            primes[cnt ++] = i;
            euler[i] = i - 1;
        }
        for (int j = 0; primes[j] <= n / i; j ++) {
            int t = primes[j] * i;
            st[t] = 1;
            if (i % primes[j] == 0) {
                euler[t] = primes[j] * euler[i];
                break;
            }
            euler[t] = (primes[j] - 1) * euler[i];
        }
    }
    for (int i = 1; i <= n; i ++) sum[i] = (sum[i - 1] + i * i % mod * euler[i] % mod) % mod;
}
int s2(int n) {
    n %= mod;
    return n * (n + 1) % mod * (2 * n + 1) % mod * inv % mod;
}
int s3(int n) {
    n %= mod;
    return (n * (n + 1) / 2) % mod * ((n * (n + 1) / 2) % mod) % mod;
}
int Sum(int n) {
    if (n < N) return sum[n];
    if (mp[n]) return mp[n];
    int res = s3(n);
    for (int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res -= (s2(r) - s2(l - 1) + mod) % mod * Sum(n / l) % mod;
        res = (res + mod) % mod;
    }
    return mp[n] = res;
}
signed main() {
    cin >> mod >> n;
    get_eulers(N - 1);
    inv = qmi(6, mod - 2);
    for (int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res += (Sum(r) - Sum(l - 1) + mod) % mod * s3(n / l) % mod;
        res = (res % mod + mod) % mod;
    }
    cout << res << endl;
}