若直线的斜率为0或者不存在斜率，则有$nC(m,3)+mC(n,3)$种方案。若直线的斜率不为0，只需考虑斜率为正的情况，最后答案再乘以2即可。枚举两个点的坐标差，设$t=\min(n,m)$，则有：
\[\begin{eqnarray*}
ans&=&\sum_{i=1}^n\sum_{j=1}^m(n-i)(m-j)(\gcd(i,j)-1)\\
&=&\sum_{d=1}^t\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d](n-i)(m-j)(\gcd(i,j)-1)\\
&=&\sum_{d=1}^t\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1](n-di)(m-dj)(d-1)\\
&=&\sum_{d=1}^t(d-1)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}(n-di)(m-dj)\sum_{k|\gcd(i,j)}\mu(k)\\
&=&\sum_{d=1}^t(d-1)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}(n-di)(m-dj)\sum_{k|i,k|j}\mu(k)\\
&=&\sum_{d=1}^t(d-1)\sum_{k}\mu(k)\sum_{k|i}\sum_{k|j}(n-di)(m-dj)\\
&=&\sum_{d=1}^t(d-1)\sum_{k}\mu(k)\sum_{k|i}(n-di)\sum_{k|j}(m-dj)\\
&=&\sum_{d=1}^t(d-1)\sum_{k}\mu(k)cal(n,d,k)cal(m,d,k)
\end{eqnarray*}\]
其中：
\[\begin{eqnarray*}
cal(n,d,k)&=&\sum_{k|i}(n-di)\\
&=&n\sum_{k|i}1-d\sum_{k|i}i\\
&=&n\lfloor\frac{n}{dk}\rfloor-\frac{dk(1+\lfloor\frac{n}{dk}\rfloor)\lfloor\frac{n}{dk}\rfloor}{2}
\end{eqnarray*}\]
设$D=dk$，则有：
\[\begin{eqnarray*}
F(n,D)=cal(n,d,k)&=&n\lfloor\frac{n}{dk}\rfloor-\frac{dk(1+\lfloor\frac{n}{dk}\rfloor)\lfloor\frac{n}{dk}\rfloor}{2}\\
&=&n\lfloor\frac{n}{D}\rfloor-\frac{D(1+\lfloor\frac{n}{D}\rfloor)\lfloor\frac{n}{D}\rfloor}{2}\\
ans&=&\sum_{d=1}^t(d-1)\sum_{k}\mu(k)cal(n,d,k)cal(m,d,k)\\
&=&\sum_{k=1}^t\mu(k)\sum_{d}(d-1)cal(n,d,k)cal(m,d,k)\\
&=&\sum_{k=1}^t\mu(k)\sum_{k|D}(\frac{D}{k}-1)F(n,D)F(m,D)\\
&=&\sum_{D=1}^t F(n,D)F(m,D)\sum_{k|D}\mu(k)(\frac{D}{k}-1)\\
&=&\sum_{D=1}^t F(n,D)F(m,D)(\sum_{k|D}\mu(k)\frac{D}{k}-\sum_{k|D}\mu(k))\\
&=&\sum_{D=1}^t F(n,D)F(m,D)(\varphi(D)-[D=1])
\end{eqnarray*}\]
综上所述，只需要一遍线性筛求出所有欧拉函数然后计算即可，时间复杂度$O(n)$。
#include<cstdio>
typedef long long ll;
const int N=50010,P=1000000007;
int n,m,t,i,j,ans,phi[N],v[N],p[N],tot;
inline ll cal(ll n,ll d){return n/d*n-d*(n/d)*(n/d+1)/2;}
ll C(ll n){return n*(n-1)*(n-2)/6;}
int main(){
scanf("%d%d",&n,&m);t=n<m?n:m;
for(phi[1]=1,i=2;i<=t;i++){
if(!v[i])p[++tot]=i,phi[i]=i-1;
for(j=1;j<=tot;j++){
if(i*p[j]>t)break;
v[i*p[j]]=1;
if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}else phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
for(i=1;i<=t;i++)ans=(ans+cal(n,i)%P*cal(m,i)%P*(phi[i]-(i==1))%P)%P;
return printf("%d",(ans*2%P+C(m)%P*n%P+C(n)%P*m%P)%P),0;
}