Эффективная генерация числа в заданном интервале

В подавляющем большинстве моих постов о генерации случайных чисел рассматривались в основном свойства различных схем генерации. Это может оказаться неожиданным, но производительность алгоритма рандомизации может зависеть не от выбранной схемы генерации, а от других факторов. В этом посте (на который меня вдохновила превосходная статья Дэниела Лемира) мы исследует основные причины снижения производительности генерации случайных чисел, которые часто перевешивают производительность движка ГПСЧ.

Представьте такую ситуацию:

В качестве домашнего задания Хуан и Саша реализуют одинаковый рандомизированный алгоритм на C++, который будет выполняться на одном университетском компьютере и с одним набором данных. Их код почти идентичен и отличается только в генерации случайных чисел. Хуан торопится на свои занятия по музыке, поэтому просто выбрал вихрь Мерсенна. Саша, с другой стороны, потратил несколько лишних часов на исследования. Саша провёл бенчмарки нескольких самых быстрых ГПСЧ, о которых недавно узнал из соцсетей, и выбрал наиболее быстрый. При встрече Саше не терпелось похвастаться, и он спросил Хуана: «Какой ГПСЧ ты использовал?»

«Лично я просто взял вихрь Мерсенна — он встроен в язык и вроде неплохо работает».

«Ха!», — ответил Саша. «Я использовал jsf32. Он намного быстрее, чем старый и медленный вихрь Мерсенна! Моя программа выполняется за 3 минуты 15 секунд!».

«Хм, неплохо, а моя справляется меньше, чем за минуту», — говорит Хуан и пожимает плечами. «Ну ладно, мне пора на концерт. Пойдёшь со мной?»

«Нет», — отвечает Саша. «Мне… эээ… нужно снова взглянуть на свой код».

Эта неловкая вымышленная ситуация не особо и вымышлена; она основана на реальных результатах. Если ваш рандомизированный алгоритм выполняется не так быстро, как хотелось бы, и узким местом похоже является генерация случайных чисел, то, как это ни странно, проблема может быть и не в генераторе случайных чисел!

Введение: случайные числа на практике

Большинство современных качественных генераторов случайных чисел создаёт машинные слова, заполненные случайными битами, то есть обычно генерирует числа в интервале [0..2³²) или [0..2⁶⁴). Но во многих случаях использования пользователям нужны числа в определённом интервале — например, для броска кубика или выбора случайной игральной карты нужны числа в небольших постоянных интервалах. Однако многим алгоритмам, от перемешивания и reservoir sampling до рандомизированных двоичных деревьев поиска требуются числа, берущиеся из других интервалов.

Методы

Мы рассмотрим множество различных методов. Чтобы упростить обсуждение, вместо генерации чисел в интервале [i..j) или [i..j] мы будем генерировать числа в интервале [0..k). Имея такую схему, мы сможем, например, генерировать числа в интервале [i..j), задав k = j — i, сгенерировав число в интервале [0..k), а затем прибавив к нему i.

Встроенные средства C++

Во многих языках есть встроенные средства для получения случайного числа в указанном интервале. Например, чтобы вытащить карту из колоды с 52 картами на таких скриптовых языках, как Perl и Python, мы можем написать соответственно int(rand(52)) и random.randint(0,52). В C++ мы аналогично можем использовать uniform_int_distribution.

Код на C++ для реализации такого подхода прост:

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    std::uniform_int_distribution<uint32_t> dist(0, range-1);

    return dist(rng);
}

Обычно во встроенных средствах используется одна из описанных ниже техник, но большинство пользователей просто использует эти средства, не задумываясь о том, что происходит «под капотом», считая, что эти средства правильно спроектированы и достаточно эффективны. В C++ встроенные средства более сложны, потому что они должны иметь возможность работать с довольно произвольными движками генерации — генератор, выдающий значения в интервале от -3 до 17, может быть вполне допустим и применяться с std::uniform_int_distribution для создания чисел в любом интервале, например [0..1000). То есть встроенные средства C++ слишком переусложнены для большинства случаев, в которых они применяются.

Классический остаток от деления (с перекосом)

Давайте перейдём от переусложнённого подхода к слишком упрощённому.

