Орел или решка. C++ Edition / forpes.ru

Главная
Орел или решка. C++ Edition

Орел или решка. C++ Edition +6

27.09.2025 06:11

AskePit 5 787 Источник

Джон Эдмунд Керрич известен, как человек, подкинувший монетку 10 000 раз, чтобы проверить закон больших чисел на практике. Достоверно неизвестно, сколько времени у него занял подобный утомительный эксперимент, но обстоятельства ему так или иначе помогли: он был в немецком плену, и — как иронично или издевательски это не прозвучит — у него было много свободного времени. К тому же у Джона был помощник Эрик Кристенсен, который, вероятно, удвоил скорость проведения эксперимента.

Эта история заставила меня задуматься о том, какой тяжелой ценой и какими временными затратами в докомпьютерную эпоху ученым давались некоторые эксперименты. Этой статьей я хочу отдать дань уважения таким упорным мужам, как Керрич, и показать, насколько далеко мы продвинулись в прогрессе. И как это можно круто замесить с C++, оптимизацией, многопоточностью и щепоткой ненормального программирования.

Чего я хочу

В результате эксперимента Керрича орел и решка выпали 5067 и 4933 раза соответственно. Нетрудные подсчеты показывают, что это 50.67% против 49.33%. То есть результаты приближаются к теоретическим 50%, хотя все еще имеют видимые отклонения. Я бы сказал, очень видимые отклонения. Все, конечно, зависит от вашей жадности и от точности, которую вы хотите получить. Кому-то и 50.67% может показаться достаточно близким к 50%, а эксперимент исчерпывающим и полностью законченным.

Я же глядя на результаты эксперимента, задумался, до какого предела можно его довести, пусть и не в виде манипуляций с реальной монеткой. Меня больше интересовал вопрос, в какой же момент мы действительно опасно приблизимся к 50%, и сколько времени нам на это понадобится? Вопрос, конечно, скорее провокационный, ну, или в лучшем случае наивный, нежели серьезный, потому что в нем видны зияющие дыры:

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

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

Первичная идея

Первое, что мне захотелось — подбросить монетку столько раз, сколько мне позволит софт и железо, но чтобы не сильно упороться в процессе. Например, подбросить монетку столько раз, сколько умещается в 32-битом беззнаковом целом. А еще лучше — 64-битном. В тупую проитерироваться по всему диапазону чисел, в процессе кидать монетку, записывать результат, считать затраченное время.

Я буду проводить эксперимент в следующей среде: Windows, x64, Visual Studio 2022, C++20.

Поехали.

Тщательно готовимся к эксперименту

Будем заходить издалека.

Для начала мне нужно было написать игрушечный фреймворк, на базе которого я смог бы гибко проводить все расчеты:

Многократно выполнять кусок полезной работы
Замерять время
Выводить его на экран в удобоворимом виде

Все три вида работ я в идеале видел, как матрешку вида:

measure( spin(doWork) )
// осторожно, псевдокод! здесь не уважаются настоящий синтаксис и семантика языка

Иными словами: мы хотим измерить и отобразить (measure) многократно выполненную (spin) работу (doWork). Бонусом будет тотальное обобщение подхода, чтобы я при необходимости смог раздеребанить этот код на другие свои повседневные нужды.

Замеряем время

Я хочу функцию measure , которая примет в качестве параметра интересующий нас коллбек, запустит его и замерит время его исполнения. В самом базовом виде я вижу ее так:

using namespace std::chrono;
using Clock = high_resolution_clock;

template<typename Func>
void measure(Func&& f)
{
    auto start = Clock::now();

    f();

    auto dt = Clock::now() - start;
    std::cout << dt << std::endl;
}

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

void doWork() {
    std::cout << "FFFUUU " << std::endl;
}

measure(doWork);

В результате в консоль выведет:

FFFUUU
36100ns

Это действительно работает, но у текущей версии measure есть проблемы:

Сигнатура вызываемой функции не подразумевает никаких аргументов, а возвращаемое функцией значение не извлечь наружу из вызова measure(...)
Точность часов high_resolution_clock в моем окружении — наносекунды, поэтому когда мы пытаемся распечатать любой замер, он выведется в наносекундах, сколь большим не получился бы результат. Это не читаемо

Сначала разберемся с первой проблемой. Ее можно решить, приправив measure неплохой такой щепоткой всех самых современных (а может и не очень) фишек C++:

template<typename Func, typename... Args>
auto measure(Func&& f, Args... args) -> decltype(f(std::forward<Args>(args)...))
{
    auto start = Clock::now();

    if constexpr (std::is_same_v<decltype(f()), void>) {
        f(std::forward<Args>(args)...);
        auto dt = Clock::now() - start;
        std::cout << dt << std::endl;
    } else {
        auto res = f(std::forward<Args>(args)...);
        auto dt = Clock::now() - start;
        std::cout << dt << std::endl;
        return res;
    }
}

void echo(std::string_view s) {
    std::cout << s << std::endl;
}

int sum(int a, int b) {
    return a + b;
}

int main()
{
    measure(std::bind(echo, "hey!"));
    std::cout << std::endl;
    
    int ab = measure(std::bind(sum, 12, 2000));
    std::cout << ab << std::endl;
    return 0;
}

Вывод на экран:

hey!
151100ns

200ns
2012

Вариадик шаблоны, decltype и constexpr слились в едином экстазе на крошечной площади, чтобы measure превратилась в многоликого джокера, который может возвращать значение, а может и не возвращать; может принимать одни аргументы, а может другие.

Фактически мы сделали из measure подобие декоратора из Python. Можно еще завернуть в лямбду повторяющийся кусок кода:

template<typename Func, typename... Args>
auto measure(Func&& f, Args... args) -> decltype(f(std::forward<Args>(args)...))
{
    auto start = Clock::now();

    auto summarize = [&start]() {
        auto dt = Clock::now() - start;
        std::cout << dt << std::endl;
    };

    if constexpr (std::is_same_v<decltype(f()), void>) {
        f(std::forward<Args>(args)...);
        summarize();
    } else {
        auto res = f(std::forward<Args>(args)...);
        summarize();
        return res;
    }
}

Я думаю, я остановлюсь на таком варианте.

Теперь к проблеме вывода времени. У STL нет единого решения для "умного" вывода duration. Если ваш duration<> — это seconds, он вам и выведет количество целочисленных секунд, а хвоста попросту не будет, т.е. 1.5 секунды превратятся в 1s. Если ваш duration<> — это milliseconds, покажет ms и т.д. Смотрите сами, что умеет стандартная библиотека:

using namespace std::chrono_literals;

auto d = 987'654'321'978ns;

std::cout << d << std::endl;
std::cout << duration_cast<microseconds>(d) << std::endl;
std::cout << duration_cast<milliseconds>(d) << std::endl;
std::cout << duration_cast<seconds>(d) << std::endl;
std::cout << duration_cast<minutes>(d) << std::endl;

Вывод на экран:

200000000000ns
200000000us
200000ms
200s
3min

Примечание: Такой вывод времени, который вы можете наблюдать выше, доступен только с C++20. Еще в C++17 стандарте std::cout << d вообще бы не скомпилировался — только так: std::cout << d.count(). И вывод будет без суффиксов ms, s, min и пр. — только голое число.

Я призадумался, а как бы я вообще хотел видеть время замеров. Как визуально адекватно представить 200 000 000 000 наносекунд? Поначалу я склонялся к варианту, где мы поочередно кастим duration, в часы, минуты, секунды, миллисекунды и т.п. до тех пор, пока не получим результат, больший нуля. Тогда 200 000 000 000 наносекунд превратились бы в 3min. И вроде бы ничего, и хвост вроде бы не особо важен при таком автоскейле.. или важен? Я колебался, а потом решил, что если разложить длительность на все-все составляющие вплоть до наносекунд — например, 2min 36s 649ms 551us 558ns — то я получу ультимативность, читаемость, точность и отсутствие сомнений. Поэтому реализовал именно этот вариант:

template <typename Dur>
void prettyPrintDuration(Dur dur)
{
    auto h = duration_cast<hours>(dur);
    if (h.count()) { std::cout << h << " "; dur -= h; }

    auto m = duration_cast<minutes>(dur);
    if (m.count()) { std::cout << m << " "; dur -= m; }

    auto s = duration_cast<seconds>(dur);
    if (s.count()) { std::cout << s << " "; dur -= s; }

    auto ms = duration_cast<milliseconds>(dur);
    if (ms.count()) { std::cout << ms << " "; dur -= ms; }

    auto us = duration_cast<microseconds>(dur);
    if (us.count()) { std::cout << us << " "; dur -= us; }

    auto ns = duration_cast<nanoseconds>(dur);
    if (ns.count()) { std::cout << ns << " "; dur -= ns; }

    std::cout << std::endl;
}

Тестируем:

using namespace std::chrono_literals;

auto d = 987'654'321'978ns;
prettyPrintDuration(d);

Получаем:

16min 27s 654ms 321us 978ns

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

Контролируемое кручение

Еще одна функция общего назначения, которая нам понадобится: spin. Она должна выполнять одну и ту же работу заданное количество раз.

В простейшем виде это может быть так:

template<typename Func>
void spin(Func&& f, size_t count)
{
    for (size_t i = 0; i<count; ++i) {
        f();
    }
}

Все очень просто. Но так как моя идея базировалась на том, что я хочу пройтись по всему диапазону беззнакового целого, я хотел отдать все заботы по вычислению границ диапазона так же на откуп spin. При этом на текущем этапе я не знаю, буду ли использовать uint8_t, uint16_t, uint32_t или uint64_t в качестве типа. Знаю, что uint8_t точно понадобится для быстрых тестов — 256 итераций имеют шансы целиком уместиться на экран и выполнить мгновенно любую разумно потребляющую работу.

Окей, достаточно пробросить желаемый тип в качестве шаблонного аргумента:

template<std::unsigned_integral T, typename Func>
void spin(Func&& f)
{
    for (T i = 0; i < std::numeric_limits<T>::max(); ++i) {
        f();
    }
}

Но это ловушка! Если вы запустите spin<uint8_t>(...), он прокрутит числа от 0 до 254 — на одну итерацию меньше, чем возможный диапазон [0; 255]. Исправляем < на <=:

template<std::unsigned_integral T, typename Func>
void spin(Func&& f)
{
    for (T i = 0; i <= std::numeric_limits<T>::max(); ++i) {
        f();
    }
}

Но это ловушка! Поздравляю, вы в бесконечном цикле — ведь любое число всегда <= своего максимума. Вот так и оказалось, что красивой записи для итерации по всему диапазону как бы и не существует. Нет, конечно, если вы дадите переменной i тип uintmax_t, то для 8, 16 и (возможно) 32 бит все сработает, как задумано. Но когда мы дойдем до 64 бита, чем uintmax_t на многих платформах и является, нас ждет все та же западня.

В итоге я пришел к раскоряке:

