Подсчёт слов

Привет! Хотел бы рассказать о процессе решения одной задачи. Она имеет очень простую формулировку и некоторые дополнительные условия, которые существенно влияют на подход к решению.

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

Уже за рамками "вступительного испытания" ради спортивного интереса можно было посоревноваться с авторским решением в скорости. Спустя примерно год после упомянутых событий у меня появилось свободное время, пришли новые идеи и я попытался найти предельно быстрое решение, о чём и пойдёт речь в статье.

Условие задачи

Исходное условие задачи:

Напишите считалку частотного словаря.
Словом считается набор латинских букв, [a-zA-Z]+. Все символы, кроме [a-zA-Z], считаются пробелами. Слова должны быть приведены к нижнему регистру. Выводить нужно частоту, затем само слово, причём результат должен быть отсортирован по частоте по убыванию, а слова, имеющие одинаковую частоту, должны быть отсортированы по алфавиту. Пример:

 $ cat in.txt
I agree with you
I say what you do about the world
It's poison
And it's sick
And you want to get out of it
 $ freq in.txt out.txt
 $ cat out.txt
3 it
3 you
2 and
2 i
2 s
1 about
1 agree
1 do
1 get
1 of
1 out
1 poison
1 say
1 sick
1 the
1 to
1 want
1 what
1 with
1 world

Решить нужно в один поток.

Как потом оказалось (и это существенный момент), считалку требовалось написать для конкретного файла: pg.txt (по ссылке — архив, содержащий файл). Файл содержит корпус из сконкатенированных текстов на английском, немецком и, возможно, некоторых других языках.

 $ md5sum *.txt
850944413ba9fd1dbf2b9694abaa930d *lf.txt
d1362c6c9de664c34c3066ec1d01078f *crlf.txt

соотв две эталонных md5

Также известно, что лучшее существующее решение отрабатывает за секунду на "среднем" современном процессоре архитектуры x86-64. Требуется дать как можно более быстрое решение.

Где всё найти

Код проекта находится в ветке master репозитория.

Окружение

Далее в тексте результаты замеров времени выполнения приведены для Intel(R) Xeon(R) CPU E5-2667 v2 @ 3.30GHz. Это Ivy Bridge. Память - DDR3 1866MHz ECC REG, 4 channels.

Мой выбор языка — C++. Компиляторы — GCC/clang актуальных версий и компилятор из Visual Studio 2019 и MinGW-w64 (в случае Windows).

На Windows я в итоге забил, так как результаты замеров времени на нём систематически оставляли желать лучшего, несмотря на настройку энергопотребления Performance.

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

Характеристики текста

 $ du -b pg.txt
336183276       pg.txt
 $ LC_ALL="C" wc -m pg.txt
336183276 pg.txt

  $ LC_ALL="C" awk -F '[^A-Za-z]+' '
 {
     for (i = 1; i <= NF; ++i)
         if ($i)
             ++j
 }
 END {
     print j
 }
 ' pg.txt
 58801281

 $ LC_ALL="C" awk -F '[^A-Za-z]+' '
{
    for (i = 1; i <= NF; ++i)
        if ($i)
            ++w[tolower($i)]
}
END {
    print length(w)
}
' pg.txt
213637

 $ LC_ALL="C" awk -F '[^A-Za-z]+' '
{
    for (i = 1; i <= NF; ++i)
        if ($i)
            ++l[length($i)]
}
END {
    for (j in l)
        print j, l[j]
}
' pg.txt | sort -k1g | column -t
1   3018380
2   10389242
3   13784482
4   10605768
5   6407052
6   4713535
7   3845034
8   2404953
9   1742863
10  1020733
11  455583
12  249480
13  108981
14  38231
15  12161
16  3546
17  881
18  219
19  73
20  20
21  15
22  8
23  2
24  7
25  1
27  2
28  2
29  5
32  3
34  1
35  1
37  3
42  1
44  2
48  1
51  1
53  1
58  1
61  2
64  1
66  2
68  1
70  1

 $ LC_ALL=C awk -F '[^A-Za-z]+' '
{
    for (i = 1; i <= NF; ++i)
        if ($i)
            print length($i)
}
' pg.txt | sort -n | jq -s '
{
    count: length,
    mininum: min,
    average: (add / length),
    maximum: max,
    median: .[length / 2  floor],
    p90: .[length * 0.90 | floor],
    p95: .[length * 0.95 | floor],
    p99: .[length * 0.99 | floor],
    p999: .[length * 0.999 | floor]
}
'
{
  "count": 58801281,
  "mininum": 1,
  "average": 4.25427611007318,
  "maximum": 70,
  "median": 4,
  "p90": 8,
  "p95": 9,
  "p99": 11,
  "p999": 13
}

 $ LC_ALL=C awk -F '[^A-Za-z]+' '
{
    for (i = 1; i <= NF; ++i)
        if ($i)
            print length($i)
}
' pg.txt | sort -n | awk '
/17/ { if (!i) i = NR }
END { print i / NR }
'
0.999979

