整除分块

整除分块

对于一个类似于 \(\sum_{i=1}^n \lfloor \frac{n}{i}\rfloor\) 的式子，可以用整除分块在 \(\mathcal{O}(\sqrt{n})\) 的时间内快速求解。

首先发现让 \(\lfloor \frac{n}{i}\rfloor\) 值不同的 \(i\) 很少，因为 \(n\bmod i\) 会有很多种取值。

然后定义决策点。令 \(d_i\) 表示第 \(i\) 个决策点。特别地，\(d_1=1\)。

若 \(\lfloor \frac{n}{j}\rfloor\neq \lfloor \frac{n}{j+1}\rfloor\)，则称 \(j+1\) 为决策点。

那么知道了定义该怎么求解呢？

考虑递推。

考虑 \(nk\le a\) 这个式子的意思，其中 \(n\in \Z,a,k>0\)。

令 \(\max\{n\}=m\)，则 \(m=\lfloor \frac{a}{k}\rfloor\)。

即 \(\lfloor \frac{a}{k}\rfloor\) 表示满足 \(nk\le a\) 的最大的 \(n\)。

那么已知 \(d_i\)，如何推出 \(d_{i+1}\)？

令 \(v=\lfloor \frac{n}{d_i}\rfloor\)，考虑求出极长区间 \(\left[l,r\right]\) 满足 \(\forall l\le i\le r,\lfloor \frac{n}{i}\rfloor=v\)。

那么 \(r\) 的意思就是满足 \(\lfloor \frac{n}{r}\rfloor=v\) 的最大值。

注意到 \(\lfloor \frac{n}{r}\rfloor=v\Rightarrow vr\le n\)，联系上文可知 \(r=\lfloor\frac{n}{v}\rfloor\)，也就是 \(r=\lfloor\frac{n}{\lfloor\frac{n}{d_i}\rfloor}\rfloor\)。

那么 \(d_{i+1}=r+1=\lfloor\frac{n}{\lfloor\frac{n}{d_i}\rfloor}\rfloor+1\)。

所以递推式为：

\[d_{i+1}=\lfloor\frac{n}{\lfloor\frac{n}{d_i}\rfloor}\rfloor+1 \]

那么我们就将 \(n\) 划分为了若干区间 \(\left[d_i,d_{i+1}-1\right]\)，值全是 \(\lfloor\frac{n}{d_i}\rfloor\)。

这显然是很好求的。

给个代码：

void init()
{
	d[1]=1;
	for(int i=2;;++i)
	{
		d[i]=n/(n/d[i-1])+1;
		if(d[i]>n)break;
	}
	return;
}

但是在实际运用上可以不用求出 \(d\) 数组，可以动态计算。具体看下文。

例题：P2261 [CQOI2007] 余数求和

首先拆式子：

\[\sum_{i=1}^n k\bmod i\\ =\sum_{i=1}^n k-i\lfloor\frac{k}{i}\rfloor\\ =kn-\sum_{i=1}^n i\lfloor\frac{k}{i}\rfloor \]

所以对 \(k\) 进行整除分块。

我们枚举 \(i\)，但表示的是当前决策点（即 \(d_i\)）

令 \(j=d_{i+1}-1=\lfloor\frac{k}{\lfloor\frac{k}{i}\rfloor}\rfloor\)，即右端点。

注意到当 \(i>k\) 时，以后都有 \(\lfloor\frac{k}{i}\rfloor=0\)，所以可以不用考虑了。

对于区间 \(\left[i,j\right]\)，令 \(v=\lfloor\frac{k}{i}\rfloor\)，对答案的贡献是：

\[\sum_{l=i}^j l\lfloor\frac{k}{l}\rfloor\\ =\sum_{l=i}^j lv\\ =v\sum_{l=i}^jl\\ =\frac{v(i+j)(j-i+1)}{2} \]

然后注意下边界就做完了。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ljl;
const int N=1e9+5,KN=1e5+5;
ljl n,k,d[KN],td,ans;

int main(){
	ios::sync_with_stdio(0);
	cin>>n>>k;ans=n*k;
	for(ljl i=1,r=1;i<=n;i=r+1/*更新端点*/)
	{
		if(!(k/i))r=n;//以后都是无意义的
		else r=min(n,k/(k/i));
		ans=ans-((k/i)*(i+r)*(r-i+1)/2);
	}cout<<ans<<'\n';
	return 0;
}

完结撒花。