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)` constructs an array of size `n` with indices [0, `n`), and all values
  initialized to `v`.
- `SqrtDecomposition<T>(lo, hi)` constructs an array from two random-access iterators, initialized
  to the elements of the range in the same order.
- `size()` returns the size of the array.
- `at(i)` returns the value at index `i`.
- `query(lo, hi)` returns the result of `combine()` applied to all indices from `lo` to `hi`,
  inclusive.
- `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()`.
- O(sqrt n) per call to `at()`, `update()`, and `query()`.

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

*/

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

template<class 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>(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<class It>
  SqrtDecomposition(It lo, It hi) : len(hi - lo), value(lo, hi) {
    init();
  }

  int size() const { return len; }
  T at(int i) const { return query(i, i); }

  T query(int lo, int hi) const {
    int blocklo = ceil(static_cast<double>(lo) / 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) {
    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 i = b * blocklen + 1; i < blockhi; i++) {
      block[b] = combine(block[b], value[i]);
    }
  }
};

/*** Example Usage and Output:

Values: 6 -2 4 8 10

***/

#include <cassert>
#include <iostream>
using namespace std;

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