在一张(n imes m)数表中，第i行第j列的数字为即能整除i又能整除j的自然数之和，Q组询问，询问这张数表中不大于a的数之和，1 < =n．m < =10^5 ， 1 < =Q < =2*10^4。
解
转换为约数计数问题，不难有
[ans=sum_{i=1}^nsum_{j=1}^msigma(gcd(i,j))
]
因为限制数表中的数不大于a，于是我们得设的式子中要有(sigma(gcd(i,j)))，后面推式也要保留其，于是设(nleq m),有
[ans=sum_{d=1}^nsigma(d)sum_{i=1}^nsum_{j=1}^m(gcd(i,j)==d)
]
设
[f(d)=sum_{i=1}^nsum_{j=1}^m(gcd(i,j)==d)
]
[F(d)=[n/d][m/d]
]
由Mobius反演定理有
[f(d)=sum_{d|x}[n/x][m/x]mu(x/d)
]
代入有
[ans=sum_{d=1}^nsigma(d)sum_{d|x}[n/x][m/x]mu(x/d)=
]
[sum_{x=1}^n[n/x][m/x]sum_{d|x}sigma(d)mu(x/d)
]
注意到(sigma)为积性函数，我们可以线性递推，而后式显然可以维护，前式又可以整除分块，于是如果没有a的限制，我们已经可以做了，
而此时能贡献的只有(sigma(d)leq a)，维护的东西要求小于某个数，而询问又可以离线，于是考虑给a排序，给(sigma)排序，按照a的增加，慢慢地加入我们需要维护的后式的前缀和，所以我们需要支持快速修改，而树状数组就是一个不错的选择，所以不难得知时间复杂度应为(O(nlog(n)+Qsqrt{n}))。
参考代码：
#include <iostream>
#include <cstdio>
#include <algorithm>
#define il inline
#define ri register
#define ll long long
#define sxr 100000
#define yyb 2147483647
#define swap(x,y) x^=y^=x^=y
using namespace std;
struct query{
int n,m,a,id;
il bool operator<(const query&x){
return a<x.a;
}
il void sort(){
if(n>m)swap(n,m);
}
}q[20001];
template<class f1,class f2,f2 key>
struct lowbit{
f1 s[key+1];
il void add(int p,int x){
while(p<=key)
s[p]+=x,p+=(-p)&p;
}
il f1 sum(int p){
f1 ans(0);
while(p)ans+=s[p],p-=(-p)&p;
return ans;
}
};
struct ysh{
int id,dat;
il bool operator<(ysh&x){
return dat<x.dat;
}
}ds[sxr+1];
bool check[sxr+1];
int prime[sxr+1],pt,mb[sxr+1],
df[sxr+1],ans[20001];
lowbit<int,int,sxr+1>czf;
il ll ask(int);
il int min(int,int);
il void prepare(int),read(int&);
int main(){
int lsy,i,j,k,l;prepare(sxr),read(lsy);
for(i=1;i<=lsy;++i)
read(q[i].n),read(q[i].m),read(q[i].a),
q[i].id=i,q[i].sort();
sort(q+1,q+lsy+1);
for(i=1,j=1;i<=lsy;++i){
while(ds[j].dat<=q[i].a&&j<=sxr){
for(k=ds[j].id;k<=sxr;k+=ds[j].id)
czf.add(k,ds[j].dat*mb[k/ds[j].id]);
++j;
}
for(k=1;k<=q[i].n;k=l+1)
l=min(q[i].n/(q[i].n/k),q[i].m/(q[i].m/k)),
ans[q[i].id]+=(q[i].n/k)*(q[i].m/k)*(czf.sum(l)-czf.sum(k-1));
}for(i=1;i<=lsy;++i)printf("%d
",ans[i]&yyb);
return 0;
}
il void read(int &x){
x&=0;ri char c;while(c=getchar(),c<'0'||c>'9');
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
}
il int min(int a,int b){
return a<b?a:b;
}
il void prepare(int n){
int i,j;mb[1]|=ds[1].id|=ds[1].dat|=df[1]|=true;
for(i=2;i<=n;++i){
if(!check[i])mb[i]=-1,ds[i].dat=i+1,
df[i]=i+1,prime[++pt]=i;
for(j=1;j<=pt&&prime[j]<=n/i;++j){
check[i*prime[j]]|=true;
if(!(i%prime[j])){
df[i*prime[j]]=df[i]*prime[j]+1;
ds[i*prime[j]].dat=ds[i].dat/df[i]*df[i*prime[j]];
break;
}mb[i*prime[j]]=~mb[i]+1,
ds[i*prime[j]].dat=ds[i].dat*ds[prime[j]].dat,
df[i*prime[j]]=df[prime[j]];
}ds[i].id=i;
}sort(ds+1,ds+n+1);
}