莫比乌斯反演习题：bzoj 2956 模积和

题目描述

求$\sum_{i=1}^{n}\sum_{j=1}^{m}(n \text{ mod } i) \times (m \text{ mod } j)$

其中$1 \le i \le n,1 \le j \le m,i \not = j$

题解

$$
\begin{aligned}
A
&=\sum_{i=1}^{n}\sum_{j=1}^{m}(n \text{ mod } i) \times (m \text{ mod } j) \\
&=(\sum_{i=1}^{n}n \text{ mod } i)(\sum_{j=1}^{m}m \text{ mod } j) \\
&=f(n)f(m)
\end{aligned}
$$

$$
\begin{aligned}
f(n)
&=\sum_{i=1}^{n}n \text{ mod } i \\
&=\sum_{i=1}^{n}n - i \lfloor \frac{n}{i} \rfloor \\
&=n^2-\sum_{i=1}^{n}i \lfloor \frac{n}{i} \rfloor
\end{aligned}
$$

$$
\begin{aligned}
g(k)
&=\sum_{i=1}^{n}i \lfloor \frac{k}{i} \rfloor
\end{aligned}
$$

$$
\begin{aligned}
B
&=\sum_{i=1}^{n}(n \text{ mod } i)(m \text{ mod } i) \\
&=\sum_{i=1}^{n}(n - i \lfloor \frac{n}{i} \rfloor)(m - i \lfloor \frac{m}{i} \rfloor) \\
&=\sum_{i=1}^{n}nm+i^2\lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor-im\lfloor \frac{n}{i} \rfloor-in\lfloor \frac{m}{i} \rfloor \\
&=n^2m+(\sum_{i=1}^{n}i^2\lfloor \frac{n}{i} \rfloor\lfloor \frac{m}{i} \rfloor)-(m\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor)-(n\sum_{i=1}^{n}i\lfloor \frac{m}{i} \rfloor) \\
&=n^2m+(\sum_{i=1}^{n}i^2\lfloor \frac{n}{i} \rfloor\lfloor \frac{m}{i} \rfloor)-mg(n)-ng(m)
\end{aligned}
$$

$$
\begin{aligned}
ans
&=A-B \\
&=f(n)f(m)+mg(n)+ng(m)-n^2m-(\sum_{i=1}^{n}i^2\lfloor \frac{n}{i} \rfloor\lfloor \frac{m}{i} \rfloor)
\end{aligned}
$$

注意$19940417=2848631 \times 7$，而且计算平方和的时候会爆long long，所以需要手动算逆元

代码

#include "bits/stdc++.h"
using namespace std;
typedef long long ll;
const int p = 19940417;

ll n, m;

ll g(ll k) {
    ll r = 0;
    for(ll i = 1, j ; i <= n ; i = j + 1) {
        if(k / i == 0) break;
        j = min(n, k / (k / i));
        r = (r + (i + j) * (j - i + 1) / 2 % p * (k / i) % p) % p;
    }
    return r;
}

ll f(ll n) {
    ll r = 0;
    for(ll i = 1, j ; i <= n ; i = j + 1) {
        j = n / (n / i);
        r = (r + (i + j) * (j - i + 1) / 2 % p * (n / i)) % p;
    }
    return (n * n % p - r % p) % p;
}

ll sum(ll n) {
    return n * (n + 1) % p * (2 * n + 1) % p * 3323403 % p;
}

ll w() {
    ll r = 0;
    for(int i = 1, j ; i <= n ; i = j + 1) {
        if(n / i == 0 || m / i == 0) break;
        j = min(n / (n / i), m / (m / i));
        r = (r + (sum(j) - sum(i - 1)) % p * (n / i) % p * (m / i) % p) % p;
    }
    return r;
}

int main() {
    scanf("%lld%lld", &n, &m); if(n > m) swap(n, m);
    ll ans = f(n) * f(m) % p + m * g(n) % p + n * g(m) % p - n * n % p * m % p - w() % p;
    printf("%lld\n", (ans % p + p) % p);
}