这题不要用莫比乌斯反演，用欧拉反演更快
【分析】
设点 \((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;
}
最后安利一下 本蒟蒻的博客