template<std::unsigned_integral T, typename Func>
void spin(Func&& f)
{
    for (T i = 0; ; ++i) {
        f();
        if (i == std::numeric_limits<T>::max()) break;
    }
}

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

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

void doWork() {
    std::cout << "F";
}

int main()
{
    measure([](){ spin<uint8_t>(doWork); });
    return 0;
}

И получаем на экран 256 () отборных F:

FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF6ms 342us 500ns

На этом подготовительное отклонение от главной темы статьи закончено — начинаем кидать монетку!

Бросаем монетку

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

std::mt19937 gen(std::random_device{}());
std::uniform_int_distribution uni(0, 1);

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

using BigInt = uint8_t;
std::mt19937 gen(std::random_device{}());
std::uniform_int_distribution uni(0, 1);
BigInt heads = 0;
BigInt tails = 0;

void doWork() {
    int res = uni(gen);
    if (res == 1) {
        ++heads;
    } else if (res == 0) {
        ++tails;
    } else {
        assert(false);
    }
}

int main()
{
    measure([](){ spin<BigInt>(doWork); });

    const double headsPercent = heads / static_cast<double>(std::numeric_limits<BigInt>::max()) * 100.0;
    const double tailsPercent = tails / static_cast<double>(std::numeric_limits<BigInt>::max()) * 100.0;

    std::cout << "Heads: " << +heads << ", " << headsPercent << "%" << std::endl;
    std::cout << "Tails: " << +tails << ", " << tailsPercent << "%" << std::endl;

    return 0;
}

Примечание: заметьте — я написал +heads и +tails при выводе на экран. Это не случайно, и нужно для типа uint8_t, который иначе интерпретируется стримом std::cout как тип char, и на экран выводятся ASCII-кракозябры вместо чисел. + позволяет скастить символ к числу и избежать этого.

В качестве нашего "большого числа" — BigInt — тестово берем uint8_t. Подбрасываем монетку нашей measure/spin крутилкой, подсчитываем проценты, и получаем:

2us 700ns
Heads: 122, 47.8431%
Tails: 134, 52.549%

Ура, вот мы и провели первый эксперимент по подбрасыванию монетки, и уложились всего в 2.7 микросекунды! Да, 256 это вам не 10 000 Керрича, поэтому проценты разбросало сильно. Но в целом все выглядит нормально. По крайней мере я так думал. Однако же стоит быть бдительным — пока я игрался с кодом, в какой-то момент случайно получил вот такой идеальный результат, где орлы и решки совпали:

2us 700ns
Heads: 128, 50.1961%
Tails: 128, 50.1961%

и понял, что мои вычисления грешат. Как вы можете увидеть, сумма вероятностей 50.1961% и 50.1961% равна чему-то чуть большему, чем 100%.

Причиной было деление на std::numeric_limits<BigInt>::max(). А фактически делить надо на одно число больше. Да только где ж его взять, если это максимум? По итогу я решил не морочиться, а приводить это число к double и уже потом к нему прибавлять единицу и делить. Благо double резиновый, хоть и не гарантирует нам возможности точно представить желаемое число на больших величинах, так что нужно держать в голове, что решение не точное. Но на первых порах, нам сойдет и такое:

const double headsPercent = heads / (static_cast<double>(std::numeric_limits<BigInt>::max()) + 1.0) * 100.0;
const double tailsPercent = tails / (static_cast<double>(std::numeric_limits<BigInt>::max()) + 1.0) * 100.0;

Громоздко, но у нас тут C++. Зато через десяток попыток я снова попал на 128/128:

2us 304ns
Heads: 128, 50%
Tails: 128, 50%

И проценты уже посчитаны верно.

Кстати, вопрос к математикам: а какова вероятность получить такую вероятность? :) Я в короткий промежуток времени получил 128 орлов против 128 решек целых два раза! Это при том, что при бóльших величинах такое совпадение, наверное, едва ли возможно, и не встречалось мне более нигде за пределами 256 подбрасываний монетки.

Теперь давайте побьем рекорд ученого, подбросив монетку более 10 000 раз. Тип uint16_t дает нам возможность провести уже 65 536 итераций — как раз то, что нужно. Меняем BigInt на uint16_t и получаем:

569us 400ns
Heads: 32685, 49.8734%
Tails: 32851, 50.1266%

Что же — наши 50.1266% уже имеют меньшую погрешность, чем 50.67% Керрича.

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

За полвека до того, как 10 000 раз монету подбросил Джон Эдмунд Керрич, другой ученый Карл Пирсон сделал то же самое вместе со своими подопечными. И они провели еще больше итераций — сумасшедшие 24 000 подкидывания монетки с результатами 12 012 против 11 988, или 50.05% против 49.95%.

Мой генератор 65 536 подкидываний монетки совершает в два с половиной раза больше итераций подкидывания, но при этом процент чаще находится в районе [50.1%; 50.2%], что "хуже" результатов Пирсона, хотя иногда и спускается и к его результатам, а иногда и несущественно ниже. Из этого можно сделать осторожный вывод, что 24 000 или 65 000 подбрасываний не существенно отличаются друг от друга в плане статистического результата.

Либо Пирсон врал.

Либо мой генератор врет. Либо я вру.

Изолируем эксперимент

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

Я не хотел каждый раз менять тип BigInt, перекомпилировать программу и запускать эксперимент. Ведь мы за автоматизацию. К тому же меня посетила мысль, что было бы неплохо видеть сводную таблицу с прогрессией: чтобы последовательно были произведены эксперименты для 8-, 16-, 32- и 64-битных диапазонов, все было посчитано и запротоколировано.

Но для такого множественного эксперимента наш код совершенно не готов. Чего нам не хватает?:

Нам нужно уметь итерироваться по [uint8_t, uint16_t, uint32_t, uint64_t]. Смекаете? — нам нужно уметь итерироваться по типам. Уже невесело
Нам нужно инкапсулировать всю логику и данные эксперимента в класс, чтобы этот код можно было переиспользовать повторно

Хорошо, начнем с того, что полегче, и объявим давно напрашивающийся класс Experiment:

template<typename BigInt>
struct Experiment
{
    std::mt19937 gen {std::random_device{}()};
    std::uniform_int_distribution<int> uni {0, 1};
    BigInt heads = 0;
    BigInt tails = 0;

    void tossCoin() {
        int res = uni(gen);
        if (res == 1) {
            ++heads;
        }
        else if (res == 0) {
            ++tails;
        }
        else {
            assert(false);
        }
    }
};

Ну, это было слишком просто. Поэтому перейдем к сложной части.

Итерация по типам — это что-то на шаблонном, на метапрограммировании, да? Вроде того. Давайте мыслить поступательно: мы хотим что-то вроде

for (type : Types) {
    f<type>();
}
// осторожно, псевдокод! здесь не уважаются настоящий синтаксис и семантика языка

где Types — список типов, а f — шаблонная функция. Вот только это псевдокод.

А как получить C++? Есть variadic templates, распаковка которых может как бы "генерировать" код в виде виртуального списка подобно старым добрым макросам. Т.е. задачу можно свести к подобной:

template<typename ...Types, typename Func>
void forEachType(Types..., Func f)
{
    (f<Types>(), ...);
}
// осторожно, псевдокод! здесь не уважаются настоящий синтаксис и семантика языка

Чтобы вы не ломали голову, сразу подскажу, что запись (f<Types>(), ...); должна теоретически раскрыться в

(f<uint8_t>(), f<uint16_t>(), f<uint32_t>(), f<uint64_t>());

т.е. в comma-operator, который перечислит запись для каждого типа и последовательно выполнит каждую. Только не забываем, что мы все еще работаем с псевдокодом, и пока что у нас не компилируется!

Но наш последний псевдокод был недалек от правды, и если его немного обработать напильником под C++, то оно заведется:

template<typename... Ts, typename F>
void forEachType(F&& f, Ts...)
{
    (f.template operator()<Ts> (), ...);
}

int main()
{
	forEachType([]<typename T>() {
	    std::cout << +std::numeric_limits<T>::max() << std::endl;
	}, uint8_t{}, uint16_t{}, uint32_t{}, uint64_t{});
}

Только не спрашивайте у меня ничего про template operator(), умоляю, я сам едва понимаю его необходимость, но работает только с ним.

Вывод на экран:

255
65535
4294967295
18446744073709551615

Скажу сразу, что мне эта реализация не нравится по ряду причин:

Список типов можно писать только в конце списка аргументов, и никак иначе. Вы не сможете объявить forEachType как void forEachType(Ts..., F&& f) — попытка вызвать такую функцию вызовет ошибку компиляции
Приходится делать фейковые переменные uint32_t{}, чтобы вызов функции был синтаксически корректен. На каждый тип по фейковому объекту

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

template<typename... Ts>
struct TypeList {};

Теперь variadic template спрячется за объект TypeList, а снаружи — это цельный единичный объект, и его можно передвинуть вперед в forEachType. Вот так работает:

template<typename... Ts>
struct TypeList {};

template<typename... Ts, typename F>
void forEachType(TypeList<Ts...>, F&& f)
{
    (f.template operator()<Ts> (), ...);
}

int main()
{
	forEachType(TypeList<uint8_t, uint16_t, uint32_t, uint64_t>{}, []<typename T>() {
	    std::cout << +std::numeric_limits<T>::max() << std::endl;
	});
}

Да, мы все еще создаем фейковый объект — TypeList<>{}, но он всего один, так что я считаю это хорошим компромиссом. Получился еще один неплохой трюк, который можно сохранить себе в другие проекты.

Первые результаты

Теперь мы можем исполнить задуманное и запустить эксперимент с итерациями разной размерности:

8 bits
8us 200ns
Heads: 120, 46.875%
Tails: 136, 53.125%

16 bits
571us 600ns
Heads: 32735, 49.9496%
Tails: 32801, 50.0504%

32 bits
35s 839ms 225us 500ns
Heads: 2147474281, 49.9998%
Tails: 2147493015, 50.0002%

64 bits

Что же, 4 миллиарда бросков из 32-битного домена спустя 35 секунд уже дают неплохой результат: 50.0002% против 49.9998%. Это существенно лучше экспериментов господ Пирсена и Керрича.

Вам может показаться, что я недокопипастил часть текста и по недосмотру не вставил сюда результат выполнения 64 bits, но нет. Те, кто еще главу назад почуял неладное, можете ликовать, вы оказались правы — выполнение 64-битного диапазона может занять очень большое время.

Я не стал сдаваться и решил подождать подольше — оставил программу выполняться на ночь. По утру я застал консоль за тем же не изменившимся экраном и решил подождать еще чуточку дольше. И только когда это не помогло, и к вечеру экран все так же показывал пустоту ниже строки 64 bits — вот только в этот момент я стал что-то подозревать.

Посмотрим на количество чисел, умещаемых в разных битовых разрядностях:

 8 bit: 256
16 bit: 65536
32 bit: 4294967296
64 bit: 18446744073709551616

