骗分带师: 莫队Ⅰ

比利♂海灵顿 · 2021-05-12 21:44:54 · 个人记录

莫队算法入门

离线算法, 解决多次区间询问问题. 算法的条件是一个问题区间 [L, R] 询问的答案可以由 [L - 1, R], [L, R - 1], [L + 1, R], [L, R + 1] 中的任意某个区间的答案以 O(1) 的复杂度求出. 通过将区间排序, 然后暴力地用上一个答案求下一个答案. 复杂度为 O(n\sqrt n).

结合莫队的经典模板题小Z的袜子来讲解.

题面简述

一行 n 只袜子, 求区间 [L, R] 中随机挑出两只后, 这两只袜子颜色相同的几率, m 次询问.

算法框架

开一个数组 Cnt 存储当前区间每种袜子数量, Cnt_i 表示当前区间第 i 个颜色的袜子个数.

开一个数组 Ans 存每个询问的答案, 用于按询问的输入顺序输出答案. 由于要输出最简分数, 所以每个答案是两个整数, 分别代表分子和分母.

分块, 以 L 的块数为第一关键字, R 为第二关键字排序. 依次枚举所有的询问, 暴力地从上一个答案推得本次答案.

增量

增量就是用相邻区间 (即只有一个端点相差为 1, 另一个端点相同的两个区间) 的答案求目标区间答案的操作. 莫队算法复杂度的正确性要求增量操作复杂度为 O(1).

设第 i 只袜子颜色为 a_i, 区间长度为 Len.

暴力维护 Cnt 数组, 用 [L - 1, R] 的 Cnt 求 [L, R] 的 Cnt 时, 需要将 Cnt_{a_i} 减少 1. 复杂度 O(1). 已知 [L \pm 1, R] 或 [L, R \pm 1] 求 [L, R] 的情况同理.

保证了 Cnt 数组的 up to date, 接下来求概率. 考虑某种颜色的袜子被选中一双的概率. 第一只被选中的袜子颜色为 Col 的几率为 \frac{Cnt_{Col}}{Len}, 第二只仍为 Col 的几率为

