[BZOJ3529] [Sdoi2014]数表
Description
有一张 n×m 的数表，其第 i 行第 j 列（1 <= i <= n, 1 <= j <= m）的数值为能同时整除 i 和 j 的所有自然数之和。给定 a , 计算数表中不大于 a 的数之和。
Input
输入包含多组数据。输入的第一行一个整数Q表示测试点内的数据组数接下来Q行，每行三个整数n，m，a(|a| < =10^9)描述一组数据。1 < =N．m < =10^5 ， 1 < =Q < =2×10^4
Output
对每组数据，输出一行一个整数，表示答案模2^31的值。
Sample Input
2
4 4 3
10 10 5
Sample Output
20
148
试题分析
说人话就是求：$$sum_{i=1}nsum_{j=1}m sum_{d|i,d|j}[dleq a] d$$
条件很烦人，先把它略去好了，转而变成：$$sum_{i=1}nsum_{j=1}m f(gcd(i,j))$$
其中(f(d))为d的约数和。
枚举(gcd(i,j))得到：$$sum_{d=1}^{min(n,m)} f(d) sum_{i=1}^{lfloor frac{n}{d}
floor} sum_{j=1}^{lfloor frac{m}{d}
floor} [gcd(i,j)=1]$$。
后面的已经是套路，就直接写出：
[sum_{d=1}^{min(n,m)} f(d) sum_{t=1}^{min(lfloor frac{n}{d}
floor,lfloor frac{m}{d}
floor)} mu(t) lfloor frac{n}{dt}
floor lfloor frac{m}{dt}
floor$$。
然后令$T=dt$，枚举$T$得到：$${sum_{T=1}^{min(n,m)}lfloorfrac{n}{T}
floorlfloorfrac{m}{T}
floorsum_{d|T}f(d)μ(frac{T}{d})}$$。
至此，基本上大功告成，前面可以用数论分块在$sqrt{n}$时间内做完，如果加上a的限制直接将$f$排序然后树状数组加入维护前缀和即可。
```c++
#include<iostream>
#include<cstring>
#include<cstdio>
#include<vector>
#include<algorithm>
using namespace std;
#define LL long long
inline int read(){
int x=0,f=1; char c=getchar();
for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
for(;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
const int MAXN = 100010;
const int INF = 2147483647;
int mu[MAXN+1],g[MAXN+1],c[MAXN+1];
int pri[MAXN+1]; bool vis[MAXN+1]; int N,cnt;
struct data{int d,s;}f[MAXN+1];
struct query{int n,m,a,id;}q[MAXN+1];
int t[MAXN+1];
bool cmp(data a,data b){return a.s<b.s;}
bool cmp2(query a,query b){return a.a<b.a;}
inline int lowbit(int x){return x&(-x);}
inline int Query(int x){
int res=0; for(;x;x-=lowbit(x)) res+=c[x]; return res;
}
inline void Add(int x,int k){
for(;x<=MAXN;x+=lowbit(x)) c[x]+=k; return ;
}
inline int Pow(int a,int b){
int res=1;
while(b){
if(b&1) res=res*a;
b>>=1; a=a*a;
} return res;
}
inline void init(){
mu[1]=f[1].d=f[1].s=1;
for(int i=2;i<MAXN;i++){
f[i].d=i;
if(!vis[i]){mu[i]=-1; f[i].s=t[i]=i+1; g[i]=1; pri[++cnt]=i;}
for(int j=1;j<=cnt&&pri[j]*i<MAXN;j++){
vis[i*pri[j]]=true;
if(i%pri[j]==0){
mu[i*pri[j]]=0;
g[i*pri[j]]=g[i]+1; t[i*pri[j]]=t[i]+Pow(pri[j],g[i]+1);
f[i*pri[j]].s=f[i].s/t[i]*t[i*pri[j]];
break;
}
mu[i*pri[j]]=-mu[i]; f[i*pri[j]].s=f[i].s*f[pri[j]].s;
g[i*pri[j]]=1; t[i*pri[j]]=pri[j]+1;
}
}
sort(f+1,f+MAXN,cmp);
return ;
}
int ans[MAXN+1];
int main(){
//freopen(".in","r",stdin);
//freopen(".out","w",stdout);
N=read(); init();
for(int i=1;i<=N;i++) q[i].n=read(),q[i].m=read(),q[i].a=read(),q[i].id=i;
sort(q+1,q+N+1,cmp2); int j=1;
for(int i=1;i<=N;i++){
while(f[j].s<=q[i].a&&j<MAXN){
for(int t=1;t*f[j].d<MAXN;t++) Add(f[j].d*t,f[j].s*mu[t]); ++j;
} 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]+=(Query(r)-Query(l-1))*(n/l)*(m/l)%INF;
} ans[q[i].id]&=INF;
}
for(int i=1;i<=N;i++) printf("%d
",ans[i]);
return 0;
}
```]