题解 JSOISC 模拟赛【A.小K的算式】（数学，整除分块）

纪念一下这场膜你赛中唯一AC的题。

---

题面：

小 K 翻看高中时期留下的数学笔记，发现了一道题：

给定四个正整数 $x_1, x_2, y_1, y_2$，计算

$$\sum_{i=x_1}^{x_2}\sum_{j=y_1}^{y_2}\left(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}{i}\rfloor+\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}{j}\rfloor\right)^2$$

可是笔记上没有解答，由于上了大学之后就再也没碰过数学，小 K 现在不会做了，于是他来请教你。

由于答案可能很大，你只要输出答案对 $10^9+7$ 取模后的值。

数据范围：

对于 $30\%$ 的数据，$T = 10, x_2 - x_1 \leq 10^3, y_2 - y_1 \leq 10^3$

对于 $60\%$ 的数据，$T = 10, x_2 - x_1 \leq 10^6, y_2 - y_1 \leq 10^6$

对于 $100\%$的数据，$1 \leq T \leq 100, 1 \leq x_1 \leq x_2 \leq 10^9, 1 \leq y_1 \leq y_2 \leq 10^9$。

---

solution：

对于 30pts ，暴力 $O(n^2)$ 枚举

无脑推柿子。

先把字母换得好看点：

$$ans=\sum_{i=a}^{b}\sum_{j=c}^d (\lfloor\dfrac{i}{a}\rfloor+\lfloor\dfrac{b}{i}\rfloor+\lfloor\dfrac{j}{c}\rfloor+\lfloor\dfrac{d}{j}\rfloor)^2$$

经典套路，裂项求和：

$$ans1=\sum_{i=a}^{b}\sum_{j=c}^d 2(\lfloor\dfrac{i}{a}\rfloor+\lfloor\dfrac{b}{i}\rfloor)(\lfloor\dfrac{j}{c}\rfloor+\lfloor\dfrac{d}{j}\rfloor)$$

$$ans2=\sum_{i=a}^{b}\sum_{j=c}^d(\lfloor\dfrac{i}{a}\rfloor+\lfloor\dfrac{b}{i}\rfloor)^2$$

$$ans3=\sum_{i=a}^{b}\sum_{j=c}^d(\lfloor\dfrac{i}{c}\rfloor+\lfloor\dfrac{d}{i}\rfloor)^2$$

则有：

$$ans=ans1+ans2+ans3$$

对于三项分别进行处理：

对于 $ans1$ ，适当变换求和顺序：

$$ans1=2\sum_{i=a}^{b}(\lfloor\dfrac{i}{a}\rfloor+\lfloor\dfrac{b}{i}\rfloor)\sum_{j=c}^d (\lfloor\dfrac{j}{c}\rfloor+\lfloor\dfrac{d}{j}\rfloor)$$

对于 $ans2$ 和 $ans3$，发现有一层和式可以直接扔掉：

$$ans2=\sum_{i=a}^{b}(\lfloor\dfrac{i}{a}\rfloor+\lfloor\dfrac{b}{i}\rfloor)^2\times(d-c+1)$$

$$ans3=\sum_{i=c}^{d}(\lfloor\dfrac{i}{c}\rfloor+\lfloor\dfrac{d}{i}\rfloor)^2\times(a-b+1)$$

处理到这步时，我们发现已经可以 $O(n)$ 进行求解，60pts 到手。

但如果我们追求更高的分数，还要继续优化：

---

$$ans1=2\sum_{i=a}^{b}(\lfloor\dfrac{i}{a}\rfloor+\lfloor\dfrac{b}{i}\rfloor)\sum_{j=c}^d (\lfloor\dfrac{j}{c}\rfloor+\lfloor\dfrac{d}{j}\rfloor)$$

$$ans2=\sum_{i=a}^{b}(\lfloor\dfrac{i}{a}\rfloor^2+2\lfloor\dfrac{i}{a}\rfloor\lfloor\dfrac{b}{i}\rfloor+\lfloor\dfrac{b}{i}\rfloor^2)\times(d-c+1)$$

