简要题意
给出三个整数 \(T,n,m\)，\(T\) 组询问，每组询问给出四个整数 \(i_1,j_1,i_2,j_2\)（数据保证 \(i_1,j_1\leq n\ \ i_2,j_2\leq m\)），计算：
\[\sum_{i=i_1}^{i_2}\sum_{j=j_1}^{j_2}{\gcd(i,j)}\bmod{10^9+7} \]\(1\leq T\leq500,1\leq n,m\leq5\times10^4\)
思路
这是一道欧拉反演题。
首先我们可以转换为求下面的式子：
\[F(x,y)=\sum_{i=1}^{x}\sum_{j=1}^{y}{\gcd(i,j)}\bmod{10^9+7} \]为什么，因为可以通过差分的思想将题目中的式子改为 \(F(i_2,j_2)-F(i_1-1,j_2)-F(i_2,j_1-1)-F(i_1-1,j_1-1)\)。而且从 \(1\) 开始求和的式子更容易变形。
下面设 \(F(x,y)\) 中 \(x<y\)。
定理 1（欧拉反演公式）：\(i=\sum_{d\mid i}{\varphi(d)}\)
证明：\(i=\sum_{d\mid i}\sum_{j=1}^{i}{[\gcd(i,j)=d]}=\sum_{d\mid i}\sum_{j=1}^{\lfloor\frac{i}{d}\rfloor}{[\gcd(i,j)=1]}=\sum_{d\mid i}\varphi(\frac{i}{d})=\sum_{d\mid i}{\varphi(d)}\)
推论 1：\(\gcd(i,j)=\sum_{p\mid i,p\mid j}\varphi(p)\)
证明：由欧拉反演公式 \(i=\sum_{d\mid i}{\varphi(d)}\) 得 \(\gcd(i,j)=\sum_{p\mid\gcd(i,j)}\varphi(p)\)。根据最大公约数的性质即可变形成 \(\gcd(i,j)=\sum_{p\mid i,p\mid j}\varphi(p)\)。
然后代入原式（以下过程省略模数）：
\[F(x,y)=\sum_{i=1}^{x}\sum_{j=1}^{y}{\sum_{p\mid i,p\mid j}\varphi(p)} \]改为先枚举 \(p\)，得：
\[F(x,y)=\sum_{p=1}^{x}\varphi(p)(\sum_{i=1}^{x}\sum_{j=1}^{y}{[p\mid i][p\mid j]}) \]将右边的式子化简得：
\[F(x,y)=\sum_{p=1}^{x}\varphi(p)\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor \]预处理 \(\varphi\) 函数前缀和，就可以使用数论分块计算了。
时间复杂度预处理 \(O(m)\)，单次询问 \(O(\sqrt{n})\)。
代码
#include <bits/stdc++.h>
#define int long long
using namespace std;
int phi[10000005],vis[10000005],prime[10000005],sum[10000005],tot;
const int mod = 1e9+7;
int M(const int x){
return (x%mod+mod)%mod;
}
void sieve(){
phi[1]=1;
for(int i=2;i<=50000;i++){
if(!vis[i]){
prime[++tot]=i;
phi[i]=i-1;
}
for(int j=1;j<=tot&&i*prime[j]<=50000;j++){
vis[i*prime[j]]=1;
if(!(i%prime[j])){
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}
for(int i=1;i<=50000;i++){
sum[i]=M(sum[i-1]+phi[i]);
}
}
int solve(int a,int b) {
int ans=0;
if (a>b) swap(a,b);
for (int l=1,r=0;l<=a;l=r+1) {
r=min(a/(a/l),b/(b/l));
ans=M(ans+M(M(M(sum[r]-sum[l-1])*(a/l))*(b/l)));
}
return ans;
}
int n,m,T;
signed main() {
sieve();
cin>>T>>n>>m;
while (T--) {
int i1,j1,i2,j2;
cin>>i1>>j1>>i2>>j2;
cout<<M(solve(i2,j2)-solve(i1-1,j2)-solve(i2,j1-1)+solve(i1-1,j1-1))<<'\n';
}
return 0;
}