Референсное решение

Несмотря на то, что есть md5 hash результата, для отладки разрабатываемого алгоритма потребуется реализация пусть медленная, но корректная и ответ которой можно будет посмотреть глазами и сравнить с ответом более продвинутой, но ещё пока неверной, реализации.

Самое простое (для меня) решение — это написать считалку "на bash" (что подразумевает использование sed, awk, sort, uniq и т.п.). Выглядит она так:

LC_ALL="C" awk -F '[^A-Za-z]+' '
{
    for (i = 1; i <= NF; ++i)
        if ($i)
            ++w[tolower($i)]
}
END {
    for (i in w)
        print w[i], i
}
' | sort -k1gr,2

Помещаю этот код в файл freq.sh и запускаю, передавая на вход pg.txt, а (эталонный) вывод сохраняя в ref.txt.

 $ time bash freq.sh <pg.txt | tee ref.txt | md5sum -
real    0m24.834s
user    0m24.759s
sys     0m0.090s
850944413ba9fd1dbf2b9694abaa930d  -

Криптографический хэш результата верный, а следовательно верным является и содержимое (коллизии md5 практически невероятны). Время сильно зависит от версии awk. Год назад у меня это решение отрабатывало за 120 секунд в том же окружении.

Характеристики результата

 $ head -16 ref.txt | column -t
3343241  the
1852717  and
1715705  of
1560152  to
1324244  a
956926   in
933954   i
781286   he
713514   that
690876   was
665710   it
608451   you
580460   his
498427   with
442288   for
433687   had

 $ LC_ALL="C" awk -F '[^A-Za-z]+' '
{
    for (i = 1; i <= NF; ++i)
        if ($i)
            w[tolower($i)]
}
END {
    for (i in w)
        ++l[length(i)]
    for (j in l)
        print j, l[j]
}
' pg.txt | sort -k1g | column -t
1   26
2   442
3   2958
4   10757
5   21406
6   31287
7   36105
8   33983
9   28101
10  20362
11  12875
12  7684
13  4045
14  1955
15  921
16  384
17  181
18  60
19  33
20  13
21  12
22  7
23  2
24  7
25  1
27  2
28  2
29  5
32  3
34  1
35  1
37  2
42  1
44  2
48  1
51  1
53  1
58  1
61  2
64  1
66  2
68  1
70  1

 $ LC_ALL=C awk -F '[^A-Za-z]+' '
{
    for (i = 1; i <= NF; ++i)
        if ($i)
            w[tolower($i)]
}
END {
    for (i in w)
        print length(w)
}
' pg.txt | sort -n | tee /tmp/1 | jq -s '
{
    count: length,
    sum: add,
    mininum: min,
    average: (add / length),
    maximum: max,
    median: .[length / 2 | floor],
    p90: .[length * 0.90 | floor],
    p95: .[length * 0.95 | floor],
    p99: .[length * 0.99 | floor],
    p999: .[length * 0.999 | floor]
}
'
{
  "count": 213637,
  "sum": 1669418,
  "mininum": 1,
  "average": 7.814273744716505,
  "maximum": 70,
  "median": 8,
  "p90": 11,
  "p95": 12,
  "p99": 14,
  "p999": 17
}

Методика измерений

Здесь стоит сразу оговорить методику измерений.

Важно создать среду, в которой воспроизводимость результатов измерений будет сохраняться как можно более стабильной.

Джиттер неизбежен, но можно свести его к минимуму. Источником джиттера могут быть прерывания потока вследствие обычной работы планировщика ОС, изменения частоты работы ядра, на котором выполняется поток, вследствие изменения режима работы процессора.

sudo find /sys/devices/system/cpu -name scaling_governor \
    -exec sh -c 'echo performance >{}' ';'
sudo sh -c 'echo off >/sys/devices/system/cpu/smt/control'
sudo tuna --cpus=1-7 --isolate

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

 $ grep MHz /proc/cpuinfo
cpu MHz         : 3043.824
cpu MHz         : 3722.037
cpu MHz         : 3721.075
cpu MHz         : 3674.557
cpu MHz         : 3592.679
cpu MHz         : 3550.398
cpu MHz         : 3513.916
cpu MHz         : 3612.262
cpu MHz         : 3468.005
cpu MHz         : 3764.145
cpu MHz         : 3586.930
cpu MHz         : 3576.853
cpu MHz         : 3564.029
cpu MHz         : 3587.192
cpu MHz         : 3244.385
cpu MHz         : 2761.846

Вторая команда отключает HyperThreading (это интелловское название SMT). В результате половина ядер "исчезает" (например, из вывода htop) — остаются ядра с номерами 0-7. В некоторых случаях это даёт 30% ускорение для однопоточных задач.

