题解 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 )

即 ${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) $\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;
}

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