笔记：整除分块

· · 算法·理论

核心思想

整除分块是这样一个算法：

求 \displaystyle \sum^{n}_{i=1}\lfloor \dfrac{n}{i}\rfloor 的值，时间复杂度为 O(\sqrt{n})。

是怎么实现的呢？

首先引入一个重要的概念——决策点。

这个概念是怎么来的？对于每一个 i，可能会有很多的 \lfloor \dfrac{n}{i} \rfloor 是相等的。比如当 n=10 时：

i \lfloor\dfrac{n}{i}\rfloor
1 10
2 5
3 3
4 2
5 2
6 1
7 1
8 1
9 1
10 1

容易发现，出现了很多的 2 和 1。我们将一段连续且相等的数字定义为一个连续段。那么，对于每一个连续段，其开头位置就被定义为决策点。有一个直观的图形表示如下：

其中被蓝色圈起来了的叫做决策点。

利用决策点，我们就可以对一段连续且相等的数进行统一处理。接下来，为了表述简洁，我们使用 d_i 表示从小到大第 i 个决策点。

容易得出一个式子：

\sum^{n}_{i=1}\lfloor\dfrac{n}{i}\rfloor=\sum(d_{i+1}-d_i)\lfloor\dfrac{n}{d_i}\rfloor

其中 d_{i+1}-d_i 表示 d_i 这个决策点所代表连续段的长度。然后 \lfloor\dfrac{n}{d_i}\rfloor 代表的就是 d_i 这个决策点所代表连续段中的值。因此等式右边的意思就是把每个连续段的和给统计出来。

所以，整除分块的关键在于从 d_i 推出 d_{i+1}。具体应该怎么做？有如下公式：

d_{i+1}=\lfloor\dfrac{n}{\lfloor\dfrac{n}{d_i}\rfloor}\rfloor+1

如何理解？首先，\lfloor\dfrac{n}{d_i}\rfloor 就是当前连续段中的值。接着我们思考 \lfloor\dfrac{n}{k}\rfloor 代表着什么。其实这个就代表满足 \lfloor\dfrac{n}{i}\rfloor\geq k 的 i 的最大值。又由于这里 \lfloor\dfrac{n}{d_i}\rfloor 是一个实际可以取到的数，所以 \lfloor\dfrac{n}{\lfloor\dfrac{n}{d_i}\rfloor}\rfloor 代表的就是 d_i 所代表连续段的右端点。+1 当然就是下一个连续段的决策点啦~

于是乎，代码就很好写了。用一道例题给出代码：

例题 1

H(n)。

分析一下代码发现就是求 \displaystyle\sum^{n}_{i=1}\lfloor\dfrac{n}{i}\rfloor。于是直接用之前的分析即可。代码非常简洁，比暴力长不了多少：

#include <bits/stdc++.h>
#define endl '\n'
using namespace std;
long long H(int n) {
    long long res = 0;
    for(int i = 1, j; i <= n; i = j) {
        j = n / (n / i) + 1;
        res += 1LL * (j - i) * (n / i);
    }
    return res;
}
int main() {
    ios :: sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    int T;
    cin >> T;
    while(T--) {
        int n;
        cin >> n;
        cout << H(n) << endl;
    }
    return 0;
}

同时，我们分析一下这个代码的时间复杂度。

显然，每一次循环会处理一个连续段，所以时间复杂度是 O(\text{连续段数量}) 的。接下来证明连续段数量为 O(\sqrt{n})。

首先，对于 i\in[1,\sqrt{n}]，i 只有 \sqrt{n} 个，所以这一部分的连续段数量不超过 \sqrt{n} 个。对于 i\in[\sqrt{n},n] 这一部分，\lfloor\dfrac{n}{i}\rfloor\leq\sqrt{n}，因此连续段数量也不超过 \sqrt{n} 个。综上，连续段数量不超过 2\sqrt{n}。

所以这个算法是 O(\sqrt{n}) 的。

例题 2

余数求和。

首先发现题目要求的式子是 \displaystyle\sum^{n}_{i=1}k\bmod i，这个式子好像和整除分块没有什么关系啊？其实，我们可以使用一些巧妙的变形。

\sum^{n}_{i=1}k\bmod i=\sum^{n}_{i=1}(k-i\lfloor\dfrac{k}{i}\rfloor)=nk-\sum^{n}_{i=1}i\lfloor\dfrac{k}{i}\rfloor

其中 k\bmod i=k-i\lfloor\dfrac{n}{i}\rfloor 是一个很重要的性质，化简式子的时候经常会用！

好的，现在就是一个整除分块的标准形式了，只不过多乘了一个 i，所以连续段内的式子要重写：

\sum^{n}_{i=1}i\lfloor\dfrac{k}{i}\rfloor=\sum\dfrac{(d_{i+1}-d_i)(d_{i+1}+d_i-1)\lfloor\dfrac{k}{d_i}\rfloor}{2}

这个变形就是等差数列求和公式。

所以代码就很好写啦~

#include <bits/stdc++.h>
#define endl '\n'
using namespace std;
int n, k;
int main() {
    ios :: sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    cin >> n >> k;
    long long ans = 1LL * n * k;
    for(int i = 1, j; i <= n; i = j) {
        if(k / i != 0) j = min(k / (k / i) + 1, n + 1);
        else j = n + 1;
        ans -= 1LL * (k / i) * (i + j - 1) * (j - i) / 2; 
    }
    cout << ans << endl;
    return 0;
}

例题 3

模积和。

这一道题的式子非常毒瘤，看似十分简单，实则暗藏五个整除分块。

