[莫比乌斯反演+数状数组] BZOJ3529: [Sdoi2014]数表

题意

有一张N×m的数表，其第i行第j列（1 <=i <=n，1 <=j <=m）的数值为

能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。

多次询问，输入的第一行一个整数Q表示测试点内的数据组数，接下来Q行，每行三个整数n，m，a(|a| < =10^9)描述一组数据。

n,m,Q <=10^5

题解

设g(i)为gcd(x,y)等于i的数对个数(x<=n,y<=m)，显然

g(i)=∑i|dμ(di)⌊nd⌋⌊md⌋

设F(i)表示i的约数和，则答案为

∑i=1F(i)g(i)=∑i=1F(i)∗∑i|dμ(di)⌊nd⌋⌊md⌋=∑d=1⌊nd⌋⌊md⌋∑i|dF(i)∗μ(di)

这个F(i) 可以用线性筛搞出来。

还有个条件就是F(i)<=a才对答案有贡献。

思路还是办法把∑i|dF(i)∗μ(di) 的前缀和搞出来。

怎么满足a的约束条件呢？离线搞就行了，把询问按a排序，把F(i)的值从小到大考察，然后用树状数组实现不断的插入F(i)的贡献，维护前缀和。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define Fir first
#define Sec second
using namespace std;
const int maxn=100005, N=100000;
int Q,mu[maxn],p[maxn],ans[maxn];
bool vis[maxn];
pair<int,int> F[maxn];
struct data{
int n,m,w,id;
bool operator < (const data &b)const{
return w<b.w;
}
} q[maxn];
void get_mu(){
memset(vis,1,sizeof(vis));
mu[1]=1;
for(int i=2;i<=N;i++){
if(vis[i]) p[++p[0]]=i, mu[i]=-1;
for(int j=1,k;j<=p[0]&&(k=p[j]*i)<=N;j++){
vis[k]=false;
if(i%p[j]==0){ mu[k]=0; break; }
mu[k]=-mu[i];
}
}
}
int t1[maxn],t2[maxn]; // t1=(1+p1^1+p1^2+p1^3+...p1^k1)  t2=p1^k1
void get_F(){
memset(vis,1,sizeof(vis)); p[0]=0;
F[1].Fir=1;
for(int i=2;i<=N;i++){
if(vis[i]) p[++p[0]]=i, F[i].Fir=1+i, t1[i]=1+i, t2[i]=i;
for(int j=1,k;j<=p[0]&&(k=p[j]*i)<=N;j++){
vis[k]=false;
if(i%p[j]==0){
t2[k]=t2[i]*p[j]; t1[k]=t1[i]+t2[k];
F[k].Fir=F[i].Fir/t1[i]*t1[k];
break;
}
F[k].Fir=F[p[j]].Fir*F[i].Fir;
t2[k]=p[j]; t1[k]=1+p[j];
}
}
for(int i=1;i<=N;i++) F[i].Sec=i;
sort(F+1,F+1+N);
}
int bit[maxn];
void Updata(int x,int val){
for(;x<=N;x+=(x&(-x))) bit[x]+=val;
}
int Query(int x){
int res=0;
for(;x;x-=(x&(-x))) res+=bit[x];
return res;
}
void Put(int i){
for(int j=1;F[i].Sec*j<=N;j++) Updata(F[i].Sec*j,F[i].Fir*mu[j]);
}
int main(){
freopen("bzoj3529.in","r",stdin);
freopen("bzoj3529.out","w",stdout);
get_mu(); get_F();
scanf("%d",&Q);
for(int i=1;i<=Q;i++) scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].w), q[i].id=i;
sort(q+1,q+1+Q);
for(int i=1,last=0;i<=Q;i++){
while(last<N&&F[last+1].Fir<=q[i].w) Put(++last);
int res=0, n=q[i].n, m=q[i].m;
if(n>m) swap(n,m);
for(int d=1,nxt;d<=n;d=nxt+1){
nxt=min(n/(n/d),m/(m/d));
res+=(Query(nxt)-Query(d-1))*(n/d)*(m/d);
}
ans[q[i].id]=res&0x7fffffff;
}
for(int i=1;i<=Q;i++) printf("%d\n",ans[i]);
return 0;
}