整除分块入门小记

· · 个人记录

如果你是从莫比乌斯反演与数论函数这里来,那就对了。

否则建议你去那里看看。

一维整除分块

先看题目。

求\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})种取值:

##################################################  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块。

################################################## 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)。

明显可见i'的最大值为N/(N/i)。

也就是说我们知道了一个区间的左端点之后,我们能知道右端点,然后左端点++。

写出代码来是这样子的。

#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

有一道简单题可以用上述代码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

(有细节需要注意!)

#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]模积和

好像是上面那道题的升级版？

题意:求\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)

后半部分: $\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这个东西好像并不是很好办？

令n=10,m=15,得到。 
********************************************************************************(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:

#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;
}

恭喜你,又学习到了整除分块!