RadixSort и резкая просадка производительности на больших значениях

Итак, есть следующая, немного сырая реализация поразрядной сортировки radix.hpp:

#include <string.h>
#include <stdint.h>
#include <functional>

template <typename T>
struct element_t {
    unsigned char* key;
    T* element;
};

/**
 * @brief Сортирует сложные объекты ( структуры и классы ), которые можно представить в виде числогого
 *        ключа. Имеет сложность O(nw), где n - количество элементов и w - количество байт на ключ.
 * @param begin - указатель на начало массива
 * @param end - указатель на следующий за конечным элемент
 * @param out - куда положить результат. Указатель НЕ должен указывать начало сортируемого элемента
 *              ( т.е. не должен быть равен `begin` ).
 * @param keyFunc - функция, возвращающая числовой ключ для переданного значения. В параметр outkey
 *                  записывается выходной ключ, а в keylen длина ключа в байтах.
 * @param freeKeys - если `true`, то все ключи из keyFunc высвобождаются из памяти
*/
template <typename T>
void radixSort(T* begin, T* end, T* out, const std::function<void(T& element, void** outkey, uint32_t* keylen)> &keyFunc, bool freeKeys = false) {
    uint64_t arraySize = end - begin;
    uint32_t keySize = 0;

    // *************************************************
    // *                   ПОДГОТОВКА                  *
    // *************************************************

    // Одним new аллоцируется место для обоих массивов сразу
    // Работает немного быстрее, да и в кэш должно попадать охотнее
    element_t<T>* elements = new element_t<T>[arraySize * 2];
    element_t<T>* elementBuffer = elements + arraySize;

    uint32_t currentKeySize;

    // Генерируем пары элемент-ключ с помощью keyFunc

    for(uint64_t i = 0; i < arraySize; i++) {
        element_t<T> &element = elements[i];

        keyFunc(begin[i], (void**)&element.key, &currentKeySize);
        keySize = std::max(keySize, currentKeySize);

        element.element = begin + i;
    }

    uint64_t counts[256];

    // *************************************************
    // *                  СОРТИРОВКА                   *
    // *************************************************

    for(uint32_t key = 0; key < keySize; key++) {
        memset(counts, 0, 256 * sizeof(uint64_t));

        // Конкрентно в этой реализации я сортирую значения по ключу ( текущему байту )
        // записью в `counts[i]` количества раз, сколько мне попался байт со значением
        // `i` в текущем положении ключа. Затем, я прохожусь по всем `counts` и для каждого
        // i-того элемента записываю в `counts[i]` сумму значений с 0 по i-ый элемент.
        // Теперь, отняв от каждого единицу, мы получим индекс крайнего правого индекса
        // для каждого из значений, а значит, мы можем пройтись по каждому значению и
        // записывать это значение в индекс, на который указывает `counts[i]`, параллельно
        // декрементируя `counts[i]` ( т.к. мы уже заняли один из индексов ).

        // Реализацию с похожей идеей можно нагуглить, но никакого адекватного объяснения
        // там нет, и понять это можно только доперев самому :/

        for(uint64_t i = 0; i < arraySize; i++) {
            counts[elements[i].key[key]]++;
        }

        for(uint64_t i = 1; i < 256; i++) {
            counts[i] += counts[i - 1];
        }

        for(uint64_t i = arraySize; i > 0; i--) {
            uint64_t idx = elements[i - 1].key[key];

            // Декрементируем сначала, а не в конце, чтобы получить индексы,
            // начинающиеся с 0
            counts[idx]--;

            elementBuffer[counts[idx]] = elements[i - 1];
        }

        // После swap'а массив `elements` всегда будет содержать отсортированные значения по текущему
        // байту ключа.
        std::swap(elements, elementBuffer);
    }

    // Теперь, выводим отсортированный результат
    for(uint64_t i = 0; i < arraySize; i++) {
        out[i] = *elements[i].element;
    }

    if(freeKeys) {
        for(uint64_t i = 0; i < arraySize; i++) {
            delete[] elements[i].key;
        }
    }

    delete[] elements;
}

