题目链接：数表
我们一起来膜PoPoQQQ大爷的题解吧Orz
首先我们来考虑没有(a)的限制该怎么做。显然交换(n)，(m)答案不变，所以后面默认(n le m)。
我们定义两个函数：
[f(x)=sum_{d|x}d]
[g(x)=sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=x]]
那么显然有：[ans=sum_{i=1}^nf(i)g(i)]
(g)函数我们可以考虑化简一下：
egin{aligned}
g(x)&=sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=x] \
&=sum_{i=1}^{lfloor frac{n}{x}
floor}sum_{j=1}^{lfloor frac{m}{x}
floor}sum_{d|i,d|j}mu(d) \
&=sum_{d=1}^{lfloor frac{n}{x}
floor}mu(d)lfloor frac{n}{dx}
floorlfloor frac{m}{dx}
floor \
&=sum_{x|d}mu(frac{d}{x})lfloor frac{n}{d}
floorlfloor frac{m}{d}
floor
end{aligned}
于是可以得到：
egin{aligned}
ans=&sum_{i=1}^nf(i)sum_{i|d}mu(frac{d}{i})lfloor frac{n}{d}
floorlfloor frac{m}{d}
floor \
=&sum_{d=1}^nlfloor frac{n}{d}
floorlfloor frac{m}{d}
floorsum_{i|d}f(i)mu(frac{d}{i})
end{aligned}
我们令(h(x)=sum_{i|x}f(i)mu(frac{x}{i}))，那么我们只要有(h(x))的前缀和，就可以在(O(sqrt{n}))的时间内求出(ans)了。
然后我们来考虑(a)的限制。显然，只有(f(x)le a)的(f(x))才会对答案有影响。所以，我们把所有的询问按(a)从小到大排好序，并且把所有的(f(x))从小到大依次加进来，使用树状数组维护前缀和，依次处理每个询问，就可以在(O(Qsqrt{n}log{n}+nlog^2 n))的时间内解决所有询问了。
下面贴代码：
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
#define maxn 100010
using namespace std;
typedef long long llg;
struct data{
int n,m,a,b;
bool operator < (const data &h)const{return a<h.a;}
}s[maxn];
int n,m,nm,a[maxn],T,c[maxn],ans[maxn];
int mu[maxn],f[maxn],pr[maxn],lp;
const int mod=2147483647;
bool vis[maxn];
int getint(){
int w=0;bool q=0;
char c=getchar();
while((c>'9'||c<'0')&&c!='-') c=getchar();
if(c=='-') c=getchar(),q=1;
while(c>='0'&&c<='9') w=w*10+c-'0',c=getchar();
return q?-w:w;
}
bool cmp(int x,int y){return f[x]<f[y];}
void add(int x,int y){if(y) while(x<=nm) c[x]+=y,x+=x&(-x);}
int sum(int x){
int t=0;
while(x) t+=c[x],x-=x&(-x);
return t;
}
int main(){
File("a");
mu[1]=f[1]=a[1]=1; T=getint();
for(int i=1;s[i].b=i,i<=T;i++){
s[i].n=getint(),s[i].m=getint();
if(s[i].n>s[i].m) swap(s[i].n,s[i].m);
s[i].a=getint(); nm=max(nm,s[i].n);
}
for(int i=2;a[i]=i,i<=nm;i++){
if(!vis[i]) pr[++lp]=i,mu[i]=-1,f[i]=i+1;
for(int j=1;pr[j]*i<=nm;j++){
vis[pr[j]*i]=1;
if(i%pr[j]) mu[pr[j]*i]=-mu[i],f[pr[j]*i]=f[i]*f[pr[j]];
else{ f[pr[j]*i]=f[i]+(f[i]-f[i/pr[j]])*pr[j]; break;}
}
}
sort(s+1,s+T+1); sort(a+1,a+nm+1,cmp);
for(int i=1,now=1,b,la=0;i<=T;i++){
n=s[i].n; m=s[i].m; b=s[i].b;
while(now<=nm && f[a[now]]<=s[i].a){
for(int j=a[now];j<=nm;j+=a[now])
add(j,f[a[now]]*mu[j/a[now]]);
now++;
}
for(int j=1,nt,na;j<=n;j=nt+1){
nt=min(n/(n/j),m/(m/j)); na=sum(nt);
ans[b]+=(na-la)*(n/j)*(m/j); la=na;
}
}
for(int i=1;i<=T;i++) printf("%d
",ans[i]&mod);
return 0;
}