bzoj2005（洛谷P1447）: [Noi2010]能量采集

题意

(1,1)到(n,m)每个点，与(0,0)点连线上的整点个数若为 $k$ ，其权值即为 $2k+1$，求权值和。

---

题解

画个样例我们能发现，整点个数也就是 $gcd(i,j)-1$

那么结果就是 $ans=\sum\limits_{i=1}^n\sum\limits_{j=1}^m2(gcd(i,j)-1)+1=\sum\limits_{i=1}^n\sum\limits_{j=1}^m(2gcd(i,j)-1)=2\sum\limits_{i=1}^n\sum\limits_{j=1}^mgcd(i,j)-nm$

问题就转化为求解 $\sum\limits_{i=1}^n\sum\limits_{j=1}^mgcd(i,j)$

这个可用容斥 $O(nlogn)$ 求

但用 $mobius$ 化简的话，也就是

$$\begin{aligned} \sum\limits_{i=1}^n\sum\limits_{j=1}^mgcd(i,j)&=\sum_{d=1}^{n}d\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{m}{d}}\varepsilon(gcd(i,j))\\ &=\sum_{d=1}^{n}d\sum\limits_{t=1}^{\frac{n}{d}}\mu(t)\frac{n}{dt}\frac{m}{dt}\\ &=\sum_{k=1}^n\sum_{t|k}\mu(t)\frac{k}{t}\frac{n}{k}\frac{m}{k}\qquad (令k=dt) \\ &=\sum_{k=1}^n\varphi(k)\frac{n}{k}\frac{m}{k}\qquad (\varphi(k)=\sum_{t|k}\mu(t)\frac{k}{t}) \\ \end{aligned}$$

---

代码

$mobius$ 版

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define ll long long
#define N 100010
int n,m,mu[N],prime[N],tot=0;
bool notprime[N];
void mobius(){
    mu[1]=1;
    memset(notprime,false,sizeof(notprime));
    for(int i=2;i<=100000;i++){
        if(!notprime[i]){
            prime[++tot]=i;
            mu[i]=-1;
        }for(int j=1;j<=tot,prime[j]*i<=100000;j++){
            notprime[prime[j]*i]=true;
            if(i%prime[j]==0){
                mu[prime[j]*i]=0;
                break;
            }mu[prime[j]*i]=-mu[i];
        }
    }for(int i=2;i<=100000;i++) mu[i]+=mu[i-1];
}
ll f(int x,int y){
    if(x>y) swap(x,y);
    ll ans=0;
    for(int i=1,last=1;i<=x;i=last+1){
        last=min(x/(x/i),y/(y/i));
        ans+=(ll)(mu[last]-mu[i-1])*(x/i)*(y/i);
    }return ans;
}
int main(){
    mobius();ll ans=0;
    scanf("%d%d",&n,&m);
    if(n>m) swap(n,m);
    for(int i=1;i<=n;i++) ans+=(ll)i*f(n/i,m/i);
    printf("%lld\n",2*ans-(ll)n*m);
    return 0;
}

$euler$ 版

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define ll long long
#define N 100010
int n,m,prime[N],tot=0;
bool notprime[N];
ll ans=0,phi[N];
void euler(){
    phi[1]=1;
    memset(notprime,false,sizeof(notprime));
    for(int i=2;i<=100000;i++){
        if(!notprime[i]) prime[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot,prime[j]*i<=100000;j++){
            notprime[prime[j]*i]=true;
            if(i%prime[j]==0){
                phi[prime[j]*i]=phi[i]*prime[j];
                break;
            }phi[prime[j]*i]=phi[prime[j]]*phi[i];
        }
    }for(int i=2;i<=100000;i++) phi[i]+=phi[i-1];
}
int main(){
    euler();
    scanf("%d%d",&n,&m);
    if(n>m) swap(n,m);
    for(int i=1,last=1;i<=n;i=last+1){
        last=min(n/(n/i),m/(m/i));
        ans+=(phi[last]-phi[i-1])*(n/i)*(m/i);
    }printf("%lld\n",2*ans-(ll)m*n);
    return 0;
}