题目

这题不要用莫比乌斯反演,用欧拉反演更快


【分析】

设点 \((x,y)\) 的能量损失为 \(f(x,y)\)

则 \(\displaystyle Ans=\sum_{i=1}^n\sum_{j=1}^m f(i,j)\)


我们先解决 \(f(x,y)\) 具体是多少:

显然 是 \(2gcd(x,y)-1\)

证明如下 (不想了解的小伙伴可以跳过)

设点 \((x,y)\) 与 \((0,0)\) 之间有且仅有 \(k\) 个整数点,它们分别为 \(({x\over \lambda _1},{y\over \lambda_1})\) ~ \(({x\over \lambda _k},{y\over \lambda_k})\)

( \(\forall i\in[1,k],\lambda_i>1\) )

\(\because \forall i\in[1,k]\) 都有 \(({x\over \lambda_i},{y\over \lambda_i})\) 为整数点

即 \({x\over \lambda_i},{y\over \lambda_i}\in Z\)

因为不保证 \(\lambda\in Z\) 所以不能写成 \(\lambda_i\mid x\bigwedge\lambda_i\mid y\)

那么,我们令 \(p={x\over \lambda_i},q={y\over \lambda_i},p,q\in Z\)

\(\therefore gcd(x,y)=gcd(p\lambda_i,q\lambda_i)\)

\(\therefore {gcd(x,y)\over \lambda_i}=gcd(p,q)\)

\(\because\) \(\lambda_i\) 有且仅有 \(k\) 个取值

\(\therefore\) \(gcd(p,q)\) 有且仅有 \(k\) 个取值,且必定为整数

$\therefore $ 对 \(\forall n>1,{gcd(x,y)\over n}\) 有且仅有 \(k\) 个整数取值

注:\(n\) 可为分数

\(\therefore\) 对 \(\forall n>0,{gcd(x,y)\over n}\) 有且仅有 \((k+1)\) 个整数取值

这不就说明了 \(gcd(x,y)=k+1\) 吗?

\(\therefore f(x,y)=2k+1=2(k+1)-1=2gcd(x,y)-1\)


好了,题目很显然了,求 \(\displaystyle Ans=\sum_{i=1}^n\sum_{j=1}^m[2gcd(x,y)-1]\)

考虑 \(-1\) 的贡献: \(\displaystyle Ans=\sum_{i=1}^n\sum_{j=1}^m2gcd(x,y)-nm\)

求和符号中提取公因式 \(2\) :\(\displaystyle Ans=2\sum_{i=1}^n\sum_{j=1}^mgcd(x,y)-nm\)

接下来,有请 欧拉 繁衍 反演 闪亮登场:

它的根本是一个公式 \(\displaystyle \sum_{d\mid n}\boldsymbol \varphi(d)=n\)

这个公式怎么来的?我上次看到某个 dalao 给出了形象的证明:

假设我们有 \({1\over n},{2\over n},{3\over n}\dots {n-1\over n},{n\over n}\) 共 \(n\) 个分数。将它们能约分的全部约分,则约分后的分母 \(d\) ,必定满足 \(d\mid n\) 。且对 \(\forall d\mid n\) ,分母为 \(d\) 的分数一共出现 \(\boldsymbol \varphi(d)\) 次。所以 \(\displaystyle \sum_{d\mid n}\boldsymbol \varphi(d)=n\) 。

还是不了解的可以试着用莫比乌斯反演证明,也是可行的

所以原式又能继续化简: \(\displaystyle Ans=2\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid gcd(i,j)}\boldsymbol \varphi(d)-nm\)

因为 \(n,m\) 对答案的贡献是对称的,我们不妨设 \(n<m\)

调换顺序,枚举 \(d\)

\(\therefore\displaystyle Ans=2\sum_{d=1}^n\boldsymbol \varphi(d)\sum_{i=1}^n\sum_{j=1}^m[d\mid gcd(i,j)]-nm\)

\(\displaystyle \qquad\quad=2\sum_{d=1}^n\boldsymbol \varphi(d)\sum_{i=1}^n\sum_{j=1}^m[d\mid i\bigwedge d\mid j]-nm\)

\(\displaystyle \qquad\quad=2\sum_{d=1}^n\boldsymbol \varphi(d)\sum_{i=1}^n[d\mid i]\sum_{j=1}^m[d\mid j]-nm\)

分别考虑 \(i,j\) 的贡献:

\(\displaystyle Ans=2\sum_{d=1}^n\boldsymbol \varphi(d)\sum_{i=1}^{\lfloor{n\over d}\rfloor}[1\mid i]\sum_{j=1}^{\lfloor{m\over d}\rfloor}[1\mid j]-nm\)

\(\displaystyle \qquad=2\sum_{d=1}^n\boldsymbol \varphi(d){\lfloor{n\over d}\rfloor}{\lfloor{m\over d}\rfloor}-nm\)

剩下的整除分块搞一搞就出来了


【代码】

那本蒟蒻就放 我码风极丑的 代码了

