题解 P1447 【[NOI2010]能量采集】

题目

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

---

【分析】

设点 \((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;
}

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