简要题意

给出三个整数 \(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;
}

最新文章

  1. LINUX运维实战案例之文件已删除但空间不释放问题的分析与解决办法
  2. django 模版语法及使用
  3. spring读取properties的方法
  4. 《APUE》第七章笔记
  5. mysql locktables
  6. terminal命令
  7. show_space/get_alert_log/get_trace_file
  8. 3D人脸识别预处理,3D face recognition preprocess
  9. .net项目svn项目管理文件清单
  10. Mego开发文档 - 快速概述
  11. 题解:[JSOI2004]平衡点 / 吊打XXX
  12. ubuntu配置小飞机
  13. 对配置文件 xml 进行操作
  14. dispatch_barrier_async--屏障是一个同步点
  15. 正则表达式中test,match,exec区别
  16. linux存储管理之磁盘阵列
  17. js中的setTimeout第三个参数
  18. Spring使用fastjson处理json数据
  19. 【Luogu4512】多项式除法(FFT)
  20. 关于Java类和包的那些事

热门文章

  1. 【算法】浅学 LCA
  2. 网络工程知识(二)VLAN的基础和配置:802.1q帧;Access、Trunk、Hybrid接口工作模式过程与配置;VLANIF的小实验
  3. Chrony时间同步服务
  4. 一篇文章带你了解服务器操作系统——Linux简单入门
  5. 利用递归的方式在JSON 数据中找到某个节点的多有父节点
  6. python(27)反射机制
  7. &lt;四&gt;构造函数初始化列表
  8. 2022春每日一题:Day 15
  9. Thinkphp6使用腾讯云发送短信步骤
  10. cJson 学习笔记