「SDOI 2014」数表

传送门

---

problem

有一张 $n\times m$ 的数表，其第 $i$ 行第 $j$ 列的数值为能同时整除 $i$ 和 $j$ 的所有自然数之和。

给定 $a$，计算数表中 $\le a$ 的数之和。

数据范围：$1\le n,m\le 10^5$，$1\le T\le 2\times 10^4$，$a\le 10^9$，其中 $T$ 是数据组数。

---

solution

设 $\sigma (n)$ 表示 $n$ 的所有约数之和，即 $\sigma (n)={\sum }_{d|n}d$。

我们先不考虑 $a$ 的限制来推一下式子：

$$ans=\sum _{i=1}^n\sum _{j=1}^m\sigma (\gcd (i,j))$$

由于 $a\mid i,a\mid j⟹a\mid \gcd (i,j)$，所以后面那个直接是 $\sigma (\gcd (i,j))$。

那么枚举 $\gcd$（设 $up=\min (n,m)$）：

$$ans=\sum _{k=1}^{up}\sigma (k)\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }[\gcd (i,j)=1]$$

然后就是莫反：

$$\begin{array}{rl}ans& =\sum _{k=1}^{up}\sigma (k)\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\sum _{d|i,d|j}\mu (d)\\ & =\sum _{k=1}^{up}\sigma (k)\sum _{d=1}^{\lfloor \frac{up}{k}\rfloor }\mu (d)\sum _{i=1}^{\lfloor \frac{n}{kd}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{kd}\rfloor }1\\ & =\sum _{k=1}^{up}\sigma (k)\sum _{d=1}^{\lfloor \frac{up}{k}\rfloor }\mu (d)\lfloor \frac{n}{kd}\rfloor \lfloor \frac{m}{kd}\rfloor \end{array}$$

设 $T=kd$，那么有：

$$ans=\sum _{T=1}^{up}\sum _{k|T}\sigma (k)\mu (\frac{T}{k})\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor$$

所以如果没有 $a$ 的限制，这道题直接就做完了。

我们设 $f(n)={\sum }_{k|n}\sigma (k)\mu (\frac{n}{k})$，那么只有 $\sigma (k)\le a$ 的部分才会对答案产生贡献。

把询问离线下来，对 $a$ 排序，然后把 $\sigma$ 的函数值也记下来排序，那么就用双指针，当进行下一个询问时，就先把 $\le a$ 的所有 $\sigma$ 值全部加进 $f$ 函数里（就是枚举倍数加就可以了）再处理。

发现维护 $f$ 时，需要单点加，区间查值（在整除分块的时候），因此我们用树状数组维护即可。

时间复杂度 $O(T\sqrt{n}\log n+n{\log }^{2}n)$。

PS：由于要对 $2^{31}$ 取模，我们可以先 unsigned int 自然溢出（对 $2^{32}$ 取模），最后再对 $2^{31}$ 取模。

---

code

#include<bits/stdc++.h>
#define N 100005
#define uint unsigned int
#define lowbit(x) (x&-x)
using namespace std;
struct Query{int n,m,a,id;}Q[N],S[N];
uint mu[N],sigma[N],ans[N],bit[N];
bool operator<(const Query &p,const Query &q)  {return p.a<q.a;}
void add(int x,uint val)  {for(;x<N;x+=lowbit(x))bit[x]+=val;}
uint Query(int x,uint ans=0)  {for(;x;x-=lowbit(x))ans+=bit[x];return ans;}
void modify(int x)  {for(int i=1;i*x<N;++i)add(i*x,sigma[x]*mu[i]);}
int sum,prime[N],mark[N];
void linear_sieves(){
	mu[1]=1;
	for(int i=2;i<N;++i){
		if(!mark[i])  prime[++sum]=i,mu[i]=-1;
		for(int j=1;j<=sum&&i*prime[j]<N;++j){
			mark[i*prime[j]]=1;
			if(i%prime[j])  mu[i*prime[j]]=-mu[i];
			else  {mu[i*prime[j]]=0;break;}
		}
	}
	for(int i=1;i<N;++i)
		for(int j=1;i*j<N;++j)  sigma[i*j]+=i;
	for(int i=1;i<N;++i)  S[i].id=i,S[i].a=sigma[i];
	sort(S+1,S+N);
}
uint solve(int n,int m){
	if(n>m)  swap(n,m);
	uint ans=0;
	for(int i=1,j;i<=n;i=j+1){
		j=min(n/(n/i),m/(m/i));
		ans+=(n/i)*(m/i)*(Query(j)-Query(i-1));
	}
	return ans;
}
int main(){
	int T;
	scanf("%d",&T);
	for(int i=1;i<=T;++i){
		scanf("%d%d%d",&Q[i].n,&Q[i].m,&Q[i].a),Q[i].id=i;
	}
	sort(Q+1,Q+T+1),linear_sieves();
	int now=0;
	for(int i=1;i<=T;++i){
		int n=Q[i].n,m=Q[i].m,a=Q[i].a;
		while(now<N-1&&S[now+1].a<=a)  modify(S[++now].id);
		ans[Q[i].id]=solve(n,m);
	}
	uint P=1ll<<31;
	for(int i=1;i<=T;++i)  printf("%d\n",ans[i]%P);
	return 0;
}