2021CCPC湘潭全国邀请赛-C-Calculate

2021 CCPC 湘潭全国邀请赛-C-Calculate

水一水，题目链接（HDU 6940）：这是链接

题面

You are given 4 positive integers x1,x2,y1,y2. Now you need to calculate

${\sum }_{i=x1}^{x2}{\sum }_{j=y1}^{y2}(\lfloor \frac{i}{x1}\rfloor +\lfloor \frac{x2}{i}\rfloor +\lfloor \frac{j}{y1}\rfloor +\lfloor \frac{y2}{j}\rfloor {)}^2$

where ⌊x⌋ denotes the biggest integer that is not bigger than x.

The answer may be too large, so you just need to output it modulo ($10^9+7$).

题意

给定$4$个数$x1,x2,y1,y2$，求上式。

思路

一、初见

首先明显可以看出上式由两种类型的部分组成:

一类是$\lfloor \frac{i}{x1}\rfloor$，另一类是$\lfloor \frac{x2}{i}\rfloor$

第一类

第一类经过打表可以发现规律，比如：$x1=4,x2=17$时：

i=4  ans=1
i=5  ans=1
i=6  ans=1
i=7  ans=1
i=8  ans=2
i=9  ans=2
i=10 ans=2
i=11 ans=2
i=12 ans=3
i=13 ans=3
i=14 ans=3
i=15 ans=3
i=16 ans=4
i=17 ans=4

不难看出每$x1$个数$\lfloor \frac{i}{x1}\rfloor$的值会加一。

总结出规律式子：

设$n$为完整的组数，$n=(x2-x2+1)/x1$，可以得到：$((x2-x1+1)-n\ast x1)$为末尾剩下的个数

则$ans=x1\ast \frac{n\ast (n+1)}{2}+((x2-x1+1)-n\ast x1)\ast (n+1)$

需要注意的是除数2需要使用逆元

第二类

很明显可以看出，第二类是一个典型的数论分块的形式。

数论分块模板：

int ans = 0,l = 0,r = 0;
for(l = 1; l <= n; l = r + 1) {
    r = n / (n / l);
    ans += (r - l + 1) * (n / l);
} 	

二、细思

我们对于整个式子已经有了初步的了解，那么如何处理这个复杂的式子呢？

我们最快能想到的就是展开，为了便于书写我们令$\lfloor \frac{i}{x1}\rfloor =A,\lfloor \frac{x2}{i}\rfloor =B,\lfloor \frac{j}{y1}\rfloor =C,\lfloor \frac{y2}{j}\rfloor =D$

展开得：${\sum }_{i=x1}^{x2}{\sum }_{j=y1}^{y2}{A}^2+B^2+C^2+D^2+2\ast (A\ast B+A\ast C+A\ast D+B\ast C+B\ast D+C\ast D)$

我们逐一分析：

1、$A^2\&C^2$

这个部分由两个相同的第一类式子相乘构成，我们展开书写$A^2$得：${\sum }_{i=x1}^{x2}{\sum }_{j=y1}^{y2}\lfloor \frac{i}{x1}{\rfloor }^2$

其实不难看出，就是第一类式子的平方，当$x1=4,x2=17$时，数列为：$1,1,1,1,4,4,4,4,9,9,9,9,16,16$。

同理可得：设$n$为完整的组数，$n=(x2-x2+1)/x1$，可以得到：$((x2-x1+1)-n\ast x1)$为末尾剩下的个数

则$ans=x1\ast \frac{n(n+1)(2n+1)}{6}+((x2-x1+1)-n\ast x1)\ast (n+1)\ast (n+1)$

2、$B^2\&D^2$

这个部分由两个相同的第二类式子相乘构成，我们展开书写$B^2$得：${\sum }_{i=x1}^{x2}{\sum }_{j=y1}^{y2}\lfloor \frac{x2}{i}{\rfloor }^2$

稍微对数论分块部分进行修改：

