题目链接
式中(sigma(d)=sumlimits_{d'|d}{d'})，这是个积性函数，可以欧拉筛(O(n))求出来。
先不考虑(a)的限制，得
(ans=sumlimits_{i=1}^n{sumlimits_{j=1}^m{sigma({gcd(i,j))}}})
(=sumlimits_{d=1}{ sigma(d) sumlimits_{i=1}^{n}{sumlimits_{j=1}^{m}{[gcd(i,j)=d]}}})
设
(f(d)=sumlimits_{i=1}^{n}{sumlimits_{j=1}^{m}{[gcd(i,j)=d]}})
(F(d)=sumlimits_{i=1}^{n}{sumlimits_{j=1}^{m}{[d mid gcd(i,j)]}})
(F(d)=sumlimits_{d'}{f(d')}=frac{n}{d} cdot frac{m}{d})
(f(d)=sumlimits_{d mid d'}{mu(frac{d'}{d})F(d')})
( herefore ans=sumlimits_{d=1}{sigma(d)sumlimits_{k=1}{mu(k)frac{n}{kd}frac{m}{kd}}})
至此我们可以(O(n))计算出(ans)，对于(a)的限制，我们可以离线，将(sigma(d) leq a)的添加入树状数组中，分块时查询区间和。
然后。。。就(T)了。
所以我们要通过某种方式预处理，优化时间复杂度。
因为(mn=10^{10})，预处理不现实，我们对(sigma)，(mu)预处理。
(m)，(n)很烦人，把他们提到前面。
(ans=sumlimits_{T=1}{frac{n}{T}frac{m}{T}sumlimits_{dmid T}{mu(frac{T}{d})sigma(d)}})
用树状树状储存(sumlimits_{dmid T}{mu(frac{T}{d})sigma(d)})这样我们每次加入(sigma(d))式就枚举(d)的倍数，加入树状数组中。
大家应该知道调和级数式(O(logn))级别的吧。。。所以复杂度为(O(nlog^2 n+qsqrt n log n))。
(frak{code:})
#include<bits/stdc++.h>
#define IL inline
#define LL unsigned int
using namespace std;
const int N=1e5+3;
const LL p=1u<<31;
struct que{
int n,m,a,id;
bool operator<(const que &b) const{
return a<b.a;}
}q[N];
struct hh{
LL val;int pos;
bool operator<(const hh &a) const{
return val<a.val;}
}b[N];
int n,m,num,t,now,mu[N],vis[N],pri[N],sum[N],Max[N];
LL s[N],ans[N],pre[N];
IL int in(){
char c;int f=1;
while((c=getchar())<'0'||c>'9')
if(c=='-') f=-1;
int x=c-'0';
while((c=getchar())>='0'&&c<='9')
x=x*10+c-'0';
return x*f;
}
struct BIT{
LL c[N];
IL int lb(int x){return x&-x;}
IL void add(int y,LL x){for(;y<N;y+=lb(y)) c[y]+=x;}
IL LL ask(int y){LL res=0;for(;y;y-=lb(y)) res+=c[y];return res;}
}T;
IL void init(int n){
mu[1]=s[1]=1;
for(int i=2;i<=n;++i){
if(!vis[i]){vis[i]=i,mu[i]=-1,pri[++pri[0]]=i,sum[i]=s[i]=i+1,Max[i]=i;}
for(int j=1,k;(k=i*pri[j])<=n;++j){
vis[k]=pri[j];
if(i%pri[j]) sum[k]=pri[j]+1,s[k]=s[i]*sum[k],mu[k]=-mu[i],Max[k]=pri[j];
else{
Max[k]=Max[i]*pri[j],sum[k]=sum[i]+Max[k],
s[k]=s[i]/sum[i]*sum[k];break;
}
}
}
for(int i=1;i<=n;++i) b[i]=(hh){s[i],i};
sort(b+1,b+n+1);
}
IL LL solve(LL n,LL m){
LL res=0;
if(n>m) swap(n,m);
for(LL l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
res+=(n/l)*(m/l)*(T.ask(r)-T.ask(l-1));
}
return res;
}
int main()
{
init(1e5),t=in();
for(int i=1;i<=t;++i) q[i]=(que){in(),in(),in(),i};
sort(q+1,q+t+1);
for(int i=1;i<=t;++i){
while(b[now+1].val<=q[i].a&&now<1e5)
for(int j=b[++now].pos,k=1;j<=1e5;j+=b[now].pos,++k) T.add(j,b[now].val*mu[k]);
ans[q[i].id]=solve(q[i].n,q[i].m)%p;
}
for(int i=1;i<=t;++i) printf("%u
",ans[i]);
return 0;
}