Для понимания этой статьи необходимо знать базовое ДО и желательно понимать, как работает обычное 2D ДО.

Решим следующую задачу:

Дана таблица $$$n \cdot m$$$, поступают 2 вида запросов:

1) Прибавить $$$x$$$ на прямоугольнике $$$[l1, r1], [l2, r2]$$$.

2) Посчитать сумму на прямоугольнике $$$[l1, r1], [l2, r2]$$$.

Почему нельзя использовать пуши? Для массовых операций с пушами должно быть выполнено свойство, что отложенная операция может складываться, но в случае 2D ДО мы должны складывать не числа, а отрезки, и их объединять уже не получиться. Но помимо пушей есть другой тип массовых операций.

Сперва научимся делать прибавление на отрезке и сумму на отрезке без пушей в одномерном ДО.

В каждой вершине будем хранить $$$sum[v]$$$ и $$$add[v]$$$ — сумму на отрезке и сколько мы прибавим каждому элементу на отрезке. Рассмотрим запрос прибавления: для всех вершин, посещенных нашей функцией (на картинке желтые и зелёные), мы корректно пересчитаем значение в вершине, для суммы это $$$sum[v] += (min(r, rx)-max(l, lx)) \cdot x$$$, но сумма в вершинах ниже останется старой (красные). Чтобы это учесть, в тех вершинах, в которых остановилась функция (зелёные) сделаем $$$add[v] += x$$$.

Номер корня равен 0, реализация на полуинтервалах

void upd(int l, int r, int lx, int rx, int v, int x) {
    if (l >= rx || r <= lx) return;
    sum[v] += (min(r, rx)-max(l, lx)) * x;

    if (l >= lx && r <= rx) {
        add[v] += x;
        return;
    }
    int m = (l + r) / 2;
    upd(l, m, lx, rx, v*2+1, x);
    upd(m, r, lx, rx, v*2+2, x);
}

Теперь запрос суммы. Чтобы узнать актуальную сумму в вершине, нужно взять сумму $$$add[p]$$$ по всем предкам, умноженную на длину отрезка, и $$$sum[v]$$$. Чтобы лучше понять, почему это так, попробуйте посчитать сумму для красной вершины на картинке. Функция запроса будет такой же, как и в обычном ДО, только нужно поддерживать сумму на пути от корня до вершины.

В реализации можно даже не поддерживать сумму на пути, а просто когда мы спускаемся добавлять к ответу $$$add[v] \cdot (min(r, rx)-max(l, lx))$$$, это равносильно тому, что описано выше. Позже именно такой способ будет использоваться для больших размерностей.

int ask(int l, int r, int lx, int rx, int v) {
    if (l >= rx || r <= lx) return 0;
    if (l >= lx && r <= rx) {
        return sum[v];
    }
    int m = (l + r) / 2;
    return ask(l, m, lx, rx, v*2+1) + ask(m, r, lx, rx, v*2+2) + (min(r, rx)-max(l, lx)) * add[v];
}

Теперь этот же способ можно использовать для двумерного ДО.

Сделаем ДО по первой координате. В вершине, отвечающей за отрезок $$$[x1,x2]$$$ будет массив, отвечающий за прямоугольник $$$[x1, x2], [0, m-1]$$$, где i-тый элемент равен $$$\displaystyle\sum_{i = x1...x2} mas[j][i]$$$, в частности, в листах ДО будут просто строки таблицы, а в корне будет массив сумм столбцов (посмотрите на рисунок для лучшего понимания). На этих массивах в каждой вершине делаем ДО $$$sum[v]$$$ и изначально заполненный нулями $$$add[v]$$$.

Sorry for my bad Paint

Выглядеть запрос прибавления будет так же, как и в одномерном случае: корректно пересчитываем $$$sum[v]$$$ для всех посещенных вершин, для этого прибавим на отрезке $$$sum[v][l2...r2] += x \cdot (min(rx, r1) - max(lx, l1))$$$. В вершинах, где остановилась функция, запомним, что для всех вершин ниже текущей для второй координаты нужно прибавить $$$x$$$ на отрезке $$$[l2, r2]$$$.

void upd(int l1, int r1, int lx, int rx, int l2, int r2, int v, int x) {
    if (l1 >= rx || r1 <= lx) return;

    sum[v].upd(l2, r2, x * (min(rx, r1) - max(lx, l1)));
    if (l1 >= lx && r1 <= rx) {
        add[v].upd(l2, r2, x);
        return;
    }
    int med = (l1 + r1) / 2;
    upd(l1, med, lx, rx, l2, r2, v*2+1, x);
    upd(med, r1, lx, rx, l2, r2, v*2+2, x);
} 

Запрос суммы тоже не будет отличаться от одномерного случая: в конечных вершинах прибавляем к ответу $$$sum[v][l2...r2]$$$, и спускаясь по ДО прибавляем к ответу $$$add[v][l2...r2] \cdot (min(rx, r1) - max(lx, l1))$$$.

int ask(int l1, int r1, int lx, int rx, int l2, int r2, int v) {
        if (l1 >= rx || r1 <= lx) return 0;
        if (l1 >= lx && r1 <= rx) {
            return sum[v].ask(l2, r2);
        }
        int med = (l1 + r1) / 2;
        return ask(l1, med, lx, rx, l2, r2, v*2+1) + ask(med, r1, lx, rx, l2, r2, v*2+2)
               + add[v].ask(l2, r2) * (min(rx, r1) - max(lx, l1));
    }

Можно заметить, что код практически не отличается от одномерного случая, на самом деле можно обобщить этот способ для k-мерного ДО, так как для этого нужно будет иметь k-1-мерное ДО. Однако на каждом слое создаётся 2 массива размера 4n, так что памяти будет $$$O(8^k * n)$$$ и один запрос $$$O(log^k n)$$$.

Что делать, если в задаче координаты до $$$10^9$$$? Разумеется можно сделать оба ДО неявными, но в таком случае потребление памяти станет неприлично большим, так что лучше сжать по первой координате, а по второй сделать неявное ДО, или же можно хитро сжать координаты отдельно для значений на отрезках каждой вершины, тогда память будет линейная.

Насколько в целом эта техника применима? Например таким ДО можно было сдать в лоб задачу A3 недавно прошедшего отбора tinkoff generation, но в целом техника довольно специфическая, разве что неплохо ускоряет задачу "сумма на отрезке и прибавить на отрезке" для одномерного случая(у меня ~2 раза быстрее пушей).

Вот полный код 2D ДО: https://pastebin.com/pPaGhizX