Mo's Algorithm: Secondary Offline Processing

Overview

Sometimes we encounter problems that look well suited to Mo's algorithm, but a single transition is not O(1)O(1). Applying Mo's algorithm directly would give the wrong complexity even after adjusting the block size.

If the contribution of each transition to the answer can be expressed by differences, we can split these transitions, take them offline, and process them in batches with another algorithm.

Let f(x,l,r)f(x, l, r) denote the contribution produced by xx with respect to [l,r][l, r].

When we extend the current interval [l,r][l, r] to [l,r+1][l, r + 1], we need f(ar+1,l,r)f(a_{r + 1}, l, r). If this can be expressed by differences, we can write it as f(ar+1,1,r)−f(ar+1,1,l−1)f(a_{r + 1}, 1, r) - f(a_{r + 1}, 1, l - 1). We can precompute the first term for every rr, while the latter terms can be taken offline and stored at their corresponding l−1l - 1, then enumerated in increasing order and processed with a scanline. The other transition directions can be handled similarly.

This algorithm, which takes the offline Mo's algorithm and processes its transitions offline once more, is called Mo's algorithm with secondary offline processing.

We explain it through concrete problems.

Example problems

Luogu P5047 [Ynoi2019 Mock Contest] Yuno loves sqrt technology II

Given a sequence of length nn called aa and mm queries, find the number of inversions in each queried interval.

Data range: 1≤n,m≤1051 \leq n,m \leq 10^5, 0≤ai≤1090 \leq a_i \leq 10^9.

Directly applying Mo's algorithm requires at least O(log⁡n)O(\log n) time per transition. Observe that each transition asks for the rank of a number within an interval; this information can be expressed by interval differences, so we consider secondary offline processing.

We extend [l,r][l, r] to [l,r+1][l, r + 1], and need to count how many numbers in [l,r][l, r] are greater than ar+1a_{r + 1}. Let f(x,r)f(x, r) denote the number of elements in [1,r][1, r] greater than axa_x, and let g(x,r)g(x, r) denote the number of elements in [1,r][1, r] less than axa_x. The change in the answer can then be written as f(r+1,r)−f(r+1,l−1)f(r + 1, r) - f(r + 1, l - 1).

Similarly, when [l,r][l, r] shrinks to [l,r−1][l, r - 1], the change is −f(r,r−1)+f(r,l−1)-f(r, r - 1) + f(r, l - 1); when [l,r][l, r] expands to [l−1,r][l - 1, r], it is g(l−1,r)−g(l−1,l−2)g(l - 1, r) - g(l - 1, l - 2); and when [l,r][l, r] shrinks to [l+1,r][l + 1, r], it is −g(l,r)+g(l,l−1)- g(l, r) + g(l, l - 1).

For f(x,x−1)f(x, x - 1) and g(x,x−1)g(x, x - 1) in these expressions, we can use a Fenwick tree to precompute them in O(nlog⁡n)O(n \log n) time. For the remaining f(x,p)f(x, p) and g(x,p)g(x, p) terms, we take xx offline at ll and process them in batches.

Here is a space optimization. When processing each query, the values taken offline at pp, denoted by xx, form a contiguous segment, so we do not need to store every move; storing only the adjusted segment reduces the space from the total number of moves O(nm)O(n\sqrt{m}) to the number of queries O(m)O(m).

We now handle the problem obtained by taking these terms offline a second time: insert numbers into a set and query the rank of a number in the set. The total number of endpoint moves in Mo's algorithm is O(nm)O(n\sqrt{m}), while the array has only O(n)O(n) elements, so we consider value-domain decomposition with O(n)O(\sqrt{n})-time insertion and O(1)O(1)-time queries.

Thus, we solve the problem in O(nm+nn)O(n \sqrt{m} + n \sqrt{n}) time and O(n+m)O(n + m) space.

Finally, note that what we obtain is the change in the answer at each step, not the answer itself; a prefix sum is needed to obtain the final answer.

Example code
#include <algorithm>
#include <cmath>
#include <cstring>
#include <iostream>
#include <vector>

constexpr int MAX = 1e5 + 5;
constexpr int MAXX = 4e2 + 5;
constexpr int MX = 1e5;
typedef long long ll;

int n, m, sz, _sz;
int a[MAX], b[MAX], bl[MAX], st[MAXX], ed[MAXX], b1[MAXX], b2[MAX], p1[MAX],
    p2[MAX];
ll ans[MAX];

struct N {
  int l, r, id, x, y;
};

std::vector<N> v[MAX];

struct Q {
  int l, r, id;
  ll ans = 0;
} q[MAX];

