要求$ans=\sum_{i=1}^n \sum_{j=1}^m (n-i)(m-j)(gcd(i,j)-1)$
可以看做枚举矩阵的大小,然后左下右上必须取的方案数。
这是斜率单增的情况
然后大力反演即可。
最后$ans=ans*2+C(n,3)*m+C(m,3)*n$
$\Theta (n \log n)$
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define F(i,j,k) for (int i=j;i<=k;++i)
#define D(i,j,k) for (int i=j;i>=k;--i)
#define ll long long
#define md 1000000007
#define inf 0x3f3f3f3f
#define maxn 50005
ll vis[maxn],mu[maxn],pr[maxn],top;
void init1()
{
mu[1]=1;
F(i,2,maxn)
{
if (!vis[i])
{
mu[i]=-1;
pr[++top]=i;
}
F(j,1,top)
{
if ((ll)i*pr[j]>=maxn) break;
vis[i*pr[j]]=1;
if (i%pr[j]==0) {mu[i*pr[j]]=0;break;}
mu[i*pr[j]]=-mu[i];
}
}
}
ll f1[maxn],f2[maxn],f3[maxn],ans=0;
ll Sum(ll n)
{
n=(((n+1)*n)>>1)%md;
return n;
}
void solve(ll n,ll m)
{
if (n>m) swap(n,m);
ll ret=0;
F(d,1,n)
{
ll tmp=0;
F(p,1,n/d)
{
tmp+=mu[p]*(n/p/d)*(m/p/d)*m*n; tmp%=md;
tmp+=mu[p]*d*d*p*p*Sum(n/p/d)*Sum(m/p/d); tmp%=md;
tmp-=mu[p]*m*d*p*(m/p/d)*Sum(n/p/d); tmp%=md;
tmp-=mu[p]*n*d*p*(n/p/d)*Sum(m/p/d); tmp%=md;
}
ret+=tmp*(d-1);
}
ans=(2*ret)%md;
}
ll n,m;
ll C(ll n)
{
n%=md;
return (n*(n-1)*(n-2)/6)%md;
}
int main()
{
init1();//init2();
scanf("%lld%lld",&n,&m);
solve(n,m);
printf("%lld\n",(ans+(n*C(m))%md+(m*C(n))%md)%md);
}