Если пристально посмотреть, можно увидеть, что последовательность растет квадратично, и 18 446 744 073 709 551 616 — это ни что иное как 4 294 967 296 * 4 294 967 296. Иными словами, если подбрасывание монетки в 32-битном домене заняло у меня 35 секунд, то 64-битный домен закончит исполнение примерно через время в 4 294 967 296 раза дольше..

35 * 4 294 967 296 секунд =
150 323 855 360 секунд =
2 505 397 589.33 минут =
41 756 626.48 часов =
1 739 859.44 суток = 
4 766.74 года

О как.

Дата пришествия Иисуса ближе к нам, чем дата завершения эксперимента.

В общем, на прохождение по полноразмерному 64-битному диапазону рассчитывать, судя по всему, нельзя. Даже если представить, что подбрасывание монетки в 32-битном домене займет всего 1 миллисекунду, то 64-битный домен должен отработать за примерно

4 294 967 296 (четыре миллиарда) миллисекунд =
4 294 967.296 (четыре миллиона) секунд =
71 582.79 минут = 
1 193.05 часов =
49.7 суток

Полтора месяца! Что конечно не четыре тысячи лет, но и 32-битное исполнение мы, наверное, не приблизим к 1 миллисекунде.

На этот моменте стало ясно, что текущий план нужно пересмотреть, поскольку 32-битный диапазон для нас слишком мал, а 64-битный непостижимо огромен, и насколько бы много мы в итоге раз не подбросили монетку, это количество будет больше 4 294 967 296 (четырех миллиардов), но существенно меньше 18 446 744 073 709 551 616 (восемнадцати квинтиллионов).

Держа в уме эти цифры, я внес некоторые корректировки в курс, который держит эксперимент:

Мы отказываемся от концепции итерирования по диапазонам 8-, 16-, 32-, 64-битных чисел — теперь наш рабочий диапазон это 64 бита, и мы можем работать в его рамках, не особо волнуясь за переполнение. А итерирование по типам, которое мы городили буквально только что, уходит на выброс
Мы берем uintmax_t, он же uint64_t на подавляющем большинстве платформ, и делаем из него универсальный тип using BigInt = uintmax_t для дальнейшей работы в рамках эксперимента
Самое важное: мы берем курс на жесточайшую оптимизацию исполнения программы, чтобы достичь максимального количества подброшенных монет на единицу времени. C++ goes hard.

Оптимизации

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

Для того, чтобы итеративно улучшать производительность, нужны объективные замеры. Если очередная правка кода выдает более быстрый результат (при условии нормальной работоспособности алгоритма), мы засчитываем ее и идем дальше.

Одним из ключевых вопросов остается, как добиться объективных замеров, и я вам сразу скажу, что, пожалуй, не стану заморачиваться с тем, чтобы бенчмарк был стерильно идеальным. Я клоню к тому, что всегда найдутся люди, которые скажут: "ты неправильно мерял". Например, потому что замеры происходили только на одной машине; или потому что я несерьезно отнесся к проблеме "(не)прогретого" кэша; или потому что в какой-то момент прирост был рамках статистической погрешности; или кому-то покажется, что итераций было слишком мало, и сам эксперимент стоит воспроизвести 10 000 раз и усреднить для уверенности в результатах. Ну и далее по списку из любых других замечаний, которые непременно всплывают, когда кто-то хочет устроить бенчмарк.

Поэтому оговорим условия замеров:

Одна машина: i5-12600K, 3700 Mhz, 16 логических процессоров, RAM 64Gb (но в память мы не упремся)
Я просто использую результаты своей функции measure в качестве замера
Я не делаю тысячи повторных прогонов эксперимента, поскольку сама работа по подбрасыванию монетки — и так сплошной многократный, многомиллиардный прогон одной и той же операции. Я рассудил, что этого достаточно
Иногда я вручную запускал бенчмарк по нескольку раз, просто чтобы убедиться, что результат замера стабильный. Эту часть я не автоматизировал
Я в курсе про Google Benchmark. Но я здесь, чтобы трогать код руками и донимать свой мозг, а не интегрировать в мой бесподобный код third party решения и убивать веселье

Еще один немаловажный момент: чтобы результаты измерений были корректно сравниваемы друг с другом, стоит изначально определиться с количеством подбросов монетки, время которого мы будем измерять. До применения оптимизаций было трудно угадать подходящее число. Можно было подумать, что 4 294 967 296 — хороший кандидат, поскольку сейчас он выполняется за 35 секунд, т.е. достаточно долго, и есть смысл начать с него. Но если мы будем слишком хороши в оптимизациях, это число может резко упасть до жалких милли- или микросекунд, а на таких величинах результаты двух бенчмарков уже могут выдавать погрешность, в рамках которой трудно будет разглядеть истинную разницу.

Я смухлюю задним числом, поскольку, как вы понимаете, статья пишется уже после доведения эксперимента до финала, и я знаю, что нас ждет в конце. Поэтому я возьму за эталон 68 719 476 736 итераций. Это 4 294 967 296 x 16, т.е. шестнадцать прогонов 32-битного диапазона чисел.

Почему я не перейду на круглые и красивые числа, в конце концов? — спросите вы. Ведь если больше не нужны все эти $2^{32}$ , $2^{64}$ , почему бы не делать 10 000 000 000, 50 000 000 000 итераций и т.д.? Первый ответ будет: я забылся, и не подумал об этом! — даже после упразднения первоначальной идеи я продолжил просто домножать 4 294 967 296 на все большие и большие величины. Деформация мозга, не иначе. Второй ответ: нам это пригодится в дальнейшем. В последующем в процессе оптимизаций у нас появится ограничение: количество подбрасываний монетки должно быть кратно 64, а потом это требование возрастет еще больше, но все так же будет крепко завязано на степени двойки. А с моими зубодробильными 4 294 967 296 это условие выполняется автоматически.

Зафиксируем, что имеем

Зафиксируем основной код до оптимизаций, чтобы у вас была более-менее цельная картина:

template<typename Func>
void spin(BigInt n, Func&& f)
{
    for (BigInt i = 0; i < n; ++i) {
        f();
    }
}

struct Experiment
{
    std::mt19937 gen {std::random_device{}()};
    std::uniform_int_distribution<int> uni {0, 1};
    BigInt heads = 0;
    BigInt tails = 0;

    void tossCoin() {
        int res = uni(gen);
        if (res == 1) {
            ++heads;
        }
        else if (res == 0) {
            ++tails;
        }
        else {
            assert(false);
        }
    }
};

int main()
{
	constexpr BigInt N = 68'719'476'736ll;

	Experiment experiment;
	measure([&experiment, N]() {
	    spin(N,
	        std::bind(&Experiment::tossCoin, &experiment)
	    );
	});
}

Время исполнения этого кода:

Time:  20min 41s 897ms 383us 600ns
Heads: 34 359 699 204, 49.9999%
Tails: 34 359 777 532, 50.0001%

Мучительно долгие 20 минут 42 секунды на подбрасывание монетки 68 миллиардов раз. Это наша стартовая точка.

Обратите внимание на вывод консоли и оцените, как я успел подсуетиться и стал отображать большие числа в отформатированном варианте: 34 359 699 204 вместо 34359699204. Повставлять пробелы в строку несложно, а читаемость на числах такого порядка значительно, значительно улучшается.

Как оптимизировать

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

Расчленяем рандом

Сначала я решил заняться вот этими ребятами:

std::mt19937 gen {std::random_device{}()};
std::uniform_int_distribution<int> uni {0, 1};

Они были слишком подозрительными, и их было слишком много.

Сначала я попробовал пересесть на сишную классику: srand + rand. Это оказалось катастрофой — было медленно.

Затем я попробовал использовать другие классы распределений. Некоторые из них давали абсолютно некорректные результаты, не сводившиеся к 50/50. А некоторые работали и даже были чуть быстрее std::uniform_int_distribution, например std::bernoulli_distribution. К тому же это распределение как раз выдавало bool. Но выхлоп был минимальный.

И я задался вопросом: а нужен ли мне вообще distribution, и что будет если я просто его выброшу? По сути нам не так много-то и надо — получить ноль или единицу. Сам генератор — std::mt19937 и множество его товарищей — вполне могут самостоятельно генерировать рандомные числа — у них есть для этого operator(). Загвоздка в том, что мы получим какое-то случайное число (его тип зависит от класса генератора), и с ним нужно будет дальше что-то делать, чтобы привести его к диапазону [0; 1]. Я наивно рассудил, что если просто взять нулевой бит этого числа, то должно получиться что-то вполне случайное. Оказалось, что не со всеми генераторами это хорошо работает, но тем не менее такие генераторы, как std::minstd_rand и std::mt19937 продолжали выдавать распределение, стремящееся к 50/50, поэтому я без сомнений выкинул std::uniform_int_distribution и внес вот такую правку в код броска монеты:

-int res = uni(gen);
+int res = gen() & 1;

Результат:

Time:  10min 50s 93ms 886us 900ns
Heads: 34 359 738 368, 50%
Tails: 34 359 738 368, 50%

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

Кто-то скажет, что я пожертвовал правильным рандомом ради оптимизации. Но, если честно, мне все равно на академическую "правильность" рандома до тех пор, пока я вижу адекватную картину при отрабатывании программы: увеличивающееся количество подбрасываний монет приближает соотношение орлов и решек к 50/50. В конце концов, когда Керрич подбрасывал монету, его распределение тоже не был идеальным, поскольку не была идеальной его монета, не была идеальной и его рука — полагаю, индивидуальная манера подбрасывать монету теоретически тоже может влиять на исход.

Фиксируем прирост x1.9 от удаления uniform-распределения.

Группируем биты

Я уже говорил, что gen() генерирует случайный целочисленный результат, который мы приводим к одному биту. Но зачем же нам так расточительствовать, если можно задействовать все биты числа, которое отдал нам генератор? std::mt19937 за один вызов генерирует 32-битное число, и если его задействовать по полной, это будет равносильно 32 подбрасываниям монетки сразу.

Единственная проблема с этой идеей — класс Experiment подбрасывает одну монетку за раз, а внешняя функция spin запускает N вызовов Experiment::tossCoin, и если нам хочется оставаться в рамках этой парадигмы, нам нужно как-то извертеться и вызывать генератор раз в 32 итерации, а остальное время брать результат подбрасывания из кэша. У меня вышло неуклюже:

struct Experiment
{
    Experiment()
    {
        regenBits();
    }

    using Generator = std::mt19937;
    using BitsType = Generator::result_type;

    BigInt heads = 0;
    BigInt tails = 0;

    Generator gen {std::random_device{}()};
    BitsType bits{};
    int bitsCounter = 0;

    void tossCoin() {
        if (yieldBit()) {
            ++heads;
        } else {
            ++tails;
        }
    }

    bool yieldBit() {
        if (bitsCounter >= std::numeric_limits<BitsType>::digits) {
            regenBits();
        }

        return (bits >> bitsCounter++) & 1;
    }

    void regenBits() {
        bitsCounter = 0;
        bits = gen();
    }
};