int ans = 0,l = 0,r = 0;
for(l = 1; l <= n; l = r + 1) {
    r = n / (n / l);
    ans += (r - l + 1) * (n / l)* (n / l);
} 	

==${\sum }_{i=x1}^{x2}{\sum }_{j=y1}^{y2}\lfloor \frac{x2}{i}{\rfloor }^2$不等于${\sum }_{j=y1}^{y2}({\sum }_{i=x1}^{x2}\lfloor \frac{x2}{i}\rfloor )\ast ({\sum }_{i=x1}^{x2}\lfloor \frac{x2}{i}\rfloor )$==此处不做赘述

注意不要漏了把最后的结果乘上$(y2-y1+1)$[由于有${\sum }_{j=y1}^{y2}$]

3、$A\ast B\&C\ast D$

这两个部分是由变量相同的一个第一类式子和一个第二类式子相乘得到，展开$A\ast B$得:${\sum }_{i=x1}^{x2}{\sum }_{j=y1}^{y2}\lfloor \frac{i}{x1}\rfloor \lfloor \frac{x2}{i}\rfloor$

这个式子的处理相对比较复杂，我们列表可以更好的理解，如$x1=4,x2=17$时：

i=4  left=1 right=4
i=5  left=1 right=3
i=6  left=1 right=2
i=7  left=1 right=2
i=8  left=2 right=2
i=9  left=2 right=1
i=10 left=2 right=1
i=11 left=2 right=1
i=12 left=3 right=1
i=13 left=3 right=1
i=14 left=3 right=1
i=15 left=3 right=1
i=16 left=4 right=1
i=17 left=4 right=1

我们前面可以很好的求出第一类式子的前缀和和第二类式子的分块区间，那么我们就有了一个大胆的想法：由于某一个分块区间中的$right$值恒定，且起止位置已知，我们可以用这个恒定值乘以$left$在此分块区间中的差分值，如：$right=2$时，左侧的和$SUM=sum(9)-sum(8-1)=(1+1+1+1+2)-(1+1+1)$，$right\ast sum$即为此区间的结果

同样不要忘了乘上$(y2-y1+1)$

4、其他部分

其他部分都是形如$A\ast D$的式子，此类式子是有变量不同的一个第一类式子和一个第二类式子相乘得到，展开$A\ast D$得${\sum }_{i=x1}^{x2}{\sum }_{j=y1}^{y2}\lfloor \frac{i}{x1}\rfloor \lfloor \frac{y2}{j}\rfloor$

这类式子有一个特点就是可以分别计算第一类式子和第二类式子并相乘得到结果，证明过程略。

三、完整AC代码（有点乱）

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll mod = 1e9 + 7;

ll fenkuai1(ll lower, ll high)
{
	ll ans = 0;
	for (ll l = lower, r = 0; l <= high; l = r + 1)
	{
		r = high / (high / l) % mod;
		ans = (ans + (r - l + 1) % mod * (high / l) % mod) % mod;
	}
	return ans % mod;
}

ll fenkuai2(ll lower, ll high)
{
	ll ans = 0;
	for (ll l = lower, r = 0; l <= high; l = r + 1)
	{
		r = high / (high / l) % mod;
		ans =
			(ans % mod + (r - l + 1) % mod * (high / l) % mod * (high / l) % mod) %	mod;
	}
	return ans % mod;
}

ll solve1(ll x1, ll x2)
{
	ll n = (x2 - x1 + 1) % mod / x1 % mod;
	return (((x1 * n) % mod * (n + 1) % mod * 500000004) % mod +
			(((x2 - x1 + 1) % mod - n * x1 % mod) % mod * (n + 1) % mod) % mod) % mod;
}

