题解 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;
}
posted @ 2022-02-09 16:32  Legitimity  阅读(45)  评论(0编辑  收藏  举报