Третья команда настраивает планировщик Linux не назначать потоки (что почти всегда соблюдается ОС) на ядра 1-7. Таким образом, все процессы в системе вытесняются на ядро 0. Также желательно закрыть браузер, выключить Wi-Fi и/или Ethernet и прочие источники прерываний.

Следующей командой можно посчитать количества переключений контекста на всех ядрах за 10 секунд и убедиться, что подавляющее большинство из них происходят на ядре 0:

sudo perf stat -e 'sched:sched_switch' -a -A --timeout 10000

Делаем множество запусков программы на ядрах 1-7. В качестве результата берём минимальное значение среди всех полученных (здесь нет ошибки: именно минимальное).

for i in {1..20}
do
    time taskset --cpu-list 1-7 ./freq pg.txt out.txt
done

После проведения замеров времени выполнения требуется вернуть систему в исходное состояние:

sudo sh -c 'echo on >/sys/devices/system/cpu/smt/control'
sudo tuna --cpus=0-15 --include
sudo find /sys/devices/system/cpu -name scaling_governor \
    -exec sh -c 'echo powersave >{}' ';'

Полное руководство по настройке системы вы можете прочитать здесь: Low latency guide. Другие статьи там также крайне интересны.

Наивное решение

#include <algorithm>
#include <fstream>
#include <functional>
#include <streambuf>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <unordered_map>
#include <vector>

#include <cstdlib>
#include <cctype>

int main(int argc, char * argv[])
{
    if (argc < 3) {
        return EXIT_FAILURE;
    }

    std::ifstream i(argv[1]);
    if (!i.is_open()) {
        return EXIT_FAILURE;
    }

    i.seekg(0, std::ios::end);   
    auto size = size_t(i.tellg());
    i.seekg(0, std::ios::beg);

    std::string input;
    input.resize(size);
    i.read(input.data(), size);

    auto toLower = [](char c) {
        return std::tolower(std::make_unsigned_t<char>(c));
    };
    std::transform(std::cbegin(input), std::cend(input),
            std::begin(input), toLower);

    std::unordered_map<std::string_view, size_t> wordCounts;

    auto isAlpha = [](char c) {
        return std::isalpha(std::make_unsigned_t<char>(c));
    };
    auto end = std::end(input);
    auto beg = std::find_if(std::begin(input), end, isAlpha);
    while (beg != end) {
        auto it = std::find_if(beg, end, std::not_fn(isAlpha));
        ++wordCounts[{beg, it}];
        beg = std::find_if(it, end, isAlpha);
    }

    std::vector<const decltype(wordCounts)::value_type *> output;
    output.reserve(wordCounts.size());
    for (const auto & wordCount : wordCounts) {
        output.push_back(&wordCount);
    }

    auto isLess = [](auto lhs, auto rhs) -> bool
    {
        return std::tie(rhs->second, lhs->first) < 
            std::tie(lhs->second, rhs->first);
    };
    std::sort(std::begin(output), std::end(output), isLess);

    std::ofstream o(argv[2]);
    if (!o.is_open()) {
        return EXIT_FAILURE;
    }
    for (auto wordCount : output) {
        o << wordCount->second << ' ' << wordCount->first << '\n';
    }
    
    return EXIT_SUCCESS;
}

Компилирую, запускаю и убеждаюсь, что результат верный:

clang++ -std=c++20 -march=native -Ofast naive.cpp -o naive &&
    time LC_ALL="C" ./naive pg.txt out.txt &&
    md5sum out.txt ref.txt
real    0m5.399s
user    0m5.248s
sys     0m0.151s
850944413ba9fd1dbf2b9694abaa930d  out.txt
850944413ba9fd1dbf2b9694abaa930d  ref.txt

Программа отрабатывает существенно быстрее, чем freq.sh: 5.4 секунды.

Что здесь можно оптимизировать? Можно попробовать другие реализации хэш-таблиц. Та же std::unordered_map из libc++ (выше были приведены результаты для libstdc++) отрабатывает уже за 3.9 секунды.

• unordered_map.cpp

• folly::F14ValueMap - за 4.1 (folly.cpp)

• google::sparse_hash_map - за 3.6 секунды (sparse_hash_map.cpp)

• absl::flat_hash_map - за 2.6 (absl.cpp)

• google::dense_hash_map - за 2.4 (dense_hash_map.cpp)

Также можно ускорить ввод, его предварительную обработку (приведение к нижнему регистру) и вывод: с ~0.5 секунды до ~0.1 секунды.

Общий код для всех библиотечных контейнеров: common.hpp.

Лучшее время - 2.4 секунды.

Второй подход: решение через trie

trie.cpp

Результатам "наивного решения" до одной секунды далеко. Нужно искать и пробовать альтернативные способы решения.

Следующей идеей было использовать префиксное дерево. Модификация префиксного дерева, которая подходит для подсчёта слов, хранит в каждом узле количество раз, которое слово, соответствующее пути, оканчивающемуся в данном узле, встретилось в тексте. После того, как дерево построено, все листья содержат ненулевое значение счётчика. Также ненулевое значение счётчик может иметь и во внутренних узлах, так как одно слово может быть префиксом другого слова.

