Description
求
[sum_{i=1}^{n}sum_{j=1}^{m} (nmod i)(mmod j)[i
e j]
]
Solution
看到取模就会想整除分块的吧，那个([i
e j])其实也很好处理啊，就是最后再减去(i=j)的情况就行了。
我式子推的有点丑……
[sum_{i=1}^{n}sum_{j=1}^{m} (nmod i)(mmod j)[i
e j]\
= sum_{i=1}^{n}sum_{j=1}^{m} (nmod i)(mmod j) - sum_{i=1}^{min(n, m)} (nmod i)(mmod i)\
= sum_{i=1}^{n}sum_{j=1}^{m} (n - lfloordfrac{n}{i}
floor i)(m - lfloordfrac{m}{j}
floor j) - sum_{i=1}^{min(n, m)} (n - lfloordfrac{n}{i}
floor i)(m - lfloordfrac{m}{i}
floor i)\
= sum_{i=1}^{n}sum_{j=1}^{m} (nm - mlfloordfrac{n}{i}
floor i - nlfloordfrac{m}{j}
floor j + lfloordfrac{n}{i}
floor ilfloordfrac{m}{j}
floor j) - sum_{i=1}^{min(n, m)}mbox{...}\
= n^2m^2 - m^2sum_{i=1}^{n}lfloordfrac{n}{i}
floor i - mbox{...}
]
后面就是一些sb的去括号了，没啥好说的了。
Code
#include <cstdio>
#include <algorithm>
typedef long long LL;
const LL MOD = 19940417;
const LL I2 = 9970209;
const LL I6 = 3323403;
int n, m;
LL sp(LL x) {
return x * (2 * x + 1) % MOD * (x + 1) % MOD * I6 % MOD;
}
int main() {
scanf("%d%d", &n, &m);
if (n < m) std::swap(n, m);
LL ans = n % MOD * n % MOD * m % MOD * m % MOD;
LL si = 0, sj = 0, sij = 0;
for (int i = 1, j; i <= n; i = j+1) {
j = n / (n / i);
si = (si + (n / i) * (i + j) % MOD * (j - i + 1) % MOD * I2 % MOD) % MOD;
}
for (int i = 1, j; i <= m; i = j+1) {
j = m / (m / i);
sj = (sj + (m / i) * (i + j) % MOD * (j - i + 1) % MOD * I2 % MOD) % MOD;
}
ans = ((ans - m % MOD * m % MOD * si % MOD - n % MOD * n % MOD * sj % MOD) %MOD + MOD) % MOD;
sij = si * sj % MOD;
ans = (ans + sij) % MOD;
ans = (ans - m % MOD * m % MOD * n % MOD + MOD) % MOD;
si = 0;
for (int i = 1, j; i <= m; i = j+1) {
j = std::min(m, n / (n / i));
si = (si + (n / i) * (i + j) % MOD * (j - i + 1) % MOD * I2 % MOD) % MOD;
}
ans = (ans + m % MOD * si % MOD + n % MOD * sj % MOD) % MOD;
sj = 0;
for (int i = 1, j; i <= m; i = j+1) {
j = std::min(n / (n / i) , m / (m / i));
sj = (sj + (n / i) % MOD * (m / i) % MOD * (sp(j) - sp(i-1) + MOD) % MOD) % MOD;
}
ans = (ans - sj + MOD) % MOD;
printf("%d
", ans);
return 0;
}