bzoj 3529: [Sdoi2014]数表

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

解题报告:
这题唯一的难点在于不大于a的限制
如果没有限制：
设i的约数和为\(f[i],g[i]\)为gcd为i的数对的个数
即求\(\sum_{i=1}^n\sum_{j=1}^nf[gcd(i,j)]\)
\(\sum_{i=1}^ng[i]*f[i]\)
\(\sum_{i=1}^nf[i]*\sum_{i|d}\mu(d/i)*(n/d)*(m/d)\)
然后是关键的变换:
\(\sum_{d=1}^n(n/d)*(m/d)*\sum_{i|d}f[i]*\mu(d/i)\)
这样我们就可以单独拿出后面的\(\sum_{i|d}f[i]*\mu(d/i)\)然后我们离线询问：
对于每一组我们把<=a的\(f[i]\)的贡献都加入到前缀和中，然后就是基本的数论分块,不过要支持插入和询问，所以要用到树状数组

#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define RG register
#define il inline
#define iter iterator
#define Max(a,b) ((a)>(b)?(a):(b))
#define Min(a,b) ((a)<(b)?(a):(b))
using namespace std;
typedef long long ll;
const int N=1e5+5,M=2e4+5;
const ll mod=(ll)1<<31;
struct node{
	int n,m,a,id;
	bool operator <(const node &pp)const{
		return a<pp.a;
	}
}q[M];
int n=0,num=0,prime[N],mu[N];bool vis[N];
struct yut{
	ll x;int id;
	bool operator <(const yut &pp)const{
		return x<pp.x;
	}
}a[N];
ll Tree[N];
void add(int sta,ll ad){
	for(int i=sta;i<=n;i+=(i&(-i))){
		Tree[i]+=ad;
		if(Tree[i]>=mod)Tree[i]-=mod;
	}
}
ll query(int sta){
	ll ret=0;
	for(int i=sta;i>=1;i-=(i&(-i))){
		ret+=Tree[i];if(ret>=mod)ret-=mod;
	}
	return ret;
}
void solve(){
	int to;mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]){
			prime[++num]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=num && i*prime[j]<=n;j++){
			to=i*prime[j];vis[to]=true;
			if(i%prime[j])mu[to]=-mu[i];
			else{mu[to]=0;break;}
		}
	}
	for(int i=1;i<=n;i++){
		a[i].id=i;
		for(int j=i;j<=n;j+=i)a[j].x+=i;
	}
	sort(a+1,a+n+1);
}
ll ot[M];
void work()
{
	int Q;cin>>Q;
	for(int i=1;i<=Q;i++){
		scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a);
		if(q[i].n>q[i].m)swap(q[i].n,q[i].m);
		if(q[i].n>n)n=q[i].n;q[i].id=i;
	}
	sort(q+1,q+Q+1);
	solve();
	ll ans=0;RG int j=1;
	for(int i=1;i<=Q;i++){
		while(j<=n && a[j].x<=q[i].a){
			for(int k=a[j].id;k<=n;k+=a[j].id)
				if(mu[k/a[j].id])add(k,mu[k/a[j].id]*a[j].x);
			j++;
		}
		ans=0;int k;
		for(int l=1;l<=q[i].n;l=k+1){
			k=Min(q[i].n/(q[i].n/l),q[i].m/(q[i].m/l));
			ans+=((query(k)-query(l-1))%mod+mod)%mod
			*(q[i].n/l)*(q[i].m/l)%mod;
			if(ans>=mod)ans-=mod;
		}
		ot[q[i].id]=ans;
	}
	for(int i=1;i<=Q;i++)
		printf("%lld\n",ot[i]);
}

int main()
{
	work();
	return 0;
}