BZOJ 3529: [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\times 10^4$
Output
对每组数据，输出一行一个整数，表示答案模2^31的值。
Sample Input
2
4 4 3
10 10 5
Sample Output
20
148

这题出得很刁钻……

预备知识（莫比乌斯反演）
思路：
莫比乌斯反演+树状数组维护。
1.推导公式。
2.理解树状数组的用处。
3.线性筛中的细节。
4.小优化。

1.推导公式

首先，我们设两个函数：
$F(a)=\lfloor n/a\rfloor$*$\lfloor m/a\rfloor$
表示$a|gcd(i,j)$的方案数.
$f(a)$表示$a=gcd(i,j)$的方案数。
显然，$F(a)={\sum }_{a|d}f(d)$。

设n<=m，那么，我们忽略a的存在，总方案数就是
${\sum }_{d=1}^n{\sum }_{i=1}^n{\sum }_{j=1}^{n}f[d]\ast [gcd(i,j)==d](这是bool表达式，真为1，假为0）$=
${\sum }_{d=1}^n{\sum }_{i=1}^{\lfloor n/d\rfloor }{\sum }_{j=1}^{\lfloor m/d\rfloor }f[d]\ast [gcd(i,j)==1]$（把d当作部分单位元或者说只有i为d的倍数且j为d的倍数时，才满足$d|gcd（i,j）$）=
${\sum }_{d=1}^{n}f[d]\ast {\sum }_{i=1}^{\lfloor n/d\rfloor }{\sum }_{j=1}^{\lfloor m/d\rfloor }{\sum }_{k|i且k|j}\mu [k]$=
${\sum }_{d=1}^{n}f[d]\ast {\sum }_{k=1}^{\lfloor n/d\rfloor }\mu [k]\ast F[k\ast d]$=
${\sum }_{d=1}^{n}f[d]\ast {\sum }_{d|i}^n\mu [i/d]\ast F[i]$（又把1当作单位元）=
(改变枚举顺序）${\sum }_{i=1}^{n}F[i]{\sum }_{d|i}f[d]\ast \mu [i/d]$.
对于$f[i]$我们可以用整除分块来减少枚举。

2.理解树状数组的用处

此时再来考虑a的影响。我们可以用离线处理数据，并让a非降序排列。
如果$f[x]<=a$,那么就可以用x为$F[i]$贡献系数。
还是为了减少枚举，我们边线性筛，边处理f，并以约数和从小到大排序f。
对于每个$f[x]<=a$，我们令x去贡献系数，并用树状数组维护系数。

3.线性筛中的细节

对于线性筛，我们的任务是求$\mu 和f$,至于$\mu$的求法这里不细讲（不会请翻预备知识），这里提一下f的求法。
方法一：时间复杂度$O(NlnN)\approx O(NlogN)$

for(int i=1;i<=10000000;i++)//Eratosthenes筛法
	for(int j=i;j<=10000000;j+=i)
		f[j]+=i;

方法二：时间复杂度$O(N)$
根据约数和公式：若$x=p_1^{c_1}\ast p_2^{c_2}\ast \dots \dots p_k^{c_k}$,
则$f(x)={\prod }_{i=1}^k{\sum }_{j=0}^{c_i}{p}_i^j$

那么我们可以实现简单继承：
当$p_j|i$时，$设k={\prod }_{i̸=y}{p}_i^{c_i},f(i\ast p_j)=f(k)\ast p_y^{c_y+1}+f[i]=f(k)+f(i)\ast p_y$(读者可自行证明，其实就是代入公式，这里不赘述）
否则，$f[i\ast p_y]=f[i]\ast f[p_y]（约数和是积性函数，读者亦可自行证明,其实又是代入公式)$

4.一个小优化：

因为答案要$mod{2}^{31}$，我们就用int存答案，溢出其实是$mod{2}^{32}$（需要理解），所以不影响答案，只是答案为负就转为正数即可。

代码:

#include<cstdio>
#include<cstring>
#include<algorithm>
#define g getchar()
using namespace std;
typedef long long ll;
const int N=1e5+10;
const int M=2e4+10;
const ll mod=1LL<<31;

void qr(int &x)
{
	char c=g;x=0;bool v=0;
	while(!(('0'<=c&&c<='9')||c=='-'))c=g;
	if(c=='-')v=1,c=g;
	while('0'<=c&&c<='9')x=x*10+c-'0',c=g;
	if(v)x=-x;
}

void write(int x)
{
	if(x/10)write(x/10);
	putchar(x%10+'0');
}

struct node
{
	int n,m,a,id;
	bool operator <(node b)const{return a<b.a;}
}a[M];

struct pnode
{
	int s,id;
	bool operator <(pnode b)const{return s<b.s;}
}f[N];

int T,n,m,prime[N/10],tot,miu[N],mo[N],ans[M],maxn,tr[N];
//f[i]表示i的约数和，对于mo[i](mother),设i除去最小的质数的所有因子后得到k，mo[i]=f[k];
bool v[N];

void get_prime()
{
	miu[1]=f[1].s=1;
	for(int i=2;i<=maxn;i++)
	{
		if(!v[i])prime[++tot]=i,miu[i]=-1,mo[i]=1,f[i].s=i+1;
		for(int j=1,k;j<=tot&&(ll)i*prime[j]<=maxn;j++)
		{
			v[k=i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				miu[k]=0;
				f[k].s=mo[i] + f[i].s*prime[j];
				mo[k]=mo[i];
				break;
			}
			else 
			{
				miu[k]=-miu[i];
				f[k].s=f[i].s*f[prime[j]].s;
				mo[k]=f[i].s;
			}
		}
	}
	for(int i=1;i<=maxn;i++)f[i].id=i;
	sort(f+1,f+maxn+1);
}

int lowbit(int x){return x&-x;}

void add(int x,int v)
{
	while(x<=maxn)
	{
		tr[x]+=v;
		x+=lowbit(x);
	}
}

void change(int i)
{
	int d=f[i].id,s=f[i].s;
	for(int j=1;j*d<=maxn;j++)if(miu[j])
		add(d*j,miu[j]*s);
}

int solve(int x)
{
	int tot=0;
	while(x>0)
	{
		tot+=tr[x];
		x-=lowbit(x);
	}
	return tot;
}

int main()
{
	qr(T);
	for(int i=1;i<=T;i++)
	{
		qr(a[i].n);qr(a[i].m);qr(a[i].a);a[i].id=i;
		if(a[i].n>a[i].m)swap(a[i].n,a[i].m);
		if(a[i].n>maxn)maxn=a[i].n;
	}
	sort(a+1,a+T+1);
	get_prime();
	for(int i=1,j=1;i<=T;i++)
	{
		while(f[j].s<=a[i].a&&j<=N-10)change(j++);
		int &sum=ans[a[i].id],n=a[i].n,m=a[i].m;
		sum=0;
		for(int l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			sum+=(solve(r)-solve(l-1))*(n/l)*(m/l);
		}
		if(sum<0)sum+=mod;
	}
	for(int i=1;i<=T;i++)write(ans[i]),puts("");
	return 0;
}