问题分析
首先有个很厉害的结论：
[d(ij)=sigma_0(ij)=sum_{x|i}sum_{y|j}[gcd(x,y)=1]
]
Prove:
[sigma_0(ij)=sumlimits_{d|ij}1=sumlimits_{d_1mid i,d_2mid j}1=sumlimits_{d_1mid i}sumlimits_{d_2mid j}[gcd(frac {i} {d_1},d_2)=1]
]
解释：
我们拆分(d=d_1d_2)，如果直接拆为(sumlimits_{d_1mid i}sumlimits_{d_2|j}1)，那么显然会记重。那么我们令(gcd(frac{i}{d_1},d_2)=1)，也就是说让(d_1)塞的因数尽可能多，那么就不会记重了。而等式右边与原式右侧等价，即证。
然后就有一下操作：
[Ans=sum_{i=1}^Nsum_{j=1}^Msum_{x|i}sum_{y|j}[gcd(x,y)=1]=sum_{x=1}^Nsum_{y=1}^M[gcd(x,y)=1]lfloorfrac{N}{x}
floorlfloorfrac{M}{y}
floor
]
然后一波莫比乌斯反演：
[f(t)=sum_{x=1}^Nsum_{y=1}^M[gcd(x,y)=1]lfloorfrac{N}{x}
floorlfloorfrac{M}{y}
floor\
F(t)=sum_{x=1}^Nsum_{y=1}^M[t|gcd(x,y)]lfloorfrac{N}{x}
floorlfloorfrac{M}{y}
floor=sum_{x=1}^{lfloorfrac{N}{t}
floor}sum_{y=1}^{lfloorfrac{M}{t}
floor}lfloorfrac{M}{tx}
floorlfloorfrac{N}{ty}
floor=sum_{t|d}f(d)\
f(t)=sum_{t|d}mu(frac{d}{t})F(d)
]
所以
[Ans=f(1)=sum_{d}^{min(N,M)}mu(d)(sum_{x=1}^{lfloorfrac{N}{d}
floor}lfloorfrac{N}{xd}
floor)(sum_{y=1}^{lfloorfrac{M}{d}
floor}lfloorfrac{M}{yd}
floor)
]
然后我们预处理(sumlimits_{x=1}limits^{t}lfloorfrac{t}{d} floor)，然后整除分块正常操作就好了。
参考代码
#include <bits/stdc++.h>
using namespace std;
const int Maxn = 50010;
void EulerSieve();
void Cal();
int main() {
EulerSieve();
int TestCases;
scanf( "%d", &TestCases );
for( ; TestCases--; )
Cal();
return 0;
}
int Vis[ Maxn ], Prime[ Maxn ], Num, Mu[ Maxn ], Sum[ Maxn ];
long long F[ Maxn ];
void EulerSieve() {
Mu[ 1 ] = 1;
for( int i = 2; i < Maxn; ++i ) {
if( !Vis[ i ] ) {
Prime[ ++Num ] = i;
Mu[ i ] = -1;
}
for( int j = 1; j <= Num && 1ll * i * Prime[ j ] < 1ll * Maxn; ++j ) {
Vis[ i * Prime[ j ] ] = 1;
if( i % Prime[ j ] == 0 ) break;
Mu[ i * Prime[ j ] ] = -Mu[ i ];
}
}
for( int i = 1; i < Maxn; ++i ) Sum[ i ] = Sum[ i - 1 ] + Mu[ i ];
for( int i = 1; i < Maxn; ++i )
for( int l = 1, r; l <= i; l = r + 1 ) {
r = i / ( i / l );
F[ i ] += 1ll * ( r - l + 1 ) * ( i / l );
}
return;
}
void Cal() {
int N, M;
scanf( "%d%d", &N, &M );
long long Ans = 0;
for( int i = 1, j; i <= N && i <= M; i = j + 1 ) {
j = min( N / ( N / i ), M / ( M / i ) );
Ans += 1ll * ( Sum[ j ] - Sum[ i - 1 ] ) * F[ N / i ] * F[ M / i ];
}
printf( "%lld
", Ans );
return;
}