Оптимизация операций для двоичного индексированного дерева (дерево Фенвика)

1,00
р.
Я занимаюсь оптимизацией в одном проекте, который сильно не укладывается в требования по производительности. Тут широко используются деревья Фенвика для быстрого получения префиксных сумм (AKA промежуточные суммы) для массива с постоянно изменяющимся содержимым. Сама по себе структура очень полезная и сильно выручает, но здесь кроме стандартных операций инкремента произвольных элементов очень часто используются две дополнительные нестандартные операции (я не нашел аналогов в литературе). И на эти операции приходится значительная часть вычислительной нагрузки.
Получение значений (выгрузка во внешний массив) префиксных сумм для элементов из непрерывного диапазона от i до i + k. Обнуление (или, скажем, заполнение константой) элементов из непрерывного диапазона от i до i + k с корректным пересчетом всех затрагиваемых сумм.
Важно: Для обоих операций среднее значение k много меньше n (количества элементов в дереве).
На данный момент реализованы эти операции тупо и просто:
В цикле (отдельно!) для каждого элемента из диапазона выполняется нахождение его значения (стандартная операция для дерева Фенвика). Итого: сложность операции составляет O(k*log(n)) для k элементов. В цикле (отдельно!) для каждого элемента выполняется нахождение его значения, после чего выполняется его инкремент на число обратное его текущему значению (для сброса в ноль). Сложность, опять же, равна O(k*log(n)).
Я задумался над структурой деревьев и этими операциями. И мне кажется, обе они могут быть реализованы с меньшей сложностью. Ведь log(n) это глубина прохода от корня до произвольного элемента, а в этих операциях выполняется обход группы последовательно расположенных элементов. Из этого факта возможно извлечь пользу, так как ближайший общий узел дерева в лучшем случае лежит на глубине log(k) от любого из конечных узлов в непрерывном диапазоне длинной k. Значит для обхода всего диапазона можно спуститься от корня за log(n), после чего подниматься и спускаться от одного соседнего элемента к другому за log(k). Итого: сложность алгоритма должна составить O(log(n) + k*log(k)) (что намного лучше текущего варианта). Это для лучшего случая. Но и средняя оценка должна быть близка к этому, так как худший случай O(k*log(n)) маловероятен при k ≪ n.
Сейчас я залип над реализаций алгоритмов для этих двух операций. Пока что мне не совсем понятно как перейти от элемента i к элементу i+1, не возвращаясь к корню. По идее надо подниматься до ближайшего общего узла и спускаться обратно к следующему элементу, по дороге подсчитывая суммы, необходимые для получения значения следующего элемента. А для операции (2), видимо, следует на этом же проходе одновременно корректировать значения промежуточных узлов как при выполнении стандартного инкремента (если это вообще возможно сделать одновременно).
Это все общие соображения для бинарных деревьев. Как это должно реализовываться конкретно для дерева Фенвика я пока не представляю.
Буду благодарен за помощь с конкретными алгоритмами на псевдокоде или любом императивном языке.
Если мои рассуждения на счет сложности неверны, то хотелось бы узнать в чем конкретно я ошибаюсь.

Ответ
Да, можно попробовать ускорить эти операции. Приведу для начала свою реализацию:
int next2n(int x) { // only 32bit x = x - 1 x = x | (x >> 1) x = x | (x >> 2) x = x | (x >> 4) x = x | (x >> 8) x = x | (x >> 16) return x+1 } class BIT_plus { std::vector in std::vector num int n public: BIT_plus(int n) { init(n) }
BIT_plus(const std::vector& list) { init(list) }
void init(const std::vector& list) { init(list.size()) for(auto i=0u i < list.size() ++i) { inc(i, list[i]) } }
std::vector sum_dip(int start, int len) { int summ = sum(start) std::vector res for(int i=start i < start + len i++) { res.push_back(summ) summ += num[i+1] } return res }
void init(int n) { in.assign(n, 0) num.assign(n, 0) this->n = n }
int size() { return n }
void inc (int i, int delta) { num[i] += delta
for ( i < n i = (i | (i+1))) in[i] += delta }
int sum (int r) const { int result = 0 for ( r >= 0 r = (r & (r+1)) - 1) result += in[r] return result }
int sum (int l, int r) const{ return sum (r) - sum (l-1) }
void fill_const(int start, int len, int val) { int n2n = next2n(start + len + 1)-1 int sumdelta = 0 for(int i = start i < start + len i++) { int delta = val - num[i] sumdelta += delta for(int j = i j < n2n j = (j|(j+1))) { in[j] += delta } num[i] = val } for(int i = n2n i < n i = (i|(i+1))) { in[i] += sumdelta } }
auto get_num() const { return num } auto get_in() const { return in } }
Она немножко раздута лишними методами, но всё самое важное в ней есть. (За некоторую основу взял реализацию отсюда). В ней по сравнению с обычной реализацией есть два новых метода: fill_const (заполняет определенным числом диапазон исходной последовательности) и sum_range (возвращает префиксные суммы элементов из диапазона).

Основная идея реализации довольно проста. Она заключается в том, что помимо непосредственно массива для дерева Фенвика(in) мы храним и сами исходные элементы (num). 1. Тогда нахождение всех префиксных сумм из диапазона [i, i+k] происходит так: сначала мы находим префиксную сумму для i-го элемента. Это происходит как и в обычной реализации. Сложность вычисления будет O(log(n)). А дальше, чтобы найти sum(i+1) прибавить число num[i+1]. И действительно если sum(i) -- сумма первых i элементов, то sum(i) + num[i+1] -- сумма первых i+1 элемента. 2. Вторая часть с заполнением константой чуть сложнее. Сначала поймём, какие элементы нам нужно изменить в исходном массиве in. Если мы меняем элемент с индексом i, то сначала мы должны сменить элемент с индексом i, потом f(i), f(f(i)), где f(i) = i|(i+1) ( | -- побитовое или). Заметим, что если в какой-то момент мы попали в 2k-1, то дальше мы также будем попадать в 2i-1. Более того, если мы начинаем из некого числа i, то в следующее за ним число вида 2k-1 мы2 не пропустим (иначе как мы получим единицу в следующем разряде, а функция у нас монотонно возрастающая). Воспользуемся этим так: от числа j из [i, i+k) будем изменять ячейки (по стандартным правилам. Прибавляя к ним delta) пока номер ячейки, которую мы меняем будет меньше, чем следующее за i+k число вида 2k-1 (next2n(i+k)). Дальше же индексы массива, который мы будем менять для всех j совпадут. И менять мы их будем на сумарную величину изменений всех ячеек.
Вероятно, второй алгоритм можно несколько оптимизировать (сейчас он работает, кажется, за O(k*log(i)+log(N)), что асимптотически не быстрее (i может быть порядка n), но быстрее практически особенно при маленькиx i), но уже в такой его реализации он будет работать быстрее, чем независимое изменение для каждой переменной.