Mo's algorithm
(misc/mo.hpp)

• View this file on GitHub
• Last update: 2021-03-26 19:14:01+09:00
• Include: #include "misc/mo.hpp"

Mo’s algorithm

$[0,N)$上の区間に対する$Q$個のクエリを$\mathrm{O}(\alpha N \sqrt{Q})$で計算するライブラリ。(ただし区間の幅を$1$変化させる伸縮が$\mathrm{O}(\alpha)$で行えるとする。)

概要

アルゴリズムはsnukeさんのブログが詳しい。

考察ミスしていて定数倍を大きく悪くしていたので、反省を込めて定数倍について考察したものを簡単に記しておく。

区間の伸縮の計算量が(定数倍込みで)等しい時、アルゴリズムの実行時間は伸縮の回数に依存するので、伸縮の回数を最小化すれば実行時間が高速になると考えられる。そこで、ランダムなクエリに対して区間の伸縮回数を最小化する区間の分割幅$W$の値を考察する。左端と右端が動く回数は、

• 左端：1回あたり区間$[0,W]$上のランダムな$2$点の距離と等しいので$\frac{W}{3} \times Q$回
• 右端：$i$番目の区間を処理するとき$N - Wi$個動くので$\sum_{i \leq \frac{N}{W}} N-Wi = \frac{N^2}{2W}+\mathrm{O}(\frac{N}{W})$

であるので、$\frac{N^2}{2W} + \frac{WQ}{3}$が最小になる$W$を選べばよく、AM-BMより$W=\frac{\sqrt{3}N}{\sqrt{2Q}}$を得る。

(修正前のライブラリはうっかり$W=\frac{\sqrt{3}N}{\sqrt{Q}}$にしていたがこれだと定数倍が$3$倍ほど悪いようだ。これはひどい…)

Required by

• 二項係数のprefix sumの多点評価 (modulo/multipoint-binomial-sum.hpp)

Verified with

• verify/verify-unit-test/multipoint-binomial-sum.test.cpp
• verify/verify-yosupo-ds/yosupo-static-range-inversions-query.test.cpp

Code

#pragma once

struct Mo {
  int width;
  vector<int> left, right, order;

  Mo(int N, int Q) : order(Q) {
    width = max<int>(1, 1.0 * N / max<double>(1.0, sqrt(Q * 2.0 / 3.0)));
    iota(begin(order), end(order), 0);
  }

  void insert(int l, int r) { /* [l, r) */
    left.emplace_back(l);
    right.emplace_back(r);
  }

  template <typename AL, typename AR, typename DL, typename DR, typename REM>
  void run(const AL &add_left, const AR &add_right, const DL &delete_left,
           const DR &delete_right, const REM &rem) {
    assert(left.size() == order.size());
    sort(begin(order), end(order), [&](int a, int b) {
      int ablock = left[a] / width, bblock = left[b] / width;
      if (ablock != bblock) return ablock < bblock;
      if (ablock & 1) return right[a] < right[b];
      return right[a] > right[b];
    });
    int nl = 0, nr = 0;
    for (auto idx : order) {
      while (nl > left[idx]) add_left(--nl);
      while (nr < right[idx]) add_right(nr++);
      while (nl < left[idx]) delete_left(nl++);
      while (nr > right[idx]) delete_right(--nr);
      rem(idx);
    }
  }
};

/**
 * @brief Mo's algorithm
 * @docs docs/misc/mo.md
 */

#line 2 "misc/mo.hpp"

struct Mo {
  int width;
  vector<int> left, right, order;

  Mo(int N, int Q) : order(Q) {
    width = max<int>(1, 1.0 * N / max<double>(1.0, sqrt(Q * 2.0 / 3.0)));
    iota(begin(order), end(order), 0);
  }

  void insert(int l, int r) { /* [l, r) */
    left.emplace_back(l);
    right.emplace_back(r);
  }

  template <typename AL, typename AR, typename DL, typename DR, typename REM>
  void run(const AL &add_left, const AR &add_right, const DL &delete_left,
           const DR &delete_right, const REM &rem) {
    assert(left.size() == order.size());
    sort(begin(order), end(order), [&](int a, int b) {
      int ablock = left[a] / width, bblock = left[b] / width;
      if (ablock != bblock) return ablock < bblock;
      if (ablock & 1) return right[a] < right[b];
      return right[a] > right[b];
    });
    int nl = 0, nr = 0;
    for (auto idx : order) {
      while (nl > left[idx]) add_left(--nl);
      while (nr < right[idx]) add_right(nr++);
      while (nl < left[idx]) delete_left(nl++);
      while (nr > right[idx]) delete_right(--nr);
      rem(idx);
    }
  }
};

/**
 * @brief Mo's algorithm
 * @docs docs/misc/mo.md
 */

Back to top page