7396. 【2021.11.18 NOIP提高组联考】B. I love math

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;
}