NOI2010 能量采集

http://www.lydsy.com/JudgeOnline/problem.php?id=2005

莫比乌斯函数。

不妨设N>=M。

我们发现，坐标(x,y)到(0,0)的连线上有gcd(x,y)个点（包含自己）。

所以答案就是：

$$sum_{1leq xleq N,1leq yleq M}2[gcd(x,y)-1]+1$$
$$=sum_{1leq xleq N,1leq yleq M}2gcd(x,y)-1$$
$$=2sum_{1leq xleq N,1leq yleq M}gcd(x,y)-N*M$$

$其实就是求：$

$$sum_{1leq xleq N,1leq yleq M}gcd(x,y)$$

$设g=gcd(x,y)，则x=gx',y=gy',gcd(x',y')=1，原式变成：$

$$sum_{1leq gleq M}gsum_{1leq x'leq left lfloorfrac{N}{g} ight floor,1leq y'leq left lfloorfrac{M}{g} ight floor}[(x',y')==1]$$

$有一个结论：$

$$sum_{d|n}mu (d)=[n==1]$$

$所以：$

$$sum_{1leq gleq M}gsum_{1leq x'leq left lfloorfrac{N}{g} ight floor,1leq y'leq left lfloorfrac{M}{g} ight floor}sum_{d|(x',y')}mu(d)$$

$$sum_{1leq gleq M}gsum_{1leq x'leq left lfloorfrac{N}{g} ight floor,1leq y'leq left lfloorfrac{M}{g} ight floor}sum_{d|x',d|y'}mu(d)$$

$设x'=dx'',y'=dy''，原式变成：$

$$sum_{1leq gleq M}sum_{1leq dleq left lfloorfrac{M}{g} ight floor}gmu (d)sum_{1leq x''leq left lfloor frac{N}{gd} ight floor}sum_{1leq y''leq left lfloor frac{M}{gd} ight floor}1$$

$$sum_{1leq gleq M}sum_{1leq dleq left lfloorfrac{M}{g} ight floor}gmu (d)leftlfloor frac{N}{gd} ight floor leftlfloor frac{M}{gd} ight floor$$

$我们枚举g，其实对于1leq dleq leftlfloor frac{M}{g} ight floor，leftlfloor frac{N}{gd} ight floor的取值最多只有sqrt{leftlfloor frac{N}{g} ight floor}个，leftlfloor frac{M}{gd} ight floor的取值最多只有sqrt{leftlfloor frac{M}{g} ight floor}个。$

$我们可以不直接枚举d，而同时从小到大枚举leftlfloor frac{N}{gd} ight floor和leftlfloor frac{M}{gd} ight floor的取值，然后记住mu (d)的前缀和，累加即可。$

#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<fstream>
#include<algorithm>
#include<cstring>
#include<string>
#include<cmath>
#include<queue>
#include<stack>
#include<map>
#include<utility>
#include<set>
#include<bitset>
#include<vector>
#include<functional>
#include<deque>
#include<cctype>
#include<climits>
#include<complex>
//#include<bits/stdc++.h>适用于CF,UOJ，但不适用于poj
 
using namespace std;

typedef long long LL;
typedef double DB;
typedef pair<int,int> PII;
typedef complex<DB> CP;

#define mmst(a,v) memset(a,v,sizeof(a))
#define mmcy(a,b) memcpy(a,b,sizeof(a))
#define re(i,a,b)  for(i=a;i<=b;i++)
#define red(i,a,b) for(i=a;i>=b;i--)
#define fi first
#define se second
#define m_p(a,b) make_pair(a,b)
#define SF scanf
#define PF printf
#define two(k) (1<<(k))

template<class T>inline T sqr(T x){return x*x;}
template<class T>inline void upmin(T &t,T tmp){if(t>tmp)t=tmp;}
template<class T>inline void upmax(T &t,T tmp){if(t<tmp)t=tmp;}

const DB EPS=1e-9;
inline int sgn(DB x){if(abs(x)<EPS)return 0;return(x>0)?1:-1;}
const DB Pi=acos(-1.0);

inline int gint()
  {
        int res=0;bool neg=0;char z;
        for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar());
        if(z==EOF)return 0;
        if(z=='-'){neg=1;z=getchar();}
        for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar());
        return (neg)?-res:res; 
    }
inline LL gll()
  {
      LL res=0;bool neg=0;char z;
        for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar());
        if(z==EOF)return 0;
        if(z=='-'){neg=1;z=getchar();}
        for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar());
        return (neg)?-res:res; 
  }

const int maxN=100000;

int N,M;
LL ans;

int flag[maxN+100],cnt,prime[maxN+100];
int mo[maxN+100],sum[maxN+100];

inline void build()
  {
      int i,j;
      mo[1]=1;
      re(i,2,M)
        {
            if(!flag[i])prime[++cnt]=i,mo[i]=-1;
            for(j=1;j<=cnt && prime[j]*i<=M;j++)
              {
                  flag[i*prime[j]]=1;
                  if(i%prime[j]==0)
                    {
                        mo[i*prime[j]]=0;
                        break;
                    }
                  else mo[i*prime[j]]=-mo[i];
              }
        }
      re(i,1,M)sum[i]=sum[i-1]+mo[i];
  }

int main()
  {
      freopen("energy.in","r",stdin);
        freopen("energy.out","w",stdout);
        int i;
        N=gint();M=gint();
        if(N<M)swap(N,M);
        build();
        int g;
        re(g,1,M)
          {
              int n=N/g,m=M/g,l,r;
              for(l=1;l<=m;l=r+1)
                {
                    int v1=n/l,v2=m/l;
                    r=min(n/v1,m/v2);
                    upmin(r,m);
                    ans+=LL(g)*LL(sum[r]-sum[l-1])*LL(v1)*LL(v2);
                }
          }
        ans=2*ans-LL(N)*LL(M);
        cout<<ans<<endl;
        return 0;
    }