bool cmp(Q a, Q b) {
  if (a.l / _sz != b.l / _sz) {
    return a.l < b.l;
  }
  return ((a.l / _sz) & 1) ? a.r < b.r : a.r > b.r;
}

void init() {
  sz = sqrt(MX);
  for (int i = 1; i <= sz; ++i) {
    st[i] = (MX / sz) * (i - 1) + 1;
    ed[i] = (MX / sz) * i;
  }
  ed[sz] = MX;
  for (int i = 1; i <= sz; ++i) {
    for (int j = st[i]; j <= ed[i]; ++j) {
      bl[j] = i;
    }
  }
}

void clear() {
  memset(b1, 0, sizeof(b1));
  memset(b2, 0, sizeof(b2));
}

void mdf(int x) {
  for (int i = bl[x]; i <= sz; ++i) {
    ++b1[i];
  }
  for (int i = x; i <= ed[bl[x]]; ++i) {
    ++b2[i];
  }
}

int qry(int x) {
  if (!x) {
    return 0;
  }
  return b1[bl[x] - 1] + b2[x];
}

int main() {
  scanf("%d%d", &n, &m);
  init();
  _sz = sqrt(m);
  for (int i = 1; i <= n; ++i) {
    scanf("%d", &a[i]);
    b[i] = a[i];
  }
  std::sort(b + 1, b + n + 1);
  int len = std::unique(b + 1, b + n + 1) - b - 1;
  for (int i = 1; i <= n; ++i) {
    a[i] = std::lower_bound(b + 1, b + len + 1, a[i]) - b;
  }
  for (int i = 1; i <= n; ++i) {
    p1[i] = qry(a[i] - 1);
    p2[i] = qry(n) - qry(a[i]);
    mdf(a[i]);
  }
  for (int i = 1; i <= m; ++i) {
    scanf("%d%d", &q[i].l, &q[i].r);
    q[i].id = i;
  }
  std::sort(q + 1, q + m + 1, cmp);
  clear();
  for (int i = 1, L = 1, R = 0; i <= m; ++i) {
    int l = q[i].l, r = q[i].r;
    if (L > l) {
      v[R].push_back({l, L - 1, i, 1, 0});
    }
    while (L > l) {
      --L;
      q[i].ans -= p1[L];
    }
    if (R < r) {
      v[L - 1].push_back({R + 1, r, i, -1, 1});
    }
    while (R < r) {
      ++R;
      q[i].ans += p2[R];
    }
    if (L < l) {
      v[R].push_back({L, l - 1, i, -1, 0});
    }
    while (L < l) {
      q[i].ans += p1[L];
      ++L;
    }
    if (R > r) {
      v[L - 1].push_back({r + 1, R, i, 1, 1});
    }
    while (R > r) {
      q[i].ans -= p2[R];
      --R;
    }
  }
  for (int i = 1; i <= n; ++i) {
    mdf(a[i]);
    for (N j : v[i]) {
      for (int k = j.l; k <= j.r; ++k) {
        if (j.y == 0) {
          q[j.id].ans += j.x * qry(a[k] - 1);
        } else {
          q[j.id].ans += j.x * (i - qry(a[k]));
        }
      }
    }
  }
  for (int i = 2; i <= m; ++i) {
    q[i].ans += q[i - 1].ans;
  }
  for (int i = 1; i <= m; ++i) {
    ans[q[i].id] = q[i].ans;
  }
  for (int i = 1; i <= m; ++i) {
    printf("%lld\n", ans[i]);
  }
}
Luogu P5501 [LnOI2019] Coming and Going

Repeatedly query the sum of the “Abbi values” of all numbers in [l,r][l, r].

Definition of the Abbi value: if aia_i is the [l,r][l,r] interval's kk-th smallest element, its “Abbi value” is kaika_i.

Let f(x,r)f(x,r) be the sum of the numbers in [1,r][1,r] greater than axa_x, and let g(x,r)g(x,r) be the number of numbers in [1,r][1,r] greater than axa_x. Then, when moving the right endpoint to the right, the contribution is f(r,r−1)−f(r,l−1)+ar(r−l+1−(g(r,r−1)−g(r,l−1)))f(r,r-1)-f(r,l-1) + a_r(r-l+ 1-(g(r,r-1)-g(r,l-1))); the other directions can be written similarly and are omitted here.

The f(r,r−1)f(r, r - 1) and g(r,r−1)g(r, r - 1) terms above can still be precomputed, while the remaining terms are taken offline at the other endpoint and processed with a scanline. We still need O(n)O(n) insertions and O(nm)O(n\sqrt{m}) rank queries, so value-domain decomposition solves this problem as well.

