整除分块入门小记

如果你是从[莫比乌斯反演与数论函数](https://www.luogu.org/blog/command-block/mu-bi-wu-si-fan-yan-ji-ji-ying-yong)这里来,那就对了。 否则建议你去那里看看。 ## 一维整除分块 先看[题目](https://www.luogu.org/problemnew/show/P2261)。 求$\sum\limits_i^nk\%i$ 别看到$\%$一脸懵逼,老司机会把$k\%i$变成$k-i\left\lfloor\dfrac{k}{i}\right\rfloor$。 然后就变为$\sum\limits_i^n(k-i\left\lfloor\dfrac{k}{i}\right\rfloor)=kn-\sum\limits_i^ni\left\lfloor\dfrac{k}{i}\right\rfloor$ 问题就在于后面那个$\sum\limits_i^ni\left\lfloor\dfrac{k}{i}\right\rfloor$怎么求。 **弱化版:** 求$\sum\limits_i^n\left\lfloor\dfrac{n}{i}\right\rfloor$ 具体怎么做呢?观察到$n/i$只有$O(\sqrt{n})$(约两倍$\sqrt{n}$)种取值: ```cpp ################################################## 25/1=25 * ######################## 25/2=12 * ################ 25/3=8 * ############ 25/4=6 * ########## 25/5=5 * ######## 25/6=4 * ###### 25/7=3 * ###### 25/8=3 #### 25/9=2 * #### 25/10=2 #### 25/11=2 #### 25/12=2 ## 25/13=1 * ## 25/14=1 ## 25/15=1 ## 25/16=1 ## 25/17=1 ## 25/18=1 ## 25/19=1 ## 25/20=1 ## 25/21=1 ## 25/22=1 ## 25/23=1 ## 25/24=1 ## 25/25=1 ``` 得到$25/i$只有9种取值,也就是9块。 ```cpp ################################################## 25 ---------- ######################## 12 ---------- ################ 8 ---------- ############ 6 ---------- ########## 5 ---------- ######## 4 ---------- ###### 3*2 ###### ---------- #### 2*4 ... #### ---------- ## 1*12 ... ## ---------- ``` 如果能算出每块的“长度”与“高度”,就可以算出“面积总和”,也就是**总值**。 假如说我们求出了某一块,即$i=[l,r]$时,$\left\lfloor\dfrac{n}{i}\right\rfloor$都相同。 对于每一组$[l,r]$对应着一个长方形,长为$r-l+1$宽为$n/l$。 于是答案加上$(r-l+1)(n/l)$。 考虑我们知道了$[l,r]$这一块之后,如何求出下一块。 设$N/i'$与$N/i$相等,则$i'$的最大值为$N/(N/i)$。 - 这个没有沉下去之前,都有整除和非整除的概念 $\color{blue}\text{证明:}$ **必要性**:$\left\lfloor\dfrac{N}{i'}\right\rfloor$与$\left\lfloor\dfrac{N}{i}\right\rfloor$相等,也就是$\left\lfloor\dfrac{N}{i'}\right\rfloor≥\left\lfloor\dfrac{N}{i}\right\rfloor$(废话) ,得到$\dfrac{N}{i'}≥\left\lfloor\dfrac{N}{i}\right\rfloor$ 两边取倒数,不等号变向$\dfrac{i'}{N}≤1/(\left\lfloor\dfrac{N}{i}\right\rfloor)$ 两边同乘$N$,得到$i'≤N/(\left\lfloor\dfrac{N}{i}\right\rfloor)$ 也就是$i'≤\left\lfloor\dfrac{N}{\left\lfloor\dfrac{N}{i}\right\rfloor}\right\rfloor$ **充分性**:令$i'=\left\lfloor\dfrac{N}{\left\lfloor\dfrac{N}{i}\right\rfloor}\right\rfloor$得到$\left\lfloor\dfrac{N}{i'}\right\rfloor=\left\lfloor\dfrac{N}{\left\lfloor\dfrac{N}{\left\lfloor\dfrac{N}{i}\right\rfloor}\right\rfloor}\right\rfloor=\left\lfloor\dfrac{N}{i}\right\rfloor$ 明显可见$i'$的最大值为$N/(N/i)$。 也就是说我们知道了一个区间的左端点之后,我们能知道右端点,然后左端点++。 写出代码来是这样子的。 ```cpp #include<iostream> #define mod 998244353 using namespace std; long long n; long long calc(long long n) { long long ans=0,l=1,r; for (;l<=n;){ r=n/(n/l);//得到右端点 ans=(ans+(r-l+1)*(n/l))%mod; l=r+1;//前往下一块 }return ans; } int main() { cin>>n; cout<<calc(n)<<endl; } ``` 重复一遍,这个程序求的是$\sum\limits_i^n\left\lfloor\dfrac{n}{i}\right\rfloor$ 有一道[简单题](https://www.luogu.org/problemnew/show/P1403)可以用上述代码AC。 (不过建议大家还是要推一下,这是个经典结论。) 然而,这只是弱化版,我们真正的目的是~~蓝题~~求$\sum\limits_i^ni\left\lfloor\dfrac{k}{i}\right\rfloor$。 这里有一个统一的技巧。 假设我们要求$\sum\limits_i^nf(i)\left\lfloor\dfrac{k}{i}\right\rfloor$,($f(x)是某数论函数$) 我们看$\left\lfloor\dfrac{k}{i}\right\rfloor$值相同的某一块$[l,r]$ 得到$f(l)*\left\lfloor\dfrac{k}{l}\right\rfloor+f(l+1)*\left\lfloor\dfrac{k}{l}\right\rfloor+...+f(r)*\left\lfloor\dfrac{k}{l}\right\rfloor$ $=(f(l)+f(l+1)+...+f(r))\left\lfloor\dfrac{k}{l}\right\rfloor$ 这里我们明显看到:对于一块,$ANS+=(k/l)*\text{函数对应区间和}$ 我们维护一个前缀和即可。 **附:** 前面的$\sum\limits_i^n\left\lfloor\dfrac{n}{i}\right\rfloor$其实就是$\sum\limits_i^nI(i)\left\lfloor\dfrac{n}{i}\right\rfloor$,$I(i)[l,r]$的和就是(r-l+1),和上面推导的结果相同。 那么你也看出来了,$\sum\limits_i^ni\left\lfloor\dfrac{k}{i}\right\rfloor$就是$\sum\limits_i^nid(i)\left\lfloor\dfrac{k}{i}\right\rfloor$,那么对$id(i)$做一个区间和,不就是等差数列求和嘛。 **Code** **(有细节需要注意!)** ```cpp #include<iostream> using namespace std; long long n,k,ans; inline long long idsum(long long l,long long r) {return (r-l+1)*(l+r)>>1;} long long calc(long long n) { long long ans=0,l=1,r; for (;l<=n&&k/l;){ //这里要注意,假如k<n的话,k/l有可能出现0,直接跳过 r=min(k/(k/l),n); //这里要注意min一下n ans+=(k/l)*idsum(l,r); l=r+1; }return ans; } int main() { cin>>n>>k; cout<<n*k-calc(n); return 0; } ``` ## 二维整除分块 [P2260 [清华集训2012]模积和](https://www.luogu.org/problemnew/show/P2260) 好像是上面那道题的升级版？ 题意:求$\sum\limits_{i=1}^n\sum\limits_{j=1}^m[i≠j](n\%i)(m\%j)$,对某个数取模。 默认$n<m$，下面一定要注意大小关系！ $=\sum\limits_{i=1}^n\sum\limits_{j=1}^m(n-\lfloor n/i\rfloor*i)(m-\lfloor m/j\rfloor*j)[i≠j]$ $=\sum\limits_{i=1}^n\sum\limits_{j=1}^m(n-\lfloor n/i\rfloor*i)(m-\lfloor m/j\rfloor*j)-\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)(m-\lfloor m/i\rfloor*i)$ 前半部分: $=\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)\sum\limits_{j=1}^m(m-\lfloor m/j\rfloor*j)$ 设$F(n)=\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)$ 则前半部分$=F(n)F(m)$ $F(n)=\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)$ $=n^2-\sum\limits_{i=1}^n\lfloor n/i\rfloor*i$ 设$G(n,k)=\sum\limits_{i=1}^n\lfloor k/i\rfloor*i$,上式=$n^2-G(n,n)$ $G(n,k)$直接如上题整除分块即可。 后半部分: $\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)(m-\lfloor m/i\rfloor*i)$ $=\sum\limits_{i=1}^n(nm-m\lfloor n/i\rfloor*i-n\lfloor m/i\rfloor*i+\lfloor n/i\rfloor\lfloor m/i\rfloor*i^2)$ $=n^2m-mG(n,n)-nG(n,m)+\sum\limits_{i=1}^n\lfloor n/i\rfloor\lfloor m/i\rfloor*i^2$ What?$\sum\limits_{i=1}^n\lfloor n/i\rfloor\lfloor m/i\rfloor*i^2$这个东西好像并不是很好办？ > $\lfloor n/i\rfloor\lfloor m/i\rfloor$也有类似$\lfloor n/i\rfloor$的规律,也会形成块状。 > 令n=10,m=15,得到。 ```cpp ********************************************************************************(10/1)*(15/1)=150 ***********************************(10/2)*(15/2)=35 ***************(10/3)*(15/3)=15 ******(10/ 4)*(15/ 4)=6 ******(10/ 5)*(15/ 5)=6 **(10/ 6)*(15/ 6)=2 **(10/ 7)*(15/ 7)=2 *(10/ 8)*(15/ 8)=1 *(10/ 9)*(15/ 9)=1 *(10/10)*(15/10)=1 ``` > 可以看到也形成块状了，但是具体怎么分块? > 其实就是两个块状楼梯叠合在一起(形象理解!~~MC中掉落沙子~~)。 > 所以拐角的个数还是$O(\sqrt{n})$,(约4倍$\sqrt{n}$) > 我们找到了一个左端点, > 那么右端点=($(n/i)$楼梯的右端点,与$(m/i)$楼梯的右端点中,最近的一个)。 > 因为上文:$N/i'$与$N/i$相等,则$i'$的最大值为$N/(N/i)$。 > 所以右端点=$min(n/(n/i),m/(m/i))$; 还有一个$i^2$怎么办？，我们通过小学奥数（代码里面的S函数）可以得到$i^2$的区间和，其余同上。 综上，$ANS=F(n)F(m)-n^2m+mG(n,n)+nG(n,m)-\sum\limits_{i=1}^n\lfloor n/i\rfloor\lfloor m/i\rfloor*i^2$ 复杂度$O(\sqrt{n})$ **Code:** ```cpp #include<algorithm> #include<cstdio> #define mod 19940417 #define inv6 3323403 #define inv2 9970209 using namespace std; int n,m; int G(int n,int lim) { int l=1,r,ans=0; for (;l<=lim;l=r+1){ r=min(n/(n/l),lim); ans=(ans+1ll*(n/l)*(r+l)%mod*(r-l+1)%mod*inv2)%mod; }return ans; } int F(int n) {return (1ll*n*n-G(n,n))%mod;} int S(int n) {return 1ll*n*(n+1)%mod*(n+n+1)%mod*inv6%mod;} int main() { scanf("%d%d",&n,&m); if (n>m)swap(n,m); int ans=1ll*F(n)*F(m)%mod; ans=(ans+1ll*n*G(m,n)+1ll*m*G(n,n))%mod; ans=(ans-1ll*n*n%mod*m)%mod; int l=1,r; for (;l<=n;l=r+1){ r=min(n/(n/l),m/(m/l)); ans=(ans-1ll*(n/l)*(m/l)%mod*(S(r)-S(l-1)))%mod; }printf("%d",(ans+mod)%mod); return 0; } ``` 恭喜你,又学习到了**整除分块**!