【BZOJ3529】[Sdoi2014]数表 莫比乌斯反演+树状数组

【BZOJ3529】[Sdoi2014]数表

Description

    有一张N×m的数表，其第i行第j列（1 < =i < =礼，1 < =j < =m）的数值为能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。

Input

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

Output

    对每组数据，输出一行一个整数，表示答案模2^31的值。

Sample Input

2
4 4 3
10 10 5

Sample Output

20
148

HINT

1 < =N．m < =10^5  ， 1 < =Q < =2×10^4

题解：首先我们要知道一个数的约数和是nloglogn级别的，所以先不考虑a的限制，我们还是采用熟悉的莫比乌斯反演。

$ans=sumlimits_{i=1}^nsumlimits_{j=1}^nf(gcd(i,j))\=sumlimits_{d=1}^nf(d)sumlimits_{i=1}^{lfloor frac n d floor}sumlimits_{j=1}^{lfloor frac m d floor } [gcd(i,j)==1]\=sumlimits_{d=1}^nf(d)sumlimits_{e=1}^{lfloor frac n d floor}mu(e) lfloor frac n {de} floor lfloor frac m {de} floor\=sumlimits_{D=1}^{n}sumlimits_{d|D}f(d)mu(frac D d) lfloor frac n {D} floor lfloor frac m {D} floor$

那么如果考虑a的限制呢？我们将所有询问离线，按a排序，然后从小到大处理所有询问，每处理到一个询问，就将所有f(d)<=a的d的贡献都统计出来。也就是说我们需要用某个数据结构来维护f*mu的前缀和，树状数组即可。

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
const int N=100000;
typedef unsigned int ui;
int pri[N/10];
ui mu[N+10],f[N+10],s[N+10],ans[N];
int q[N+10];
int num,Q;
bool np[N+10];
struct node
{
	int n,m,a,org;
}p[N];
inline int rd()
{
	int ret=0,f=1;	char gc=getchar();
	while(gc<'0'||gc>'9')	{if(gc=='-')	f=-f;	gc=getchar();}
	while(gc>='0'&&gc<='9')	ret=ret*10+(gc^'0'),gc=getchar();
	return ret*f;
}
bool cmp(const node &a,const node &b)
{
	return a.a<b.a;
}
bool cmpq(const int &a,const int &b)
{
	return f[a]<f[b];
}
inline void updata(int x,ui val)
{
	for(int i=x;i<=N;i+=i&-i)	s[i]+=val;
}
inline ui query(int x)
{
	ui ret=0;
	for(int i=x;i;i-=i&-i)	ret+=s[i];
	return ret;
}
int main()
{
	int i,j,k,last;
	mu[1]=1;
	for(i=2;i<=N;i++)
	{
		if(!np[i])	pri[++num]=i,mu[i]=-1;
		for(j=1;j<=num&&i*pri[j]<=N;j++)
		{
			np[i*pri[j]]=1;
			if(i%pri[j]==0)	break;
			mu[i*pri[j]]=-mu[i];
		}
	}
	for(i=1;i<=N;i++)	for(j=i;j<=N;j+=i)	f[j]+=i;
	for(i=1;i<=N;i++)	q[i]=i;
	Q=rd();
	for(i=1;i<=Q;i++)
	{
		p[i].n=rd(),p[i].m=rd(),p[i].org=i;
		if(p[i].n>p[i].m)	swap(p[i].n,p[i].m);
		p[i].a=max(rd(),0);
	}
	sort(q+1,q+N+1,cmpq);
	sort(p+1,p+Q+1,cmp);
	for(i=j=1;i<=Q;i++)
	{
		for(;f[q[j]]<=p[i].a;j++)	for(k=q[j];k<=N;k+=q[j])	updata(k,f[q[j]]*mu[k/q[j]]);
		for(k=1;k<=p[i].n;k=last+1)
		{
			last=min(p[i].n/(p[i].n/k),p[i].m/(p[i].m/k));
			ans[p[i].org]+=(p[i].n/k)*(p[i].m/k)*(query(last)-query(k-1));
		}
	}
	for(i=1;i<=Q;i++)	printf("%u
",ans[i]&0x7fffffff);
	return 0;
}