$$ans3=\sum_{i=c}^{d}(\lfloor\dfrac{i}{c}\rfloor^2+2\lfloor\dfrac{i}{c}\rfloor\lfloor\dfrac{d}{i}\rfloor+\lfloor\dfrac{d}{i}\rfloor^2)\times(a-b+1)$$

令

$$f1(n,m)=\sum_{i=n}^m\lfloor\frac{i}{n}\rfloor$$

$$f2(n,m)=\sum_{i=n}^m\lfloor\frac{m}{i}\rfloor$$

$$f3(n,m)=\sum_{i=n}^m\lfloor\frac{i}{n}\rfloor^2$$

$$f4(n,m)=\sum_{i=n}^m\lfloor\frac{m}{i}\rfloor^2$$

则有：

$$ans1=2\times(f1(a,b)+f2(a,b))\times(f1(c,d)+f2(c,d))$$

$$ans2=(d-c+1)f3(a,b)+(d-c+1)f4(a,b)+2(d-c+1)\sum_{i=a}^{b}\lfloor\dfrac{i}{a}\rfloor\lfloor\dfrac{b}{i}\rfloor$$

$$ans3=(b-a+1)f3(c,d)+(b-a+1)f4(c,d)+2(b-a+1)\sum_{i=c}^{d}\lfloor\dfrac{i}{c}\rfloor\lfloor\dfrac{d}{i}\rfloor$$

对于 $ans2$ 和 $ans3$ 的最后一项，回顾一下上次讲数论分块时提到的（不知道的看这里）：

有这样一个和式：

$$\sum_{i=1}^{n}f_i\lfloor\dfrac{n}{i}\rfloor$$

我们就可以通过数论分块来求解她。

那么在数论分块的过程中，对于块 $[l,r]$ ，她的贡献是：

$$ans=\lfloor\dfrac{n}{i}\rfloor\times(S_{r}-S_{l-1})$$

带入到 $\sum_{i=c}^{d}\lfloor\frac{i}{c}\rfloor\lfloor\frac{d}{i}\rfloor$ 中，前缀和 $S$ 是什么？

我们惊喜的发现 $S_{r}-S_{l-1}$ 就是 $f1(l-1,r)$

那么只要能快速求出 $f1,f2,f3,f4$ ，这个问题就解决了。

---

对于 $f2$ 和 $f4$，是个显然的数论分块b有坑，因为边界不等于分子，不能省略对边界的特判）。

对于 $f1$ 和 $f3$，乍一看不可做，但是我们考虑这个东西的本质：

$$\lfloor\dfrac{i}{x}\rfloor(a\leq i\leq b)$$

$i$ 原本是一个等差数列，除以了一个 $x$ 后：

$$1,1,1,2,2,2,3,3,3,4,4,4\ldots$$

变成了每个数连续出现 $x$ 遍的等差数列！

不妨就处理这个等差数列，最后乘上 $x$。

序列 $\lfloor\dfrac{i}{x}\rfloor(a\leq i\leq b)$，我们根据值进行分类，会有 $\lfloor\dfrac{b-a+1}{x}\rfloor$ 个整块，每块 $x$ 个数值相同；最后还有 $(b-a+1)\bmod x$ 个散数，值为 $\lfloor\dfrac{b}{x}\rfloor$。

那么接下来就好处理了：

对于 $f1$，直接套等差数列求和公式：

$$f1(a,b)=\lfloor\dfrac{b-a+1}{x}\rfloor(\lfloor\dfrac{b-a+1}{x}\rfloor+1)\div 2+((b-a+1)\bmod x)\lfloor\dfrac{b}{x}\rfloor$$

这个柿子可以 $O(1)$ 求出。

对于 $f3$，学过小学奥数的我们知道平方和公式：

$$f3(a,b)=\lfloor\dfrac{b-a+1}{x}\rfloor(\lfloor\dfrac{b-a+1}{x}\rfloor+1)(2\lfloor\dfrac{b-a+1}{x}\rfloor+1)\div 6+((b-a+1)\bmod x)\lfloor\dfrac{b}{x}\rfloor^2$$

这个柿子也可以 $O(1)$ 求出。

大功告成！代码也很好打，没有啥高级的数论算法。（就是细节挺多的。）

复杂度为 $O(T\sqrt n \times\texttt{巨大常数})$，瓶颈在于多次数论分块。