но сильно заморачиваться я не стал и замерял, как есть:

Time:  7min 8s 600ms 933us
Heads: 34 359 585 551, 49.9998%
Tails: 34 359 891 185, 50.0002%

Еще минус три минуты и прирост скорости в x1.5. Это тоже очень хороший результат. Особенно радует, что этот код есть куда оптимизировать по части эффективности "кручения" функции подброса монетки.

Реструктуризируем кручение

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

Еще нам нужно как-то эффективно "расчленить" число на биты и посчитать количество нулей и единиц. На деле задачу можно свести к подсчету только единиц или только нулей, а вторую половину подсчитать вычитанием из 32. Пока что произведем подсчет "руками" через битовые операции:

template <std::unsigned_integral T>
unsigned countOnes(T val) {
    unsigned n = 0;
    for (int i = 0; i < std::numeric_limits<T>::digits; ++i) {
        n += val & 1;
        val >>= 1;
    }
    return n;
}

Кто-то заметил бы, что алгоритм можно оптимизировать, если выходить из цикла, как только val достигает нуля. Например, сделать вот так:

while (val) {
	n += val & 1;
	val >>= 1;
}

И это звучит логично с точки зрения алгоритмов и здравой логики, но работает в два раза медленнее на моей машине. Ставлю на то, что ~~гигачад~~-for пользуется благами SIMD и loop unrolling, а чимс while в машинном коде так и остается щупать val на ноль, и бранчинг мешается CPU под ногами снова и снова.

В своем новом обличии и с дополнительными полномочиями класс Experiment стал выглядеть так:

struct Experiment
{
    using Generator = std::mt19937;
    using BitsType = Generator::result_type;

    BigInt n = 0;
    BigInt heads = 0;
    BigInt tails = 0;

    Generator gen {std::random_device{}()};
    static constexpr size_t BITS_COUNT = std::numeric_limits<BitsType>::digits;

    Experiment(BigInt n_)
        : n(n_) {
    }

    void spin() {
        for (BigInt i = 0; i < n; i += BITS_COUNT) {
            BitsType bits = gen();
            heads += countOnes(bits);
        }

        tails = n - heads;
    }
};

Его вызов в бенчмарке:

Experiment experiment(n);
measure(std::bind(&Experiment::spin, &experiment));

Заметим BITS_COUNT, который неявно требует, чтобы количество подбрасывания монетки было кратно 32 битам.

Замеряем время:

Time:  19s 826ms 601us
Heads: 34 359 613 380, 49.9998%
Tails: 34 359 863 356, 50.0002%

Сначала я думал, что меня приглючило, но нет — вместо минут ожидания, то же количество подбрасываний теперь уложилось в двадцать секунд! Это колоссальный прирост скорости в x21.6 раз — нечто феноменальное.

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

Считаем биты правильно

На самом деле мне сразу было ясно, что подсчет битов-единиц в числе можно выполнить оптимальнее: например, использовать интринсики или платформозависимые решения. Я же, как истинный Cpp-enjoyer, решил наслаждаться благами C++20 и его функцией std::popcount, которая выдает количество бит-единичек для любого беззнакового целого.

-heads += countOnes(bits);
+heads += std::popcount(bits);

Time:  5s 460ms 958us 600ns
Heads: 34 359 743 765, 50%
Tails: 34 359 732 971, 50%

x3.63. И это неудивительно: вызов std::popcount() разворачивается в одну машинную инструкцию POPCNT, и с этим решением по эффективности потягаться просто невозможно.

Расширяем генератор

Я стал внимательно рассматривать всю линейку генераторов, которую нам предоставляет STL. Вот они все:

std::minstd_rand0
std::minstd_rand
std::mt19937
std::mt19937_64
std::ranlux24_base
std::ranlux48_base
std::ranlux24
std::ranlux48
std::knuth_b

Оказалось, что некоторые генераторы выдают 64-битные числа, а это потенциальная заявка на новую двухкратную оптимизацию! Недолго думая, делаем однострочную замену:

-using Generator = std::mt19937;
+using Generator = std::mt19937_64;

и получаем улучшенный результат:

Time:  2s 554ms 51us 200ns
Heads: 34 359 573 276, 49.9998%
Tails: 34 359 903 460, 50.0002%

Теперь мы одновременно подбрасываем по 64 монетки вместо 32 и закономерно получаем соизмеримый выхлоп: x2.14. Так же заметим, что теперь количество подбрасываний должно быть кратно 64.

Так же ради приличия давайте пропишем этот инвариант в явном виде, чтобы мы наверняка знали, если нарушим его:

Experiment(BigInt n_)
	: n(n_) {
	if (n % (sizeof(BigInt) * CHAR_BIT) != 0) {
		std::cerr << "Error: n must be a multiple of 64, got " << n << std::endl;
		std::abort();
	}
}

Автоматическая многопоточность

В C++17 язык обзавелся обновлением заголовочника <algorithm>: появилась возможность запускать функции std::for_each, std::sort, std::find и пр. в параллельном режиме. Нужно всего лишь передать алгоритмической функции первым параметром std::execution::par_unseq (остальные политики для нашей задачи имеют мало смысла):

std::sort(std::execution::par_unseq, v.begin(), v.end());

И вектор волшебным образом отсортируется с задействованием всех ваших ядер. На самом деле это не точно — STL сама решает, каким именно образом она производит распараллеливание, сколько потоков будет использовано и будет ли вообще — вы лишь можете указать рекомендацию в виде execution policy.

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

std::vector<int> v(1000, 1);
int sum = 0;

std::for_each(std::execution::par_unseq, v.begin(), v.end(),
	[&](int x) {
		sum += x; // data race!
	}
);

Представьте, что STL создала 1000 потоков (мы не знаем, сколько в действительности), которые одновременно стали инкрементить sum. Эта переменная не является атомиком, поэтому у вас проблемы, и результат едва ли будет равен 1000.

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

Сперва я попробовал ворваться с ноги в этот параллельный мир написанием наивного std::for_each закрыв глаза на гонки данных. Зачем? Просто чтобы быстро пощупать скорость выполнения алгоритма без примитивов синхронизации или атомиков, чтобы знать, на какую скорость выполнения я могу рассчитывать в пределе (спойлер: все не настолько просто, как я предполагал, и так не работает). Код:

void spin() {
	BigInt steps = static_cast<size_t>(n / BITS_COUNT);
	std::vector<BigInt> indices(steps);
	std::iota(indices.begin(), indices.end(), 0);

	std::for_each(std::execution::par_unseq,
		indices.begin(), indices.end(),
		[this](size_t idx) {
			BitsType bits = gen();        // BUG: data race on random generator
			heads += std::popcount(bits); // BUG: non-atomic access from several threads
		}
	);

	tails = n - heads;
}

Результаты его выполнения:

Time:  11s 199ms 665us 700ns
Heads: 6 069 669 426, 8.83253%
Tails: 62 649 807 310, 91.1675%

Сразу же обращаем внимание на очевидную некорректность подсчета — тотальная несходимость к 50/50. Я списываю это на потоконебезопасное изменение heads. Что меня расстроило больше: я получил 11 секунд, что в четыре раза хуже однопоточного исполнения. Но как же так?

Чисто ради интереса я продолжил копаться в параллельных алгоритмах и рассудил, что могу переписать мой std::for_each на std::transform + std::reduce чтобы убрать гонку за heads, а гонку за генератор устранить за счет того, что на каждый поток будем иметь свою копию:

void spin() {
	BigInt steps = static_cast<size_t>(n / BITS_COUNT);
	std::vector<BigInt> res(steps);

	std::transform(std::execution::par_unseq,
		res.begin(), res.end(),
		res.begin(),
		[this](BitsType) {
			thread_local Generator gen{std::random_device{}()};
			BitsType bits = gen();
			return std::popcount(bits);
		}); 

	heads = std::reduce(std::execution::par_unseq, res.begin(), res.end());
	tails = n - heads;
}

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

Результаты выполнения:

Time:  2s 250ms 10us 700ns
Heads: 34 359 867 682, 50.0002%
Tails: 34 359 609 054, 49.9998%

Т.е. время выполнения все же не улучшилось в сравнении с нашим однопоточным вариантом. Но и не ухудшилось. А должно было улучшиться.

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

Ручная многопоточность

Разочаровавшись в параллельных алгоритмах STL, мы переходим к задаче разбить подбросы монеток на вручную созданные потоки. Благо, наша задача достаточно изолирована и примитивна, и не требует зависимостей или общих ресурсов. Кроме генератора, конечно, но мы уже клонировали генератор для разных потоков в предыдущей главе и видели, что рандом при этом не ломается. Мне снова могут сказать, что я занимаюсь откровенной хулой, и генератор по-хорошему должен быть один на всю задачу; а я снова отвечу, что покуда рандом ведет себя, как должен, то я считаю его работающим.

Вкратце алгоритм действий:

Определяемся с количеством потоков. Обычно выбирают std::thread::hardware_concurrency() потоков, и это дает оптимальную загрузку логических ядер
Разбиваем количество подбрасываний монеты, которые предстоит сделать, на количество потоков
Одновременно запускаем все потоки. Каждый поток локально посчитает свою часть результата и вернет ее
Результаты всех потоков складываются вместе и получается итоговый результат переменной heads

Вот реализация:

void spin() {
	const BigInt steps = n / BITS_COUNT;
	const size_t threadsCount = std::max(std::thread::hardware_concurrency(), 1u);
	const BigInt chunkSize = steps / threadsCount;

	const auto thread = [this, chunkSize]() -> BigInt {
		Generator gen {std::random_device{}()};
		BigInt localHeads = 0;
		for (BigInt i = 0; i < chunkSize; ++i) {
			localHeads += std::popcount(gen());
		}
		return localHeads;
	};

	std::vector<std::future<BigInt>> threadHeads(threadsCount - 1);
	for (auto& f : threadHeads) {
		f = std::async(thread);
	}

	heads += thread();
	for (auto&& fut : threadHeads) {
		heads += fut.get();
	}

	tails = n - heads;
}

Из интересного:

Я запустил N-1 потоков и отдельно загрузил главный поток приложения. Не, ну а че ему простаивать
Я использовал std::async, а не std::thread, поскольку std::async умеет возвращать результат, к тому же вызов std::async() без параметров std::launch::async или std::launch::defered работает в автоматическом режиме — создает потоки, пока их количество остается приемлемым; перестает создавать новые потоки и отрабатывает в главном, если мы переборщили с числом запрашиваемых потоков.

Время выполнения:

Time:  126ms 479us 900ns
Heads: 12 884 816 423, 49.9997%
Tails: 12 884 987 353, 50.0003%

Вот он — выхлоп от многопоточности. Фиксируем еще один гигантский прирост x20.2!

