这个题观察一下之后发现，答案就是求
[sum_{i=1}^nsum_{j=1}^mgcd(i,j) *2 - 1
]
那我们的目标就是求(sum_{i=1}^nsum_{j=1}^mgcd(i,j))，先老套路转化成枚举gcd：
[sum_{i=1}^nsum_{j=1}^msum_{d=1}^nd[gcd(i,j)=d]
]
继续套路，把d除进去，然后因为莫比乌斯函数(sum_{d|n}mu(d) = [n=1])的性质，我们转化式子
[sum_{d=1}^ndsum_{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}^ndsum_{p=1}^{frac{n}{d}}mu(p)leftlfloorfrac{n}{dp}
ight
floorleftlfloorfrac{m}{dp}
ight
floor
]
这个式子乱七八糟的。有大佬在题解中说：看着这堆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}^nsum_{d|T}dmu(frac{T}{d})leftlfloorfrac{n}{T}
ight
floorleftlfloorfrac{m}{T}
ight
floor
]
然后根据莫比乌斯反演，我们发现(sum_{d|T}dmu(frac{T}{d}) = varphi(T))
所以我们最后要求的就是
[sum_{T=1}^nvarphi(T)leftlfloorfrac{n}{T}
ight
floorleftlfloorfrac{m}{T}
ight
floor
]
这个我们直接把欧拉函数筛出来求前缀和，用整除分块做一下就好了。
或者还有一种更简便的推导。
直接用欧拉函数的性质，我们对于每个gcd枚举他的因子，改为用这些因子的欧拉函数之和表示。
[sum_{i=1}^nsum_{j=1}^msum_{d|i,d|j} varphi(d)
]
还是把各项(sum)之间的限制条件进行变换，直接可以得到
[sum_{d=1}^nvarphi(d)leftlfloorfrac{n}{d}
ight
floorleftlfloorfrac{m}{d}
ight
floor
]
#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('
')
#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
",solve(n,m));
return 0;
}