Префиксные суммы за логарифм шагов: Хиллис - Стил и Блеллох

Если ты знаешь, чем скан Хиллиса - Стила отличается от скана Блеллоха по работе и глубине и как устроены подъём и спуск по дереву, листай до «Теперь сам».

Как думаешь, получится ли

Есть восемь чисел 3, 1, 4, 1, 5, 9, 2, 6, нужны все префиксные суммы: 3, 4, 8, 9, 14, 23, 25, 31. Обычный цикл делает семь шагов, и каждый ждёт предыдущий. Можно ли получить все восемь сумм за три шага, если исполнителей восемь?

Ответь себе до чтения. Сумму ВСЕХ восьми чисел дерево даёт за три шага - это было в прошлой главе. Но префиксов восемь, и кажется, что для каждого нужно своё дерево. Ответ - да, за три шага можно, а если согласиться на шесть шагов, то ещё и почти без лишних сложений (так выходит исключающий скан; включающий - ещё одним проходом).

На что это похоже

Колонна из восьми солдат, каждый держит карточку с числом, и каждый должен узнать сумму карточек от первого до себя. Можно передавать сумму по цепочке - семь передач подряд. А можно так: на первом шаге каждый узнаёт число соседа впереди и прибавляет. Теперь у каждого сумма двух. На втором шаге каждый спрашивает того, кто на два места впереди, - и у него сумма четырёх. На третьем - на четыре места. Через три шага у каждого сумма всех стоящих впереди: каждый шаг удваивает охват.

Механизм: две схемы

Хиллис - Стил (1986). На шаге с расстоянием d каждый элемент с номером i ≥ d прибавляет к себе значение x[i − d], взятое ДО этого шага. После шага с d в ячейке i лежит сумма последних 2d элементов до неё включительно (или всех, если их меньше). Шагов ⌈log₂ n⌉, но на каждом работает почти n сложений: при n, равном степени двойки, всего n · log₂ n − (n − 1). Ячейку x[i − d] на шаге читают двое - её владелец и сосед справа, - так что это шаг модели CREW; на EREW ячейку пришлось бы сначала скопировать для соседа - лишний шаг на каждом уровне.

Блеллох (1990) делает две фазы над тем же деревом, что и сумма из прошлой главы. Подъём: на шаге с расстоянием d ячейка i = 2d − 1, 4d − 1, … прибавляет к себе ячейку i − d. После подъёма в последней ячейке - сумма всего массива, а в промежуточных - суммы блоков. Спуск: в последнюю ячейку кладут ноль и идут вниз с d = n/2, n/4, …, 1: каждый узел отдаёт своё значение левому ребёнку, а правому - своё плюс старое значение левого. На выходе - исключающий скан: y[i] = x[0] + … + x[i − 1], без самого x[i]. Обе фазы - по log₂ n шагов и по n − 1 сложению.

Правило

Скан ассоциативной операции над n элементами считается за O(log n) шагов. Хиллис - Стил: ⌈log₂ n⌉ шагов, работа n · log₂ n - не работоэффективен. Блеллох: подъём и спуск по дереву, 2 log₂ n шагов, 2(n − 1) сложений - работоэффективен, но даёт исключающий скан, а включающий получается из него прибавлением x[i].

Разбор: Хиллис - Стил на восьми и на 1024 числах

Каждый шаг строится в новом массиве: в PRAM все читают старые значения и пишут новые разом. Если прибавлять на месте, слева направо, ячейка i прочитает уже обновлённую i − d, и суммы удвоятся.

def hillis_steele(x, pokazat):
    # включающий скан: на шаге с d элемент i >= d прибавляет x[i - d]
    x = x[:]
    n = len(x)
    d = 1
    shagov = 0
    slozhenij = 0
    while d < n:
        novyj = x[:]
        for i in range(d, n):
            novyj[i] = x[i] + x[i - d]
            slozhenij = slozhenij + 1
        x = novyj
        shagov = shagov + 1
        if pokazat:
            print("d =", d, x)
        d = d * 2
    return [shagov, slozhenij]

print(hillis_steele([3, 1, 4, 1, 5, 9, 2, 6], True))
print(hillis_steele([1] * 1024, False))

d = 1 [3, 4, 5, 5, 6, 14, 11, 8]
d = 2 [3, 4, 8, 9, 11, 19, 17, 22]
d = 4 [3, 4, 8, 9, 14, 23, 25, 31]
[3, 17]
[10, 9217]

Три шага и все восемь префиксов - ответ на вопрос из начала. Но сложений 17 вместо 7. На 1024 числах - 10 шагов и 9217 сложений против 1023 у цикла: сложений в девять раз больше (копирование массива на каждом шаге - приём программы, изображающей синхронный шаг, а не часть алгоритма). На восьми ядрах это хуже цикла; на видеокарте с тысячами потоков и коротким массивом - лучше: глубина минимальна.

