题解 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
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