2.4.3 Square Root Decomposition

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

/*

Maintain a fixed-size array while supporting point updates and contiguous range aggregate queries.
Square root decomposition partitions the array into contiguous blocks of about $\sqrt{n}$ elements,
each caching the aggregate of its block. A range query combines the cached aggregates of the whole
blocks contained in the range with the individual elements of the partial blocks at either end,
while a point update refreshes only the single block containing the changed index.

The query operation is defined by an associative aggregate function `combine(a, b)`. The default
definition below assumes a numerical array type, supporting queries for the "min" of the target
range. Another possible query operation is "sum", in which case `combine(a, b)` should return
`a + b`.

The point update operation is defined by `apply_delta(v, d)`, which returns the new value at a
single updated index. The default definition below supports updates that "set" the chosen array
index to a new value. Another possible update operation is "increment", in which case
`apply_delta(v, d)` should return `v + d`.

- `SqrtDecomposition<T>(n, v = T{})` constructs an array of size `n` with indices $[0, `n`)$, and
  all values initialized to `v`.
- `SqrtDecomposition<T>(lo, hi)` constructs an array from the half-open random-access iterator range
  $[`lo`, `hi`)$.
- `size()` returns the size of the array.
- `at(i)` returns the value at index `i`.
- `query(lo, hi)` returns the aggregate of the values at indices in $[`lo`, `hi`]$.
- `update(i, d)` assigns the value `v` at index `i` to `apply_delta(v, d)`.

The supported operations are identical to those of the point-update segment tree in this section.

Time Complexity:
- O(n) per call to both constructors, where $n$ is the size of the array.
- O(1) per call to `size()` and `at()`.
- O(sqrt(n)) per call to `query()` and `update()`.

Space Complexity:
- O(n) for storage of the array elements.
- O(1) auxiliary for all operations.

*/

#include <algorithm>
#include <cassert>
#include <cmath>
#include <vector>

template<typename T>
class SqrtDecomposition {
  static T combine(const T &a, const T &b) { return std::min(a, b); }
  static T apply_delta(const T &v, const T &d) { return d; }

  int len, blocklen;
  std::vector<T> value, block;

  void init() {
    blocklen = std::max(1, static_cast<int>(std::sqrt(len)));
    int nblocks = (len + blocklen - 1) / blocklen;
    for (int i = 0; i < nblocks; i++) {
      T blockval = value[i * blocklen];
      int blockhi = std::min(len, (i + 1) * blocklen);
      for (int j = i * blocklen + 1; j < blockhi; j++) {
        blockval = combine(blockval, value[j]);
      }
      block.push_back(blockval);
    }
  }

 public:
  explicit SqrtDecomposition(int n, const T &v = T{}) : len(n), value(n, v) { init(); }

  template<typename It>
  SqrtDecomposition(It lo, It hi) : len(hi - lo), value(lo, hi) {
    init();
  }

  int size() const { return len; }

  T at(int i) const {
    assert(0 <= i && i < len);
    return value[i];
  }

  T query(int lo, int hi) const {
    assert(0 <= lo && lo <= hi && hi < len);
    int blocklo = (lo + blocklen - 1) / blocklen, blockhi = (hi + 1) / blocklen - 1;
    if (blocklo > blockhi) {
      T res = value[lo];
      for (int i = lo + 1; i <= hi; i++) {
        res = combine(res, value[i]);
      }
      return res;
    }
    T res = block[blocklo];
    for (int i = blocklo + 1; i <= blockhi; i++) {
      res = combine(res, block[i]);
    }
    for (int i = lo; i < blocklo * blocklen; i++) {
      res = combine(res, value[i]);
    }
    for (int i = (blockhi + 1) * blocklen; i <= hi; i++) {
      res = combine(res, value[i]);
    }
    return res;
  }

  void update(int i, const T &d) {
    assert(0 <= i && i < len);
    value[i] = apply_delta(value[i], d);
    int b = i / blocklen;
    int blockhi = std::min(len, (b + 1) * blocklen);
    block[b] = value[b * blocklen];
    for (int j = b * blocklen + 1; j < blockhi; j++) {
      block[b] = combine(block[b], value[j]);
    }
  }
};

/*** Example Usage ***/

#include <cassert>
using namespace std;

int main() {
  vector<int> a{6, -2, 1, 8, 10};
  SqrtDecomposition<int> sd(a.begin(), a.end());
  sd.update(2, 4);
  vector<int> expected{6, -2, 4, 8, 10};
  for (int i = 0; i < sd.size(); i++) {
    assert(sd.at(i) == expected[i]);
  }
  assert(sd.query(0, 3) == -2);
  return 0;
}