题意
求
[sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} ext{lcm}(i,j)(mod 20101009)
]
思路
容易想到原式等价于
[sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}frac{i* j}{gcd(i,j)}
]
枚举(i,j)的最大公约数(d)，显然(gcd(frac id,frac jd)=1)，即(frac id)和(frac jd)互质
[sumlimits_{i=1}^{n}sumlimits_{j=1}^msumlimits_{d|i,d|j,gcd(frac id,frac jd)=1}frac{i*j}d
]
变换求和顺序
[sumlimits_{d=1}^{n}dsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor}[gcd(i,j)=1]i*j
]
记(sum(n,m)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=1]i*j)
对其进行化简，用(varepsilon(gcd(i,j)))替换([gcd(i,j)=1])
[sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}sumlimits_{d|gcd(i,j)}mu(d)*i*j
]
转化为首先枚举约数
[sumlimits_{d=1}^{min(n,m)}sumlimits_{d|i}^{n}sumlimits_{d|j}^{m}mu(d)*i*j
]
设(i=i'*d,j=j'*d)，则可以进一步转化
[sumlimits_{d=1}^{min(n,m)}mu(d)*d^2*sumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor}i*j
]
前半段可以处理前缀和，后半段可以(O(1))求，设
[Q(n,m)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}i*j=frac{n*(n+1)}{2}*frac{m*(m+1)}{2}
]
显然可以(O(1))求解
到现在
[sum(n,m)=sumlimits_{d=1}^{min(n,m)}mu(d)*d^2*Q(lfloorfrac nd
floor,lfloorfrac md
floor)
]
可以用数论分块求解
回带到原式中
[sumlimits_{d=1}^{min(n, m)}d*sum(lfloorfrac nd
floor,lfloorfrac md
floor)
]
又可以数论分块求解了
然后就做完啦
代码
/*
Author:loceaner
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define int long long
using namespace std;
const int A = 1e7 + 11;
const int B = 1e6 + 11;
const int mod = 20101009;
const int inf = 0x3f3f3f3f;
inline int read() {
char c = getchar(); int x = 0, f = 1;
for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
return x * f;
}
bool vis[A];
int n, m, mu[A], p[B], sum[A], cnt;
void getmu() {
mu[1] = 1;
int k = min(n, m);
for (int i = 2; i <= k; i++) {
if (!vis[i]) p[++cnt] = i, mu[i] = -1;
for (int j = 1; j <= cnt && i * p[j] <= k; ++j) {
vis[i * p[j]] = 1;
if (i % p[j] == 0) break;
mu[i * p[j]] = -mu[i];
}
}
for (int i = 1; i <= k; i++) sum[i] = (sum[i - 1] + i * i % mod * mu[i]) % mod;
}
int Sum(int x, int y) {
return (x * (x + 1) / 2 % mod) * (y * (y + 1) / 2 % mod) % mod;
}
int solve2(int x, int y) {
int res = 0;
for (int i = 1, j; i <= min(x, y); i = j + 1) {
j = min(x / (x / i), y / (y / i));
res = (res + 1LL * (sum[j] - sum[i - 1] + mod) * Sum(x / i, y / i) % mod) % mod;
}
return res;
}
int solve(int x, int y) {
int res = 0;
for (int i = 1, j; i <= min(x, y); i = j + 1) {
j = min(x / (x / i), y / (y / i));
res = (res + 1LL * (j - i + 1) * (i + j) / 2 % mod * solve2(x / i, y / i) % mod) % mod;
}
return res;
}
signed main() {
n = read(), m = read();
getmu();
cout << solve(n, m) << '
';
}