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

#алгоритм #любой_язык #дерево #структуры_данных #сложность

Я занимаюсь оптимизацией в одном проекте, который сильно не укладывается в требования
по производительности. Тут широко используются деревья Фенвика для быстрого получения
префиксных сумм (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), видимо, следует на этом же проходе одновременно корректировать
значения промежуточных узлов как при выполнении стандартного инкремента (если это вообще
возможно сделать одновременно).

Это все общие соображения для бинарных деревьев. Как это должно реализовываться конкретно
для дерева Фенвика я пока не представляю.

Буду благодарен за помощь с конкретными алгоритмами на псевдокоде или любом императивном
языке.

Если мои рассуждения на счет сложности неверны, то хотелось бы узнать в чем конкретно
я ошибаюсь.
    

Ответы

Ответ 1

Да, можно попробовать ускорить эти операции. Приведу для начала свою реализацию: 

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), но уже в такой его реализации он будет работать
быстрее, чем независимое изменение для каждой переменной.


Ответ 2

Операция выгрузки ответов во внешний массив делается за O(k + log n) следующим образом:


Отдельно храним исходный массив
Выгружаем за O(log n) сумму до первого элемента (i).
За O(1) прибавляем к ней элемент i+1, получаем сумму до элемента i+1.
Аналогично считаем суммы до всех элементов до i+k, итого O(k) операций.


С изменением сложнее - дерево Фенвика такую операцию так просто не поддерживает.
Однако вам может помочь дерево отрезков с групповым изменением - оно весьма популярно
в спортивном программировании. По этому поводу есть, например, руководство на e-maxx,
но никаким разумным стандартам читаемости и поддерживаемости кода оно не удовлетворит,
поэтому надо вникать и писать свою реализацию (возможно, есть библиотечная, но я не
изучал).

Дерево отрезков позволит сделать подсчёт суммы и массовое присваивание на любом отрезке
за O(log n) (независимо от k), а выгрузку операций можно либо сделать за O(k * log
n) "в лоб", либо за O(k + log n) примерно той же оптимизацией, что и с Фенвиком - сначала
считаем сумму до элемента i, а далее по одному находим значения следующих элементов.
Это можно делать массово, и там получится O(k) операций (придётся полностью обойти
несколько поддеревьев суммарного размера O(k)).


Ответ 3

Алгоритмы с использованием массива исходных данных (из которых строилось дерево)
будут, конечно, работать быстрее. Так что, если проект позволяет хранить копию исходных
данных параллельно с деревом Фенвика, рекомендую обратиться в первую очередь к ответам
других участников.

Мне же было интересно повозиться только с самим деревом Фенвика. Поэтому мои решения
не столь быстры, зато позволяют не хранить копию исходных данных.

Пункт 1 - получение значений (выгрузка во внешний массив) префиксных сумм для элементов
из непрерывного диапазона от i до (i+k)

Перейти от префиксной суммы i-го элемента к префиксной сумме (i+1)-го на самом деле
несложно.
У меня получился следующий алгоритм:
1) собираем частичные суммы, составляющие i-й элемент, кладём их в структуру данных
"стек", в порядке их индексов в дереве от наименьших к наибольшим;
   в целях оптимизации можно хранить в стеке не сами частичные суммы а нарастающий
итог по ним;
   вершина стека при этом будет хранить значение i-й префиксной суммы (возвращаем
её в ответ);
например, при i = 10, стек будет таким:
   stack[0] = tree_item[7]
   stack[1] = tree_item[7] + tree_item[9]
   stack[2] = tree_item[7] + tree_item[9] + tree_item[10]  // вершина стека, равная
10-му префиксу
   кстати, максимальный размер этого стека совсем невелик - log2(max_size), т.е.
когда количество исходных элементов описывается 32-разрядным числом, размер данного
стека будет всего 32 элемента;
2) рассматриваем двоичную запись числа (i+1), записывая биты слева направо от старших
к младшим;
   если двоичная запись заканчивается нулём, переходим к следующему шагу;
   если двоичная запись заканчивается некоторой непрерывной последовательностью единиц,
то подсчитываем количество этих единиц, и удаляем с вершины стека столько элементов,
сколько обнаружили завершающих единиц;
например, при i = 10, имеем:
   (i+1) = 11 (десятичная запись) = 0001011 (двоичная запись);
   двоичная запись числа заканчивается группой из 2-х единиц;
   следовательно, снимаем со стека 2 элемента, и содержимое стека становится таким:
   stack[0] = tree_item[7]
3) добавляем в стек (i+1)-й элемент дерева;
   при этом вершина стека станет содержать значение (i+1)-й префиксной суммы (возвращаем
её в ответ);
например, при i = 10, имеем:
   stack[0] = tree_item[7]
   stack[1] = tree_item[7] + tree_item[11]  // вершина стека, равная 11-му префиксу
4) повторяем шаги 2-3 для индексов (i+2),(i+3),...(i+k).

для большей ясности выпишу ещё несколько шагов:
   (i+2), двоичная запись 0001100; из стека ничего не удаляется;
   stack[0] = tree_item[7]
   stack[1] = tree_item[7] + tree_item[11]
   stack[2] = tree_item[7] + tree_item[11] + tree_item[12]  // 12-й префикс
   (i+3), двоичная запись 0001101; из стека удаляется только 1 элемент;
   stack[0] = tree_item[7]
   stack[1] = tree_item[7] + tree_item[11]
   stack[2] = tree_item[7] + tree_item[11] + tree_item[13]  // 13-й префикс
   (i+4), двоичная запись 0001110; из стека ничего не удаляется;
   stack[0] = tree_item[7]
   stack[1] = tree_item[7] + tree_item[11]
   stack[2] = tree_item[7] + tree_item[11] + tree_item[13]
   stack[3] = tree_item[7] + tree_item[11] + tree_item[13] + tree_item[14]  // 14-й
префикс
   ...
все эти операции со стеком очень легко понять если смотреть на картинку, изображающую
дерево Фенвика в 2-мерном виде, примерно как здесь

По части владения Big-O-notation у меня не очень, поэтому не могу точно определить
сложность данного алгоритма. Могу только сказать, что это точно будет намного лучше,
чем O(k*log(i)), когда i достаточно велико.

Пункт 2 - заполнение константным значением элементов из непрерывного диапазона от
i до (i+k) с корректным пересчетом всех затрагиваемых сумм

Здесь у меня получилось очень похоже на то, что описано в ответе участника retorta,
поэтому не буду зря дублировать.
Упомяну только о том, что тут приходилось дополнительно вычислять исходные значения.
value[i] = prefix(i) - prefix(i-1)
Тривиальная реализация требует 2*k обращений к функции вычисления префиксной суммы,
и это легко оптимизировать до (k+1) обращений.
Однако, можно зайти и с другой стороны: неплохо оптимизируется функция, вычисляющая
i-е значение исходного массива, если основываться на том, что каждый чётный элемент
дерева есть значение исходное, а из нечётных каждый второй - равен сумме i-го и (i-1)-го
исходных значений, и т.д...
В общем же случае - опять "выстреливает" подсчет количества единичных битов в младших
разрядах числа, как и по пункту 1 :)

Реализация
Текста получилось уже и так довольно много, поэтом код я решил вынести на github,
надеюсь это не запрещается :)
Закодил в итоге на C++ в MS Visual Studio Community 2015.
В коде отсутствуют некоторые проверки, зато есть простенькие тесты на валидность
и производительность кода.
Отдельно хочу напомнить, что результаты замеров сильно различаются в Debug- и Release-конфигурациях,
в последней выигрыш обычно меньше.
Ссылка на код