跳转至

数论分块

一个简单且有趣的小技巧。

我以前似乎想到过类似的东西，并且过了一个题，但是我忘记是啥了。

Easyψ(｀∇´)ψ

给定 \(n, k\)，求 \(\sum\limits_{i = 1}^n \lfloor\dfrac{k}{i}\rfloor\) 的值。

\(1\le n, k \le 10^5\)。

有一个观察是，\(\lfloor\dfrac{k}{i}\rfloor\) 这个东西其实是一块一块的。

比如 \(k = 10 \to 10, 5, 3, 2, 2, 1, 1, 1, 1, 1\)。

所以我们有一个想法是考虑算出每一段的长度，那么复杂度就可以降到 \(O(B)\)，其中 \(B\) 是段数。

我们考虑能整除 \(k\) 的数在 \(1 \sim k\) 当中有多少个，这个是 \(O(\sqrt k)\) 级别的，如果在 \(i > k\) 的部分，那都是 \(0\) 了。

所以 \(O(B) = O(\sqrt k)\)，现在要做的就是，怎么分块，换句话说怎么快速知道每一段是哪一些。

结论：如果 \(i < j\)，并且 \(\lfloor\dfrac{k}{i}\rfloor = \lfloor\dfrac{k}{j}\rfloor\)，那么 \(j \in \mathbb{Z}\) 的 \(\max j = \left\lfloor\dfrac{k}{\lfloor\dfrac{k}{i}\rfloor}\right\rfloor\)。

结论的证明很简单：

首先 \(\lfloor\dfrac{x}{y}\rfloor \le \dfrac{x}{y}\) 这个是显然的，然后我们带入上面的条件：

\[ \lfloor\dfrac{k}{i}\rfloor =\lfloor\dfrac{k}{j}\rfloor \le \dfrac{k}{j} \]

所以：

\[ \lfloor\dfrac{k}{i}\rfloor \le \dfrac{k}{j} \]

有：

\[ j \le \dfrac{k}{\lfloor\dfrac{k}{i}\rfloor} \]

然后因为 \(j \in \mathbb{Z}\) 所以：

\[ j \le \left\lfloor\dfrac{k}{\lfloor\dfrac{k}{i}\rfloor}\right\rfloor \]

这里证明了充分性，必要性是显然的。

那么每次找到段首 \(l\)，令 \(r =\left\lfloor\dfrac{k}{\lfloor\dfrac{k}{l}\rfloor}\right\rfloor\) 即可。

Code 
1
2
3
4
5
6
7
8
9
int calc(int n, int k) {
    int l = 1, r = 0, ret = 0;
    for(; l <= n && k / l; ) {
        r = (k / (k / l)), r = min(n, r);
        ret += (r - l + 1) * (k / l);
        l = r + 1;
    }
    return ret;
}

注意代码里要判 \(i > k\) 和 \(i > n\) 的情况。

Normalψ(｀∇´)ψ

如果求 \(G(n, k) = \sum\limits_{i = 1}^n i\lfloor\dfrac{k}{i}\rfloor\) 呢？

这里是个 Trick，注意到分段之后，段中每一个位置的 \(\lfloor\dfrac{k}{i}\rfloor\) 是相等的。

这部分可以乘法结合律，所以我们只需要求一个 \(\sum\limits_{i = l}^r i\) 即可。

如果是一个函数 \(f(x)\)，那么我们要多求的就是 \(\sum\limits_{i = l}^r f(i)\)。

code 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
int sum(int l, int r) { return (r - l + 1) * (r + l) / 2ll; }
int calc(int n, int k) {
    int l = 1, r = 0, ret = 0;
    for(; l <= n && k / l; ) {
        r = (k / (k / l)), r = min(n, r);
        ret += sum(l, r) * (k / l);
        l = r + 1;
    }
    return ret;
}

Hardψ(｀∇´)ψ

Source

求 \(F(n, k) = \sum\limits_{i = 1}^n (k\bmod i)\)。

数据范围不变。

拆成 \(k\bmod i = k - i\lfloor\dfrac{k}{i}\rfloor\) 即可。

Lunaticψ(｀∇´)ψ

Source

求 \(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m (n\bmod i) \cdot (m \bmod j), i \not= j\)。

不慌，我们先拆一下：

