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×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)=leftlfloor n/a ight floor$ * $leftlfloor m/a ight floor$
表示 $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]*[gcd(i,j)==d](这是bool表达式，真为1，假为0）$ =
$sum_{d=1}^n sum_{i=1}^{leftlfloor n/d ight floor} sum_{j=1}^{leftlfloor m/d ight floor} f[d]*[gcd(i,j)==1]$ （把d当作部分单位元或者说只有i为d的倍数且j为d的倍数时，才满足 $d|gcd（i,j）$ ）=
$sum_{d=1}^n f[d]*sum_{i=1}^{leftlfloor n/d ight floor} sum_{j=1}^{leftlfloor m/d ight floor}sum_{k|i且k|j}mu[k]$ =
$sum_{d=1}^n f[d]*sum_{k=1}^{leftlfloor n/d ight floor}mu[k]*F[k*d]$ =
$sum_{d=1}^n f[d]*sum_{d|i}^nmu[i/d]*F[i]$ （又把1当作单位元）=
(改变枚举顺序） $sum_{i=1}^n F[i] sum_{d|i} f[d]*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}*p_2^{c_2}*……p_k^{c_k}$ ,
则 $f(x)=prod_{i=1}^ksum_{j=0}^{c_i}p_i^j$

那么我们可以实现简单继承：
当 $p_j|i$ 时， $设k=prod_{i e y}p_i^{c_i},f(i*p_j)=f(k)*p_y^{c_y+1}+f[i]=f(k)+f(i)*p_y$ (读者可自行证明，其实就是代入公式，这里不赘述）
否则， $f[i*p_y]=f[i]*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;
}