完全不会的数学神题，正解留着以后填坑
将一个口胡的部分分做法，我们考虑计算格点多边形（包括三角形）面积的皮克公式：
[S=a+frac{1}{2}b-1 ext({a为图形内部节点个数，b为边界上的点数})
]
那么我们枚举每一个点，考虑算出它作为内部节点的总方案数以及作为边界上的点的方案数
然后考虑还有一个(-1)的常数，应该减去的是三角形的个数
所以我们大力组合容斥算出三角形个数就得到了一个优秀的(O(nm))做法
正解也许是推式子+容斥，放个CODE先坑了
Python3的：
n,m=sorted(map(int,input().split()))
mu=[i*i for i in range(0,n+1)]
for i in range(1,n+1):
for j in range(i+i,n+1,i):
mu[j]-=mu[i]
ans=n*(n-1)*m*(m-1)*(11*n*(n+1)*m*(m+1)+6*(n*(n+1)+m*(m+1)))//144
for i in range(1,n+1):
ans-=mu[i]*((n-i+n%i)*(n//i)//2)*((m-i+m%i)*(m//i)//2)
print(ans//3%1004535809)
C++的：
#include<cstdio>
#define RI register int
#define CI const int&
using namespace std;
const int N=3005,mod=1004535809;
int n,m,mu[N],ans;
inline void swap(int& x,int& y)
{
int t=x; x=y; y=t;
}
inline int sum(CI x,CI y)
{
int t=x+y; return t>=mod?t-mod:t;
}
inline void dec(int& x,CI y)
{
if ((x-=y)<0) x+=mod;
}
inline int inv(int x,int p=mod-2,int mul=1)
{
for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
int main()
{
RI i,j; scanf("%d%d",&n,&m); if (n>m) swap(n,m);
for (i=1;i<=n;++i) mu[i]=i*i; for (i=1;i<=n;++i)
for (j=i<<1;j<=n;j+=i) mu[j]-=mu[i];
ans=1LL*n*(n-1)%mod*m%mod*(m-1)%mod*sum(11LL*n*(n+1)%mod*m%mod*(m+1)%mod,
6LL*sum(1LL*n*(n+1)%mod,1LL*m*(m+1)%mod)%mod)%mod*inv(144)%mod;
for (i=1;i<=n;++i) dec(ans,1LL*mu[i]*(n-i+n%i)%mod*(n/i)%mod*
inv(2)%mod*(m-i+m%i)%mod*(m/i)%mod*inv(2)%mod);
return printf("%d",1LL*ans*inv(3)%mod),0;
}