Далее следует фрагмент решения. Класс Timer (для замера времени) и классы InputStream и OutputStream (кастомный ввод/вывод большими блоками) не приведены для краткости.

#include "helpers.hpp"
#include "io.hpp"
#include "timer.hpp"

#include <fmt/color.h>
#include <fmt/format.h>

#include <algorithm>
#include <chrono>
#include <iterator>
#include <memory>
#include <string_view>
#include <utility>
#include <vector>

#include <cassert>
#include <cstdint>
#include <cstdio>
#include <cstdlib>

namespace
{

constexpr std::size_t kAlphabetSize = 'z' - 'a' + 1;

alignas(__m128i) char input[1 << 29];
auto inputEnd = input;

struct TrieNode
{
    uint32_t count = 0;
    uint32_t children[kAlphabetSize] = {};
};

}  // namespace

int main(int argc, char * argv[])
{
    Timer timer{fmt::format(fg(fmt::color::dark_green), "total")};

    if (argc != 3) {
        fmt::print(stderr, "usage: {} in.txt out.txt\n", argv[0]);
        return EXIT_FAILURE;
    }

    using namespace std::string_view_literals;

    auto inputFile =
        (argv[1] == "-"sv) ? wrapFile(stdin) : openFile(argv[1], "rb");
    if (!inputFile) {
        fmt::print(stderr, "failed to open '{}' file to read\n", argv[1]);
        return EXIT_FAILURE;
    }

    auto outputFile =
        (argv[2] == "-"sv) ? wrapFile(stdout) : openFile(argv[2], "wb");
    if (!outputFile) {
        fmt::print(stderr, "failed to open '{}' file to write\n", argv[2]);
        return EXIT_FAILURE;
    }

    timer.report("open files");

    std::size_t readSize =
        readInput(std::begin(input), std::size(input), inputFile);
    if (readSize == 0) {
        return EXIT_SUCCESS;
    }
    inputEnd += readSize;
    timer.report("read input");

    toLower(input, inputEnd);
    timer.report("make input lowercase");

    std::vector<TrieNode> trie(1);
    {
        uint32_t index = 0;
        for (auto i = input; i != inputEnd; ++i) {
            if (*i != '\0') {
                uint32_t & child = trie[index].children[*i - 'a'];
                if (child == 0) {
                    child = uint32_t(trie.size());
                    index = child;
                    trie.emplace_back();
                } else {
                    index = child;
                }
            } else if (index != 0) {
                ++trie[index].count;
                index = 0;
            }
        }
    }
    fmt::print(stderr, "trie size = {}\n", trie.size());

    timer.report(fmt::format(fg(fmt::color::dark_blue), "count words"));

    std::vector<std::pair<uint32_t, uint32_t>> rank;
    std::vector<char> words;

    std::vector<char> word;
    auto traverseTrie = [&](const auto & traverseTrie,
                            const auto & children) -> void {
        size_t c = 0;
        for (uint32_t child : children) {
            if (child != 0) {
                const TrieNode & node = trie[child];
                word.push_back('a' + c);
                if (node.count != 0) {
                    rank.emplace_back(node.count, uint32_t(words.size()));
                    words.insert(std::cend(words), std::cbegin(word),
                                 std::cend(word));
                    words.push_back('\0');
                }
                traverseTrie(traverseTrie, node.children);
                word.pop_back();
            }
            ++c;
        }
    };
    traverseTrie(traverseTrie, trie.front().children);
    assert(word.empty());
    fmt::print(stderr, "word count = {}, length = {}\n", rank.size(),
               words.size());

    timer.report("recover words from trie");

    std::stable_sort(std::begin(rank), std::end(rank),
                     [](auto && l, auto && r) { return r.first < l.first; });

    timer.report(fmt::format(fg(fmt::color::dark_orange), "sort words"));

    OutputStream<> outputStream{outputFile};
    for (const auto & [count, word] : rank) {
        if (!outputStream.print(count)) {
            fmt::print(stderr, "output failure");
            return EXIT_FAILURE;
        }
        if (!outputStream.putChar(' ')) {
            fmt::print(stderr, "output failure");
            return EXIT_FAILURE;
        }
        if (!outputStream.print(std::next(words.data(), word))) {
            fmt::print(stderr, "output failure");
            return EXIT_FAILURE;
        }
        if (!outputStream.putChar('\n')) {
            fmt::print(stderr, "output failure");
            return EXIT_FAILURE;
        }
    }
    timer.report("write output");

    return EXIT_SUCCESS;
}

1.48 секунды. Уже существенно ближе к секунде, чем любой из вариантов с хэш-таблицами из библиотек. К сожалению, решение через trie не удастся ощутимо ускорить. Это связано с тем как на низком уровне функционирует работа с памятью в системе.

 $ getconf -a | grep -P '[123]_D?CACHE_'
