Description
作为数论之神的你相当喜欢数论，小A打算考考你，给你四个参数 \(x_{1}, x_{2}, y_{1}, y_{2}\) 。
请你求出 \(\sum_{i=x_{1}}^{x_{2}} \sum_{j=y_{1}}^{y_{2}}\left(\left[\frac{i}{x_{1}}\right]+\left[\frac{x_{2}}{i}\right]+\left[\frac{j}{y_{1}}\right]+\left[\frac{y_{2}}{j}\right]\right)^{2}\) 的计算结果。
其中 \([x]\) 表示高斯函数，意思是保留 \(x\) 的整数部分。
请对计算结果取模 \(10^{9}+7\) 。
\(T\le 100,1\le x1,x2,y1,y2\le 10^9\)。
Solution
\[\sum_{i=x_{1}}^{x_{2}} \sum_{j=y_{1}}^{y_{2}}\left(\left[\frac{i}{x_{1}}\right]+\left[\frac{x_{2}}{i}\right]+\left[\frac{j}{y_{1}}\right]+\left[\frac{y_{2}}{j}\right]\right)^{2}=\sum_{i=x_{1}}^{x_{2}} \sum_{j=y_{1}}^{y_{2}}\left(\left[\frac{i}{x_{1}}\right]+\left[\frac{x_{2}}{i}\right]\right)^{2}+\sum_{i=x1}^{x2}\sum_{j=y1}^{y2}\left(\left[\frac{j}{y1}\right]+\left[\frac{y2}{j}\right]\right)^2+2\sum_{i=x1}^{x2}\sum_{j=y1}^{y2}\left(\left[\frac{i}{x_{1}}\right]+\left[\frac{x_{2}}{i}\right]\right)\left(\left[\frac{j}{y1}\right]+\left[\frac{y2}{j}\right]\right)\\=(y2-y1+1)\sum_{i=x_{1}}^{x_{2}} \left(\left[\frac{i}{x_{1}}\right]+\left[\frac{x_{2}}{i}\right]\right)^{2}+(x2-x1+1)\sum_{j=y1}^{y2}\left(\left[\frac{j}{y1}\right]+\left[\frac{y2}{j}\right]\right)^2+2\sum_{i=x1}^{x2}\left(\left[\frac{i}{x_{1}}\right]+\left[\frac{x_{2}}{i}\right]\right)\sum_{j=y1}^{y2}\left(\left[\frac{j}{y1}\right]+\left[\frac{y2}{j}\right]\right)
\]
然后对于 \(\sum_{i=x1}^{x2}\left(\left[\frac{i}{x_{1}}\right]+\left[\frac{x_{2}}{i}\right]\right)^{2}\) 和 \(\sum_{i=x1}^{x2}\left(\left[\frac{i}{x_{1}}\right]+\left[\frac{x_{2}}{i}\right]\right)\) 就可以愉快的整除分块了（对于 \(y1,y2\) 也同理）
细节很多，请注意。
Code
#include<cstdio>
#define mod 1000000007
#define inv_6 166666668
#define ll long long
using namespace std;
ll T,x1,x2,y1,y2,ans;
ll gs(ll x) {return x*(x+1)%mod*(2*x+1)%mod*inv_6%mod;}
ll sum1(ll l,ll r) {return ((l+r)*(r-l+1)/2)%mod;}
ll sum2(ll l,ll r) {return (gs(r)-gs(l-1)+mod)%mod;}
ll f1(ll l,ll r,ll p)//sum_{i=l}^{r} (i/p)
{
if (l/p==r/p) return (r-l+1)*(l/p)%mod;
return ((((((l/p+1)*p-1)-l+1)*(l/p)%mod)+((r-(r/p*p)+1)*(r/p)%mod))%mod+(p*sum1(l/p+1,r/p-1)%mod))%mod;
}
ll f2(ll l,ll r,ll p)//sum_{i=l}^{r} (i/p)^2
{
if (l/p==r/p) return (r-l+1)*(l/p)%mod*(l/p)%mod;
return ((((((l/p+1)*p-1)-l+1)*(l/p)%mod*(l/p)%mod)+((r-(r/p*p)+1)*(r/p)%mod*(r/p)%mod)%mod)+p*sum2(l/p+1,r/p-1)%mod)%mod;
}
ll g(ll l,ll r)
{
ll res=0;
for (ll i=l,j;i<=r;i=j+1)
{
j=r/(r/i);
res=(res+(j-i+1)*(r/i)%mod)%mod;
}
return res;
}
ll calc1(ll l,ll r) {return (f1(l,r,l)+g(l,r))%mod;}//一次
ll calc2(ll l,ll r)//二次
{
ll res=0;
for (ll i=l,j;i<=r;i=j+1)
{
ll k=r/i;j=r/k;
res=(res+(j-i+1)*k%mod*k%mod)%mod;
res=(res+f2(i,j,l))%mod;
res=(res+2*k*f1(i,j,l)%mod);
}
return res;
}
int main()
{
freopen("B.in","r",stdin);
freopen("B.out","w",stdout);
scanf("%lld",&T);
while (T--)
{
scanf("%lld%lld%lld%lld",&x1,&x2,&y1,&y2);
ans=0;
ans=(ll)2*calc1(x1,x2)%mod*calc1(y1,y2)%mod;
ans=(ans+(y2-y1+1)*calc2(x1,x2)%mod)%mod;
ans=(ans+(x2-x1+1)*calc2(y1,y2)%mod)%mod;
printf("%lld\n",ans);
}
return 0;
}