BZOJ 3529 【SDOI2014】 数表

题目链接：数表

　　我们一起来膜PoPoQQQ大爷的题解吧Orz

　　首先我们来考虑没有(a)的限制该怎么做。显然交换(n)，(m)答案不变，所以后面默认(n le m)。

　　我们定义两个函数：

[f(x)=sum_{d|x}d]

[g(x)=sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=x]]

　　那么显然有：[ans=sum_{i=1}^nf(i)g(i)]

　　(g)函数我们可以考虑化简一下：

egin{aligned}
g(x)&=sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=x] \
&=sum_{i=1}^{lfloor frac{n}{x} floor}sum_{j=1}^{lfloor frac{m}{x} floor}sum_{d|i,d|j}mu(d) \
&=sum_{d=1}^{lfloor frac{n}{x} floor}mu(d)lfloor frac{n}{dx} floorlfloor frac{m}{dx} floor \
&=sum_{x|d}mu(frac{d}{x})lfloor frac{n}{d} floorlfloor frac{m}{d} floor
end{aligned}

　　于是可以得到：

egin{aligned}
ans=&sum_{i=1}^nf(i)sum_{i|d}mu(frac{d}{i})lfloor frac{n}{d} floorlfloor frac{m}{d} floor \
=&sum_{d=1}^nlfloor frac{n}{d} floorlfloor frac{m}{d} floorsum_{i|d}f(i)mu(frac{d}{i})
end{aligned}

　　我们令(h(x)=sum_{i|x}f(i)mu(frac{x}{i}))，那么我们只要有(h(x))的前缀和，就可以在(O(sqrt{n}))的时间内求出(ans)了。

　　然后我们来考虑(a)的限制。显然，只有(f(x)le a)的(f(x))才会对答案有影响。所以，我们把所有的询问按(a)从小到大排好序，并且把所有的(f(x))从小到大依次加进来，使用树状数组维护前缀和，依次处理每个询问，就可以在(O(Qsqrt{n}log{n}+nlog^2 n))的时间内解决所有询问了。

　　下面贴代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
#define maxn 100010

using namespace std;
typedef long long llg;

struct data{
	int n,m,a,b;
	bool operator < (const data &h)const{return a<h.a;}
}s[maxn];
int n,m,nm,a[maxn],T,c[maxn],ans[maxn];
int mu[maxn],f[maxn],pr[maxn],lp;
const int mod=2147483647;
bool vis[maxn];

int getint(){
	int w=0;bool q=0;
	char c=getchar();
	while((c>'9'||c<'0')&&c!='-') c=getchar();
	if(c=='-') c=getchar(),q=1;
	while(c>='0'&&c<='9') w=w*10+c-'0',c=getchar();
	return q?-w:w;
}

bool cmp(int x,int y){return f[x]<f[y];}
void add(int x,int y){if(y) while(x<=nm) c[x]+=y,x+=x&(-x);}
int sum(int x){
	int t=0;
	while(x) t+=c[x],x-=x&(-x);
	return t;
}

int main(){
	File("a");
	mu[1]=f[1]=a[1]=1; T=getint();
	for(int i=1;s[i].b=i,i<=T;i++){
		s[i].n=getint(),s[i].m=getint();
		if(s[i].n>s[i].m) swap(s[i].n,s[i].m);
		s[i].a=getint(); nm=max(nm,s[i].n);
	}
	for(int i=2;a[i]=i,i<=nm;i++){
		if(!vis[i]) pr[++lp]=i,mu[i]=-1,f[i]=i+1;
		for(int j=1;pr[j]*i<=nm;j++){
			vis[pr[j]*i]=1;
			if(i%pr[j]) mu[pr[j]*i]=-mu[i],f[pr[j]*i]=f[i]*f[pr[j]];
			else{ f[pr[j]*i]=f[i]+(f[i]-f[i/pr[j]])*pr[j]; break;}
		}
	}
	sort(s+1,s+T+1); sort(a+1,a+nm+1,cmp);
	for(int i=1,now=1,b,la=0;i<=T;i++){
		n=s[i].n; m=s[i].m; b=s[i].b;
		while(now<=nm && f[a[now]]<=s[i].a){
			for(int j=a[now];j<=nm;j+=a[now])
				add(j,f[a[now]]*mu[j/a[now]]);
			now++;
		}
		for(int j=1,nt,na;j<=n;j=nt+1){
			nt=min(n/(n/j),m/(m/j)); na=sum(nt);
			ans[b]+=(na-la)*(n/j)*(m/j); la=na;
		}
	}
	for(int i=1;i<=T;i++) printf("%d
",ans[i]&mod);
	return 0;
}