Portal

Description

进行\(T(T\leq10^5)\)次询问,每次给出\(x_1,x_2,y_1,y_2\)和\(d\)(均不超过\(10^5\)),求\(\sum_{i=x_1}^{x_2} \sum_{j=y_1}^{y_2} [gcd(i,j)=d]\)。

Solution

莫比乌斯反演入门题。

设\(calc(n,m)\)表示\(i\in[1,n],j\in[1,m]\)且\(gcd(i,j)=d\)的数对\((i,j)\)的个数。那么简单地进行容斥,可知\(ans=calc(x_2,y_2)-calc(x_1-1,y_2)-calc(x_2,y_1-1)+calc(x_1-1,x_2-1)\)。

于是考虑如何计算\(calc(n,m)\)。

\[f(d) = \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=d]
\]

\[\begin{align*}
F(x) &= \sum_{x|d} f(d) \\
&= \sum_{x|d} \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=d] \\
&= \sum_{k=1}^{⌊\frac{n}{x}⌋} \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=kx] \\
&= ⌊\frac{n}{x}⌋⌊\frac{m}{x}⌋
\end{align*}$$ $gcd(i,j)=kx \Leftrightarrow x|i$且$x|j$,那满足条件的$(i,j)$就有$⌊\frac{n}{x}⌋⌊\frac{m}{x}⌋$对。再进行莫比乌斯反演:
$$ f(x)= \sum_{x|d} \mu(\frac{d}{x}) F(d) = \sum_{x|d} \mu(\frac{d}{x})⌊\frac{n}{d}⌋⌊\frac{m}{d}⌋ = \sum_{k=1}^{⌊\frac{n}{x}⌋} \mu(k)⌊\frac{n}{kx}⌋⌊\frac{m}{kx}⌋ $$这个做法看起来是$O(\dfrac{n}{x})$的。不过由于$⌊\dfrac{n}{i}⌋$最多只有$\sqrt n$种取值,所以我们可以以$O(\sqrt n)$的复杂度进行计算。

|i| 1| 2| 3|4|5|6|7|8|9|10|11|12|13|14|15|
|-|-|-|-|-|-|-|-|-|-|--|--|--|--|--|--|
|15/i|15|7|5|3|3|2|2|1|1|1|1|1|1|1|1|

观察发现,一个取值为$v$的区间是以$⌊\frac{n}{v}⌋$结尾的,下一个区间是从$⌊\frac{n}{v}⌋+1$开始的,模拟这一性质去计算即可。若对于区间$k\in[L,R]$有$⌊\frac{n}{kx}⌋=v_1,⌊\frac{m}{kx}⌋=v_2$,那么该区间对答案的贡献为$v_1v_2\sum_{k=L}^R \mu(k)$,预处理出$\mu(x)$的前缀和即可。
> 时间复杂度$O(T\sqrt {10^5})$。

##Code
```cpp
//[HAOI2011]Problem b
#include <algorithm>
#include <cstdio>
using std::min; using std::swap;
typedef long long lint;
inline char gc()
{
static char now[1<<16],*s,*t;
if(s==t) {t=(s=now)+fread(now,1,1<<16,stdin); if(s==t) return EOF;}
return *s++;
}
inline int read()
{
int x=0; char ch=gc();
while(ch<'0'||'9'<ch) ch=gc();
while('0'<=ch&&ch<='9') x=x*10+ch-'0',ch=gc();
return x;
}
const int N=5e4+10;
int mu[N],pre[N];
int cntP,pr[N]; bool notP[N];
void getMu(int n)
{
mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!notP[i]) pr[++cntP]=i,mu[i]=-1;
for(int j=1;j<=cntP;j++)
{
if((lint)i*pr[j]>n) break;
int x=i*pr[j]; notP[x]=true;
if(i%pr[j]) mu[x]=-mu[i]; else {mu[x]=0; break;}
}
}
for(int i=1;i<=n;i++) pre[i]=pre[i-1]+mu[i];
}
int k;
lint calc(int x,int y)
{
x/=k,y/=k; if(x>y) swap(x,y);
lint res=0;
for(int L=1,R;L<=x;L=R+1)
{
int v1=x/L,v2=y/L; R=min(x/v1,y/v2);
res+=1LL*(pre[R]-pre[L-1])*v1*v2;
}
return res;
}
int main()
{
getMu(5e4);
int Q=read();
while(Q--)
{
int fr1=read(),to1=read(),fr2=read(),to2=read(); k=read();
printf("%lld\n",calc(to1,to2)-calc(fr1-1,to2)-calc(to1,fr2-1)+calc(fr1-1,fr2-1));
}
return 0;
}
```
##P.S.
同样的题[洛谷P2257](https://www.luogu.org/problemnew/show/P2257)。\]

最新文章

  1. linux环境下安装jdk
  2. gulp插件gulp-ruby-sass和livereload插件
  3. HTML标签CSS默认值研究
  4. Win 2003下IIS6+Mysql+php5.2  isapi搭建 升级php5.2到5.3测试 借助fastcgi实现
  5. JQuery获取浏览器窗口的高度和宽度
  6. lucas求组合数C(n,k)%p
  7. Java反射结合JDBC写的一个通用DAO
  8. 在Xcode中如何屏蔽某个源文件的编译警告信息
  9. (转)苹果推送通知服务教程 Apple Push Notification Services Tutorial
  10. EasyUI在MVC4中需要部分刷新页面时load()后页面变形问题!
  11. PHPStrom使用SASS,SCSS和Compass
  12. JSP中的九大隐式对象及四个作用域
  13. cryptojs的使用
  14. 20190320_head first pyhton学习笔记之构建发布
  15. 执行python文件报错SyntaxError: Non-ASCII character &#39;\xe8&#39; in file, but no encoding declared
  16. VS2013中Python学习笔记[基础入门]
  17. Bleve代码阅读(二)——Index Mapping
  18. Stacking调参总结
  19. linux调试工具glibc的演示分析
  20. 迷你MVVM框架 avalonjs 0.96发布

热门文章

  1. codevs 1262 不要把球传我 2012年CCC加拿大高中生信息学奥赛
  2. Fedora19添加和设置YUM源
  3. elasticsearch最全详细使用教程:入门、索引管理、映射详解、索引别名、分词器、文档管理、路由、搜索详解
  4. 浮动清楚浮动及position的用法
  5. 判断一个链表是否为回文结构 【题目】 给定一个链表的头节点head,请判断该链表是否为回 文结构。 例如: 1-&gt;2-&gt;1,返回true。 1-&gt;2-&gt;2-&gt;1,返回true。 15-&gt;6-&gt;15,返回true。 1-&gt;2-&gt;3,返回false。 进阶: 如果链表长度为N,时间复杂度达到O(N),额外空间复杂 度达到O(1)。
  6. 【整理】C#文件操作大全
  7. bootstrap历练实例: 导航元素中禁用的链接
  8. Bootstrap历练实例:垂直的按钮组
  9. shell脚本,计算1+2+3+....100等于多少?
  10. mongodb测试类