傳送門

這個題觀察一下之後發現，答案就是求

\[\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\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;
}