题解 JSOISC 模拟赛【A.小K的算式】(数学,整除分块)
纪念一下这场膜你赛中唯一AC的题。
题面:
小 K 翻看高中时期留下的数学笔记,发现了一道题:
给定四个正整数 \(x_1, x_2, y_1, y_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)\) 枚举
无脑推柿子。
先把字母换得好看点:
经典套路,裂项求和:
则有:
对于三项分别进行处理:
对于 \(ans1\) ,适当变换求和顺序:
对于 \(ans2\) 和 \(ans3\),发现有一层和式可以直接扔掉:
处理到这步时,我们发现已经可以 \(O(n)\) 进行求解,60pts 到手。
但如果我们追求更高的分数,还要继续优化:
令
则有:
对于 \(ans2\) 和 \(ans3\) 的最后一项,回顾一下上次讲数论分块时提到的(不知道的看这里):
有这样一个和式:
我们就可以通过数论分块来求解她。
那么在数论分块的过程中,对于块 \([l,r]\) ,她的贡献是:
带入到 \(\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\),乍一看不可做,但是我们考虑这个东西的本质:
\(i\) 原本是一个等差数列,除以了一个 \(x\) 后:
变成了每个数连续出现 \(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\),直接套等差数列求和公式:
这个柿子可以 \(O(1)\) 求出。
对于 \(f3\),学过小学奥数的我们知道平方和公式:
这个柿子也可以 \(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;
}