LEVEL1_DCACHE_SIZE                 32768
LEVEL1_DCACHE_ASSOC                8
LEVEL1_DCACHE_LINESIZE             64
LEVEL2_CACHE_SIZE                  262144
LEVEL2_CACHE_ASSOC                 8
LEVEL2_CACHE_LINESIZE              64
LEVEL3_CACHE_SIZE                  26214400
LEVEL3_CACHE_ASSOC                 20
LEVEL3_CACHE_LINESIZE              64

libstdc++'s policy based data structures

pb_ds.cpp

Среди расширений libstdc++ есть trie-based containers. Они кастомизируются и умеют делать то, что требуется: быть счётным префиксным деревом.

Так выглядит нужный контейнер:

#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tag_and_trait.hpp>
#include <ext/pb_ds/trie_policy.hpp>


template<typename Key, typename Value>
using Trie =
    __gnu_pbds::trie<Key, Value,
        __gnu_pbds::trie_string_access_traits<Key, 'a', 'z'>,
        __gnu_pbds::pat_trie_tag,
        __gnu_pbds::null_node_update>;

Работает он медленнее, чем ручная реализация, представленная выше: 6.57 секунды.

Дальнейшие рассуждения

Что можно существенно улучшить в решении с префиксным деревом? Я затрудняюсь ответить: я пробовал path compression, пробовал близкие к корню пары букв хранить в массиве, пробовал ещё что-то, - ничего не давало выигрыша. До секунды ещё далеко, а значит нужно двигаться другим путём.

Как могло бы выглядеть решение, которое приводило бы к меньшему количеству случайных доступов к памяти, имело бы меньшее количество непредсказуемых условных ветвлений? В таком решении входная строка читается последовательно и для каждого очередного слова делается ровно один поход в память за значением соответствующего ему счётчика, которое затем инкрементируется и записывается обратно. Этот поход - единственное место, где случайного доступа не избежать.

Знание того, что слову соответствует некоторый конкретный счётчик, расположенный по некоторому определённому адресу в памяти - это так называемое взаимно однозначное отображение. Нам подойдёт односторонне отображение из строки в целое число - некоторая функция. Каким бы могло быть это отображение? К примеру, можно построить minimal perfect hash function. Для этого существует утилита gperf. Она строит minimal perfect hash function для некоторого набора строк. Однако, вычисление значения функции, которую генерирует gperf, является дорогим. Да и, как оказалось, построить её для 213637 значений практически нереально. Если ослабить требования и не требовать минимальности, т.е. отображения в диапазон целых чисел без пропусков, то найти perfect hash function для словаря из pg.txt вполне реально.

Генерация хэш-функций

Процессоры x86 и ARM имеют инструкции для расчёта CRC32C. Значение CRC32, посчитанное для строки, можно использовать как хэш строки. Такая хэш-функция имеет свои недостатки (см. секцию EDIT здесь), но при решении данной задачи они не важны. Плюсом использования хардварного расчёта значения хэш-функции на основе CRC32-инструкций является то, что выполняется оно очень быстро: reciprocal throughput - от 1 до 3 тактов для однобайтового варианта инструкции. Мощность множества принимаемых значений также обладает оптимальными свойствами: минимальный нативный тип на x86, позволяющий индексировать массив из 213637 элементов, - это как раз double word или uint32_t.

Как же сгенерировать целое семейство хэш-функций, имея только одну инструкцию? Инструкции, рассчитывающие CRC32, аккумулируют значение CRC32. То есть для расчёта значения CRC32 строки необходимо "скармливать" инструкции два операнда: текущее значение CRC32 и очередную порцию данных из строки, и делать следует это последовательно. Полученное в результате выполнения инструкции значение будет являться новым "текущим значением". Завершается процесс расчёта CRC32, когда строка пройдена полностью. На первом шаге "текущее значение" инициализируется некоторой константой. Варьируя значение этой константы можно получать разные хэш-функции.

Семейства хэш-функций надобятся, к примеру, в реализации фильтра Блума.

Количество коллизий

Посчитаем ожидаемое количество коллизий. Есть k = 213637 мячей и m = 2^32 корзин. Раскидываем мячи так, чтобы вероятность попасть в любую корзину была одинаковой и равной 1/m. Вероятность, что i-ый мяч попадёт в уже занятую корзину (i - 1)/m. Ожидаемое количество корзин с более, чем одним мячом таково:

Что составляет примерно 5.31 в нашем случае.

Вероятность найти хэш-функцию такую, что случится 0 коллизий для хэша длиной 32 бита (m = 2^32 значений) среди k = 213637 строк :

from decimal import Decimal
from functools import reduce
from operator import mul


def prob(k, m):
	return reduce(mul, map(lambda i: 1 - Decimal(i) / m, range(k)))

print(prob(213637, 2 ** 32))  # 0.004925409327321181763062307628

