骗分带师: 莫队Ⅰ 莫队算法入门

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

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

题面简述

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

(n, m leq 50000)

算法框架

开一个数组 (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(nsqrt 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 2sqrt{mn^2} = 2nsqrt m = O(nsqrt{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
", Ans[i][0], Ans[i][1]);
  }
  return Wild_Donkey;
}