莫比乌斯反演习题: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^2lfloor frac{n}{i} rfloor lfloor frac{m}{i} rfloor-imlfloor frac{n}{i} rfloor-inlfloor frac{m}{i} rfloor \
&=n^2m+(sum_{i=1}^{n}i^2lfloor frac{n}{i} rfloorlfloor frac{m}{i} rfloor)-(msum_{i=1}^{n}ilfloor frac{n}{i} rfloor)-(nsum_{i=1}^{n}ilfloor frac{m}{i} rfloor) \
&=n^2m+(sum_{i=1}^{n}i^2lfloor frac{n}{i} rfloorlfloor 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^2lfloor frac{n}{i} rfloorlfloor frac{m}{i} rfloor)
end{aligned}
$$

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

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

using namespace std;
typedef long long ll;
const int p = 19940417;

ll n, m;

ll (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("%lldn", (ans % p + p) % p);
}