能量采集（HYSBZ-2005）

Problem Description

栋栋有一块长方形的地，他在地上种了一种能量植物，这种植物可以采集太阳光的能量。在这些植物采集能量后，栋栋再使用一个能量汇集机器把这些植物采集到的能量汇集到一起。 栋栋的植物种得非常整齐，一共有n列，每列有m棵，植物的横竖间距都一样，因此对于每一棵植物，栋栋可以用一个坐标(x, y)来表示，其中x的范围是1至n，表示是在第x列，y的范围是1至m，表示是在第x列的第y棵。 由于能量汇集机器较大，不便移动，栋栋将它放在了一个角上，坐标正好是(0, 0)。 能量汇集机器在汇集的过程中有一定的能量损失。如果一棵植物与能量汇集机器连接而成的线段上有k棵植物，则能量的损失为2k + 1。例如，当能量汇集机器收集坐标为(2, 4)的植物时，由于连接线段上存在一棵植物(1, 2)，会产生3的能量损失。注意，如果一棵植物与能量汇集机器连接的线段上没有植物，则能量损失为1。现在要计算总的能量损失。 下面给出了一个能量采集的例子，其中n = 5，m = 4，一共有20棵植物，在每棵植物上标明了能量汇集机器收集它的能量时产生的能量损失。 在这个例子中，总共产生了36的能量损失。

Input

仅包含一行，为两个整数n和m。

Output

仅包含一个整数，表示总共产生的能量损失。

Examples

Input1

5 4

Output1

36

Input2

3 4

Output2

20

思路：

已知从一个点 (x,y) 到原点的能量损失为两点连成一条直线后，直线上的点数*2+1

那么考虑横纵坐标的关系，画个图算几组数据可以发现，将 (0,0)、(x,y) 两点连成一条线后，出去这条线上的两个端点有 GCD(x,y)-1 个点，那么每个点的贡献就是

设 f(d) 为 gcd(x,y)=d 的点的数量，那么问题就转化为计算：

由于 ，因此最终答案为：

因此，问题的关键是快速计算：

设 g(d) 为 满足 d|gcd(x,y) 的个数，即：

可以看出，g(k)、f(d) 符合莫比乌斯反演的形式：

于是有：

设 ，那么：

所以：

此时的复杂度是 O(n^2) 的，发现式子中有整除，利用整除分块两次求和，每次时间复杂度为 O(√n)，总复杂度为 O(√n*√n) 即 O(n)

Source Program

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<string>
#include<cstring>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<utility>
#include<stack>
#include<queue>
#include<vector>
#include<set>
#include<map>
#include<bitset>
#define PI acos(-1.0)
#define INF 0x3f3f3f3f
#define LL long long
#define Pair pair<int,int>
LL quickPow(LL a,LL b){ LL res=1; while(b){if(b&1)res*=a; a*=a; b>>=1;} return res; }
LL quickModPow(LL a,LL b,LL mod){ LL res=1; a=a%mod; while(b){if(b&1)res=(a*res)%mod; a=(a*a)%mod; b>>=1;} return res; }
LL getInv(LL a,LL mod){ return quickModPow(a,mod-2,mod); }
const double EPS = 1E-10;
const int MOD = 1E9+7;
const int N = 100000+5;
const int dx[] = {-1,1,0,0,-1,-1,1,1};
const int dy[] = {0,0,-1,1,-1,1,-1,1};
using namespace std;

int mu[N];
int prime[N];
bool bprime[N];
int cnt;
LL sum[N];
void getMu(int n){//线性筛求莫比乌斯函数
    cnt=0;
    mu[1]=1;//根据定义，μ(1)=1
    memset(bprime,false,sizeof(bprime));
 
    for(int i=2;i<=n;i++){//求2~n的莫比乌斯函数
        if(!bprime[i]){
            prime[++cnt]=i;//存储质数
            mu[i]=-1;//i为质数时，μ(1)=-1
        }
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){//枚举i之前的素数个数
            bprime[i*prime[j]]=true;//不是质数
            if(i%prime[j])//i不是prime[j]的整数倍时，i*prime[j]就不会包含相同质因子
                mu[i*prime[j]]=-mu[i];//mu[k]=mu[i]*mu[prime[j]]，因为prime[j]是质数，mu值为-1
            else{
                mu[i*prime[j]]=0;
                break;//留到后面再筛
            }
        }
    }
 
    for(int i=1;i<=n;i++)
        sum[i]=sum[i-1]+mu[i];
}
LL cal(LL n,LL m,LL minn){
    LL res=0;
    for(LL left=1,right;left<=minn;left=right+1){
        right=min(n/(n/left),m/(m/left));
        res+=(sum[right]-sum[left-1])*(n/left)*(m/left);
    } 
    return res;
}
int main(){
    getMu(100000);
    LL n,m;
    scanf("%lld%lld",&n,&m);

    LL res=0;
    LL minn=min(n,m);
    for (LL left=1,right;left<=minn;left=right+1){
        right=min(n/(n/left),m/(m/left));
        LL temp=(LL)(right+left)*(right-left+1)/2;
        res+=cal(n/left,m/left,minn/left)*temp;
    }
    res=res*2-n*m;
    printf("%lld\n",res);

    // system("pause");
    return 0; 
}