Example code
#include <algorithm>
#include <cmath>
#include <cstring>
#include <iostream>
#include <vector>

constexpr int MAX = 1e5 + 5;
constexpr int MAXX = 4e2 + 5;
constexpr int MX = 1e5;

typedef long long ll;

ll n, m, sz, _sz, sm;
ll st[MAXX], ed[MAXX], bl[MAX], a[MAX], ans[MAX], sm1[MAXX], sm3[MAXX],
    sm2[MAX], sm4[MAX], s[MAX], pre[MAX];

void init() {
  sz = sqrt(MX);
  for (int i = 1; i <= sz; ++i) {
    st[i] = (MX / sz) * (i - 1) + 1;
    ed[i] = (MX / sz) * i;
  }
  ed[sz] = MX;
  for (int i = 1; i <= sz; ++i) {
    for (int j = st[i]; j <= ed[i]; ++j) {
      bl[j] = i;
    }
  }
}

void mdf(ll x) {
  sm += x;
  for (int i = bl[x]; i <= sz; ++i) {
    ++sm1[i];
    sm3[i] += x;
  }
  for (int i = x; i <= ed[bl[x]]; ++i) {
    ++sm2[i];
    sm4[i] += x;
  }
}

ll qry(ll x) {
  if (!x) {
    return 0;
  }
  return sm1[bl[x] - 1] + sm2[x];
}

ll _qry(ll x) {
  if (!x) {
    return 0;
  }
  return sm3[bl[x] - 1] + sm4[x];
}

void clear() {
  sm = 0;
  memset(sm1, 0, sizeof(sm1));
  memset(sm2, 0, sizeof(sm2));
  memset(sm3, 0, sizeof(sm3));
  memset(sm4, 0, sizeof(sm4));
}

struct Q {
  ll l, r, id, ans;
} q[MAX];

struct N {
  ll l, r, id, x;
};

std::vector<N> v[MAX];

bool cmp(Q a, Q b) {
  if (a.l / _sz != b.l / _sz) {
    return a.l < b.l;
  } else {
    return ((a.l / _sz) & 1) ? a.r < b.r : a.r > b.r;
  }
}

int main() {
  scanf("%lld%lld", &n, &m);
  _sz = n / sqrt(m);
  init();
  for (int i = 1; i <= n; ++i) {
    scanf("%lld", &a[i]);
    s[i] = s[i - 1] + a[i];
  }
  for (int i = 1; i <= n; ++i) {
    pre[i] = qry(a[i] - 1) * a[i] + sm - _qry(a[i]);
    mdf(a[i]);
  }
  clear();
  for (int i = 1; i <= m; ++i) {
    scanf("%lld%lld", &q[i].l, &q[i].r);
    q[i].id = i;
  }
  std::sort(q + 1, q + m + 1, cmp);
  for (int i = 1, L = 1, R = 0; i <= m; ++i) {
    int l = q[i].l, r = q[i].r;
    if (L > l) {
      v[R].push_back({l, L - 1, i, 1});
    }
    while (L > l) {
      --L;
      q[i].ans -= pre[L];
    }
    if (R < r) {
      v[L - 1].push_back({R + 1, r, i, -1});
    }
    while (R < r) {
      ++R;
      q[i].ans += pre[R];
    }
    if (L < l) {
      v[R].push_back({L, l - 1, i, -1});
    }
    while (L < l) {
      q[i].ans += pre[L];
      ++L;
    }
    if (R > r) {
      v[L - 1].push_back({r + 1, R, i, 1});
    }
    while (R > r) {
      q[i].ans -= pre[R];
      --R;
    }
  }
  for (int i = 1; i <= n; ++i) {
    mdf(a[i]);
    for (N j : v[i]) {
      for (int k = j.l; k <= j.r; ++k) {
        q[j.id].ans += 1ll * j.x * (qry(a[k] - 1) * a[k] + sm - _qry(a[k]));
      }
    }
  }
  for (int i = 2; i <= m; ++i) {
    q[i].ans += q[i - 1].ans;
  }
  for (int i = 1; i <= m; ++i) {
    ans[q[i].id] = q[i].ans + s[q[i].r] - s[q[i].l - 1];
  }
  for (int i = 1; i <= m; ++i) {
    printf("%lld\n", ans[i]);
  }
  return 0;
}

This content comes from OI Wiki: Mo's Algorithm: Secondary Offline Processing. See the original for authors and revision history. AIOJ has adapted the formatting, in-site links, and resource paths.

Content is shared under CC BY-SA 4.0 and SATA license; additional terms are governed by the original notice. Thanks to the OI Wiki authors and contributors; per the SATA terms, please visit the OI Wiki repository and star it to show your support.