Тут стоит оговориться, что прирост будет сильно зависеть от машины. У меня 16 логических ядер. К тому же большую роль играет локальный рандом-генератор. На первый взгляд 16 копий std::mt19937_64 может показаться расточительством — тем более, что каждый такой генератор весит по 4.9Кб памяти! Да, объекты этого класса такие огромные. Но если бы генератор был общим, мы не выжали бы из многопоточности так много — даже если использовать общий генератор без защиты (или притвориться, что он потокобезопасный), потоки бы постоянно инвалидировали друг другу кэш-линии, где лежит генератор. Я полагаю, именно потому реализация с std::for_each из предыдущей главы выдавала 11 секунд, в то время как std::transform + std::reduce выдавала 2 секунды. Последний алгоритм не имел общих данных, и потоки не "сливали" друг другу кэш.

Так же добавлю, что помимо std::async я пробовал использовать и std::thread в связке с std::promise + std::future и в связке с std::packaged_task, но результат получался примерно таким же по производительности, а бойлерплейта выходило больше, поэтому остановимся на std::async, как на самом интерфейсно удобном решении для нашей задачи.

Свой генератор

Со времен, когда я учился писать шейдеры, я помнил, как часто в них приходилось на коленке писать псевдо-рандомный генератор при помощи примитивных математических операций умножения, сдвига, сложения. И это отлично работало, поскольку рандом ощущался вполне себе рандомным — а что еще было нужно?

Я задумался, насколько реалистично заменить тяжелый, но качественный std::mt19937_64 на рукописное поделие. Через некоторое время я имел на руках вот такого кандидата на замену:

struct LCG64 {
    using result_type = BigInt;
    uint64_t state;

    LCG64(uint64_t seed = 1) : state(seed) {}

    uint64_t operator()() {
        state = state * 6364136223846793005 + 1;
        return state;
    }
};

Пускай вас не смущает его примитивность! Я попробовал погенерировать им некоторое количество чисел и отобразить их в двоичном виде. Для отображения, кстати, отлично подходит трюк с std::bitset:

LCG64 gen{std::random_device{}()};
uint64_t bits = gen();
std::cout << std::bitset<64>(bits) << std::endl;

Вот вывод некоторого количества 64-битных чисел, которые выдал мне генератор:

0101100010010111001110101000000011111000110110011000010111110100
0100001110101001110010100110001100001110101101001001011111100101
1101001000111110110110110010101010010001011000100100111001000010
0000101001110111000110111010100001110111100001000111111110011011
1111111000100110100000000010011100000110110011110101001101000000
1011110101100110111111101011110110011110000101110000100111100110
1110010000001000101010001000100101000010110100111101011101101111
0011001110001110001011000000110000101110101101111110111110000100
1000100000110000011111111010111100001101111110111001011000110101
0100100000010001000010010000110001101101100101101011001001010010
1001010010100101110001111001110100100111101011100000011001101011
1101110011111010111011000100111111001111000011010011010111010000
0001001010101111010001011111100100000101000101001010010110010001
1011100010110000100011100011000101001111111010100111011000100101

Не знаю, как вам, а мне распределение нулей и единиц кажется абсолютно случайным, без паттернов, и вроде бы даже сводимым к 50/50. В любом случае это легко проверить. Мы всего лишь меняем одну строку нашего кода:

-using Generator = std::mt19937_64;
+using Generator = LCG64;

и запускаем бенчмарк:

68 719 476 736 rounds
Time:  93ms 529us 800ns
Heads: 34 359 798 049, 50.0001%
Tails: 34 359 678 687, 49.9999%

И ничего не сломалось. Судя по раскладу 50.0001% / 49.9999% наш рукописный генератор работает, как должен. К тому же он работает быстрее предыдущего, и мы получаем еще x1.4 в копилку. Я, правда, ожидал немного больше прироста, но, видимо, все-таки std::mt19937_64 достаточно эффективно выполняет свою работу.

Ещё ядрёнее

Я был достаточно твердо уверен, что std::thread::hardware_concurrency() потоков на задачу — это максимально комфортный режим, а дальше — лишь деградация и увядание. Я обнаружил, что не прав, чисто случайно, решив от безделия и безысходности увеличить количество потоков в два раза. Потом в четыре.. потом в восемь.. мои брови поднимались кратно множителю и зафиксировались на std::thread::hardware_concurrency() * 64 — здесь был sweet spot с результатом:

Time:  75ms 91us 300ns
Heads: 34 359 767 603, 50%
Tails: 34 359 709 133, 50%

А это еще приятные x1.25 к скорости.

Как это вообще может работать? На машине с 8 физическими ядрами, 16 логическими ядрами лучше всего отрабатывает 1024 потока! Странно. А как же постоянные переключения контекста, битвы за кэш и прочее? Я могу предположить здесь несколько факторов:

Особенности работы операционной системы с потоками. В данном случае мы на Windows. Какие особенности, понятия не имею, но смею предположить, что у ОС могут быть какие-то хитрые алгоритмы для оптимальной работы множества потоков
Особенности железа. Возможно, что-то подшаманивает и сама архитектура
Особенности работы std::async, который строго говоря может и не создать поток, если мы не просим его об этом явно (а мы не просим), а может просто поставить асинхронное выполнение в главном потоке. К тому же наши потоки очень быстрые и короткоживущие. Возможно, std::async как-то грамотно их ставит в очередь или хз, что он там делает

Это лишь мои предположения, в каком направлении может скрываться истина. Прошу помощи зала в объяснении.

SIMD

Может показаться, что SIMD ничем не поможет в нашей задаче: SIMD-инструкции плохо работают в "горизонтальной" плоскости — когда нужно что-то сделать внутри одного SIMD-объекта. SIMD создан, чтобы очень эффективно проводить операции векторов друг с другом. К тому же, мы по сути уперлись в производительность буквально одной машинной инструкции POPCNT, которую достаточно трудно переплюнуть. Горячий путь нашего алгоритма находится вот здесь:

    for (BigInt i = 0; i < chunkSize; ++i) {
        BitsType bits = gen();
->      localHeads += std::popcount(bits);  // пичот
    }

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

Запрос мой был лаконичен (только так с железными мешками и надо): у меня есть горстка uint64_t-чисел, я хочу провести на них массовый SIMD-attack, чтобы посчитать все их биты-единички; результат выдать в виде итоговой суммы бит-единиц во всех числах.

Во-первых, меня сразу обрадовали, что у AVX-512 есть буквально SIMD-popcount — инструкция _mm512_popcnt_epi64. И даже предложили функцию-обертку, которую можно удобно использовать:

uint64_t avx512_popcount_sum(const uint64_t* data, size_t n) {
    __m512i acc = _mm512_setzero_si512();
    size_t i = 0;
    for (; i + 8 <= n; i += 8) {
        __m512i v = _mm512_loadu_si512(&data[i]);
        __m512i pc = _mm512_popcnt_epi64(v);
        acc = _mm512_add_epi64(acc, pc);
    }
    // horizontal reduce
    uint64_t result = _mm512_reduce_add_epi64(acc);
    // tail
    for (; i < n; i++)
        result += std::popcount(data[i]);
    return result;
}

Код даже не особо страшный и длинный, как часто бывает с SIMD.

Я очень радостный пошел включать AVX-512 в проекте Visual Studio, машинально пересобрал старую версию программы (еще без использования SIMD-инструкций, и без предложенной функции), запустился. Почти с порога получил краш с формулировкой "illegal instruction". Это означало две вещи:

Компилятор умудрился нагенерить AVX-512 инструкций там, где я даже хз, что они вообще были возможны — мое почтение. Я же ведь еще даже не успел скопировать код от ChatGPT! Компилятор молодец. Упало тут:

const double headsPercent = static_cast<double>(experiment.heads) / n * 100.0;
const double tailsPercent = static_cast<double>(experiment.tails) / n * 100.0;
  
// Вот где он тут планировал применить AVX-512? В ассемблер я не полез, сорян

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

Я, конечно расстроился от такого поворота, ведь этот элитный SIMD мог фактически выполнить восемь POPCNT параллельно. Вряд ли бы это был честный x8 прирост, но какой-то бы да был.

В итоге я остался наедине с AVX2 — он-то мне доступен, это я знаю, пользовался. Проблема в том, что у него даже нет какого-то аналога _mm512_popcnt_epi64. Что есть — так это кем-то придуманный алгоритм, как можно, располагая ~~нехитрыми приспособлениями и буханкой белого (или черного) хлеба~~ лишь доступными нам инструкциями, зашаффлить биты до смерти так, чтобы на выходе получилось количество бит-единиц со всех четырех поданных на вход 64-битных чисел. То есть, мало того, что мы можем рассчитывать только на четырех-кратный прирост, так мы еще и не можем на него рассчитывать, поскольку решение будет заведомо не оптимальным с точки зрения SIMD-операций. Добавим сюда накладные расходы на перегонку чисел в SIMD и обратно, и получаем что-то на грани "то на то и выйдет". Но не проверишь, не узнаешь. Я решил попробовать.

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

uint64_t avx2Popcount(const uint64_t* data, size_t n) {
    static const __m256i lut = _mm256_setr_epi8(
        0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
        0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4
    );

    __m256i total = _mm256_setzero_si256();

    size_t i = 0;
    for (; i + 4 <= n; i += 4) {
        __m256i v = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(&data[i]));

        __m256i lo = _mm256_and_si256(v, _mm256_set1_epi8(0x0F));
        __m256i hi = _mm256_and_si256(_mm256_srli_epi16(v, 4), _mm256_set1_epi8(0x0F));

        __m256i cnt = _mm256_add_epi8(
            _mm256_shuffle_epi8(lut, lo),
            _mm256_shuffle_epi8(lut, hi)
        );

        total = _mm256_add_epi64(total, _mm256_sad_epu8(cnt, _mm256_setzero_si256()));
    }

    // Reduce total (4×64 → 1×64)
    __m128i low128  = _mm256_castsi256_si128(total);
    __m128i high128 = _mm256_extracti128_si256(total, 1);
    __m128i sum128  = _mm_add_epi64(low128, high128);
    uint64_t result = _mm_cvtsi128_si64(sum128) + static_cast<uint64_t>(_mm_extract_epi64(sum128, 1));

    // Tail
    for (; i < n; i++)
        result += std::popcount(data[i]);

    return result;
}

Отвратительно.

Но нам еще только предстоит подготовить наш код, чтобы хотя бы иметь возможность воспользоваться этой божественной функцией. Нам нужно уметь нагенерировать определенное количество uint64_t-чисел, чтобы подать их на вход функции. Вопрос — сколько? На самом деле на этот вопрос никто не сможет ответить однозначно, придется это количество подбирать/настраивать экспериментально. Объявляем для этого глобальную константу:

static constexpr size_t SIMD_BATCH = 4;

Пускай пока побудет 4. Фактически это всего на один раунд __m256, и скорее всего этого мало, но начнем хоть с чего-то. Константу мы еще успеем порегулировать, когда все напишем.

Вспоминаем, как каждый поток генерирует и обрабатывает числа сейчас:

Generator gen {std::random_device{}()};
BigInt localHeads = 0;
for (BigInt i = 0; i < chunkSize; ++i) {
	localHeads += std::popcount(gen());
}