但是如果每个颜色单独算, 单次查询时间复杂度 O(n), 一定会超时. 所以考虑一起计算. 如果一次增量中, 区间长度从 Len' 变成了 Len (当然, Len' = Len \pm 1). Cnt 没有变化的某种颜色 Col 的袜子被选中一双的概率从 \frac{{Cnt_{Col}}^2 - Cnt_{Col}}{{Len'}^2 - Len'} 变成 \frac{{Cnt_{Col}}^2 - Cnt_{Col}}{{Len}^2 - Len}.

所以对于每个 Cnt 没有变化的颜色, 被选中一双的概率都会变成原来的 \frac{{Len'}^2 - Len'}{{Len}^2 - Len}.

那么对于被加入或删除了一个袜子的颜色 Col, 选中一双颜色为 Col 的袜子的概率从原来的 \frac{{{Cnt_{Col}}'}^2 - {Cnt_{Col}}'}{{Len'}^2 - Len'} 变成 \frac{{Cnt_{Col}}^2 - Cnt_{Col}}{{Len}^2 - Len}.

因此, 设原答案为 \frac{Ans0'}{Ans1'} 新的答案 \frac{Ans0}{Ans1} 就是:

分类讨论, 推式子

• 当区间增长了 1, 即 Len = Len' + 1

• 当区间缩短了 1, 即 Len = Len' - 1

考虑约分的操作. 因为分数的计算需要大量的乘积, 所以需要随时约分. 答案都是真分数, 分母最多是 Len^2 - Len, 分子不大于分母, 所以答案的分子分母不超过 2.5*10^{9}, 不超过 unsigned 的范围.

使用 unsigned long long 进行中间计算, 在每次计算后约分. 使用欧几里得算法求分子分母的 gcd, 然后分子分母一起除以这个 gcd.

特别地, 在 Len = 1 的情况下, Len^2 - Len = 0, 分母为零无意义, 也存在概率为零的情况, 样例表明, 概率为 0 时, 输出 0/1.

当概率为 1 时, 答案为整数, 同样有样例表明, 这时输出 1/1.

但是这样的复杂度外面就会乘上一个 O(log(n^2)), 变成 O(n\sqrt n log(n^2)), 因为每次增量要求 GCD. 所以考虑更优的解法. 考虑从一开始单个颜色袜子选中一双的概率入手, 推得一个式子, 整理.

这样, 只要维护区间袜子数量平方和即可 O(1) 求答案. 由于平方和, 区间长度的平方不大于 2.5 * 10^9, 所以无需每次增量约分, 只要约分 O(m) 次即可.

维护平方和, 对于某个 Cnt 的变化对平方和的影响, 如果 Cnt 在原来的基础上增加了 1, 即 Cnt = Cnt' + 1, 则

对于 Cnt = Cnt' - 1 的情况, 有

复杂度证明

因为莫队的前提是一次增量时间复杂度 O(1), 所以只要证明莫队对长度为 n 的序列的 m 次区间查询的复杂度是 O(n \sqrt m).

在分块的块长选择 \frac{n}{\sqrt m} 的前提下, 块数就是 \sqrt{m}. 使用 algorithm 中自带的排序, 复杂度是 O(logm).

接下来对当前答案进行增量, 对于左端点和上一个询问同块的情况, 左端点最多增量 \frac{n}{\sqrt m} 次. 对于左端点在新的块的询问, 最多增量 \frac{2n}{\sqrt m} 次, 也是 O(\frac{n}{\sqrt m}). 所有左端点增量次数复杂度为 O(\frac{nm}{\sqrt m}) = O(n \sqrt m). 左端点这个块内的所有询问, 右端点递增, 所以这些询问右端点总共增量最多 n 次. 所以右端点总共的增量数应该是 O(n \sqrt m). 综上, 左右端点增量数复杂度为 O(n \sqrt m)

分析更一般的情况, 如果块长选 x, 块数就是 \frac {n}{x}. 用上面的方式分析, 左端点的增量数复杂度是 O(xm), 右端点的增量数复杂度是 O(\frac{n^2}{x}), 总复杂度为 O(xm + \frac{n^2}{x}). 使用均值不等式, 最优复杂度为 xm + \frac{n^2}{x} \geq 2\sqrt{mn^2} = 2n\sqrt m = O(n\sqrt{m}), 当 xm = \frac{n^2}{x} 时取到等号. 整理求块长:

所以, 一般莫队的最优块长是 \frac{n}{\sqrt m}

代码实现

细节不多, 码量短小, 把推的式子写成代码就好了, 是骗分的好工具.

unsigned m, n, Cnt[50005], BlockLen, BlockCnt;
long long a[50005], Ans[50005][2], Tmp0(0), Tmp1(1), TmpG, TmpSquare(1);
struct Query{
  unsigned L, R, Num, BelongToBlocks;
  inline const char operator <(const Query &x) {  // 按左端点所在块排序 
    return (this->BelongToBlocks ^ x.BelongToBlocks) ? this->BelongToBlocks < x.BelongToBlocks : this->R < x.R;
  }
}Q[50005];
inline unsigned GCD(register unsigned x, register unsigned y) {
  register unsigned tmp;
  while(y) tmp = x, x = y, y = tmp % y;
  return x;
}
int main() {
  n = RD(), m = RD();
  BlockLen = (m ^ 0) ? (n / sqrt(m)) + 1 : sqrt(n) + 1;
  for (register unsigned i(1); i <= n; ++i) {
    a[i] = RD();
  }
  for (register unsigned i(1); i <= m; ++i) {
    Q[i].L = RD(), Q[i].R = RD(), Q[i].Num = i, Q[i].BelongToBlocks = (Q[i].L + BlockLen - 1) / BlockLen;
  }
  sort(Q + 1, Q + m + 1);
  Tmp0 = 0, Tmp1 = 1, Q[0].L = 1, Q[0].R = 1, Cnt[a[1]] = 1; // 初始化当前区间为 [1, 1] 
  for (register unsigned i(1); i <= m; ++i) {
    if(Q[i].L == Q[i].R) {      // 特判, 单点查询 
      Ans[Q[i].Num][0] = 0, Ans[Q[i].Num][1] = 1;
      continue;
    }
    register unsigned Col, Len(Q[i].R - Q[i].L + 1);
    while (Q[0].L > Q[i].L) {   // 左端点左移 
      ++Cnt[Col = a[--Q[0].L]];
      TmpSquare += (Cnt[Col] << 1) - 1; // Cnt[Col] 增加, 维护方差, 下同 
    }
    while (Q[0].R < Q[i].R) {   // 右端点右移 
      ++Cnt[Col = a[++Q[0].R]];
      TmpSquare += (Cnt[Col] << 1) - 1;
    }
    while (Q[0].L < Q[i].L) {   // 左端点右移
      --Cnt[Col = a[Q[0].L++]];
      TmpSquare -= (Cnt[Col] << 1) + 1;
    }
    while (Q[0].R > Q[i].R) {   // 右端点左移
      --Cnt[Col = a[Q[0].R--]];
      TmpSquare -= (Cnt[Col] << 1) + 1;
    }
    Ans[Q[i].Num][0] = TmpSquare - Len;
    Ans[Q[i].Num][1] = Len * Len - Len;
    TmpG = GCD(Ans[Q[i].Num][0], Ans[Q[i].Num][1]);
    Ans[Q[i].Num][0] /= TmpG;
    Ans[Q[i].Num][1] /= TmpG;
  }
  for (register unsigned i(1); i <= m; ++i) {
    printf("%u/%u\n", Ans[i][0], Ans[i][1]);
  }
  return Wild_Donkey;
}