ll solve2(ll x1, ll x2)
{
	ll n = (x2 - x1 + 1) / x1 % mod;
	//	cout<<n<<' '<<n%mod*(n+1)%mod*(2*n+1)%mod*166666668%mod<<' '<<n * (n +
	//1) % mod * ((2 * n + 1) % mod) % mod * 166666668 %mod<<endl;
	return (((x1 * n) % mod * (n + 1) % mod * (2 * n % mod + 1) % mod *
			 166666668) %
				mod +
			(((x2 - x1 + 1) % mod - n * x1 % mod) % mod * (n + 1) % mod *
			 (n + 1) % mod) %
				mod) %
		   mod;
}
//((x2-x1+1)%mod-n*x1%mod)%mod*(n+1)%mod)%mod

ll solve3(ll lower, ll high)
{
	ll ans = 0;
	for (ll l = lower, r = 0; l <= high; l = r + 1)
	{
		r = high / (high / l) % mod;
		//	    cout<<l<<' '<<r<<' '<<(high / l)<<'
		//'<<solve1(lower,r)-solve1(lower,l-1)<<' '<<solve1(lower,r)<<endl;
		if (l > lower)
			ans = (ans % mod +
				   (high / l) % mod *
					   (solve1(lower, r) % mod - solve1(lower, l - 1) % mod) % mod) %
				  mod;
		else
			ans =
				(ans % mod + (high / l) % mod * (solve1(lower, r) % mod) % mod) % mod;
		//	    cout<<(r - l + 1) *(high /
		//l)*(solve1(lower,r)-solve1(lower,l-1))<<endl;
		while (ans < 0)
			ans += mod;
	}

	return ans % mod;
}
int main()
{
	ll x1, x2, y1, y2, A, B, C, D;
	int t;
	cin >> t;
	while (cin >> x1 >> x2 >> y1 >> y2)
	{
		ll ans = 0;
		A = solve1(x1, x2) % mod;
		B = fenkuai1(x1, x2) % mod;
		C = solve1(y1, y2) % mod;
		D = fenkuai1(y1, y2) % mod;
		//		cout<<solve2(x1,x2)*(y2-y1+1)%mod<<endl;
		//		cout<<fenkuai2(x1,x2)*(y2-y1+1)%mod<<endl;
		//		cout<<solve2(y1,y2)*(x2-x1+1)%mod<<endl;
		//		cout<<fenkuai2(y1,y2)*(x2-x1+1)%mod<<endl;
		//
		//		cout<<solve3(x1,x2)*(y2-y1+1)%mod<<endl;
		//		cout<<A*C%mod<<endl;
		//		cout<<A*D%mod<<endl;
		//		cout<<B*C%mod<<endl;
		//		cout<<B*D%mod<<endl;
		//		cout<<solve3(y1,y2)*(x2-x1+1)%mod<<endl;
		ans = (ans + fenkuai2(x1, x2) % mod * (y2 - y1 + 1) % mod) % mod;
		ans = (ans + solve2(x1, x2) % mod * (y2 - y1 + 1) % mod) % mod;
		ans = (ans + fenkuai2(y1, y2) % mod * (x2 - x1 + 1) % mod) % mod;
		ans = (ans + solve2(y1, y2) % mod * (x2 - x1 + 1) % mod) % mod;
		ans = (ans + solve3(x1, x2) % mod * (y2 - y1 + 1) % mod) % mod;
		ans = (ans + A * C % mod) % mod;
		ans = (ans + A * D % mod) % mod;
		ans = (ans + B * C % mod) % mod;
		ans = (ans + B * D % mod) % mod;
		ans = (ans + solve3(y1, y2) % mod * (x2 - x1 + 1) % mod) % mod;
		ans = (ans + solve3(x1, x2) % mod * (y2 - y1 + 1) % mod) % mod;
		ans = (ans + A * C % mod) % mod;
		ans = (ans + A * D % mod) % mod;
		ans = (ans + B * C % mod) % mod;
		ans = (ans + B * D % mod) % mod;
		ans = (ans + solve3(y1, y2) % mod * (x2 - x1 + 1) % mod) % mod;
		cout << ans << endl;
	}
}