bzoj 2005: [Noi2010] 能量采集

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

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

从题意可以得出一个式子：∑（i=1->n）∑(j=1->m) 2*（i，j）+1;的（i，j）意思是i，j的最大公约数，因为我们可以发现每个点与（0，0）之间相隔的点数即为该点横纵坐标的最大公约数，这也是很显然的，它们的公约数将这两个数分为相等的几份，而每增加一份就是一个新点，，举个例子：（8，4），它们的最大公约数是4，也就是说横纵坐标各可以分成相等的四份，（2，1），而每累加一次就是一个点的坐标，而该点显然在（8，4）和（0，0）之间；比如（2，1）是相隔的一个点，（4，2）是另一个点，（6，3）又是一个点。。。。。。所以说我们只要求出上面的那个式子即可；

现在让我们来看这个式子，时间复杂度是n^3的，很明显不可以，那么我们变形一下：

∑（i=1->n）∑(j=1->m) 2*（i，j）+1 =∑（i=1->n）∑(j=1->m)∑（d=1->min(i,j)）[d | i][d | j] φ(d);

这个式子的正确性也是显然的，对于任意两个数，它们的公约数也是最大公约数的约数，而最大公约数的各因子的欧拉函数之和就等于最大公约数。所以上面两个式子是等价的；

但是时间复杂度好像没有变化，那我们再变一下：

∑（i=1->n）∑(j=1->m)∑（d=1->min(i,j)）[d | i][d | j] φ(d)=∑(d=1->min(n,m)) φ(d) ∑(i=1->n) ∑(j=1->m) [d | i][d | j] ;

这个变形应该没什么问题，只是一个交换律的应用；下面再来一次变形：

∑(d=1->min(n,m)) φ(d) ∑(i=1->n) ∑(j=1->m) [d | i][d | j] = ∑(d=1->min(n,m)) φ(d) [ n/d ] [ m/d ] ;

这个式子我来解释一下，很显然我们上一个式子中的∑(i=1->n) ∑(j=1->m) [d | i][d | j] 这一项 所求的就是在n和m中能够整除d的数的个数，，而这也就等于[ n/d ][ m/d ];（[]是下取整符号）

那么我们对于欧拉函数线性预处理出来，然后计算即可；

代码如下：

#include<cstdio>
#include<iostream>
#define N 1000010
using namespace std;
int n,m,to,phi
,pri[10000];
long long ans; bool v
;
void work() {
phi[0]=phi[1]=phi[2]=1;
for(long long i=2;i<=to;i++) {
if(!v[i]) pri[++pri[0]]=i,phi[i]=i-1;
for(long long j=1;j<=pri[0]&&i*pri[j]<=to;j++) {
v[i*pri[j]]=true;
if(i%pri[j]) phi[i*pri[j]]=phi[i]*(pri[j]-1);
else { phi[i*pri[j]]=phi[i]*pri[j];break; }
}
}
}
int main () {
scanf("%d%d",&n,&m);to=min(n,m);work();
for(int d=1;d<=to;d++)  ans+=(long long) phi[d]*(n/d)*(m/d);
ans=ans*2;ans-=(long long) m*n;
printf("%I64d",ans);
}