题目描述
求(对 (20101009) 取模，(n,mle10^7) )
[sum_{i=1}^nsum_{j=1}^moperatorname{lcm}(i,j)
]
大体思路
推式子:
[egin{aligned}sum_{i=1}^nsum_{j=1}^moperatorname{lcm}(i,j)
&=sum_{i=1}^nsum_{j=1}^mfrac{i imes j}{gcd(i,j)} \
&=sum_{i=1}^nsum_{j=1}^msum_{d|i,d|j,gcd(i/d,j/d)=1}frac{i imes j}{d} \
&=sum_{d=1}^{min(n,m)} imes d imessum_{i=1}^{lfloor frac{n}{d}
floor}sum_{j=1}^{lfloor frac{m}{d}
floor}[gcd(i,j)=1] imes i imes jend{aligned}]
把式子后面那一大堆设为 (sum(n,m)) ：
[sum(n,m)=sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=1] imes i imes j
]
考虑化简一下 (sum) :
[egin{aligned}sum(n,m) &= sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=1] imes i imes j\
&=sum_{i=1}^nsum_{j=1}^msum_{d|gcd(i,j)}mu(d) imes i imes j \
&=sum_{d=1}^{min(n,m)}mu(d) imes d^2sum_{i=1}^{lfloor frac{n}{d}
floor}sum_{j=1}^{lfloor frac{m}{d}
floor}i imes jend{aligned}]
可以发现 (sum) 后面那一大堆(设为 (g(n,m)) )可以 (O(1)) 求：
[egin{aligned}g(n,m)&=sum_{i=1}^nsum_{j=1}^m i imes j \
&=frac{n imes(n+1)}{2} imes frac{m imes(m+1)}{2}end{aligned}]
那么 (sum(n,m)) 可以化为：
[sum(n,m)=sum_{d=1}^{min(n,m)}mu(d) imes d^2 imes g(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)
]
这个可以数论分块 (lfloorfrac{n}{lfloorfrac{n}{d}
floor}
floor) 求。
再回到定义 (sum) 的地方，那么：
[Ans=sum_{d=1}^{min(n,m)} imes d imes sum(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)
]
好像这个还是可以数论分块 (QwQ)
至此这道题就解决了。
细节注意事项
- (long long)一定要开呀。
- 不要写挂呀！！！
参考代码
/*--------------------------------
Code name: crash.cpp
Author: The Ace Bee
This code is made by The Ace Bee
--------------------------------*/
#include <cstdio>
#define rg register
#define int long long
#define fileopen(x)
freopen(x".in", "r", stdin);
freopen(x".out", "w", stdout);
#define fileclose
fclose(stdin);
fclose(stdout);
const int mod = 20101009;
const int MAXN = 10000010;
inline int min(int a, int b) { return a < b ? a : b; }
inline int read() {
int s = 0; bool f = false; char c = getchar();
while (c < '0' || c > '9') f |= (c == '-'), c = getchar();
while (c >= '0' && c <= '9') s = (s << 3) + (s << 1) + (c ^ 48), c = getchar();
return f ? -s : s;
}
int vis[MAXN], mu[MAXN];
int num, pri[MAXN], sum[MAXN];
inline void seive() {
mu[1] = 1;
for (rg int i = 2; i < MAXN; ++i) {
if (!vis[i]) mu[i] = -1, pri[++num] = i;
for (rg int j = 1; j <= num && i * pri[j] < MAXN; ++j) {
vis[i * pri[j]] = 1;
if (i % pri[j]) mu[i * pri[j]] = - mu[i];
else { mu[i * pri[j]] = 0; break; }
}
}
for (rg int i = 1; i < MAXN; ++i)
sum[i] = (sum[i - 1] + 1ll * i * i % mod * (mu[i] + mod) % mod) % mod;
}
inline int g(int n, int m)
{ return 1ll * n * (n + 1) / 2 % mod * (m * (m + 1) / 2 % mod) % mod; }
inline int f(int n, int m) {
int res = 0;
for (rg int i = 1, j; i <= min(n, m); i = j + 1) {
j = min(n / (n / i), m / (m / i));
res = (res + 1ll * (sum[j] - sum[i - 1] + mod) * g(n / i, m / i) % mod) % mod;
}
return res;
}
inline int solve(int n, int m) {
int res = 0;
for (rg int i = 1, j; i <= min(n, m); i = j + 1) {
j = min(n / (n / i), m / (m / i));
res = (res + 1ll * (j - i + 1) * (i + j) / 2 % mod * f(n / i, m / i) % mod) % mod;
}
return res;
}
signed main() {
// fileopen("crash");
seive();
int n = read(), m = read();
printf("%lld
", solve(n, m));
// fileclose;
return 0;
}
完结撒花(qwq)