BZOJ 2005 能量采集 (莫比乌斯反演)

题意：一个NM的矩形方格，对每个整点（i，j），每次的消耗为点（0,0）到（i，j）距离上整点（不包括(i,j)）个数2+1。
分析：设\(N\le M\)，不满足交换即可。则所有点被扫到的次数会在有1-N共N种可能。先假设每个点只要被扫到，那贡献就是2。
枚举点被扫到的次数p，设对于被扫到p次的点，对最终结果的贡献为\(f(p)\),有

\[f(p) = 2*p*\sum_{i=1}^{\lfloor \frac{N}{p} \rfloor} \sum_{j=1}^{\lfloor \frac{M}{p} \rfloor} [gcd(i,j)=1] \]

\([gcd(i,j)=1]指gcd(i,j)=1的对数\)
则最终答案为

\[ans = \sum_{p=1}^{N}f(p)-N*M \]

求\(f(p)\)的方法就是最基础的莫比乌斯反演
最后减去\(N*M\)是因为每个点如果是第一个被扫到，则对结果贡献只有1。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxn=2e5+5;
bool vis[maxn];
int prime[maxn],mu[maxn];
int sum[maxn];
void init(){
    memset(vis,false,sizeof(vis));
    mu[1] = 1;
    prime[0] = 0;
    int cnt=0;
    for(int i=2;i<maxn;++i){
        if(!vis[i]){
            mu[i] = -1;
            sum[i] = 1;
            prime[++cnt] = i;
        }
        for(int j=1;j<=cnt;++j){
            if(i*prime[j] >= maxn)  break;
            vis[i*prime[j]] = true;
            if(i % prime[j]){
                mu[i*prime[j]] = -mu[i];
            }
            else{
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
    for(int i=1;i<maxn;++i) sum[i] = sum[i-1] + mu[i];
}

void prepare(){
    int i,j,cnt=0;
    mu[1]=sum[1]=1;
    for(i=2;i<maxn;i++){
        if(!vis[i])
            prime[++cnt]=i,mu[i]=-1;
        for(j=1;prime[j]*i<maxn;j++){
            vis[prime[j]*i]=1;
            if(i%prime[j]==0){
                mu[prime[j]*i]=0;
                break;
            }
            mu[prime[j]*i]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];              //前缀和
    }
}

LL gao(LL n,LL m,LL p)      //枚举p
{
    LL ans = 0;
    LL a = n/p,b = m/p;
    for(int i=1,j;i<=a;i=j+1){
        j = min(a/(a/i),b/(b/i));
        ans += (sum[j]-sum[i-1])*(a/i)*(b/i);
    }
    return ans; 
}

int main()
{
    #ifndef ONLINE_JUDGE
        freopen("in.txt","r",stdin);
        freopen("out.txt","w",stdout);
    #endif
    init();
    LL a,b,c,d,k,N,M;
    while(~scanf("%lld %lld ",&N, &M)){
        LL res=0;
        if(N>M) swap(N,M);
        for(int i=1;i<=N;++i){
            res += 2*i*gao(N,M,i);
        }
        printf("%lld\n",res-N*M);
    } 
    return 0;
}