\[ \begin{aligned} &= \sum\limits_{i = 1}^n\sum\limits_{j = 1}^m (n - i\lfloor\dfrac{n}{i}\rfloor) \cdot (m - j\lfloor\dfrac{m}{j}\rfloor), i \not= j \\ &= \sum\limits_{i = 1}^n\sum\limits_{j = 1}^m (n - i\lfloor\dfrac{n}{i}\rfloor) \cdot (m - j\lfloor\dfrac{m}{j}\rfloor) - \sum\limits_{i = 1}^n (n - i \lfloor\dfrac{n}{i}\rfloor)\cdot(m - i\lfloor\dfrac{m}{i}\rfloor) \\ &= F(n, n)\cdot F(m, m) - n^2m + mG(n, n) + nG(n, m) - \sum\limits_{i = 1}^n i^2 \cdot\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{i}\rfloor \end{aligned} \]

\(F, G\) 就是前面说过的，直接做就行。

我们考虑最后这一坨怎么办，看起来像一个二维的东西，可以发现 \(\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{i}\rfloor\) 的值也是分块的。

你可以理解成，有两个楼梯接在了一起，然后就像 minecraft 里面的沙子一样（我觉得 cmd 大爷的这个比喻很形象就搬过来了）。

然后图形化理解一下，可以发现 \(r = \min(\left\lfloor\dfrac{n}{\lfloor\dfrac{n}{l}\rfloor}\right\rfloor, \left\lfloor\dfrac{m}{\lfloor\dfrac{m}{l}\rfloor}\right\rfloor)\)，之后做起来就很方便了。

然后乘上的 \(f(i) =\sum\limits_{j = 1}^ij^2\)，这个用高中数学教的数列知识推一下可以发现 \(f(i) = \dfrac{i(i + 1)(2i + 1)}{6}\)。

Code

代码里的 F, G 和描述里不太一样，注意一下就好了。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
// author : black_trees

#include <cmath>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <iostream>
#include <algorithm>

#define endl '\n'
#define int __int128_t

using namespace std;
using i64 = long long;

template <typename __Tp> void read(__Tp &x) {
    int f = x = 0; char ch = getchar();
    for (; !isdigit(ch); ch = getchar()) if (ch == '-') f = 1;
    for (; isdigit(ch); ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ 48);
    if (f) x = -x;
}
void read(char &ch) { for (ch = getchar(); isspace(ch); ch = getchar()); }
template <typename __Tp1, typename ...__Tp2> void read(__Tp1 &x, __Tp2 &... y) { read(x), read(y...); }
template <typename __Tp> void write(__Tp x) {
    if (x < 0) putchar('-'), x = -x;
    if (x > 9) write(x / 10);
    putchar(x % 10 + 48);
}
void write(char ch) { putchar(ch); }
void write(const char *s) { for (; *s; ++s) putchar(*s); }
template <typename __Tp1, typename ...__Tp2> void write(__Tp1 x, __Tp2 ... y) { write(x), write(y...); }

const int mod = 19940417ll;
const int inv2 = 9970209ll;
const int inv6 = 3323403ll;

int n, m;
int Sum1(int l, int r) { return (r - l + 1 + mod) % mod * (r + l) % mod * inv2 % mod; }
int sum(int n) { return (n % mod) * ((n + 1) % mod) % mod * ((2ll * n + 1) % mod) % mod * inv6 % mod; }
int Sum2(int l, int r) { return (sum(r) - sum(l - 1) + mod) % mod; } 
int F(int n, int k) {
    int ret = 0;
    int l = 1, r = 0;
    for(; l <= n && k / l; ) {
        r = (k / (k / l)), r = min(r, n);
        (ret += Sum1(l, r) * (k / l) % mod) %= mod;
        l = r + 1;
    }
    return ret % mod;
}
int G(int n, int m) {
    int ret = 0;
    int l = 1, r = 0;
    for(; l <= n; ) {
        r = min((n / (n / l)), (m / (m / l)));
        (ret += Sum2(l, r) * (n / l) % mod * (m / l) % mod) %= mod;
        l = r + 1;
    }
    return ret % mod;
}

signed main() {

    cin.tie(0) -> sync_with_stdio(false);
    cin.exceptions(cin.failbit | cin.badbit);

    read(n, m);
    if(n > m) swap(n, m);

    int ans = (mod - n * n % mod * m % mod + mod) % mod;
    (ans += (n * n % mod - F(n, n) % mod) % mod * (m * m - F(m, m) % mod) % mod) %= mod;
    (ans += m * F(n, n) % mod + n * F(n, m) % mod) %= mod;
    (ans += mod - G(n, m) % mod) %= mod;

    write(ans % mod, endl);

    return 0;
}

Extraψ(｀∇´)ψ

现在你可以去学莫比乌斯反演了！

最后更新: June 18, 2023