显然任意一个顶点是关键点的正方形都可以用两个量来刻画：以其任意一条边对应的直角边平行于坐标轴的直角三角形的两直角边的长度 \(i,j\)（在下文中记这样的正方形为正方形 \((i,j)\)）。对于 \(i+j\) 相同的两种正方形，显然一个 \(n\times m\) 的点阵中这样的两正方形个数是相同。因此考虑预处理出 \(f_x\) 表示对于所有 \(i+j=x\) 的正方形 \((i,j)\)，该正方形中小方格个数之和。处理出 \(f_x\) 后，答案就是 \(\sum\limits_{i=1}^{\min(n,m)}(n-i+1)(m-i+1)f_i\)，可以 \(O(1)\) 回答询问，于是这样任务变为求 \(f_x\)。
考虑怎么等价地表达正方形 \((i,j)\) 中包含的方格个数，将正方形用类似于弦图的方式进行划分，里面有个边长为 \(|i-j|\) 的正方形，其中包含 \((i-j)^2\) 个小方格，四周的四个直角三角形包含的方格数是一样的，那么如果我们记 \(g(i,j)\) 表示两边长分别为 \(i,j\) 的直角三角形中包含的小方格的数量。那么正方形 \((i,j)\) 中包含的方格数就是 \((i-j)^2+4g(i,j)\)。
直接计算 \(g(i,j)\) 看着需要什么类欧几里得之类的东西处理下取整，如果你往这个方向想就大错特错了。一种等价的表达是，三角形中除了两直角边的格点个数（即，内部格点数量 \(+\) 斜边上的格点数量）。而内部格点数量是可以用 picks 定理直接计算即可，可以表达为 \(i,j,\gcd(i,j)\) 的线性组合。斜边上格点数量为 \(\gcd(i,j)-1\)，因此 \(g(i,j)\) 就可以表示为 \(i,j,\gcd(i,j)\) 的线性组合。这样求 \(f_x\) 只需要求一个类似于 \(\sum\limits_{i+j=x}Ai+Bj+C\gcd(i,j)+D+(i-j)^2\) 的东西，其中 \(A,B,C,D\) 为可以手算求得的常数。求 \(\sum\limits_{i+j=x}\gcd(i,j)\) 直接莫反即可。这部分是 trivial 了，时间复杂度 \(O(n\log n)\)。
式子就不细推了，相信如果看懂了自然可以推出来。
const int MAXN=1e6;
const int MOD=1e9+7;
int f[MAXN+5],sf2[MAXN+5],sf1[MAXN+5],sf0[MAXN+5];
int getsum1(int x){return 1ll*x*(x+1)/2%MOD;}
int getsum2(int x){return 1ll*x*(x+1)*(2*x+1)/6%MOD;}
int pr[MAXN+5],prcnt,phi[MAXN+5],vis[MAXN+5];
void init(int n){
phi[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i])pr[++prcnt]=i,phi[i]=i-1;
for(int j=1;j<=prcnt&&pr[j]*i<=n;j++){
vis[pr[j]*i]=1;
if(i%pr[j]==0){
phi[pr[j]*i]=phi[i]*pr[j];
break;
}
phi[pr[j]*i]=phi[pr[j]]*phi[i];
}
}
for(int i=1;i<=n;i++){
f[i]=1ll*i*i%MOD;
f[i]=(f[i]+2ll*getsum2(i-1))%MOD;
f[i]=(f[i]-2ll*i*(i-1)%MOD+MOD)%MOD;
}
for(int i=1;i<=n;i++)for(int j=i+i;j<=n;j+=i)
f[j]=(f[j]+2ll*i*phi[j/i])%MOD;
for(int i=1;i<=n;i++){
sf2[i]=(sf2[i-1]+1ll*f[i]*i%MOD*i)%MOD;
sf1[i]=(sf1[i-1]+1ll*f[i]*i)%MOD;
sf0[i]=(sf0[i-1]+f[i])%MOD;
}
}
int main(){
init(MAXN);int qu;scanf("%d",&qu);
while(qu--){
int n,m,res=0;scanf("%d%d",&n,&m);
printf("%d\n",(1ll*(n+1)*(m+1)%MOD*sf0[min(n,m)]+1ll*(MOD-sf1[min(n,m)])*(n+1+m+1)+sf2[min(n,m)])%MOD);
}
return 0;
}