Эти расчёты говорят, что в среднем каждая двухсотая (~1/203) хэш-функция, генерируемая гипотетическим "генератором хэш-функций", не будет давать коллизий, т.е. будет совершенной.

Код для поиска значений seed, соответствующих совершенным хэш-функциям:

std::unordered_set<uint32_t> hashes;
for (uint32_t seed = 0;
        seed < std::numeric_limits<uint32_t>::max();  // Ctrl+C
        ++seed) {
    bool bad = false;
    for (const auto & word : words) {
        auto hash = seed;
        for (char c : word) {
            hash = _mm_crc32_u8(hash, uint8_t(c));
        }
        if (!hashes.insert(hash).second) {
            bad = true;
            break;
        }
    }
    hashes.clear();
    if (!bad) {
        fmt::print("{}\n", seed);
    }
}

На практике для pg.txt совершенные хэш-функции на базе CRC32C встречаются раз в 5 чаще. Объяснения этому я не имею (кстати, был бы рад увидеть в коментариях обсуждение причин такого расхождения теории с практикой).

Реализация: sparsest hash map

sparsest.cpp

Простейший вариант хэш-таблицы:

uint32_t counts[std::size_t(std::numeric_limits<uint32_t>::max()) + 1];

Для такой хэш-таблицы алгоритм подсчёта слов выглядит следующим образом. Бежим по строке и считаем хэш каждого очередного слова. Инкрементируем элемент массива counts с индексом, равным значению хэша. Если предыдущее значение счётчика было 0, то сохраняем указатель на слово в соответствующем элементе массива words. words имеет тот же размер, что и counts. После окончания подсчёта бежим по counts и для всех ненулевых элементов собираем слова из words. Сортируем полученные пары и выводим.

Проблемой такого подхода является то, что значения разбросаны по памяти и в одну страницу виртуальной памяти (4096 байт на x86) попадает в среднем примерно 1.026 элементов хэш-таблицы (либо 0). То есть буквально один элемент на одну страницу. Страницы физической памяти выделяются лениво, стоимость выделения ненулевая и происходит, если мы "трогаем" память - пишем в неё или читаем. В случае столь разреженной хэш-таблицы, как данная, получается, что почти на каждый элемент ОС выделяет отдельную страницу физической памяти.

Вторая проблема: для сбора результатов необходимо прочитать все 16 ГБ counts и найти все счётчики, не содержащие нули. При этом все незатронутые при подсчёте страницы опять же будут задеты (прочитаны) и для них будет выделена физическая память, что даст дополнительное замедление (и рост потребления памяти). Всего при подсчёте задето примерно 4.9% всех страниц, в которых располагается counts, но ОС придётся аллоцировать страницы физической памяти для всех 100% виртуальных.

В операционных системах есть функции, позволяющие узнать состояние страниц виртуальной памяти. В Linux для этих целей есть файл /proc/PID/pagemap для каждого процесса в системе. Для текущего процесса к нему можно обратиться как к /proc/self/pagemap. Файл состоит из 64-битных блоков, соответствующих страницам памяти процесса, старший бит которых установлен для тех страниц, для которых выделена физическая память. Для некоторой страницы, имеющей адрес первого байта p, вычислить смещение offset соответствующего блока в /proc/self/pagemap можно так: offset = std::uintptr_t(p) / getpagesize() - это смещение в файле в 8-байтных блоках. Используя знание о том, какие блоки виртуальной памяти не были затронуты, можно пропускать их целиком: в них наверняка нули.

Для того, чтобы код, содержащий массив из 2^32 элементов, скомпилировался, необходимо указать флаг -mcmodel=medium. Также для выставления смещения в адресном пространстве 64 - log2(4096) + log2(8) = 55 бит в файле /proc/self/pagemap потребуется включить поддержку 64-битных смещений в <stdio.h>. Делается это указанием -D_FILE_OFFSET_BITS=64 компилятору, в результате чего становится доступна нестандартная функция fseeko64.

При запуске данного варианта решения у меня возникла следующая проблема. В системе есть всего 32GB RAM (swap отключен) и аллоцировать виртуальной памяти больше, чем есть физической, ОС не даёт. systemd-run --scope -p MemoryMax=33G --user ./sparsest также не решает проблему. В итоге как минимум один из двух массивов counts и words пришлось урезать: вместо uint32_t для words (напомню, тип значения его элемента должен вмещать 1883056) я использовал 3-байтовый тип uint24:

#pragma pack(push, 1)
struct uint24 {
    unsigned value : 24;
};
#pragma pack(pop)

Результаты:

• 6.07 секунды, если читать все 16ГБ массива counts.

• 2.07 секунды, если использовать /proc/self/pagemap для отсева заведомо пустых страниц.

Такой подход не лучше, чем подход с trie. Двигаемся дальше.

Хэш-таблица поплотнее

Для размещения элементов будем использовать массив, размер которого меньше, чем размер области принимаемых значений хэш-функции.

