若直线的斜率为0或者不存在斜率，则有$nC(m,3)+mC(n,3)$种方案。若直线的斜率不为0，只需考虑斜率为正的情况，最后答案再乘以2即可。枚举两个点的坐标差，设$t=min(n,m)$，则有：
[egin{eqnarray*}
ans&=&sum_{i=1}^nsum_{j=1}^m(n-i)(m-j)(gcd(i,j)-1)\
&=&sum_{d=1}^tsum_{i=1}^nsum_{j=1}^m[gcd(i,j)=d](n-i)(m-j)(gcd(i,j)-1)\
&=&sum_{d=1}^tsum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{m}{d}
floor}[gcd(i,j)=1](n-di)(m-dj)(d-1)\
&=&sum_{d=1}^t(d-1)sum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{m}{d}
floor}(n-di)(m-dj)sum_{k|gcd(i,j)}mu(k)\
&=&sum_{d=1}^t(d-1)sum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{m}{d}
floor}(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*}]
其中：
[egin{eqnarray*}
cal(n,d,k)&=&sum_{k|i}(n-di)\
&=&nsum_{k|i}1-dsum_{k|i}i\
&=&nlfloorfrac{n}{dk}
floor-frac{dk(1+lfloorfrac{n}{dk}
floor)lfloorfrac{n}{dk}
floor}{2}
end{eqnarray*}]
设$D=dk$，则有：
[egin{eqnarray*}
F(n,D)=cal(n,d,k)&=&nlfloorfrac{n}{dk}
floor-frac{dk(1+lfloorfrac{n}{dk}
floor)lfloorfrac{n}{dk}
floor}{2}\
&=&nlfloorfrac{n}{D}
floor-frac{D(1+lfloorfrac{n}{D}
floor)lfloorfrac{n}{D}
floor}{2}\
ans&=&sum_{d=1}^t(d-1)sum_{k}mu(k)cal(n,d,k)cal(m,d,k)\
&=&sum_{k=1}^tmu(k)sum_{d}(d-1)cal(n,d,k)cal(m,d,k)\
&=&sum_{k=1}^tmu(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;
}