#include<cstdio>
#include<algorithm>
using namespace std;
#define f(a,b,c,d) for(register int a=b,c=d;a<=c;a++)
#define g(a,b,c,d) for(register int a=b,c=d;a>=c;a--)
//#define LOCAL
typedef int i32;
typedef unsigned int u32;
typedef long long int i64;
typedef unsigned long long int u64;
const i32 MAXN=1e5;
typedef i64 ar[MAXN+10]; namespace HABIT{
template<typename T> inline T Max(T a) { return a; }
template<typename T,typename... Args> inline T Max(T a,Args... args){
T b=Max(args...);
return (a>b)?a:b;
}
template<typename T> inline T Min(T a) { return a; }
template<typename T,typename... Args> inline T Min(T a,Args... args){
T b=Min(args...);
return (a<b)?a:b;
} #ifdef LOCAL
inline char gc() { return getchar(); }
#else
inline char gc() {
static char s[1<<20|1]={0},*p1=s,*p2=s;
return (p1==p2)&&(p2=(p1=s)+fread(s,1,1<<20,stdin),p1==p2)?EOF:*(p1++);
}
#endif
inline i32 read(){
register i32 ans=0;register char c=gc();register bool neg=0;
while(c<48||c>57) neg^=!(c^'-'),c=gc();
while(c>=48&&c<=57) ans=(ans<<3)+(ans<<1)+(c^48),c=gc();
return neg?-ans:ans;
} char Output_Ans[1<<20|1],*Output_Cur=Output_Ans;
inline bool output() { Output_Cur-=fwrite(Output_Ans,1,Output_Cur-Output_Ans,stdout); }
inline void print(char c) { (Output_Cur-Output_Ans+1>>20)&&output(),*(Output_Cur++)=c; }
inline void print(i64 x){
char buf[21]={0}; (Output_Cur-Output_Ans+sprintf(buf,"%lld",x)>>20)&&output();
Output_Cur+=sprintf(Output_Cur,"%lld",x);
}
}
using namespace HABIT; ar ar_d_Fc,ar_d_Prime,ar_d_Phi;
inline void pre(i32 d_Lim){
i64 *ptr_d_Prime=ar_d_Prime;
ar_d_Phi[1]=1;
f(i,2,I,d_Lim){
if(!ar_d_Fc[i]) ar_d_Phi[i]=(*(ptr_d_Prime++)=ar_d_Fc[i]=i)-1;
for(register i64 *p=ar_d_Prime;p<ptr_d_Prime&&*p*i<=d_Lim;p++){
ar_d_Fc[*p*i]=*p;
ar_d_Phi[*p*i]=*p*ar_d_Phi[i];
if(*p<ar_d_Fc[i]) ar_d_Phi[*p*i]-=ar_d_Phi[i];
else break;
}
ar_d_Phi[i]+=ar_d_Phi[i-1];
}
} int main(){
i32 d_N=read(),d_M=read();
i64 lld_Ans=-1ll*d_N*d_M;
if(d_N>d_M) d_N^=d_M^=d_N^=d_M;
pre(d_N);
for(register i32 l=1,r;l<=d_N;l=r+1){
r=Min( d_N/(d_N/l) , d_M/(d_M/l) );
lld_Ans+=(ar_d_Phi[r]-ar_d_Phi[l-1])*(d_N/r)*(d_M/r)*2;
}
print(lld_Ans);
output();
return 0;
}

最后安利一下 本蒟蒻的博客

最新文章

  1. Android笔记——permission权限大全
  2. WARN util.NativeCodeLoader: Unable to load native-hadoop library for your platform... using builtin-java classes where applicable
  3. BZOJ2408 混乱的置换
  4. web发布 将各个文件夹输出合并到其自己的程序集 注意事项
  5. silverlight 4 tools for vs2010无法在vs2010 SP1上安装的解决办法
  6. bzoj4404: [Neerc2015]Binary vs Decimal
  7. docker管理shipyard中文版v3.0.2更新
  8. Java 操作符
  9. 04---Net基础加强
  10. PHP程序员如何突破成长瓶颈
  11. 在windows上配置jdk环境
  12. Emmet快速开发
  13. ArcEngine 9.3 学习笔记(六):图层符号化(COlorRamp,MarkerSymbol,LineSymbol,FillSymbol,TextSymbol,3DChartSymbol,使用ServerStyle符号库,FeatureRender,RasterRender)
  14. js文件中调用另一个js文件:
  15. .Net开发之旅(一个年少轻狂的程序员的感慨)
  16. ElasticSearch入门 附.Net Core例子
  17. 写前端页面步骤----vue+iview
  18. Redis入门指南之二(安装及配置)
  19. 半吊子的STM32 — SPI通信
  20. Android -- EventBus解析

热门文章

  1. Maven的安装和创建项目的过程
  2. 大二暑假第一周总结--初次安装配置Hadoop
  3. 一本通1166 求f(x,n)
  4. centos7+nginx+php+mysql环境搭建
  5. HDU_1059 多重背包问题
  6. tf.summary可视化参数
  7. C#高级编程(第9版) 第08章 委托、lambda表达式和事件 笔记
  8. python期末复习—列表
  9. Linux学习-课后练习(第二章命令)20200216
  10. 面向对象变成(OOP)-创建类和使用类