Разбор: Блеллох - подъём и спуск

Тот же массив, теперь работоэффективно. Программа требует, чтобы n было степенью двойки, - как и дерево в стандартном изложении; иначе массив дополняют нейтральным элементом операции: нулями для суммы, единицами для произведения, очень малыми числами для максимума.

def blelloch(x, pokazat):
    # исключающий скан: подъём строит дерево сумм, спуск раздаёт префиксы
    x = x[:]
    n = len(x)
    shagov = 0
    slozhenij = 0
    d = 1
    while d < n:
        i = 2 * d - 1
        while i < n:
            x[i] = x[i] + x[i - d]
            slozhenij = slozhenij + 1
            i = i + 2 * d
        shagov = shagov + 1
        d = d * 2
    if pokazat:
        print("после подъёма", x)
    x[n - 1] = 0
    d = n // 2
    while d >= 1:
        i = 2 * d - 1
        while i < n:
            levyj = x[i - d]
            x[i - d] = x[i]
            x[i] = x[i] + levyj
            slozhenij = slozhenij + 1
            i = i + 2 * d
        shagov = shagov + 1
        if pokazat:
            print("спуск d =", d, x)
        d = d // 2
    return [x, shagov, slozhenij]

r = blelloch([3, 1, 4, 1, 5, 9, 2, 6], True)
print(r)
r = blelloch([1] * 1024, False)
print(r[0][:5], r[0][-1], r[1], r[2])

после подъёма [3, 4, 4, 9, 5, 14, 2, 31]
спуск d = 4 [3, 4, 4, 0, 5, 14, 2, 9]
спуск d = 2 [3, 0, 4, 4, 5, 9, 2, 23]
спуск d = 1 [0, 3, 4, 8, 9, 14, 23, 25]
[[0, 3, 4, 8, 9, 14, 23, 25], 6, 14]
[0, 1, 2, 3, 4] 1023 20 2046

После подъёма в ячейке 7 лежит 31 - сумма всех, в ячейке 3 - сумма первой половины, 9. Спуск начинается с нуля в корне. На шаге d = 4 корень отдаёт левому ребёнку (ячейке 3) свой ноль, а себе оставляет 0 + 9 = 9 - префикс правой половины. Дальше то же на каждом уровне. Результат - исключающий скан: 0, 3, 4, 8, … - каждая сумма без своего элемента. Прибавь x[i] - получишь включающий.

На 1024 числах: 20 шагов, 2046 сложений - вдвое больше цикла, но в четыре с половиной раза меньше, чем у Хиллиса - Стила. Платой стали вдвое большая глубина и вторая фаза.

Что ломается без этого

Сложение на месте в Хиллисе - Стиле. Цикл for i in range(d, n): x[i] += x[i - d] без копии читает уже обновлённые соседние ячейки - ответ неверен с первого же шага. Исправление - новый массив на каждом шаге или проход справа налево.

Скан неассоциативной операции. Перестановка скобок законна только для ассоциативной операции: сложение, максимум, умножение матриц. Вычитание или усреднение «скользящим окном» так не считаются. Сложение чисел с плавающей точкой ассоциативно лишь приближённо - параллельный скан даёт немного другой ответ, чем цикл, и это не ошибка.

Путаница включающего и исключающего. Библиотеки дают оба: inclusive_scan и exclusive_scan в C++, cumsum в numpy - включающий. Взять не тот - значит сдвинуть все индексы на единицу.

Зачем это в работе

На видеокарте скан - стандартная операция библиотек (CUB в CUDA, thrust::exclusive_scan). Через него делают упаковку массива, радиксную сортировку, построение разреженных структур - следующий урок о применениях. В C++17 есть std::inclusive_scan и std::exclusive_scan с параллельной политикой исполнения, и оговорка о точности в них та же: операция обязана быть ассоциативной, иначе результат зависит от того, как реализация расставила скобки.

Сама идея «перестроить скобки ассоциативной операции в дерево» работает далеко за пределами сложения: так считают суммы с переносом в сумматорах процессора (схема Когге - Стоуна устроена как Хиллис - Стил, схема Брента - Кунга - как подъём и спуск Блеллоха) и решают линейные рекуррентные соотношения.

Теперь сам

1. Сколько шагов и сложений сделает Хиллис - Стил на 16 числах? А Блеллох?
2. Исключающий скан Блеллоха дал 0, 2, 7, 7, 10. Каким был исходный массив, если его последний элемент равен 4?
3. Почему Хиллис - Стил на 1024 числах делает 9217 сложений, а не 1024 · 10 = 10240?