Бенчмируется она топорным кодом, который для каждого размера массива просто прогоняет функцию 30000 раз, замеряет время каджого раза с помощью std::chrono::steady_clock::now() и выводит среднее время:

#include <iostream>
#include <algorithm>
#include <chrono>
#include <random>
#include "radix.hpp" // тут реализация

const uint64_t samplesCount = 1000;

void keyConverter(uint32_t& element, void** key, uint32_t *keylen) {
    (*key) = (void*)&element;
    (*keylen) = sizeof(uint32_t);
}

int main() {
    std::function<void(uint32_t& element, void** outkey, uint32_t* keylen)> convFunc(keyConverter);

    uint64_t arraySize = 10;

    uint32_t *benchArray = new uint32_t[20000000];
    uint32_t *sortedArray = benchArray + 10000000;

    // От 1 до 7 ( не включительно ) степени десятки
    for(uint32_t sp = 1; sp < 7; sp++) {
        uint64_t totalTime = 0;

        for(uint32_t i = 0; i < samplesCount; i++) {
            // Вот этот цикл проходит 30 тактов
            for(uint32_t valueLimiter = 2; valueLimiter < 0x80000000; valueLimiter *= 2) {
                for(uint32_t j = 0; j < arraySize; j++) {
                    benchArray[j] = rand() & (valueLimiter - 1);
                }

                auto start = std::chrono::steady_clock::now();

#ifndef TESTING_STD_SORT
                radixSort(benchArray, benchArray + arraySize, sortedArray, convFunc, false);
#else
                std::sort(benchArray, benchArray + arraySize);
#endif

                auto end = std::chrono::steady_clock::now();

                totalTime += std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
            }
        }

        std::cout << "O(" << arraySize << ")\t";
        
        // Для более красивого вывода
        if(arraySize < 10000) {
            std::cout << "\t";
        }

        std::cout << " = " << totalTime / (30 * samplesCount) << "us\n";

        arraySize *= 10;
    }

    delete[] benchArray;

    return 0;
}

На всякий напомню, что в случае обычных 32-разрядных чисел ( на которых этот алгоритм и тестируется ) поразрядная сортировка - линейный алгоритм, с небольшой константой.

Компилирую с помощью GCC 12.3 и флагом -O0 ( без оптимизаций ) и получаю следующую картину ( в секундах ):

На 10'000'000 значениях у radixSort огромная просадка по производительности. Я подозревал, что это может быть кэш процессора, т.к. значений уже очень много и, следовательно, я могу получать много промахов в кэше. Конкрентно мой процессор - Intel Core i7 7740X 4.30GHz, с 128КБ L1 кэша, 1МБ L2 кэша и 8МБ L3 кэша.

Однако, чтобы убедиться, что это именно кэш, я протестировал свой алгоритм на Intel Core 2 Duo E8400 3.00GHz c 64 Кб L1 кэша и 6МБ L2 кэша ( L3 отсутствует ) и получил уже следующее:

Почти все значения в 2-4 раза больше, но на 10'000'000 значениях у radixSort абсолютно такая же просадка с почти таким же временем, разве что на 0.3 секунды дольше, что, по моему мнению, странно, ведь на компьютере с этим процессором у меня DDR3 память на частоте 1333МГц ( против моей DDR4 на 2133МГц ).

Те же тесты, но с флагами -O2 -march=core2:

В данном случае вообще получилось, что на 10'000'000 Core 2 оказался лучше, чем Core i7.

У меня алгоритм разделён на подготовку и саму сортировку ( в самом коде они визуально отделены специальными комментариями ), поэтому я решил замерить их отдельно, чтобы убедится, что обе части имеют просадки в производительности:

-O2 -march=core2

Просадка как в подготовке, так и в сортировке, как и предполагалось.

Исходя из замеров я получаю, что на 10'000'000 Core 2 Duo если не быстрее Core i7, то уступает ему не сильно. Однако это не совсем вяжется с моим предположением о том, что просадка вызвана кэшем, ведь если это так, то, по идее, разница в скорости уже должна зависеть напрямую от частоты процессора и ОЗУ. В связи с этим у меня возник вопрос:

Можно ли утверждать, что это всё же дело в кэше процессора, или дело может быть в чём-то другом?

Я не знаток плюсов, и уж тем более могу чего-то не знать на уровне железа, однако у меня чистый интерес узнать хотя-бы на уровне теории ответ на этот вопрос, даже если ответ будет очень сложным.

UPD

Я решил детальнее изучить картину с построения, прежде всего, графиков.

Рассмотрим следующий кусок кода:

for(uint32_t i = 0; i < samplesCount; i++) {
    for(uint32_t valueLimiter = 2; valueLimiter < 0x80000000; valueLimiter *= 2) {

        for(uint32_t j = 0; j < arraySize; j++) {
            benchArray[j] = rand() & (valueLimiter - 1); // rand() % valueLimiter
        }

        ...
    }
}

С помощью valueLimiter я ограничиваю значения, которые могут попасться в массиве, т.е. фактически оно определяет максимум.

Если замерить время исполнения в зависимости от значения valueLimiter, то можно получить следующую картину:

На этот раз я решил взять больше процессоров, большинство из которых - современные. Видно, что все процессоры имеют +- схожий по форме график, отличающийся лишь некой "константой", на которую они домножены. Один Core 2 Duo показывает линейную зависимость, выделяясь среди прочих.

Для тех, кому интересно, таблица размеров кэшей для процессоров:

У некоторых графиков есть [RCM] на конце. Под этими буквами подразумевается модификация кода в radix.hpp:

for(uint64_t i = arraySize; i > 0; i--) {
    uint64_t idx = elements[i - 1].key[key];

    // Декрементируем сначала, а не в конце, чтобы получить индексы,
    // начинающиеся с 0
    counts[idx]--;

    elementBuffer[counts[idx]] = elements[i - 1];
}

Тут алгоритм читает elements[i - 1].key[key] от бОльших индексов к меньшим. Как отметил @StanislavVolodarskiy, этот код не очень хорошо утилизирует кэш.

Обновлённая версия:

#ifdef REDUCE_CACHE_MISS
// [RCM]

// Вот тут оптимизированная часть алгоритма, который хоть и идёт
// от бОльших индексов к меньшим, но читает по четыре элемента
// в удобном, прямолинейном для кэша формате.
uint64_t i = arraySize;

for(; i > 4; i -= 4) {
    element_t<T>& element3 = elements[i - 4];
    element_t<T>& element2 = elements[i - 3];
    element_t<T>& element1 = elements[i - 2];
    element_t<T>& element0 = elements[i - 1];

    uint64_t idx3 = element3.key[key];
    uint64_t idx2 = element2.key[key];
    uint64_t idx1 = element1.key[key];
    uint64_t idx0 = element0.key[key];

    counts[idx0]--;
    elementBuffer[counts[idx0]] = element0;

    counts[idx1]--;
    elementBuffer[counts[idx1]] = element1;

    counts[idx2]--;
    elementBuffer[counts[idx2]] = element2;

    counts[idx3]--;
    elementBuffer[counts[idx3]] = element3;
}

for(; i > 0; i--) {
    uint64_t idx = elements[i - 1].key[key];

    counts[idx]--;

    elementBuffer[counts[idx]] = elements[i - 1];
}

#else
// не [RCM]

for(uint64_t i = arraySize; i > 0; i--) {
    uint64_t idx = elements[i - 1].key[key];

    // Декрементируем сначала, а не в конце, чтобы получить индексы,
    // начинающиеся с 0
    counts[idx]--;

    elementBuffer[counts[idx]] = elements[i - 1];
}