---

放一下代码：

已经把象征着血泪的调试语句删掉了

这份代码常数很大，请吸臭氧食用。

#pragma GCC optimize(3)
#pragma GCC optimize(2)
#include<bits/stdc++.h>
using namespace std;
#define maxn 1005
#define rg register
#define ll long long
#define inf 0x7f7f7f7f
#define djq 1000000007
inline void file(){
	freopen("test.in","r1",stdin);
}
char buf[1<<21],*p1=buf,*p2=buf;
inline int getc(){
    return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;
}
inline int read(){
	rg int ret=0,f=0;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')f=1;ch=getchar();}
    while(isdigit(ch)){ret=ret*10+ch-48;ch=getchar();}
    return f?-ret:ret;
}
char buffer[1<<25];
int q1=-1;	
const int q2=(1<<25)-1;
inline void flush(){ fwrite(buffer,1,q1+1,stdout),q1=-1; }
inline void putc(const char &x){ if(q1==q2) flush(); buffer[++q1]=x;}
void wrtn(int x){
  	static char buf[15];
  	static int len=-1;
  	if(x>=0){ do{ buf[++len]=x%10+48,x/=10; }while(x);
	}else{ putc('-'); do{ buf[++len]=-(x%10)+48,x/=10; }while(x);}
  	while(len>=0){ putc(buf[len]), --len; }
}
int _,a,b,c,d,inv6=166666668;
ll sum1,sum,sum3,sum2;
inline int pfh(int n){
	return 1ll*n*(n+1)%djq*(2*n+1)%djq*inv6%djq;
}
inline int work1(int n,int m){
	ll ret=0;
	for(rg int l=1,r=0;l<=n;l=r+1){
		r=m/l?min(m/(m/l),n):n;
		ret+=1ll*(r-l+1)*(m/l)%djq;
		ret%=djq;
	}
	return ret;
}
inline int work2(int a,int b){
	if(b<a) return 0;
	int tmp=(b-a+1)/a;
	ll ret=1ll*tmp*(tmp+1)/2%djq*a%djq;
	ret+=1ll*b/a*((b-a+1)%a)%djq;
	ret%=djq;
	return ret;
}
inline int work3(int n,int m){
	ll ret=0;
	for(rg int l=1,r=0;l<=n;l=r+1){
		r=m/l?min(m/(m/l),n):n;
		ret+=1ll*(r-l+1)*(m/l)%djq*(m/l)%djq;
		ret%=djq;
	}
	return ret;
}
inline int work4(int a,int b){
	if(b<a) return 0;
	ll ret=1ll*pfh((b-a+1)/a)*a%djq;
	ret+=1ll*(b/a)*(b/a)%djq*((b-a+1)%a)%djq;
	ret%=djq;
	return ret;
}
inline int work5(int a,int b){
	ll ret=0;
	for(rg int l=a,r=0;l<=b;l=r+1){
		r=b/(b/l);
		ret=ret+1ll*(work2(a,r)-work2(a,l-1)+djq)*(b/l)%djq;
		ret%=djq;
	}
	return ret;
}
int main(){
	_=read();
	while(_--){
		a=read(); b=read(); c=read(); d=read();
		sum=sum1=sum2=sum3=0;
		sum=1ll*work2(c,d)+work1(d,d)-work1(c-1,d)+djq;
		sum%=djq;
		sum1=1ll*work2(a,b)+work1(b,b)-work1(a-1,b)+djq;
		sum1%=djq;
		sum1*=sum*2; sum1%=djq;
		sum2=1ll*work4(a,b)+work3(b,b)-work3(a-1,b)+2*work5(a,b)+djq;
		sum2%=djq; sum2*=(d-c+1)%djq; sum2%=djq;
		sum3=1ll*work4(c,d)+work3(d,d)-work3(c-1,d)+2*work5(c,d)+djq;
		sum3%=djq; sum3*=(b-a+1)%djq; sum3%=djq;
		wrtn((sum1+sum2+sum3)%djq); putc('\n');
	}
   flush();
	return 0;
}

__EOF__

本文作者：Legitimity
本文链接：https://www.cnblogs.com/tiatto/p/15875636.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！