Когда я учился программированию, мы генерировали числа в интервале (например, для выбора карты в колоде из 52 карт) с помощью оператора остатка от деления. Чтобы получить число в интервале [0..52), мы писали rand() % 52.

На C++ такой подход можно реализовать следующим образом:

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    return rng() % range;
}

Несмотря на простоту такого подхода, он демонстрирует причину того, что получение чисел в нужном интервале обычно является медленной задачей — для него требуется деление (чтобы вычислить остаток, получаемый оператором %). Деление обычно как минимум на порядок величин медленнее, чем другие арифметические операции, поэтому единственная арифметическая операция занимает больше времени, чем вся работа, выполняемая быстрым ГПСЧ.

Но кроме низкой скорости он ещё и перекошен. Чтобы понять, почему rand() % 52 возвращает перекошенные числа, предположим, что rand() создаёт числа в интервале [0..2³²), и заметим что 52 не делит 2³² нацело, оно делит его 82 595 524 раз с остатком 48. То есть если мы используем rand() % 52, то у нас будет 82 595 525 способов выбрать первые 48 карты из колоды и всего 82 595 524 способов выбрать последние четыре карты. Другими словами, существует перекос 0,00000121% против этих последних четырёх карт (возможно, это короли!). Когда я был учеником и писал домашнюю работу о бросании кубиков или вытаскивании карт, никого обычно не волновали подобные крошечные перекосы, но при увеличении интервала перекос растёт линейно. Для 32-битного ГПСЧ ограниченный интервал меньше 2²⁴ имеет перекос меньше 0.5%, но выше 2³¹ перекос составляет 50% — некоторые числа будут возвращаться в два раза чаще, чем другие.

В этой статье мы в основном будем рассматривать методики, в которых используются стратегии для устранения систематической ошибки, но наверно стоит сказать, что для 64-битного ГПСЧ величина перекоса в обычных случаях применения скорее всего будет ничтожной.

Ещё одна проблема может заключаться в том, что некоторые генераторы имеют слабые младшие биты. Например, семейства ГПСЧ Xoroshiro+ и Xoshiro+ имеют младшие биты, которые не проходят статистических тестов. Когда мы выполняем % 52 (потому что 52 чётное), то передаём младший бит прямиком на выход.

Умножение чисел с плавающей запятой (с перекосом)

Ещё одна распространённая методика — использование ГПСЧ, генерирующего числа с плавающей запятой в интервале [0..1) с последующим преобразованием этих чисел в нужный интервал. Такой подход применяется в Perl, в нём рекомендуется использовать int(rand(10)) для генерации целого числа в интервале [0..10) генерацией числа с плавающей запятой с последующим округлением вниз.

На C++ этот подход записывается так:

static uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    double zeroone = 0x1.0p-32 * rng();
    return range * zeroone;
}

(Учтите, что 0x1.0p-32 — это двоичная константа с плавающей запятой для 2^-32, которую мы используем для преобразования случайного целого числа в интервале [0..2³²) в double в единичном интервале; вместо этого мы можем выполнить такое преобразование с помощью ldexp(rng(), -32), но когда я провёл бенчмарк такого подхода, он оказался гораздо медленнее.)

Этот подход столь же перекошен, как и классический остаток от деления, но перекос проявляется иначе. Например, если бы мы выбирали числа в интервале [0..52), то числа 0, 13, 26 и 39 встречались бы на один раз реже, чем другие.

Эта версия при обобщении до 64 бит ещё более неприятна, потому что для неё требуется тип с плавающей запятой, мантисса которого не менее 64 бит. На машинах x86 с Linux и macOS мы можем использовать long double, чтобы воспользоваться числами x86 с плавающей запятой повышенной точности, имеющими 64-битную мантиссу, но long double не портируется универсально на все системы — в некоторых системах long double эквивалентен double.

Есть и хорошая сторона — этот подход работает быстрее, чем решения на основе остатков для ГПСЧ со слабыми младшими битами.

Целочисленное умножение (с перекосом)

Метод умножения можно адаптировать к арифметике с фиксированной, а не с плавающей запятой. По сути, мы просто постоянно умножаем на 2³²,

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    uint32_t x = rng();
    uint64_t m = uint64_t(x) * uint64_t(range);
    return m >> 32;
}