Будем выбирать в качестве размера массива степени двойки. Почему? Потому что операцией взятия остатка от деления на степень двойки является операция побитового И с константой вида 0b111...111. Операция эта на полтора порядка (десятичных на этот раз) быстрее, чем целочисленное деление (частью которого и является на самом деле получение остатка, то есть код правых частей операторов x = a / b; y = a % b; оптимизируется в одну инструкцию, если они стоят рядом: idiv на x86). Замечу, что использование размеров из ряда степеней двойки - это существенная оптимизация не только для случая столь узкоспециализированной хэш-таблицы, которая потребуется в решении данной задачи, но и при разработке хэш-таблиц общего назначения. Операция взятия остатка нужна, чтобы адресоваться в массиве длиной n, используя значение хэш-функции, лежащее в диапазоне {0, 1, ..., m - 1}, когда m > n. Остаток от деления всегда будет в диапазоне {0, 1, ..., n - 1}.

Так как операция взятия остатка нарушает свойство взаимной однозначности отображения, то в одну ячейку массива теперь может попадать более одного элемента. Есть два подхода для того, чтобы иметь с этим дело. Один подход - это chaining, когда в каждой ячейке массива находится некоторый контейнер, в который складываются новые элементы, имеющие одинаковое значение хэш-функции по модулю размера масива. То есть такие, для которых произошла так называемая хэш-коллизия. Другой подход - open addressing. Он состоит в том, чтобы при хэш-коллизии искать по какому-либо фиксированному алгоритму следующую незанятую ячейку.

Open addressing предполагает какое-то сканирование памяти в надежде найти пустую ячейку. Это лишние чтения из памяти, которые даже в случае стратегии линейного поиска (англ. linear probing, при хэш-коллизии пробуем следующий-по-модулю-размера-массива элемент) будут негативно сказываться на производительности. Можно ли как-то красиво использовать chaining? Ведь chaining - это не обязательно односвязные списки, как обычно даётся в определении.

Можно в каждой ячейке хранить массив фиксированной длины (корзину), состоящий из элементов, и надеяться, что он не переполнится. Получается гибрид chaining и open addressing: в двумерном массиве при коллизиях по первому измерению постепенно заполняются элементы по второму.

Переполнение корзин

Ожидаемое количество элементов, попадающих в одну корзину, - 213637 / 2^b, где 2^b - количество корзин.
Вероятность того, что в одной из корзин находится не более a элементов:

Соответственно для 2^b корзин вероятность, что ни в одной из них нет более a элементов:

Какова же вероятность того, что нам посчастливится сгенерировать хэш-функцию, которая для 2^b корзин не даст переполнения ни в одной из них при размере корзины a? Будем рассматривать в качестве значений a некоторые степени двойки, которые могут быть количествами элементов хэш-таблицы, помещающимися в одну кэш-линию: 4, 8, 16, 32.

from decimal import Decimal, getcontext
from scipy.special import gammaincc


def prob(k, b, s):
    return Decimal(gammaincc(s + 1, k / b)) ** b


K = 213637
A = [4, 8, 16, 32]
B = range(13, 19)

print('|a\\b', end='|')
for b in B:
    print(b, end='|')
print()

print('|:---:', end='|')
for b in B:
    print(':---:', end='|')
print()

getcontext().prec = 3
for a in A:
    print(f'|{a}', end='|')
    for b in B:
        if a * 2 ** b >= K:
            print(prob(K, 2 ** b, a), end='|')
        else:
            print('-', end='|')
    print()

Выбор параметров хэш-таблицы

Допустим мы сгенерировали хэш-функцию не дающую коллизий на наших данных. Это значит, что для элементов попадающих в одну корзину, то есть имеющих совпадающий префикс хэша hashLow, суффикс хэша hashHigh не будет давать коллизий в рамках данной корзины. Если суффикс хэша хранить в корзине наряду с элементом, то, только лишь сравнивая одно значение сохранённого суффикса хэша, можно однозначно идентифицировать нужный элемент. При этом сами элементы не потребуется сравнивать на равенство.

Элементы, т.е. подстроки-слова, которые, напомню, могут иметь длину 70 байт, хранить в корзине контрпродуктивно. Хранение их в корзине раздуло бы её размер далеко за пределы одной кэш-линии (64 байта). При этом, данные эти в ходе подсчёта слов не читаются, а только единожды записываются, когда слово встречается впервые.

Таким образом в корзине будут храниться 3 значения: суффикс хэша, счётчик, указатель или индекс, указывающий на слово. При этом "горячими" из них будут являться только суффикс и счётчик. Указатель или индекс можно хранить в соседнем массиве words[hashLow][popcnt(hashHigh != 0)], имеющем те же размерности, что и массив с корзинами.

