solution
先(orz)一波大佬，然后扔上一个公式(233)
[d(ij)=sumlimits_{x | i}sumlimits_{y|j}[gcd(x,y)==1]
]
是不是看到这个公式瞬间就有思路了。
我们开始推柿子。
[sumlimits_{i=1}^nsumlimits_{j=1}^m d(ij)\ =sumlimits_{i=1}^nsumlimits_{j=1}^msumlimits_{x|i}sumlimits_{y|j}[gcd(x,y)==1]\ =sumlimits_{x=1}^nsumlimits_{y=1}^mlfloorfrac{n}{x}
floorlfloorfrac{m}{y}
floorsumlimits_{d|gcd(x,y)}mu(d)\ =sumlimits_{d=1}^nmu(d)sumlimits_{x=1}^{lfloorfrac{n}{d}
floor}sumlimits_{y=1}^{lfloorfrac{m}{d}
floor}lfloorfrac{n}{dx}
floorlfloorfrac{m}{dy}
floor\=sumlimits_{d=1}^nmu(d)sumlimits_{x=1}^{lfloorfrac{n}{d}
floor}lfloorfrac{n}{dx}
floorsumlimits_{y=1}^{lfloorfrac{m}{d}
floor}lfloorfrac{m}{dy}
floor
]
令(f(n)=sumlimits_{i=1}^nlfloorfrac{n}{i} floor)
我们可以利用整除分块在(O(nsqrt{n}))的复杂度内预处理所有的(f)
所以原式=
[sumlimits_{d=1}^nmu(d)f(lfloorfrac{n}{d}
floor)f(lfloorfrac{m}{d}
floor)
]
因为(f)都已经预处理出来，所以只要整除分块，就可以在(O(sqrt{n}))的复杂度内计算每次询问了。
code
/*
* @Author: wxyww
* @Date: 2020-04-23 19:36:37
* @Last Modified time: 2020-04-23 20:04:55
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
#include<ctime>
using namespace std;
typedef long long ll;
const int N = 50010;
ll read() {
ll x = 0,f = 1;char c = getchar();
while(c < '0' || c > '9') {
if(c == '-') f = -1; c = getchar();
}
while(c >= '0' && c <= '9') {
x = x * 10 + c - '0'; c = getchar();
}
return x * f;
}
int prime[N],tot,mu[N],vis[N];
void pre() {
mu[1] = 1;
for(int i = 2;i <= 50000;++i) {
if(!vis[i]) {
prime[++tot] = i;
mu[i] = -1;
}
for(int j = 1;j <= tot && prime[j] * i <= 50000;++j) {
vis[prime[j] * i] = 1;
if(i % prime[j])
mu[i * prime[j]] = -mu[i];
else break;
}
}
for(int i = 1;i <= 50000;++i) mu[i] += mu[i - 1];
}
ll f[N];
int main() {
pre();
for(int i = 1;i <= 50000;++i) {
for(int l = 1,r;l <= i;l = r + 1) {
r = i / (i / l);
f[i] += (r - l + 1) * (i / l);
}
}
int T = read();
while(T--) {
int n = read(),m = read();
if(n > m) swap(n,m);
ll ans = 0;
for(int l = 1,r;l <= n;l = r + 1) {
r = min(n / (n / l),m / (m / l));
ans += 1ll * (mu[r] - mu[l - 1]) * f[n / l] * f[m / l];
}
printf("%lld
",ans);
}
return 0;
}