Т.е. генерируем очередное число, сразу же считаем его биты. Теперь надо генерировать SIMD_BATCH чисел сразу, потом отдавать их в avx2Popcount. Придется выделить немного памяти под массив с числами, и чуть подшаманить с циклом:

Generator gen{ std::random_device{}() };
BigInt localHeads = 0;
BitsType bits[SIMD_BATCH];

for (BigInt i = 0; i < chunkSize / SIMD_BATCH; ++i) {
    for (int j = 0; j < SIMD_BATCH; ++j) {
        bits[j] = gen();
    }
    localHeads += avx2Popcount(&bits[0]);
}

Так же имеем в виду, что теперь мы должны гарантировать, что количество подбрасываний монетки не просто кратно 64, как было раньше, а кратно SIMD_BATCH * 64, поскольку именно столько бит (подбрасываний монетки) мы генерируем за один SIMD-раунд. Не кратное количество мы строго не хотим, поскольку обработка остаточных хвостов, которых не хватает на SIMD — это дополнительные траты, которые нам ни к чему. Проще уложить количество итераций в требования. Укладываем:

static constexpr size_t SIMD_BATCH_BITS = SIMD_BATCH * sizeof(BigInt) * CHAR_BIT;

...

Experiment(BigInt n_)
    : n(n_) {
    if (n % SIMD_BATCH_BITS != 0) {
        std::cerr << "Error: n must be a multiple of " << SIMD_BATCH_BITS << ", got " << n << std::endl;
        std::abort();
    }
}

А еще не забудем, что и сам SIMD_BATCH должен быть кратен четырем, поскольку мы работаем с __m256, так что где-нибудь впишем static_assert:

static_assert(SIMD_BATCH % 4 == 0, "shouldda be multiple of four!");

Теперь можно подбирать SIMD_BATCH, не боясь получить белиберду — если что, на нас поругаются, и мы исправимся. Я подобрал оптимальное для себя значение 32, с которым получил результат:

Time:  60ms 906us 200ns
Heads: 34 359 807 145, 50.0001%
Tails: 34 359 669 591, 49.9999%

x1.23. Это конечно, не в 4 раза быстрее, но по крайней мере я не вышел в ноль после переписывания кода.

Итоги оптимизаций

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

Веха	Время	Прирост
Начало замеров	`20min 41s 897ms 383us 600ns`	-
Без distribution	`10min 50s 93ms 886us 900ns`	x1.9
Ленивый генератор	`7min 8s 600ms 933us`	x1.5
Реструктуризация	`19s 826ms 601us`	x21.6
std::popcount()	`5s 460ms 958us 600ns`	x3.6
64-bit генератор	`2s 554ms 51us 200ns`	x2.1
Многопоточность	`126ms 479us 900ns`	x20.2
Свой генератор	`93ms 529us 800ns`	x1.4
Больше потоков	`75ms 91us 300ns`	x1.25
SIMD	`60ms 906us 200ns`	x1.23

Итоговое ускорение: x20 235! На самом деле вместо таблицы должен был быть показательный график, но ускорение программы получилось настолько масштабным, что график, только начавшись, просто падал на ось X из-за слишком большого диапазона по Y. Но я все же попытался выкрутиться и представляю вам свое ноу-хау — график с вложенными диапазонами:

В общем, не знаю как еще показать масштабность ускорения программы в двадцать тысяч раз — результат слишком хорош :) Читаем график так — как только график падает на ось X, переходите на следующий график, который детализирует ось Y. Если плохо видно надписи, на основном графике диапазон [0; 1400] секунд, на графике второго уровня [0; 20] секунд, третьего — [0; 0.14] секунд.

Еще полезно вычислить скорость нашего автоматизированного подкидывателя монеток. 68 719 476 736 подбрасывания в ~61 миллисекунду, что эквивалентно:

68 719 476 736 / 61ms
1 128 283 766 / ms
1 128 283 / us
1 128 / ns

Что же, 1к/наносек — это, конечно, не 300к/наносек, но все равно круто.

Что еще могло бы ускорить?

Ухх, что я только не пытался сделать еще:

Выравнивал данные потока по alignas(std::hardware_destructive_interference_size), чтобы они были выровнены по кэш-линии, и были минимизированы риски потокам мешать друг другу
Делал разные вариации ручного loop-unrolling
Упразднял генератор и вычислял рандомное число прямо внутри буффера BitsType bits[SIMD_BATCH]: каждый следующий элемент вычислялся из предыдущего, потом закольцовывался. Это вроде бы как локализовывало данные

Но это все не меняло скорость выполнения. На данном этапе я думаю, что мои полномочия по оптимизации все. Но у меня есть планы на будущее:

AVX-512 точно улучшит выполнение программы. Когда я дорвусь до машины, которая поддерживает этот набор инструкций, я обязательно поделюсь с вами результатами. Но, боюсь, это произойдет совсем нескоро. Когда-нибудь в 2030
Я хочу переписать программу на Rust. Мне любопытно, будет ли быстрее моей C++ версии. Еще у Rust есть u128, но он вряд ли что-то даст, поскольку как и __128 из Clang/GCC (в MSVC он отсутствует) — это лишь софтварные симуляции, а нативный размер в любом случае — 64 бит, и вы никуда от него не денетесь
Я хочу переписать программу на Zig. Уж сравнивать системные языки, так сравнивать
Когда я научусь писать compute shaders, я стану жечь GPU вместо CPU, и тоже надеюсь на прорывные результаты

Если что-то из этого случится, я напишу отдельную статью об этом.

Большие и малые числа

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

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

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

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

В чем его проблема? В том, что это число с плавающей точкой, и даже если бы оно было 128- или даже 256-битным, то мы не можем ручаться, какие числа оно сможет представить точно, а какие нет. Убедиться в этом достаточно легко, в свое время я написал IEEE-754 калькулятор, используя который можно увидеть, что даже 0.1 в форме числа с плавающей точкой не возможно представить как есть. float по-факту будет хранить что-то вроде 0.10000000149011611938, 64-битный double — 0.10000000000000000555 и т.д.

Чтобы считать точно, нам нужны fixed-point числа, так же часто называемые decimal. Суть их проста: под капотом хранится целое число — количество долей. Один из самых типовых случаев использования decimal-чисел — деньги. Под капотом храним все в копейках, наружу притворяемся, что имеем дело с рублями и его долями. Количество цифр после точки у таких денег — 2. И это количество фиксированное, потому и fixed-point. Вот нам нужен свой такой fixed-point, плюсовый STL таким не располагает.

Мы будем реализовывать decimal под нашу крайне узкую задачу:

Наш decimal будет представлять числа, близкие к 0.5. Например, 0.5000053 или 0.49999645
Мы хотим хранить как можно больше цифр после точки
Мы вообще не будем пользоваться цифрами перед точкой. То есть нам не нужно уметь представлять 1.0 или 2.45
отрицательные числа так же не нужны

Примечание: ранее мы использовали проценты для представления распределения. Теперь мы от этого уходим и вместо 50% будем использовать 0.5. Это избавит нас от лишней головной боли с конвертацией и унифицирует идею малого числа.

Под decimal нам логично выделить uint64_t под капотом. Так как все числа будут меньше единицы, мы можем отдать все цифры под дробную часть. Чтобы задействовать весь uint64_t по полной, можно обратиться к справочной константе std::numeric_limits<uint64_t>::digits10. Это максимальное количество десятичных цифр, которыми можно представить число и не бояться нарваться на overflow. Для uint64_t эта константа равна 19. И действительно, максимальное число для uint64_t — это 18 446 744 073 709 551 615. У него 20 цифр, но, как вы понимаете, не любое число с 20 цифрами поместится в 64-бита. Поэтому 19 — это безопасный вариант, и даже 9 999 999 999 999 999 999 вмещается без угрозы все сломать. Следовательно, точность нашего decimal мы можем сделать в 19 цифр после точки. Я думаю, нам хватит. А если не хватит, всегда можно сэмулировать 128 бит.

Реализовываем костяк класса Decimal:

class Decimal
{
public:
    Decimal(size_t digitsAfterComma, uint64_t underlying_ = 0)
       : underlying(underlying_)
    {
        if (digitsAfterComma > std::numeric_limits<uint64_t>::digits10) {
            std::cerr << "Error: you've exceeded limit of digits after fixed point" << std::endl;
            std::abort();
        }
        scale = static_cast<size_t>(std::pow(10, digitsAfterComma));
    }

    std::string toString() const {
        std::string res("0.");
        res += std::to_string(underlying);
        while (res[res.size() - 1] == '0') {
            res.resize(res.size() - 1);
        }

        return res;
    }

private:
    uint64_t underlying{};
    size_t scale{};
};

Публичный интерфейс почти отсутствует, но мы будем его дописывать по мере продвижения вперед и понимания, что конкретно нам нужно. Заметьте, что мы сделали возможность создавать числа любой точности и ограничили digitsAfterComma максимальной планкой в 19 цифр. Нам еще понадобятся числа разной точности, но об этом позже.

Трудности деления

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

const double headsPercent = static_cast<double>(experiment.heads) / n * 100.0;
const double tailsPercent = static_cast<double>(experiment.tails) / n * 100.0;

Возьмем наш случай с 68 719 476 736 подбрасываниями, которые мы так долго оптимизировали. Представьте, что у нас было 34 359 743 765 орлов и 34 359 732 971 решек. Я посчитал на виндовом калькуляторе — соотношение получается таким: 0.5000000785366864875 / 0.4999999214633135125.

Кстати говоря, это числа из главы "Считаем биты правильно", и double там уже захлебнулся и отобразил просто 50%/50%. А возможно, так сработало форматирование double стримом. В любом случае, мы желаем большего.

В общем, чтобы вычислить желаемый Decimal, нужно поделить 34 359 743 765 на 68 719 476 736. Хм. А чем делить, а как делить?

Делить через каст в double, затем как-то конвертировать в Decimal? Бессмысленно — мы растеряем всю точность и похороним суть Decimal. Сформировать Decimal из этих больших чисел? Так они не влазят в наш Decimal, у которого все цифры отданы на дробную часть.

