NOI2010 能量采集

传送门

这个题观察一下之后发现，答案就是求

\[\sum_{i=1}^n\sum_{j=1}^mgcd(i,j) *2 - 1 \]

那我们的目标就是求\(\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)\)，先老套路转化成枚举gcd：

\[\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^nd[gcd(i,j)=d] \]

继续套路，把d除进去，然后因为莫比乌斯函数\(\sum_{d|n}\mu(d) = [n=1]\)的性质，我们转化式子

\[\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\sum_{p|i,p|j}\mu(p) \]

把各项\(\sum\)之间的限制条件互相交换一下。之后发现，其实对于每一个枚举的p，i，j能取到的个数是一个可计算的定值，也就是式子可以变成如下情况：

\[\sum_{d=1}^nd\sum_{p=1}^{\frac{n}{d}}\mu(p)\left\lfloor\frac{n}{dp}\right\rfloor\left\lfloor\frac{m}{dp}\right\rfloor \]

这个式子乱七八糟的。有大佬在题解中说：看着这堆dp就很不爽。那我们用T来替代dp。

把T带入其实要做很多变化，比较显然的不说了，主要是两层\(\sum\)的嵌套改变比较让人难受。\(\sum_p^{\frac{n}{d}}\)变成\(\sum_{T=1}^n\)比较好理解，这个上下同乘d即可。至于\(\sum_{d=1}^n\)如何变成\(\sum_{d|T}\)，这个我不是完全明白，大致的理解是，因为d必然是T的一个因子，所以其实这和d直接进行枚举基本是一样的，只是把原来的外层循环变到里面，里层循环变到外面。

那么式子就变成了：

\[\sum_{T=1}^n\sum_{d|T}d\mu(\frac{T}{d})\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor \]

然后根据莫比乌斯反演，我们发现\(\sum_{d|T}d\mu(\frac{T}{d}) = \varphi(T)\)

所以我们最后要求的就是

\[\sum_{T=1}^n\varphi(T)\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor \]

这个我们直接把欧拉函数筛出来求前缀和，用整除分块做一下就好了。

或者还有一种更简便的推导。

直接用欧拉函数的性质，我们对于每个gcd枚举他的因子，改为用这些因子的欧拉函数之和表示。

\[\sum_{i=1}^n\sum_{j=1}^m\sum_{d|i,d|j} \varphi(d) \]

还是把各项\(\sum\)之间的限制条件进行变换，直接可以得到

\[\sum_{d=1}^n\varphi(d)\left\lfloor\frac{n}{d}\right\rfloor\left\lfloor\frac{m}{d}\right\rfloor \]

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<set>
#include<vector>
#include<map>
#include<queue>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('\n')
#define fr friend inline
#define y1 poj
#define mp make_pair
#define pr pair<int,int>
#define fi first
#define sc second
#define pb push_back

using namespace std;
typedef long long ll;
const int M = 200005;
const int INF = 1000000009;
const double eps = 1e-7;

int read()
{
    int ans = 0,op = 1;char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
    return ans * op;
}

int a,b,c,d,k,p[M],mu[M],tot,n,phi[M],m;
ll sum[M];
bool np[M];

void euler()
{
   np[1] = 1,phi[1] = 1;
   rep(i,2,M-2)
   {
      if(!np[i]) p[++tot] = i,phi[i] = i-1;
      for(int j = 1;i * p[j] <= M-2;j++)
      {
     np[i * p[j]] = 1;
     if(!(i % p[j])) {phi[i * p[j]] = phi[i] * p[j];break;}
     phi[i * p[j]] = phi[i] * (p[j] - 1);
      }
   }
   rep(i,1,M-2) sum[i] = sum[i-1] + phi[i];
}

ll solve(ll a,ll b)
{
   ll lim = min(a,b),ans = 0;
   for(int i = 1,j;i <= lim;i = j + 1)
   {
      j = min(a / (a / i),b / (b / i));
      ans += (ll)(a / i) * (ll)(b / i) * (sum[j] - sum[i-1]);
   }
   return ans * 2 - a * b;
}

int main()
{
   euler();
   n = read(),m = read();
   printf("%lld\n",solve(n,m));
   return 0;
}