Может показаться, что для этой версии требуется 64-битная арифметика, на процессорах x86 хороший компилятор скомпилирует этот код в 32-битную инструкцию mult (которая даёт нам два 32-битных выходных значения, одно из которых является возвращаемым значением). Можно ожидать, что эта версия будет быстрой, но она перекошена точно так же, как метод умножения чисел с плавающей запятой.

Деление с отбрасыванием (без перекоса)

Мы можем модифицировать схему умножения чисел с плавающей запятой в схему на основе деления. Вместо умножения x * range / 2**32 мы вычисляем x / (2**32 / range). Так как мы работаем с целочисленной арифметикой, округление в этой версии будет выполняться иначе, и иногда генерировать значения за пределами нужного интервала. Если мы будем отбрасывать эти значения (например избавляться от них и генерировать новые значения), то в результате получим методику без перекосов.

Например, в случае вытаскивания карты с помощью 32-битного ГПСЧ мы можем сгенерировать 32-битное число и разделить его на 2³²/52 = 82 595 524, чтобы выбрать карту. Эта методика работает, если случайное значение из 32-битного ГПСЧ меньше 52 × 82595524 = 2³²/32 – 48. Если случайное значение из ГПСЧ является одним из последних 48 значений верхней части интервала генератора, то его нужно отбросить и искать другое.

В нашем коде для этой версии используется трюк с делением 2³² на range без применения 64-битной математики. Для непосредственного вычисления 2**32 / range нам необходимо представить число 2³², которое слишком велико (на единицу!) для представления как 32-битного целого числа. Вместо этого мы учтём, что для беззнаковых целых унарная операция отрицания range вычисляет положительное значение 2³² – range; поделив это значение на range, мы получим ответ на единицу меньше, чем 2**32 / range.

Следовательно, код на C++ для генерации числе с помощью деления с отбрасыванием выглядит так:

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    // calculates divisor = 2**32 / range
    uint32_t divisor = ((-range) / range) + 1;
    if (divisor == 0) // overflow, it's really 2**32
        return 0;
    for (;;) {
        uint32_t val = rng() / divisor;
        if (val < range)
            return val;
    }
}

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

Остаток от деления (двойной) без перекосов — методика OpenBSD

Мы также можем применить подход с отбрасыванием для устранения перекоса в классическом методе остатка от деления. В примере с игральными картами нам снова нужно отбросить 48 значения. В этой версии вместо отбрасывания последних 48 значений мы (эквивалентно) отбрасываем первые 48 значений.

Вот реализация этого подхода на C++:

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    // calculates 2**32 % range
    uint32_t t = (-range) % range;
    for (;;) {
        uint32_t r = rng();
        if (r >= t)
            return r % range;
    }
}

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

В коде arc4random_uniform OpenBSD (который также используется в OS X и iOS) применяется эта стратегия.

Остаток от деления (одиночный) без перекоса — методика Java

В Java используется другой подход к генерации числа в интервале, в котором используется только одна операция деления с остатком, за исключением довольно редких случаев отбрасывания результата. Код:

static uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    uint32_t x, r;
    do {
        x = rng();
        r = x % range;
    } while (x - r > (-range));
    return r;
}

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

Целочисленное умножение без перекосов — методика Лемира

Почти так же, как мы устранили перекос из методики остатка от деления, мы можем устранить перекос из методики целочисленного умножения. Эта методика была придумана Лемиром.

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    uint32_t t = (-range) % range;
    do {
        uint32_t x = rng();
        uint64_t m = uint64_t(x) * uint64_t(range);
        uint32_t l = uint32_t(m);
    } while (l < t);
    return m >> 32;
}

Битовая маска с отбрасыванием (без перекосов) — методика Apple

В нашем последнем подходе полностью исключаются операции деления и остатка. Вместо них в нём используется простая операция маскирования для получения случайного числа в интервале [0..2^k), где k — наименьшее значение, такое, что 2^k больше интервала. Если значение слишком велико для нашего интервала, мы отбрасываем его и пробуем получить другое. Код показан ниже:

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    uint32_t mask = ~uint32_t(0);
    --range;
    mask >>= __builtin_clz(range|1);
    uint32_t x;
    do {
        x = rng() & mask;
    } while (x > range);
    return x;
}

