【算法】莫队算法粗略讲解

原文链接：https://cnblogs.com/ctjcalc/p/post8.html

莫队是由莫涛大佬研究出的一种非常优秀的离线算法，主要用来解决区间问题，甚至可以在非强制在线的情况下离线处理包括修改的操作。本文主要讲解普通莫队算法。

Copyright (C) 2020 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

前置知识

• 分块思想
• 熟练掌握 STL 算法

莫队思想

先来看一道例子，给出下面的一个序列，给出一个区间，求区间和：

下标	1	2	3	4	5	6
\(A\)	2	1	1	3	4	3

如果这道题不涉及修改操作，大家都会想到预处理前缀和。现在我强制要求使用莫队，要怎么做呢？假设我们要求 \(S_{[2,4]}\) 的和，并且知道了 \(S_{[2,3]}\)，可以想到用 \(S_{[2,3]}\)加上 \(A_4\)，这样就能够得到答案了。同理，如果我们知道了 \(S_{[3,4]}\)，\(S_{[2,4]} = S_{[3,4]} + A_2\)。这些操作实际上就是移动区间的一个端点，同时把新的端点的答案合并。莫队就是这样，如果已经知道了一个区间的答案，就试图通过移动区间边界，把原来的答案与新的边界的值合并，最终让当前区间与询问的区间重合，然后记录答案。

但是，现在又有一个问题了，如果有一个序列，长度为 \(n\)，询问 \(m\) 次，按照 \([1,2],[n-1,n],[3,4],[n-3,n-2],\cdots\) 询问，那时间复杂度就变成了 \(O(nm)\) 级别了，显然不利于解题。于是莫涛大佬又想出了一个解决方案，把序列分块，再把询问排序，按照排序后的询问区间处理答案，对于两个询问区间，如果两个区间的左端点在同一个块，就比较它们的右端点，否则比较左端点。对于上面的那些询问，排序后，时间复杂度减小至 \(\Theta(n+m)\)。对于一般情况下，这样的时间复杂度是 \(O(n\sqrt{n})\)。但我不会证明 QwQ。

Copyright (C) 2020 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

代码模板

通过上面的介绍，我们来系统地总结一下移动区间端点的 \(4\) 种情况（上面只举例了两种）。记当前区间为 \([l,r]\)，询问区间为 \([L,R]\)，用代码说明。

情况	代码实现
\(L < l\)	while (L < l) add(--l);
\(R > r\)	while (R > r) add(++r);
\(L > l\)	while (L > l) sub(l++);
\(R < r\)	while (R < r) sub(r--);

前两种情况是扩大区间，就把新端点的答案加到总答案里，后两种相反，把旧端点的答案从总答案里删除。

接下来，看看莫队算法的框架。

// ...
struct query {
    int l, r, id; // 询问区间以及它是第几个询问
};
 
constexpr int maxn = /* ... */; // 序列长度
constexpr int maxq = /* ... */; // 询问个数
query qs[maxq];
int arr[maxn], ans[maxq], n, q, res, blocksize; // res 表示临时计算的答案，会随着区间的移动不断更新
 
inline int blockid(int x) { return (x - 1) / blocksize + 1; }
 
inline void add(int x) {
    // ...
}
 
inline void sub(int x) {
    // ...
}
 
void solve() {
    // ...
    blocksize = sqrt(n); // 计算块的大小
    sort(qs + 1, qs + 1 + q, [](const query &a, const query &b) { // 排序区间
        return blockid(a.l) == blockid(b.l) ？ a.r < b.r : a.l < b.l;
    }); // 这里使用了 Lambda 表达式
    int l = 1, r = 0; // 一开始区间是空的，这样写是为了符合语义
    for (int i = 1; i <= q; ++i) { // 处理询问
        while (qs[i].l < l)
            add(--l);
        while (qs[i].r > r)
            add(++r);
        while (qs[i].l > l)
            sub(l++);
        while (qs[i].r < r)
            sub(r--);
        ans[qs[i].id] = res; // 记录答案
    }
    // ...
}