首先，化简一下式子：

\begin{aligned}&\sum^{n}_{i=1}\sum^{m}_{j=1}(n\bmod i)\times(m\bmod j),i\neq j\\=&\sum^{n}_{i=1}(n\bmod i)\times\sum^{m}_{j=1}(m\bmod j)-\sum^{\min(n,m)}_{i=1}(n\bmod i)\times(m\bmod i)\end{aligned}

然后对于前半部分和后半部分分别处理。对于前半部分我们使用对于取模的经典套路，可以得到这个形式：

\begin{aligned}&\sum^{n}_{i=1}(n\bmod i)\times\sum^{m}_{j=1}(m\bmod j)\\=&(n^2-\sum^{n}_{i=1}i\lfloor\dfrac{n}{i}\rfloor)(m^2-\sum^{m}_{i=1}i\lfloor\dfrac{m}{i}\rfloor)\end{aligned}

是 2 个整除分块。对于后半部分，我们同样使用套路，然后拆开式子，容易得到以下的形式：

\begin{aligned}&\sum^{\min(n,m)}_{i=1}(n\bmod i)\times(m\bmod i)\\=&\sum^{\min(n,m)}_{i=1}(n-i\lfloor\dfrac{n}{i}\rfloor)(m-i\lfloor\dfrac{m}{i}\rfloor)\\=&\sum^{\min(n,m)}_{i=1}nm-mi\lfloor\dfrac{n}{i}\rfloor-ni\lfloor\dfrac{m}{i}\rfloor+i^2\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{i}\rfloor\\=&nm\times\min(n,m)-m\sum^{\min(n,m)}_{i=1}i\lfloor\dfrac{n}{i}\rfloor-\sum^{\min(n,m)}_{i=1}i\lfloor\dfrac{m}{i}\rfloor+\sum^{\min(n,m)}_{i=1}i^2\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{i}\rfloor\end{aligned}

有三个整除分块。

不对啊，最后一个和式有两个向下取整啊？

这个东西叫二维整除分块。当求和里面有两个向下取整的时候，我们可以注意到决策点的数量还是 O(\sqrt{n}+\sqrt{m}) 的，因为决策点必定是之前两个向下取整的决策点。因此我们还是可以枚举决策点，使用整除分块的经典思路来做。

下一个决策点是 \min(\lfloor\dfrac{n}{\lfloor\dfrac{n}{i}\rfloor}\rfloor,\lfloor\dfrac{m}{\lfloor\dfrac{m}{i}\rfloor}\rfloor)+1，也就是两个决策点中较小的一个。容易证明在得到的下一个决策点和当前决策点之间，\lfloor\dfrac{n}{i}\rfloor 和 \lfloor\dfrac{m}{i}\rfloor 都是相同的。

接着用二维等差数列的求和公式：

i^2+(i+1)^2+\ldots+j^2=f(j)-f(i-1),f(x)=\dfrac{x(x+1)(2x+1)}{6}

就可以很方便的得到答案了。代码实现有一点点复杂，注意我代码里的 j 表示的是当前连续段的右端点，而不是决策点：

#include <bits/stdc++.h>
#define int __int128
#define endl '\n' 
using namespace std;
const int MOD = 19940417;
const int inv2 = 9970209;
const int inv6 = 3323403;
int n, m;
inline int read() {
    int f = 1, x = 0;
    char ch = getchar();
    while(!isdigit(ch)) {if(ch == '-') f = -f; ch = getchar();}
    while(isdigit(ch)) {x = x * 10 + ch - '0'; ch = getchar();}
    return f * x;
}
int get(int x) {
    return x * (x + 1) * (2 * x + 1) * inv6 % MOD; 
}
inline void write(int x) {
    char buf[20];
    int t = 0;
    if(x < 0) { putchar('-'); x = -x; }
    if(x == 0) { putchar('0'); return; }
    while(x > 0) { buf[++t] = x % 10 + '0'; x /= 10; }
    for(int i = t; i > 0; i--) putchar(buf[i]);
} 
signed main() {
    n = read(), m = read(); 
    int a1 = 0, a2 = 0, a3 = 0, a4 = 0, a5 = 0;
    for(int i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        a1 = (a1 + (j - i + 1) * (j + i) * inv2 * (n / i)) % MOD; 
    }
    a1 = (n * n - a1) % MOD;
    for(int i = 1, j; i <= m; i = j + 1) {
        j = m / (m / i);
        a2 = (a2 + (j - i + 1) * (j + i) * inv2 * (m / i)) % MOD;
    }
    a2 = (m * m - a2) % MOD;
    for(int i = 1, j; i <= min(n, m); i = j + 1) {
        j = min(n / (n / i), min(n, m));
        a3 = (a3 + (j - i + 1) * (j + i) * inv2 * (n / i)) % MOD;
    }
    a3 = a3 * m % MOD;
    for(int i = 1, j; i <= min(n, m); i = j + 1) {
        j = min(m / (m / i), min(n, m));
        a4 = (a4 + (j - i + 1) * (j + i) * inv2 * (m / i)) % MOD;
    }
    a4 = a4 * n % MOD;
    for(int i = 1, j; i <= min(n, m); i = j + 1) {
        j = min(min(n / (n / i), m / (m / i)), min(n, m));
        a5 = (a5 + (get(j) - get(i - 1)) * (n / i) * (m / i)) % MOD; 
    }
    write((a1 * a2 - (min(n, m) * n * m - a3 - a4 + a5) % MOD + MOD) % MOD);
    return 0;
}