Этот подход был принят компанией Apple, когда (в релизе macOS Sierra) она выполняла собственную ревизию кода arc4random_uniform.

Бенчмаркинг базовых методик

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

Мы используем три бенчмарка и протестируем методики со множеством различных ГПСЧ.

Бенчмарк Large-Shuffle

Вероятно, самый очевидный бенчмарк — это перемешивание. В этом бенчмарке мы имитируем выполнение крупномасштабного перемешивания. Для сортировки массива размером N мы должны сгенерировать числа в интервалах [0..N), [0..(N-1)), …, [0..1). В этом бенчмарке мы будем считать, что N — это максимально возможное число (для uint32_t это 2³²-1 ). Код:

for (uint32_t i = 0xffffffff; i > 0; --i) {
    uint32_t bval = bounded_rand(rng, i);
    assert(bval < i);
    sum += bval;
}

Заметьте, что мы «используем» каждое число, прибавляя его к sum (чтобы его не отбросила оптимизация), но не выполняем никакого перемешивания, чтобы сосредоточиться на генерации чисел.

Для тестирования 64-битной генерации у нас есть аналогичный тест, но будет непрактично выполнять тест, соответствующий перемешиванию массива размером 2⁶⁴ – 1 (потому что на выполнение этого более масштабного бенчмарка потребуется много тысяч лет). Вместо этого мы пересечём весь 64-битный интервал, но сгенерируем то же количество выходных значений, что и в 32-битном тесте. Код:

for (uint32_t i = 0xffffffff; i > 0; --i) {
    uint64_t bound = (uint64_t(i)<<32) | i;
    uint64_t bval = bounded_rand(rng, bound );
    assert(bval < bound);
    sum += bval;
}

Результаты вихря Мерсенна

Показанные ниже результаты демонстрируют производительность этого бенчмарка для каждого из рассмотренных нами методов при использовании вихря Мерсенна и тестировании на рассмотренном в статье 32-битном (при помощи std::mt19937 из libstdc++) и аналогичном 64-битном коде (при помощи std:mt19937_64 из libstdc++). Результаты — это геометрическое среднее 15 прогонов с разными значениями seed, которое затем нормализовано, чтобы классический метод остатка от деления имел единичное время выполнения.

Может показаться, что у нас есть чёткие ответы о производительности — похоже, можно выстроить методики по их совершенству, и задаться вопросом, о чём думали разработчики libstdc++, когда писали столь ужасную реализацию для 32-битных чисел. Но, как это часто бывает при бенчмаркинге, ситуация сложнее, чем кажется по этим результатам. Во-первых, есть риск того, что результаты могут быть специфичными для вихря Мерсенна, поэтому мы расширим множество тестируемых ГПСЧ. Во-вторых, может существовать малозаметная проблема с самим бенчмарком. Давайте сначала разберёмся с первым вопросом.

Результаты разных ГПСЧ

32-битные ГПСЧ мы протестируем с помощью arc4_rand32, chacha8r, gjrand32, jsf32, mt19937, pcg32, pcg32_fast, sfc32, splitmix32, xoroshiro64+, xorshift*64/32, xoshiro128+ и xoshiro128**, а 64-битные ГПСЧ — с помощью gjrand64, jsf64, mcg128, mcg128_fast, mt19937_64, pcg64, pcg64_fast, sfc64, splitmix64, xoroshiro128+, xorshift*128/64, xoshiro256+ и xoshiro256*. Эти наборы дадут нам несколько медленных ГПСЧ и большое количество очень быстрых.

Вот результаты:

Мы можем увидеть ключевые отличия от результатов с вихрем Мерсенна. Более быстрые ГПСЧ сдвигают равновесие в сторону ограничивающего кода, а потому разница между разными подходами становится более явной, особенно в случае 64-битных ГПСЧ. При более широком наборе ГПСЧ реализация libstc++ перестаёт казаться такой ужасной.

Выводы