Copyright (C) 2020 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

算法优化

排序后的区间顺序其实对处理的速度有较大的影响，实际上，我们一般采用一个叫做奇偶化排序的东西。如果块的编号是奇数，就按右端点升序排序，否则按降序排序。它的原理大致是——如果处理完奇数块，\(r\) 这个指针就可以不用再跑更远，从前开始往后扫描，可以看看下面这个询问的例子：

// blocksize = 3
1 1
2 25
3 50
4 1
5 25
6 50

用了奇偶化排序后，速度可以提升大约 \(30\%\)。

代码实现：

sort(qs + 1, qs + 1 + q, [](const query &a, const query &b) {
    return blockid(a.l) == blockid(b.l)
               ? (a.r == b.r ? 0 : !((blockid(a.l) & 1) && (a.r < b.r)))
               : a.l < b.l;
});

Copyright (C) 2020 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

例题讲解

SPOJ 3267 D-query

给出一个长度为 \(n\) 的序列，进行 \(q\) 次询问，每次询问给出区间 \([l,r]\)，问区间去重后有多少个数？

唯一需要注意的就是多维护一个数组，记为 cnt[i]，表示第值为 \(i\) 的数目前出现了多少次。很显然，直接在 add(x) 和 sub(x) 里更新，看代码。

#include <bits/stdc++.h>
using namespace std;
 
template <typename T> T read() {
    T x = 0, s = 1;
    char c = getchar();
    while (c < '0' || '9' < c) {
        if (c == '-')
            s = -1;
        c = getchar();
    }
    while ('0' <= c && c <= '9') {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }
    return x * s;
}
 
struct query {
    int l, r, id;
};
 
constexpr int maxn = 30000 + 5;
constexpr int maxv = 1e6 + 5;
constexpr int maxq = 200000 + 5;
query qs[maxq];
int arr[maxn], cnt[maxv], ans[maxq], n, q, res, blocksize;
 
inline int blockid(int x) { return (x - 1) / blocksize + 1; }
 
inline void add(int x) {
    ++cnt[arr[x]];
    if (cnt[arr[x]] == 1)
        ++res;
}
 
inline void sub(int x) {
    --cnt[arr[x]];
    if (!cnt[arr[x]])
        --res;
}
 
int main() {
#ifndef ONLINE_JUDGE
    freopen("Environment/project.in", "r", stdin);
    freopen("Environment/project.out", "w", stdout);
#endif
    n = read<int>();
    blocksize = sqrt(n);
    for (int i = 1; i <= n; ++i)
        arr[i] = read<int>();
    q = read<int>();
    for (int i = 1; i <= q; ++i) {
        qs[i].l = read<int>();
        qs[i].r = read<int>();
        qs[i].id = i;
    }
    sort(qs + 1, qs + 1 + q, [](const query &a, const query &b) {
        return blockid(a.l) == blockid(b.l)
                   ? (a.r == b.r ? 0 : (blockid(a.l) & 1) ^ (a.r < b.r))
                   : a.l < b.l;
    });
    int l = 1, r = 0;
    for (int i = 1; i <= q; ++i) {
        while (qs[i].l < l)
            add(--l);
        while (qs[i].r > r)
            add(++r);
        while (qs[i].l > l)
            sub(l++);
        while (qs[i].r < r)
            sub(r--);
        ans[qs[i].id] = res;
    }
    for (int i = 1; i <= q; ++i)
        printf("%d\n", ans[i]);
    return 0;
}

Copyright (C) 2020 ctjcalc，转载请注明URL，并给出原文链接，谢谢。

Luogu P1494 小Z的袜子

有一个长度为 \(n\) 的序列，进行 \(q\) 次询问，每次询问给出区间 \([l,r]\)，问随机从区间选两个数相等的概率是多少？使用最简分数输出。

