题面
100
[Ans=sum_{i=1}^nsum_{j=1}^mg(gcd(i,j))
]
其中，
[g(d)=sum_{i|d}i
]
我们注意到(gcd(i,j))最多有(O(n))种取值，所以我们枚举(d=gcd(i,j))；
就有，
[Ans=sum_{d=1}^ng(d)*f(d)
]
其中，(f(d))表示，有多少对((i,j))的最大公约数为(d)，可以使用莫比乌斯反演求出。
那么，
[Ans=sum_{d=1}^ng(d)*sum_{d|i}lfloorfrac{n}{i}
floorlfloorfrac{m}{i}
floormu(frac{i}{d})
]
更换主体，
[Ans=sum_{i=1}^nlfloorfrac{n}{i}
floorlfloorfrac{m}{i}
floorsum_{d|i}g(d)*mu(frac{i}{d})
]
令(h(d)=sum_{d|i}g(d)*mu(frac{i}{d}))，
对于单个询问而言，如果我们知道(h)的前缀和，那么就可以(O(sqrt n))求出答案。
如果(a)一定时，我们可以预处理出(h)。
但实际上(a)不定。
突破口则是本题允许离线；
我们可以对(a)偏序，顺便利用树状数组(c)维护(h)，具体是：
对于一个(g(d)<=a)的，我们给所有(c[i*d])加上一个(g(d)*mu(i))，然后(c)就是当前的(h)，
其总的复杂度是(O(sum_{i=1}^nfrac{n}{i}*logn)approx O(nlog^2n))。
加上询问，那么最终的复杂度是(O(Qsqrt nlogn+nlog^2n))
Code
#include<bits/stdc++.h>
#define ll long long
#define fo(i,x,y) for(int i=x;i<=y;i++)
#define fd(i,x,y) for(int i=x;i>=y;i--)
using namespace std;
const char* fin="ex3623.in";
const char* fout="ex3623.out";
const int inf=0x7fffffff;
const int maxn=100007;
const ll mo=1ll<<31;
int t,mu[maxn],wu[maxn],p[maxn],ans[maxn];
int c[maxn],g[maxn];
bool bz[maxn];
void modify(int v,int u){for(;v<maxn;v+=v&-v)c[v]=(0ll+u+c[v])%mo;}
int getsum(int v){
ll k=0;
for(;v;v-=v&-v) k+=c[v];
return k%mo;
}
int getsum(int l,int r){return (getsum(r)-getsum(l-1)+mo)%mo;}
struct query{
int n,m,a,id;
}q[maxn];
bool cmp(query a,query b){return a.a<b.a;}
bool cmp1(int a,int b){return wu[a]<wu[b];}
int main(){
freopen(fin,"r",stdin);
freopen(fout,"w",stdout);
scanf("%d",&t);
fo(i,1,t) scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a),q[i].id=i;
mu[1]=1;
wu[1]=1;
fo(i,2,maxn-1){
if (!bz[i]){
p[++p[0]]=i;
wu[i]=i+1;
mu[i]=-1;
}
fo(j,1,p[0]){
int k=p[j]*i;
if (k>=maxn) break;
bz[k]=true;
if (i%p[j]==0){
mu[k]=0;
wu[k]=wu[i]+p[j]*(wu[i]-wu[i/p[j]]);
break;
}else{
mu[k]=-mu[i];
wu[k]=p[j]*wu[i]+wu[i];
}
}
}
fo(i,1,maxn-1) g[i]=i;
sort(q+1,q+t+1,cmp);
sort(g+1,g+maxn,cmp1);
int l=1;
fo(i,1,t){
while (l<maxn && wu[g[l]]<=q[i].a){
for(int k=g[l];k<maxn;k+=g[l]) modify(k,wu[g[l]]*mu[k/g[l]]);
l++;
}
int n=q[i].n,m=q[i].m;
if (n>m) swap(n,m);
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans[q[i].id]=(ans[q[i].id]+1ll*getsum(l,r)*(n/l)*(m/l))%mo;
}
}
fo(i,1,t) printf("%d
",ans[i]);
return 0;
}