В какой-то момент я сообразил, что с точки зрения математики неважно, делите вы 34359743765 / 68719476736 или 0.34359743765 / 0.68719476736 — результат одинаковый! То есть мы можем эти числа затащить в наш целевой Decimal просто как Decimal(19, 34'359'743'765), и ничего страшного не произойдет. Останется поделить Decimal на Decimal.

Как правильно делят Decimal числа? Очевидно, вы не можете просто взять и поделить underlying переменные — в нашем случае мы таким макаром получим 34359743765 / 68719476736 = 0. Можно смухлевать, скастовать оба underlying в double и произвести деление. И снова убить весь смысл Decimal. В общем, правильное деление делается хитрее: оно целочисленное, но предварительно первое число домножают на scale, потом целочисленно делят, затем делят результат на scale обратно. Т.е.:

343597437650000000000000000000 / 68719476736 = 5000000785366864874

5000000785366864874 / 10000000000000000000 = 0.5000000785366864874

Хм. Ндаа, scale для чисел с 19 цифрами после точки это действительно $10^{19}$ , что как мы помним, само по себе является числом, очень близким к пределу возможностей uint64_t. И умножать его на другое огромное число никак нельзя — мгновенное переполнение и невозможность ничего посчитать. Обидная ситуация — нам нужно сделать непосильное умножение, чтобы тут же поделить все обратно. Еще обиднее, что результат-то в наших реалиях мы все равно получим близким к 0.5.

Я походил, поломал голову над этой задачей некоторое время. Посидел с ручкой, листочком. И знаете, вот на листочке у меня все получалось, как надо. Потому что я мог считать столбиком любые деления любых чисел, не то что эта тупая машина.. Wait, OH SHI~ Деление столбиком! Мы будем симулировать деление столбиком.

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

Я недолго думая, набросал такой генератор. Держим в уме, что он сильно заточен под наш случай — делимое всегда меньше делителя, а цифры выдаются те, которые после 0.:

class Divider
{
public:
    Divider(BigInt nom_, BigInt denom_)
        : nom(nom_)
        , denom(denom_)
    {}

    int operator()() {
        if (nom == 0) {
            return -1;
        }

        nom *= 10;
        if (nom < denom) {
            return 0;
        }

        BigInt acc = denom;
        while (acc <= nom) {
            acc += denom;
        }
        acc -= denom;

        int digit = static_cast<int>(acc / denom);
        nom -= acc;

        return digit;
    }

private:
    BigInt nom{};
    BigInt denom{};
};

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

// 1 / 12 = 0.083333...
{
	Divider d(1, 12);
	d(); // 0
	d(); // 8
	d(); // 3
	d(); // 3
	d(); // 3
	d(); // 3
}

// 1 / 2 = 0.5
{
	Divider d(1, 2);
	d(); // 5
	d(); // -1
	d(); // -1
	d(); // -1
}

Как я планирую собрать все воедино и совершить акт деления:

У нас есть два больших числа
Из них мы формируем Divider
Создаем пустой Decimal нужной нам точности
Вытаскиваем из делителя числа, пока он не истощится, либо пока Decimal не заполнится под завязку
Если делитель еще живой, выпросим у него еще одно число для правильного округления

Все это собираем в одной функции:

Decimal calcPercent(BigInt nom, BigInt denom, size_t digitsAfterComma) {
    Decimal res(digitsAfterComma);
    Divider divider(nom, denom);

    size_t scale = res.getScale() / 10;

    while (scale != 0) {
        int digit = divider();
        if (digit < 0) {
            break;
        }
        res.addFracts(digit * scale);
        scale /= 10;
    }

    // rounding
    int digit = divider();
    if (digit >= 5) {
        res.addFracts(1);
    }

    return res;
}

Видим, что классу Decimal нужно немного дополнить интерфейс: нужен метод addFracts(). Я на правах Ванги еще добавлю, что необходим метод isHalf() — он пригодится нам позже, увидите:

class Decimal
{
    ...
    void addFracts(size_t fracts) {
       underlying += fracts;
    }

    bool isHalf() const {
        return underlying == 5 * (scale/10);
    }
    ...
};

Теперь мы готовы. Не думал, что делить числа может быть так сложно. Но было весело.

Заменяем double-деление на новое, точное:

const Decimal headsPercent = calcPercent(experiment.heads, n, 19);
const Decimal tailsPercent = calcPercent(experiment.tails, n, 19);

и запускаем наш типовой оптимизационный прогон:

68 719 476 736 rounds
Time:  67ms 411us 500ns
Heads: 34 359 684 850, 0.4999992212106008083
Tails: 34 359 791 886, 0.5000007787893991917

Другое дело: это вам не жалкие 50.0001%, это точность на кончиках пальцев. Теперь можно переходить к следующему этапу.

Майним коины

Теперь, используя Decimal, я хочу получить самое точное стремление к 50/50, которое только можно получить за разумное время. Идея такая:

Стартуем с Decimal низкой точности. Например, c 1 цифры после точки.
Запускаем наш эксперимент, и крутим его бесконечно — ждем, пока не получим результат 0.5 / 0.5. С учетом (не)точности нашего Decimal, конечно же
Фиксируем результат: время, количество итераций, за которое мы освоили точность в 1 цифру после точки
Переходим к точности в 2 цифры после точки и закольцовываем эксперимент в бесконечном цикле

Получается бесконечно работающий майнер точности. Он подскажет нам, за какое время и за какое количество подкидываний монетки мы покоряем все более и более точное распределение. Когда/если замайнится точность 19, эксперимент будем считать полностью оконченным по техническим причинам.

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

struct Experiment
{
    ...
    
    BigInt heads = 0;
    BigInt tossed = 0;
    nanoseconds timePassed{};

    void spin(BigInt n) {
        ... // all the work and `heads` calculation
        
        // measurements
        tossed += n;
        auto dt = Clock::now() - start;
        timePassed += dt;
    }

    BigInt getTails() const {
        return tossed - heads;
    }
};

Пишем сырую версию майнера:

constexpr size_t MAX_PRECISION = std::numeric_limits<uint64_t>::digits10;
constexpr BigInt SPIN_STEP = 1ull << 21;

Experiment ex;

const auto reportSuccess = [&ex](size_t precision) -> bool {
	Decimal heads = calcPercent(ex.heads, ex.tossed, precision);
	Decimal tails = calcPercent(ex.getTails(), ex.tossed, precision);

	if (heads.isHalf() || tails.isHalf()) {
		std::cout << "Tossed: " << prettifyBigInt(ex.tossed) << std::endl;
		std::cout << "Heads:  " << prettifyBigInt(ex.heads) << ", " << heads.toString() << std::endl;
		std::cout << "Tails:  " << prettifyBigInt(ex.getTails()) << ", " << tails.toString() << std::endl;
		std::cout << "Time:   " << prettyDuration(ex.timePassed) << std::endl << std::endl;
		return true;
	}

	return false;
};

for (size_t precision = 1; precision <= MAX_PRECISION;	++precision) {
	std::cout << "precision: " << precision << '\n';

	while (true) {
		ex.spin(SPIN_STEP);
		if (reportSuccess(precision)) {
			break;
		}
	}
}

Первый запуск:

precision: 1
Tossed: 2 097 152
Heads:  1 048 560, 0.5
Tails:  1 048 592, 0.5
Time:   3ms 275us 200ns

precision: 2
Tossed: 4 194 304
Heads:  2 097 912, 0.5
Tails:  2 096 392, 0.5
Time:   4ms 796us 600ns

precision: 3
Tossed: 6 291 456
Heads:  3 147 050, 0.5
Tails:  3 144 406, 0.5
Time:   6ms 403us 800ns

precision: 4
Tossed: 16 777 216
Heads:  8 388 366, 0.5
Tails:  8 388 850, 0.5
Time:   13ms 117us 700ns

precision: 5
Tossed: 20 971 520
Heads:  10 485 662, 0.5
Tails:  10 485 858, 0.5
Time:   16ms 578us 700ns

precision: 6
Tossed: 1 333 788 672
Heads:  666 894 992, 0.5
Tails:  666 893 680, 0.5
Time:   711ms 941us 100ns

...

Робко предполагаем, что оно работает. Но есть проблемы.

Во-первых уж как-то непоказательно выходит, что когда мы бьем очередной рекорд и Decimal, упираясь в свой лимит точности говорит: "я — круглый 0.5!", мы выводим его на экран, и это всегда 0.5. И непонятно, зачем его тогда вообще отображать. Но вместе с тем посмотреть на "реальное" полученное число тоже хотелось бы. Я придумал такой трюк: я буду считать деление для Decimal текущей точности майнера и параллельно поделю все тоже самое еще раз, но для максимально возможной точности. Решение о покорении точности я буду принимать по "грубому" Decimal, а отображать буду "мелкозернистый" Decimal. Звучит запутанно, кодом будет проще:

Decimal headsRough = calcPercent(ex.heads, ex.tossed, precision);
Decimal tailsRough = calcPercent(ex.getTails(), ex.tossed, precision);

Decimal headsDetailed = calcPercent(ex.heads, ex.tossed, MAX_PRECISION);
Decimal tailsDetailed = calcPercent(ex.getTails(), ex.tossed, MAX_PRECISION);

if (headsRough.isHalf() || tailsRough.isHalf()) {
    // display `headsDetailed` and `tailsDetailed`
}

Вторая проблема не так хорошо видна — я ее скорее почувствовал, нежели увидел. Посмотрите на количество итераций для каждого нового рекорда — оно каждый раз разное. Я же стал подозревать, что по крайней мере на малых точностях один раунд может запросто "взять" сразу несколько рекордов за раз, перепрыгнув, например, от 0.51 до 0.5006. Решить эту, возможно, гипотетическую проблему легко — перед каждым проведением эксперимента проверять, не побил ли уже эксперимент ту точность, которую мы сейчас планируем крутить:

for (size_t precision = 1; precision <= MAX_PRECISION;	++precision) {
    std::cout << "precision: " << precision << '\n';
	
+   if (reportSuccess(precision)) {
+       continue;
+   }

    while (true) {
        ex.spin(SPIN_STEP);
        if (reportSuccess(precision)) {
            break;
        }
    }
}

На самом деле, проблема могла проявляться, а могла и не проявляться, в зависимости от выбранного шага кручения, а здесь мы уже подходим к последней проблеме. Эта проблема — подобрать SPIN_STEP оочень трудно. Я указал 1ull << 21 (это $2^{21}$ ), но это был результат утомительного поиска. $2^{20}$ уже давал черепашью скорость, и у меня не хватало тепрпения дождаться даже взятия precision: 1. $2^{22}$ уже летел слишком быстро и моментально упирался в слишком большие величины, которые было очень долго считать. Это кстати говоря показывает, что майнер не очень просто настроить на оптимальную скорость — она будет очень сильно зависеть от шага, который вы укажете. Я рассудил, что шаг должен быть адаптивным — он должен быть достаточно мелким на старте, а потом увеличиваться. Я ввел минимальный шаг и максимальный шаг:

constexpr BigInt MIN_SPIN_STEP = SIMD_BATCH_BITS;
constexpr BigInt MAX_SPIN_STEP = 4'294'967'296;
size_t step = MIN_SPIN_STEP;

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

step = experiment.tossed < step ? step : std::min(step * 2, MAX_SPIN_STEP);
experiment.spin(step);

Т.е. периодически удваиваем шаг, пока не упремся в максимум. Максимум тут чисто номинальный, на всякий случай, я взял проверенный немаленький 4 294 967 296.

Майнер стал работать, как надо. По крайней мере, как мне хочется. Вот первые 6 взятий:

precision: 1
Tossed: 33 552 384
 Heads: 15 724 557, 0.4686569216661325765
 Tails: 17 827 827, 0.5313430783338674235
  Time: 18ms 222us 700ns

precision: 2
Tossed: 268 433 408
 Heads: 133 163 234, 0.4960754884876326571
 Tails: 135 270 174, 0.5039245115123673429
  Time: 21ms 636us 200ns

precision: 3
Tossed: 2 147 481 600
 Heads: 1 072 710 984, 0.4995204540984192833
 Tails: 1 074 770 616, 0.5004795459015807167
  Time: 25ms 722us 500ns

precision: 4
Tossed: 21 474 834 432
 Heads: 10 736 419 663, 0.4999535478141562046
 Tails: 10 738 414 769, 0.5000464521858437954
  Time: 46ms 516us 900ns

precision: 5
Tossed: 227 633 264 640
 Heads: 113 815 517 221, 0.4999951013354671009
 Tails: 113 817 747 419, 0.5000048986645328991
  Time: 253ms 918us 100ns

precision: 6
Tossed: 1 026 497 181 696
 Heads: 513 248 105 522, 0.4999995272018192996
 Tails: 513 249 076 174, 0.5000004727981807004
  Time: 1s 38ms 463us 500ns

...

Работает отлично, а главное — теперь прогресс точности стал хорошо наблюдаем по растущим нулям и девяткам. Приятное зрелище.

Энергонезависимость

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

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

А потому нужно озаботиться сохранением данных в файл. Но и в консоль мы по-прежнему хотим выводить. Потому что это удобно. Выводить в оба источника нужно одни и те же записи, т.е. мы либо будем дублировать каждый вывод, либо напишем какое-то решение по автоматизированной записи сразу в два стрима. Само собой второй вариант предпочтительнее, поэтому я озаботился написанием класса StreamOverdub

class StreambufOverdub : public std::streambuf {
public:
    StreambufOverdub(std::streambuf* sb1, std::streambuf* sb2)
        : sb1(sb1), sb2(sb2) {
    }

protected:
    int overflow(int c) override {
        if (c == EOF) return !EOF;
        const int r1 = sb1->sputc(static_cast<char>(c));
        const int r2 = sb2->sputc(static_cast<char>(c));
        return (r1 == EOF || r2 == EOF) ? EOF : c;
    }

    int sync() override {
        int const r1 = sb1->pubsync();
        int const r2 = sb2->pubsync();
        return (r1 == 0 && r2 == 0) ? 0 : -1;
    }

private:
    std::streambuf* sb1;
    std::streambuf* sb2;
};

class StreamOverdub : public std::ostream {
public:
    StreamOverdub(std::ostream& o1, std::ostream& o2)
        : std::ostream(&tbuf)
        , tbuf(o1.rdbuf(), o2.rdbuf()) {
    }

private:
    StreambufOverdub tbuf;
};

Как видите, весь секрет состоит в оборачивании двух std::streambuf в прокси. И теперь, если создать SteamOverdub вот так:

std::ofstream file("mined coins.txt", std::ios::trunc);
StreamOverdub out(file, std::cout);

то мы сможем просто везде заменить std::cout на out и получить вывод в консоль и вывод файл по цене одного вывода. Также после каждого успешного взятия стоит вызвать out.flush(), чтобы файл обновился свежими данными наверняка.

Теперь никакие катаклизмы не смогут уничтожить полученные данные, и мы можем уверенно запустить майнер на выходные.

Что я намайнил

Итак, представляю вашему вниманию результаты двухдневного майнинга:

precision: 1
Tossed: 33 552 384
 Heads: 15 724 557, 0.4686569216661325765
 Tails: 17 827 827, 0.5313430783338674235
  Time: 18ms 222us 700ns

precision: 2
Tossed: 268 433 408
 Heads: 133 163 234, 0.4960754884876326571
 Tails: 135 270 174, 0.5039245115123673429
  Time: 21ms 636us 200ns

precision: 3
Tossed: 2 147 481 600
 Heads: 1 072 710 984, 0.4995204540984192833
 Tails: 1 074 770 616, 0.5004795459015807167
  Time: 25ms 722us 500ns

precision: 4
Tossed: 21 474 834 432
 Heads: 10 736 419 663, 0.4999535478141562046
 Tails: 10 738 414 769, 0.5000464521858437954
  Time: 46ms 516us 900ns

precision: 5
Tossed: 227 633 264 640
 Heads: 113 815 517 221, 0.4999951013354671009
 Tails: 113 817 747 419, 0.5000048986645328991
  Time: 253ms 918us 100ns

precision: 6
Tossed: 1 026 497 181 696
 Heads: 513 248 105 522, 0.4999995272018192996
 Tails: 513 249 076 174, 0.5000004727981807004
  Time: 1s 38ms 463us 500ns

precision: 7
Tossed: 10 647 223 924 736
 Heads: 5 323 611 436 183, 0.4999999505800757343
 Tails: 5 323 612 488 553, 0.5000000494199242657
  Time: 10s 355ms 188us 500ns

precision: 8
Tossed: 11 583 526 795 264
 Heads: 5 791 763 344 721, 0.4999999954322201748
 Tails: 5 791 763 450 543, 0.5000000045677798252
  Time: 11s 292ms 750us 400ns

precision: 9
Tossed: 11 785 390 258 176
 Heads: 5 892 695 123 667, 0.4999999995400237174
 Tails: 5 892 695 134 509, 0.5000000004599762826
  Time: 11s 483ms 665us

precision: 10
Tossed: 13 284 333 844 480
 Heads: 6 642 166 921 632, 0.4999999999542318036
 Tails: 6 642 166 922 848, 0.5000000000457681964
  Time: 12s 903ms 693us 300ns

precision: 11
Tossed: 41 145 786 693 632
 Heads: 20 572 893 346 782, 0.4999999999991736699
 Tails: 20 572 893 346 850, 0.5000000000008263301
  Time: 1min 36s 550ms 404us 400ns

precision: 12
Tossed: 166 129 335 007 232
 Heads: 83 064 667 503 543, 0.4999999999995605833
 Tails: 83 064 667 503 689, 0.5000000000004394167
  Time: 3min 50s 949ms 149us 600ns

precision: 13
Tossed: 6 560 038 558 627 840
 Heads: 3 280 019 279 313 828, 0.4999999999999859757
 Tails: 3 280 019 279 314 012, 0.5000000000000140243
  Time: 6h 22min 42s 628ms 678us 500ns

precision: 14
Tossed: 6 568 598 428 448 768
 Heads: 3 284 299 214 224 400, 0.5000000000000024358
 Tails: 3 284 299 214 224 368, 0.4999999999999975642
  Time: 6h 22min 51s 19ms 966us 300ns

Красиво. Из примечательного — точность до 12 цифр посчиталась всего за три минуты, а вот следующие результаты пришлось ждать 6 часов. Зато сразу пачкой взятия 13 и 14 с разницей в десять секунд. О как. Остальное время выходных машина зря жгла CPU :) Ну, или нужно было еще немного подождать, кто знает.