同样是要维护 cnt[i]。对于区间 \([l,r]\)，选取两个数的方案数是 \(\binom{r-l+1}{2}\)，而选取到相同两个数的方案数是 \(\sum_{i=1}^{N}\binom{cnt[i]}{2}\)。概率就是它们相除。我们可以维护一下后面的这个方案数，前一个直接计算。每次移动区间端点都是差不多的，以扩大区间为例，要先减去之前的答案，再加上 cnt[i] + 1 的答案，那么加在一起就是 \(\binom{cnt[i]+1}{2}-\binom{cnt[i]}{2}\)。这个式子是可以化简的：

\[\begin{aligned}\binom{cnt[i]+1}{2}-\binom{cnt[i]}{2}&=\frac{(cnt[i]+1)!}{2!(cnt[i]-1)!}-\frac{cnt[i]}{2!(cnt[i]-2)}\\&=\frac{cnt[i](cnt[i]+1)}{2}-\frac{cnt[i](cnt[i]-1)}{2}\\&=\frac{2cnt[i]}{2}\\&=cnt[i]\end{aligned} \]

缩小区间也是一样的，你可以再推导一遍，也可以先把 cnt[i] 减去 \(1\)，然后让临时答案减去 cnt[i]，可以想想为什么。

#include <bits/stdc++.h>
using namespace std;
 
template <typename T> T read() {
    T x = 0, s = 1;
    char c = getchar();
    while (c < '0' || '9' < c) {
        if (c == '-')
            s = -1;
        c = getchar();
    }
    while ('0' <= c && c <= '9') {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }
    return x * s;
}
 
struct query {
    int l, r, id;
};
 
constexpr int maxn = 50000 + 5;
query qs[maxn];
int arr[maxn], cnt[maxn], n, q, blocksize;
long long ans[maxn][2], res;
 
inline int blockid(int x) { return (x - 1) / blocksize + 1; }
 
inline void add(int x) { res += cnt[arr[x]]++; }
 
inline void sub(int x) { res -= --cnt[arr[x]]; }
 
long long GCD(long long x, long long y) { return y == 0 ? x : GCD(y, x % y); }
 
int main() {
#ifndef ONLINE_JUDGE
    freopen("Environment/project.in", "r", stdin);
    freopen("Environment/project.out", "w", stdout);
#endif
    n = read<int>();
    q = read<int>();
    blocksize = sqrt(n);
    for (int i = 1; i <= n; ++i)
        arr[i] = read<int>();
    for (int i = 1; i <= q; ++i) {
        qs[i].l = read<int>();
        qs[i].r = read<int>();
        qs[i].id = i;
    }
    sort(qs + 1, qs + 1 + q, [](const query &a, const query &b) {
        return blockid(a.l) == blockid(b.l)
                   ? (a.r == b.r ? 0 : (blockid(a.l) & 1) ^ (a.r < b.r))
                   : a.l < b.l;
    });
    int l = 1, r = 0;
    for (int i = 1; i <= q; ++i) {
        if (qs[i].l == qs[i].r) {
            ans[qs[i].id][0] = 0;
            ans[qs[i].id][1] = 1;
            continue;
        }
        while (qs[i].l < l)
            add(--l);
        while (qs[i].r > r)
            add(++r);
        while (qs[i].l > l)
            sub(l++);
        while (qs[i].r < r)
            sub(r--);
        ans[qs[i].id][0] = res;
        ans[qs[i].id][1] = (1LL * (r - l) * (r - l + 1)) >> 1;
    }
    for (int i = 1; i <= q; ++i) {
        long long g = GCD(ans[i][0], ans[i][1]);
        printf("%lld/%lld\n", ans[i][0] / g, ans[i][1] / g);
    }
    return 0;
}

Copyright (C) 2020 ctjcalc，转载请注明URL，并给出原文链接，谢谢。