题意
求 (sumlimits_{i = 1}^N sumlimits_{j = 1}^M lcm (i, j))
Solution
易知,原式
[sumlimits_{i = 1}^N sumlimits_{j = 1}^M frac{ij}{gcd (i, j)}
]
枚举 (gcd (i, j)) ,且将 (d) 提出来得
[sumlimits_{d = 1}^{min (N, M)} d sumlimits_{i = 1}^{leftlfloorfrac{N}{d}
ight
floor} sumlimits_{j = 1}^{leftlfloorfrac{M}{d}
ight
floor} ij[(i, j) = 1]
]
将公式 (sumlimits_{k | n} mu(k) = [n = 1]) 代入,得
[sumlimits_{d = 1}^{min (N, M)} d sumlimits_{i = 1}^{leftlfloorfrac{N}{d}
ight
floor} sumlimits_{j = 1}^{leftlfloorfrac{M}{d}
ight
floor} ij sumlimits_{k | (i, j)} mu(k)
]
套路枚举 (k) ,得
[sumlimits_{d = 1}^{min (N, M)} d sumlimits_{k = 1}^{min (leftlfloorfrac{N}{d}
ight
floor, leftlfloorfrac{M}{d}
ight
floor)} mu(k) sumlimits_{i = 1}^{leftlfloorfrac{N}{d}
ight
floor} sumlimits_{j = 1}^{leftlfloorfrac{M}{d}
ight
floor} ij [k | (i, j)]
]
那么 (ij) 存在贡献时其必定是 (k) 的倍数,故
[sumlimits_{d = 1}^{min (N, M)} d sumlimits_{k = 1}^{min (leftlfloorfrac{N}{d}
ight
floor, leftlfloorfrac{M}{d}
ight
floor)} mu(k) sumlimits_{ki = 1}^{leftlfloorfrac{N}{d}
ight
floor} sumlimits_{kj = 1}^{leftlfloorfrac{M}{d}
ight
floor} k^2 ij
]
将 (k) 提出,得
[sumlimits_{d = 1}^{min (N, M)} d sumlimits_{k = 1}^{min (leftlfloorfrac{N}{d}
ight
floor, leftlfloorfrac{M}{d}
ight
floor)} k^2 mu(k) ( sumlimits_{i = 1}^{leftlfloorfrac{N}{kd}
ight
floor} i) (sumlimits_{j = 1}^{leftlfloorfrac{M}{kd}
ight
floor} j)
]
那么就可以预处理 (sumlimits_{k = 1}^n k^2 mu(k)) ,后面的用整除分块就好了
Code
#include <iostream>
#include <cstdio>
#include <cstring>
#define MOD 20101009
using namespace std;
typedef long long LL;
const int MAXN = 1e07 + 10;
int prime[MAXN];
int vis[MAXN]= {0};
int pcnt = 0;
int mu[MAXN]= {0};
LL sum[MAXN]= {0};
const int MAX = 1e07;
void prime_Acqu () {
mu[1] = 1;
for (int i = 2; i <= MAX; i ++) {
if (! vis[i]) {
prime[++ pcnt] = i;
mu[i] = - 1;
}
for (int j = 1; j <= pcnt && i * prime[j] <= MAX; j ++) {
vis[i * prime[j]] = 1;
if (! (i % prime[j]))
break;
mu[i * prime[j]] = - mu[i];
}
}
for (int i = 1; i <= MAX; i ++)
sum[i] = (sum[i - 1] + 1ll * i * 1ll * i % MOD * mu[i] % MOD) % MOD;
}
int N, M;
inline LL calc (int n) {
LL fn = (LL) n;
return (fn * (fn + 1) >> 1) % MOD;
}
LL Solve () {
LL ans = 0;
int limit = min (N, M);
for (int d = 1; d <= limit; d ++) {
LL total = 0;
int minlim = min (N / d, M / d);
for (int l = 1, r; l <= minlim; l = r + 1) {
r = min ((N / d) / ((N / d) / l), (M / d) / ((M / d) / l));
total = (total + (sum[r] - sum[l - 1] + MOD) % MOD * calc (N / d / l) % MOD * calc (M / d / l) % MOD) % MOD;
}
ans = (ans + (LL) (d) * total % MOD) % MOD;
}
return ans;
}
int main () {
prime_Acqu ();
scanf ("%d%d", & N, & M);
LL ans = Solve ();
cout << ans << endl;
return 0;
}
/*
4 5
*/