【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 )描述一组数据。

Output

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

Sample Input

2

4 4 3

10 10 5

Sample Output

20

148

Hint

题目分析

设\(h\)为预处理的约数之和，根据之前【BZOJ4407】于神之怒加强版的化简套路，我们可以得到

\[\begin{split} ans&=\sum_{i=1}^n\sum_{j=1}^mh(\gcd(i,j))\\ &=\sum_{T=1}^n\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\sum_{d|T}\mu(\frac Td)h(d) \end{split} \]

按正常情况，我们应该把后面的\(\sum\limits_{d|T}\mu(\frac Td)h(d)\)预处理出前缀和形式。

然而，由于现在只有\(h(d)<a\)的部分才会产生贡献，所以我们可以将询问离线下来，按\(a\)排序，把不断更新的前缀和存在树状数组中即可。

代码实现

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<cstdlib>
#include<queue>
#define MAXN 0x7fffffff
typedef long long LL;
const int N=100005;
const LL mod=1ll<<31; 
using namespace std;
inline int Getint(){register int x=0,f=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*f;}
int mu[N],prime[N];
int h[N],low[N];
bool vis[N];

int ksm(int x,int k){
	int ret=1;
	while(k){
		if(k&1)ret=ret*x;
		x=x*x,k>>=1;
	}
	return ret;
}

struct Que{int n,m,a,id;bool operator <(const Que &q)const{return a<q.a;}}s[N];

#define mp(x,y) make_pair(x,y)
#define Pr pair<int,int>
priority_queue<Pr,vector<Pr>,greater<Pr> >q;

#define lowbit(x) x&(-x)
int tr[N],ans[N];
void Add(int x,int val){for(int i=x;i<=1e5;i+=lowbit(i))tr[i]=(1ll*tr[i]+val)%mod;}
int Query(int x){LL ret=0;for(int i=x;i;i-=lowbit(i))ret=(ret+tr[i])%mod;return (ret+mod)%mod;}
int main(){
	mu[1]=h[1]=1;
	q.push(mp(1,1));
	for(int i=2;i<=1e5;i++){
		if(!vis[i])prime[++prime[0]]=i,mu[i]=-1,h[i]=1+i,low[i]=2;
		for(int j=1;j<=prime[0]&&1ll*i*prime[j]<=1e5;j++){
			vis[i*prime[j]]=1;
			if(i%prime[j]==0){
				low[i*prime[j]]=low[i]+1;
				h[i*prime[j]]=1ll*h[i]*(ksm(prime[j],low[i]+1)-1)/(ksm(prime[j],low[i])-1);
				break;
			}
			low[i*prime[j]]=2;
			h[i*prime[j]]=h[i]*h[prime[j]];
			mu[i*prime[j]]=-mu[i];
		}
		q.push(mp(h[i],i));
	}
	
	int Q=Getint();
	for(int i=1;i<=Q;i++)s[i].n=Getint(),s[i].m=Getint(),s[i].a=Getint(),s[i].id=i;
	sort(s+1,s+1+Q);
	
	for(int t=1;t<=Q;t++){
		int n=s[t].n,m=s[t].m;
		if(n>m)swap(n,m);
		while(!q.empty()&&q.top().first<=s[t].a){
			int x=q.top().first,i=q.top().second;q.pop();
			for(int j=i;j<=1e5;j+=i)Add(j,mu[j/i]*x);
		}
		for(int l=1,r,i=s[t].id;l<=n;l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans[i]=(ans[i]+1ll*(n/l)*(m/l)%mod*(Query(r)-Query(l-1))%mod)%mod;
		}
	}
	for(int i=1;i<=Q;i++)cout<<(ans[i]+mod)%mod<<'\n';
	return 0;
}