На самом деле это не единственный мой запуск в долгую — такие результаты повторяются из раза в раз. И мы снова не освоили лимиты 64-битных чисел — до точности в 19 знаков точки нам как до луны. На таких выводах и остановимся.

Выводы

Завершая этот долгий путь, мы фиксируем следующие результаты:

За 6 часов 22 минуты мы подкинули монетку шесть с половиной квадриллионов раз (квадриллионы идут сразу после триллионов) и получили точность в 14 знаков после точки
Если же нам будет достаточно скромных 12 знаков точности, то за какие-то три минуты мы можем подкинуть монетку 166 триллионов раз и получить желаемое

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

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

constexpr решение задачи

Вообще, у задачи есть compile-time решение. Весь его код укладывается в эти строки:

#include <iostream>

int main()
{
    std::cout << "precision: inf" << std::endl;
    std::cout << "Tossed: inf" << std::endl;
    std::cout << " Heads: inf, 0.5" << std::endl;
    std::cout << " Tails: inf, 0.5" << std::endl;
    std::cout << "  Time: 0" << std::endl;

    return 0;
}

Выведет на экран:

precision: inf
Tossed: inf
 Heads: inf, 0.5
 Tails: inf, 0.5
  Time: 0

Я считаю это решение оптимальным.

Комментарии (5)

kipar
27.09.2025 08:07
#28888274
Сначала автор

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

генерирует 32 бита рандома и выкидывает 31 из них, потом начинает кешировать(?) биты чтобы было не так медленно потом при переходе к simd начинает использовать ГПСЧ нормально, а дальше даже берет нормальный для данной ситуации LCG (кстати константа совпадает с используемой в PCG32, а "шифрующая" функция оттуда в данном случае и не нужна).
Так что начинал читать с рукой прижатой к лицу но дальше все оказалось не так плохо.

da-nie
27.09.2025 08:07
#28888276
- Когда я научусь писать compute shaders, я стану жечь GPU вместо CPU, и тоже надеюсь на прорывные результаты
Да CUDA просто возьмите. :)

Jijiki
27.09.2025 08:07
#28888316
дальше только монте карло. взять багосорт и стохастически вероятностно предсказать сортировку :) интересно такое бывает? для меня пока это фантастика

Bogosort

kovserg
27.09.2025 08:07
#28888446
Интересно а какая у вас вероятность того что Heads < Tails ? Почему у вас постоянно Tails доминируют?

VBDUnit
27.09.2025 08:07
#28888492
Вероятность 50/50 предсказана для неких абстрактных «идеально‑случайных» значений. А подкидывание монетки и то, что описано в статье — это псевдослучайные числа, лишь некоторая приближённая модель случайных чисел. Вы как раз об этом и пишете.

На то же «подкидывание монетки» влияет психоэмоциональная сфера, время и длительность серии подкидываний, когнитивные искажения человека, форма ногтей, состояние суставов, процент ошибок при фиксации значений, концентрация молочной кислоты в мышце и прочие взаимосвязанные факторы. Про фишки из C и манипуляциями с SIMD я вообще молчу.

Поэтому и отличие наблюдаемого результата от ожидаемых 50/50 складывается далеко не только от »недостаточно большого количества итераций». Это как отличие идеальной абстрактной окружности из геометрии и реального колеса автомобиля. Модель работает лишь в каком‑то приближении.

Орел или решка. C++ Edition +6

Чего я хочу

Первичная идея

Тщательно готовимся к эксперименту

Замеряем время

Контролируемое кручение

Бросаем монетку

Изолируем эксперимент

Первые результаты

Оптимизации

Зафиксируем, что имеем

Как оптимизировать

Расчленяем рандом

Группируем биты

Реструктуризируем кручение

Считаем биты правильно

Расширяем генератор

Автоматическая многопоточность

Ручная многопоточность

Свой генератор

Ещё ядрёнее

SIMD

Итоги оптимизаций

Что еще могло бы ускорить?

Большие и малые числа

Трудности деления

Майним коины

Энергонезависимость

Что я намайнил

Выводы

constexpr решение задачи

Комментарии (5)

kipar

da-nie

Jijiki

kovserg

VBDUnit