Далее стоит вопрос выбора AoS vs SoA ("array of structures versus structure of arrays"). Хранить ли суффиксы хэша вперемешку со счётчиками? Если хранить суффиксы хэша в отдельном массиве hashesHigh, то массив можно загружать как вектор в SSE/AVX регистр и сравнивать с новым суффиксом хэша за одну векторную операцию без необходимости бежать циклом по массиву, перебирая все сравниваемые значения по одному. Несмотря на то, что цикл этот имел бы фиксированное максимальное количество итераций, даже будучи развёрнутым полностью он скорей всего сгенерировал бы множество инструкций условного перехода, что наверняка дало бы код более медленный, чем пара векторных инструкций, оперирующих с массивом, как целым.

Определив, какой суффикс хэша совпадает с текущим, инкрементируем соответствующий элемент массива счётчиков. Если же совпадающего значения не найдено, то значение кладём в первый незанятый элемент массива, инкрементируем соответственный элемент в массиве со счётчиками, копируем новое слово в массив со словами и сохраняем его индекс в массиве words.

Оценим возможное количество элементов в корзине. Счётчик - это uint32_t (если размер кэш-линии прижмёт, то можно опять же использовать uint24), суффикс хэша - это либо uint32_t, либо, если его размер не более 16 бит, uint16_t (или даже uint8_t, но здесь не пригодился). Решим неравенство в натуральных числах: (sizeof(hashHigh) + sizeof(counter)) * x = (4 + 4) * x <= 64. Откуда x <= 8. Для hashHigh типа uint16_t решение - x <= 10. SSE/AVX регистры имеют размерами степени двойки. Отсюда - корзина будет содержать 8 элементов.

Для b = 17 вполне реально найти совершенную хэш-функцию, не дающую переполнений даже для корзины, содержащей всего 8 элементов. При таких параметрах хэш-таблицы (a = 8, b = 17) доля занятых элементов (англ. occupancy) равна ~0.203. В одной корзине будет находится в среднем ~1.63 элемента.

На самом деле есть ещё один факт, который делает b = 17 оптимальным выбором. Для пометки незанятых элементов необходим так называемый "tombstone" - такое особое значение элемента массива hashesHigh, которое не может получиться в результате вычисления hashHigh. Если hashHigh имеет 15 значащих бит, то при типе элемента массива hashesHigh uint16_t старший бит всех незанятых элементов должен содержать единицу, чтобы любое возможное значение hashHigh, полученное как hash >> 17, отличалось от него. Существуют SSE/AVX инструкции (movemask), которые позволяют выделять значения старших битов элементов векторных регистров и помещать их в регистр общего назначения. Другие инструкции (bsf, lzcnt) далее позволяют вычислить индекс младшего установленного бита в регистре общего назначения, который и является индексом первого незанятого элемента hashesHigh. Всё это без каких-либо циклов и условных переходов.

Итоговый результат для данного подхода - это 0.892 секунды. Насколько мне известно, это лучший результат на данный момент. Когда я автору сообщил, что достиг субсекундной скорости, он ответил "ну, тогда возможно перегнал и с меня будет бутылка".

oaph.cpp

Выводы

Как видно, если понимать условие задачи буквально, то можно найти решение, которое будет эксплуатировать какие-то упущения формулировки.

Мне определённо больше по душе подход с trie. Решение же с фиксированным seed-ом, заточенное под конкретные входные данные, будет работать только для них (ну почти). Хотя и для подхода, который описан в статье, можно найти применение, например, если универсум фиксированный и известен заранее. Универсум при этом должен иметь не очень большое количество элементов, так как уже для 436492 элементов вероятность найти совершенную хэш-функцию на базе CRC32 меньше 2^-32.

m = Decimal(2 ** 32)
k = 1
value = Decimal(1)
while value * m >= 1:
	value *= 1 - Decimal(k) / m
	k += 1

print(k)  # 436492

Возможны области применения, где множество слов (или элементов иной природы), на котором нашли seed, незначительно отличается от множества слов в текстах, с которыми алгоритм работает дальше. Если ошибки, вызванные коллизиями или переполнением корзин, несущественны, то подход также может оказаться полезен. При этом переполнение следует специальным образом обработать, например, игнорировать.

On genre

На самом деле, конечно же, в статье приведено описание того, "как можно было бы научным методом прийти к решению". Само решение я нашёл случайно, хоть и перепробовал перед этим почти всё, что описано в статье. Повезло, что нужный seed (10675 при практической вероятности ~10^-5) нашёлся довольно быстро. Все вероятности были посчитаны постфактум, специально для статьи, чтобы дать предлагаемым подходам теоретическое обоснование.

Идея того, как могут быть устроены корзины, была мною найдена в статье Open-sourcing F14 for faster, more memory-efficient hash tables. Идея применить в решении совершенную хэш-функцию - собственная, как и способ работы с входной строкой.

Обращение к читателю

Обо всех неточностях, помарках и стилистических ошибках прошу сообщать в личку.

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

Спасибо Николаю Кривчанскому за корректуру и редактуру. Спасибо автору задачи за задачу.