有 \(Q\) 组询问，每次询问给定 \(n,m,a\)，求
\[\sum_{i=1}^n\sum_{j=1}^m[\sigma(\gcd(i,j))\le a]\sigma(\gcd(i,j)) \]其中 \(\sigma(n)\) 为 \(n\) 的所有正因子之和。
答案对 \(2^{31}\) 取模。\(1\le Q\le 2\times10^4,1\le n,m\le 10^5,|a|\le 10^9\)
不妨设 \(n\le m\)。简单推一下式子
\[\begin{aligned}
ans&=\sum_{d=1}^n[\sigma(d)\le a]\sigma(d)\sum_{t=1}^{n/d}\mu(t)(n/dt)(m/dt)\\
&=\sum_{T=1}^n\left(\sum_{d\mid T}[\sigma(d)\le a]\sigma(d)\mu(T/d)\right)(n/T)(m/T)
\end{aligned}\]
令 \(f(T,a)=\sum\limits_{d\mid T}[\sigma(d)\le a]\sigma(d)\mu(T/d)\)。
可以枚举所有满足 \(\sigma(d)\le a\) 的 \(d\)，在 \(d\mid T\) 处累加贡献。
我们要处理多组询问，考虑将询问全部离线，按 \(a\) 从小到大排序。
则处理第 \(i\) 组询问时，只需枚举满足 \(\sigma(d)\in (a_{i-1},a_i]\) 的 \(d\)，在之前的基础上累加贡献。
然后用整除分块计算 \(ans\)。这里要对 \(f\) 进行区间求和。
也就是说我们要写一个东西，支持 \(f\) 的单点修改和区间求和。树状数组即可。
时间复杂度为 \(O(Q\log Q+n\log^2n+Q\sqrt n\log n)\)。
#include<stdio.h>
#include<algorithm>
typedef unsigned U;
const int M=100005,K=20005;
const U inf=2147483647;
const double eps=1e-12;
int N,n,m,Q,cnt,d,t,t0,p[K];
struct node { int i,n,m,a; }q[K];
bool v[M]; int s[M],g[M],mu[M],a[M],b[M];
U S,ans[M]; double I[M];
inline bool cmp(int x,int y) { return s[x]<s[y]; }
inline bool cmp2(node p,node q) { return p.a<q.a; }
void upd(int i,int x) { while (i<=N) b[i]+=x,i+=(i&-i); }
int qry(int i) { int s=0; while (i) s+=b[i],i&=i-1; return s; }
inline int min(int x,int y) { return x<y?x:y; }
int main() {
scanf("%d",&Q);
for (int i=1; i<=Q; ++i)
scanf("%d%d%d",&n,&m,&d),
n>m&&(t=n,n=m,m=t),
N<n&&(N=n),q[i]={i,n,m,d};
std::sort(q+1,q+Q+1,cmp2); //将询问离线
s[1]=1,mu[1]=1,I[1]=1;
for (int i=2; i<M; ++i) I[i]=1./i+eps;
for (int i=2; i<=N; ++i) { // 筛 mu 和 sigma
v[i]||(p[++cnt]=i,mu[i]=-1,g[i]=i+1,s[i]=i+1);
for (int j=1; j<=cnt&&(t=p[j]*i)<=N; ++j) {
v[t]=1;
if (i%p[j]==0) {
s[t]=s[i]/g[i]*(g[t]=p[j]*g[i]+1);
break;
}
mu[t]=-mu[i],g[t]=p[j]+1,s[t]=s[p[j]]*s[i];
}
}
for (int i=1; i<=N; ++i) a[i]=i;
std::sort(a+1,a+N+1,cmp);
for (int i=1,j=1; i<=Q; ++i) {
while (s[d=a[j]]<=q[i].a&&j<=N) { // 修改
for (int t=1,T=d; T<=N; ++t,T+=d)
upd(T,mu[t]*s[d]);
++j;
}
n=q[i].n,m=q[i].m,S=0,t0=0;
for (int l=1,r,t1,t2; l<=n; l=r+1) // 查询
r=min(n*I[t1=n*I[l]],m*I[t2=m*I[l]]),
t=qry(r),S+=(U)(t-t0)*t1*t2,t0=t;
ans[q[i].i]=(S&inf);
}
for (int i=1; i<=Q; ++i)
printf("%u\n",ans[i]);
return 0;
}