题目
求(sumlimits^{n}_{i=1}sumlimits^{m}_{j=1}{lcm(i, j)})
答案模100000009输出
题解
注意100000009是1e8级别,而且不是质数,不要想逆元之类的东西。
老老实实推式子,一步一步小心不要出错。注意在提出项时可能要平方之类的。
[sumlimits_{d_1}sumlimits_{i=1}^nsumlimits_{j=1}^m{frac{ij}{d_1}[gcd(i,j)=d_1]}
]
提出(d_1)
[sumlimits_{d_1}sumlimits_{i=1}^{frac{n}{d_1}}sumlimits_{j=1}^{frac{m}{d_1}}{ijd_1[gcd(i,j)=1]}
]
[sumlimits_{d_1}sumlimits_{i=1}^{frac{n}{d_1}}sumlimits_{j=1}^{frac{m}{d_1}}{ijd_1sumlimits_{d_2|gcd(i,j)}{mu(d_2)}}
]
[sumlimits_{d_1}sumlimits_{d_2}{mu(d_2)d_1}sumlimits_{i=1}^{frac{n}{d_1}}{i[d_2|i]}sumlimits_{j=1}^{frac{m}{d_1}}{j[d_2|j]}
]
提出(d_2)
[sumlimits_{d_1}sumlimits_{d_2}{mu(d_2)d_1d_2^2}sumlimits_{i=1}^{frac{n}{d_1d_2}}{i}sumlimits_{j=1}^{frac{m}{d_1d_2}}{j}
]
令(d_1d_2=T)
[sumlimits_{T=1}^{min(n,m)}(sumlimits_{d|T}{mu(d)dT})sumlimits_{i=1}^{frac{n}{T}}{i}sumlimits_{j=1}^{frac{m}{T}}{j}
]
令(F(T)=sumlimits_{d|T}{mu(d)dT}),求其前缀和,就可以数论分块了。
#include <bits/stdc++.h>
#define endl '
'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define mp make_pair
#define seteps(N) fixed << setprecision(N)
typedef long long ll;
using namespace std;
/*-----------------------------------------------------------------*/
ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f
const int N = 1e7 + 10;
const double eps = 1e-5;
const int M = 1e8 + 9;
int mu[N];
ll f[N];
int prime[N];
int cnt;
int vis[N];
void init() {
mu[1] = 1;
for(int i = 2; i < N; i++) {
if(!vis[i]) {
mu[i] = -1;
prime[cnt++] = i;
}
for(int j = 0; j < cnt; j++) {
int p = prime[j];
if(i * p > N) break;
vis[i * p] = 1;
if(i % p == 0) {
mu[i * p] = 0;
break;
}
mu[i * p] = -mu[i];
}
}
for(int i = 1; i < N; i++) {
for(int j = i; j < N; j += i) {
f[j] += 1ll * mu[i] * i % M;
f[j] %= M;
}
}
for(int i = 2; i < N; i++) {
f[i] = f[i] * i % M;
f[i] += f[i - 1];
f[i] %= M;
}
}
inline ll sum(int n) {
return 1ll * (1 + n) * n / 2 % M;
}
int main() {
IOS;
init();
int t;
cin >> t;
while(t--) {
ll ans = 0;
int n, m;
cin >> n >> m;
int cur = 1;
while(cur <= min(n, m)) {
int p1 = n / (n / cur);
int p2 = m / (m / cur);
int nt = min(min(n, m), min(p1, p2));
ans += 1ll * sum(n / cur) * sum(m / cur) % M * (f[nt] - f[cur - 1]) % M;
ans %= M;
cur = nt + 1;
}
cout << (ans + M) % M << endl;
}
}