#endif

Перейдём к анализу самой части сортировки. Замерим исполнение более маленьких участков кода. Начнём с части распределения элементов по корзинам:

for(uint64_t i = 0; i < arraySize; i++) {
    counts[elements[i].key[key]]++;
}

Как указал @MBo, этот участок кода беспощадно прыгает по памяти и, отсюда, плохо кэшируется. Замеры:

Что интересно, оба процессора по-разному отрабатывают в этих ситуациях ( причём проверенно: компилятор не делает никаких оптимизаций для кэша ). Core 2 показывает линейную зависимость, в то время как Core i7 поступает как-то по-своему и сильно ему уступает ( видимо в попытках закешировать по-умному тратит много своего времени ).

Теперь, замерим часть записи элементов в отсортированном формате:

#ifdef REDUCE_CACHE_MISS
// [RCM]

// Вот тут оптимизированная часть алгоритма, который хоть и идёт
// от бОльших индексов к меньшим, но читает по четыре элемента
// в удобном, прямолинейном для кэша формате.
uint64_t i = arraySize;

for(; i > 4; i -= 4) {
    element_t<T>& element3 = elements[i - 4];
    element_t<T>& element2 = elements[i - 3];
    element_t<T>& element1 = elements[i - 2];
    element_t<T>& element0 = elements[i - 1];

    uint64_t idx3 = element3.key[key];
    uint64_t idx2 = element2.key[key];
    uint64_t idx1 = element1.key[key];
    uint64_t idx0 = element0.key[key];

    counts[idx0]--;
    elementBuffer[counts[idx0]] = element0;

    counts[idx1]--;
    elementBuffer[counts[idx1]] = element1;

    counts[idx2]--;
    elementBuffer[counts[idx2]] = element2;

    counts[idx3]--;
    elementBuffer[counts[idx3]] = element3;
}

for(; i > 0; i--) {
    uint64_t idx = elements[i - 1].key[key];

    counts[idx]--;

    elementBuffer[counts[idx]] = elements[i - 1];
}

#else
// не [RCM]

for(uint64_t i = arraySize; i > 0; i--) {
    uint64_t idx = elements[i - 1].key[key];

    // Декрементируем сначала, а не в конце, чтобы получить индексы,
    // начинающиеся с 0
    counts[idx]--;

    elementBuffer[counts[idx]] = elements[i - 1];
}

#endif

Замеры:

Core i7 с -DREDUCE_CACHE_MISS ускорился в ~3 раза по сравнению с обычной версией, и на этот раз превосходит Core 2 Duo. Сам Core 2 с такой оптимизацией перестал показывать линейную зависимость и его график стал похож на график Core i7.

Я также провёл замеры с -DREDUCE_CACHE_MISS, чтобы узнать, стал ли Core i7 быстрее:

На этот раз во всех тестах Core i7 быстрее.

Из всех этих графиков напрашивается вывод о том, что система кэширования Core 2 отличается от Core i7, однако, по всей видимости, она это делает несколько более прямым и топорным способом, от чего растёт линейно, а Core i7 делает это умнее, рассчитывая на более частый, нежели на более общий, случай.

Получается, что если это действительно так, то ответом на мой вопрос

Можно ли утверждать, что это всё же дело в кэше процессора, или дело может быть в чём-то другом?

будет:

Да, это дело в кэше, а конкрентно просто в другом подходе к кэшированию информации.

а столь неожиданное превосходство Core 2 перед Core i7 заключается лишь в более эффективном для данной ситуации подходе к кэшированию, несмотря на, казалось бы, логичное предположение о том, что современные процессоры должны отрабатывать лучше везде.

В итоге, вопрос можно сузить до следующего:

Насколько данная теория может быть правдоподобна?