В этом бенчмарке со значительным отрывом выигрывает в скорости подход на основе умножения с перекосом. Есть множество ситуаций, в которых границы будут маленькими относительно размера ГПСЧ, а производительность абсолютно критична. В таких ситуациях незначительный перекос скорее всего не окажет заметного влияния, зато скорость ГПСЧ окажет. Одним из таких примеров является Quicksort со случайной опорной точкой. Из методов без перекосов многообещающей выглядит техника битовых масок.

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

Бенчмарк Small-Shuffle

Этот бенчмарк похож на предыдущий, но выполняет гораздо меньшее «перемешивание массива» (многократное). Код:

for (uint32_t j = 0; j < 0xffff; ++j) {
    for (uint32_t i = 0xffff; i > 0; --i) {
        uint32_t bval = bounded_rand(rng, i);
        assert(bval < i);
        sum += bval;
    }
}

Результаты вихря Мерсенна

Результаты разных ГПСЧ

Выводы

Этот бенчмарк избегает слишком большого упора на большие границы и более точно отражает реальные случаи применения, но теперь полностью отбрасывает большие границы.

Бенчмарк для всех интервалов

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

for (uint32_t bit = 1; bit != 0; bit <<= 1) {
    for (uint32_t i = 0; i < 0x1000000; ++i) {
        uint32_t bound = bit | (i & (bit - 1));
        uint32_t bval = bounded_rand(rng, bound);
        assert(bval < bound);
        sum += bval;
    }
}

Результаты вихря Мерсенна

Результаты разных ГПСЧ

Выводы

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

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

Вносим улучшения

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

Более быстрое отбрасывание на основе порогов

В некоторых наших алгоритмах есть код с использованием порогового значения, например:

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    // calculates 2**32 % range
    uint32_t t = (-range) % range;
    for (;;) {
        uint32_t r = rng();
        if (r >= t)
            return r % range;
    }
}

Когда range мал по сравнению с выходным интервалом ГПСЧ, чаще всего число будет намного больше порога. То есть если мы сможем добавить предварительную оценку порога, который может быть немного больше, то сэкономим на затратной операции взятия остатка деления.

С этой задачей справляется такой код:

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    uint32_t r = rng();
    if (r < range) {
        uint32_t t = (-range) % range;
        while (r < t)
            r = rng();
    }
    return r % range;
}

Это изменение можно применить и к «двойному Mod без перекосов» (см. выше), и к «целочисленному умножению без перекосов». Идею придумал Лемир, который применил её ко второму методу (но не к первому).

Результаты бенчмарка Large-Shuffle

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

Результаты бенчмарка Small-Shuffle

С другой стороны, это изменение значительно ускоряет бенчмарк small-shuffle и для метода умножения целых чисел, и для метода двойного остатка от деления. В обоих случаях из производительность сдвигается ближе к результатам вариантов без перекосов. Производительность метода двойного остатка (OpenBSD) теперь почти равна показателям метода одиночного остатка (Java).

Результаты бенчмарков для всех интервалов

Похожее улучшение мы видим и в бенчмарке для всех интервалов.

Похоже, что мы можем объявить нового всеобщего победителя: оптимизированный метод умножения целых чисел Лемира без перекосов.

Оптимизация остатка от деления

Обычно для вычисления a % b требуется деление, но в ситуациях, когда a < b результатом будет просто a, а деление не требуется. А когда a/2 < b, результат равен просто a - b. Следовательно, вместо вычислений

a %= b;

мы можем выполнить

if (a >= b) {
    a -= b;
    if (a >= b) 
        a %= b;
}

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

Результаты бенчмарка Large-Shuffle

Добавление этой оптимизации значительно улучшает результаты бенчмарка large-shuffle. Это снова более заметно в 64-битном коде, где операция взятия остатка затратнее. В методе с двойным остатком (в стиле OpenBSD) показываны версии с оптимизацией только для одной операции остатка и для обеих.

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

Результаты бенчмарка Small-Shuffle

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

Результаты бенчмарка для всех интервалов

В бенчмарке для всех интервалов изменения тоже малы.

Бонус: результаты сравнения ГПСЧ

Основная причина использования множества ГПСЧ для тестирования схем чисел в интервалах заключалась в том, чтобы избежать непреднамеренного искажения результатов вследствие особенностей работы отдельных схем ГПСЧ. Но мы можем использовать те же результаты внутренних тестов и для сравнения самих схем генерации.

ГПСЧ с выводом 32-битных чисел

График ниже показывает производительность разных 32-битных схем генерации, усреднённых для всех методов и пятнадцати прогонов, нормализованные по производительности 32-битного вихря Мерсенна:

С одной стороны, я рад видеть, что pcg32_fast и в самом деле быстр — его победил только небольшой вариант Xoroshiro (который не проходит статистические тесты). Но это показывает и то, почему я редко расстраиваюсь из-за производительности современных высокопроизводительных ГПСЧ общего назначения — разница между разными методами очень незначительна. В частности, самые быстрые четыре схемы отличаются по производительности меньше чем на 5%, и я считаю, что это просто вызвано «шумом».

ГПСЧ с выводом 64-битных чисел

На графике показана производительность разных 64-битных схем генерации, усреднённых среди всех техник и пятнадцати прогонов, нормализованных по производительности 32-битного вихря Мерсенна. Может показаться странным, что нормализация выполняется по 32-битному вихрю Мерсенна, но это позволяет нам увидеть дополнительные затраты от использования 64-битной генерации в случаях, когда достаточно 32-битной генерации.

Эти результаты подтверждают, что mcg128_fast невероятно быстр, но последние четыре техники снова отличаются всего примерно на 5%, поэтому из самых быстрых методов выбирать сложно. pcg64 и pcg64_fast обязаны быть медленнее mcg128_fast, потому что их базовыми генераторами являются 128-битный линейные конгруэнтные генераторы (ЛКГ, LCG) и 128-битные мультипликативные конгруэнтные генераторы (МКГ, MCG). Несмотря на то, что они не являются самыми быстрыми техниками в этом наборе, pcg64 всё равно на 20% быстрее 64-битного вихря Мерсенна.

Но возможно более важно то, что эти результаты также показывают, что если вам не нужен 64-битный вывод, то 64-битный ГПСЧ обычно медленнее, чем 32-битный.

Выводы

Из наших бенчмарков мы можем увидеть, что переход от стандартно используемых ГПСЧ (например, 32-битного вихря Мерсенна) к более быстрым ГПСЧ снизило время выполнения бенчмарков на 45%. Но переход от стандартного метода поиска числа в интервале к нашему самому быстрому методу позволило снизить время бенчмарка примерно на 66%; другими словами, до одной трети от исходного времени.

Самый быстрый метод (без перекосов) — это метод Лемира (с моей дополнительной оптимизацией). Вот он:

uint32_t bounded_rand(rng_t& rng, uint32_t range) {
    uint32_t x = rng();
    uint64_t m = uint64_t(x) * uint64_t(range);
    uint32_t l = uint32_t(m);
    if (l < range) {
        uint32_t t = -range;
        if (t >= range) {
            t -= range;
            if (t >= range) 
                t %= range;
        }
        while (l < t) {
            x = rng();
            m = uint64_t(x) * uint64_t(range);
            l = uint32_t(m);
        }
    }
    return m >> 32;
}

Использование метода Лемира улучшит производительность большинства рандомизированных алгоритмов больше, чем переход от быстрого движка генерации к работающему немного быстрее.

Приложения: примечания по тестированию

Код всех тестов выложен на GitHub. В целом я протестировал 23 метода для bounded_rand с помощью 26 разных ГПСЧ (13 32-битных ГПСЧ и 13 64-битных), в двух компиляторах (GCC 8 и LLVM 6), что дало мне 26 * 23 * 2 = 1196 исполняемых файлов, каждый из которых выполнялся с одинаковыми 15 seed, что даёт 1196 * 15 = 17 940 уникальных прогонов тестов, в каждом из которых объединено три бенчмарка. В основном я прогонял тесты на 48-ядерной машине с четырьмя процессорами Xeon E7-4830v3 с частотой 2,1 ГГц. Выполнение полного набора тестов потребовало чуть меньше месяца процессорного времени.

В конце вернёмся к ситуации из введения статьи. Представим, что Саша использовал jsf32.STD-libc++, а Хуан — mt19937.BIASED_FP_MULT_SCALE. В бенчмарке 3, последний занимает на 69,6% меньше времени. То есть время из этой вымышленной ситуации основано на данных из реальности.

Автор: PatientZero

Источник