StarPilgrim 3 июн в 16:01

Недистрибутивность деления, или Как я считал среднюю величину

Средний

8 мин

5.1K

Блог компании RUVDS.comC++*Алгоритмы*Ненормальное программирование*

Обзор

Казалось бы: сложно отыскать формулу проще, чем нахождение среднего арифметического. Однако код — не формула, вдобавок, если вы пишете на С++, то разного (и в основном неприятного) рода сюрпризы могут ожидать вас где угодно.

Постановка задачи: реализовать функцию uint32_t average(uint32_t a, uint32_t b), не используя типов шире, чем uint32_t, и затем обобщить этот подход на произвольное количество аргументов.

Небольшая преамбула: зачем я этим всем занимаюсь?

Позвольте мне привести здесь небольшой диалог (приписываемый Давиду Гилберту):

Математик Давид Гилберт отвечает на вопросы:
— Какая задача сейчас для науки наиболее важна?
— Поймать муху на обратной стороне Луны!
— А кому это надо?
— Это никому не надо. Но подумайте над тем, сколько важных сложнейших задач надо решить, чтобы это осуществить!

Два значения

Сразу оговорюсь: всюду под средним будет пониматься результат деления суммы аргументов нацело на их количество, иначе говоря, average(1, 1) это 1, и average(1, 2) это тоже 1, потому что $inline$ будет 1 — то есть округляем всегда вниз.

Итак, сначала рассмотрим вариант наивный и, очевидно, неверный:

uint32_t average(uint32_t a, uint32_t b)
{
    return (a + b) / 2;
}

К сожалению или нет, но такое буквальное переложение школьной формулы в код будет очень часто давать неверный ответ, например: average(0x80000000, 0x80000000) будет равно нулю, а должно быть равно 0x80000000 — это происходит из-за переполнения uint32_t при сложении a + b, при этом важно заметить, что и аргументы, и конечный результат во всех случаях влезают в тип uint32_t . Таким образом, задача-то вполне корректна, в отличие от кода.

Почему бы не рассмотреть всё то же самое, но с типом uint8_t

Это подход, который я исповедовал в своих предыдущих статьях об алгоритмах деления: тут и тут. В данном случае он не сработает в своей наивной форме, так как компилятор неявно использует int для промежуточного значения. Для деталей можно обратиться к секции integral promotions.

Первая пришедшая мне в голову идея была основана лишь на интуиции и заключалась в следующем исправлении:

uint32_t average(uint32_t a, uint32_t b)
{
    return a / 2 + b / 2;
}

Проблемы с переполнением тут, конечно, быть не может, потому что a /2 и b /2 не превосходят половины от numeric_limits<uint32_t >::max(), но проблема с неправильными ответами остаётся актуальной: average(3, 3) будет равен...2, а не 3.

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

uint32_t average(uint32_t a, uint32_t b)
{
    return a / 2 + b / 2 + (a & b & 1);
}

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

Для самых дотошных

Конечно, это не все способы найти среднее двух чисел, а именно: есть вариант избавиться от одного деления, но он может страдать от underflow, что потребует кода сравнения и перестановки двух чисел: $inline$ .

А ещё есть способ избавиться от обоих делений, используя немного битовой магии:

(a & b) + ((a ^ b) >> 1).

Однако обобщить эти способы на большее количество чисел не представляется возможным.

Три значения

Прежде чем идти дальше, повторим материальную часть: для пары натуральных чисел a и b (b не равно нулю) существует единственное представление вида:

$display$

Где

$q \ge 0, r\ge0\ и\ r < b$

Здесь q — частное, r — остаток, а процесс нахождения чисел q и r, собственно говоря, и называется делением.

Если мы запишем все три значения в вышеприведённом виде, учтя, что b = 3, затем перегруппируем их, а уже потом разделим, то, во-первых, очевидно, ничего не изменится (от перестановки мест слагаемых сумма не меняется), а, во-вторых, получится следующее:

uint32_t average(uint32_t x1, uint32_t x2, uint32_t x3)
{
    return x1 / 3 + x2 / 3 + x3 / 3  + (x1 % 3 + x2 % 3 + x3 % 3) / 3;
}

Проговорим алгоритм словами: сначала сложим все частные, потом посмотрим: «добежали ли» остатки до целого значения или нет.

Выглядит неплохо, но может ли здесь быть переполнение? Короткий ответ — нет.

Действительно: x1 / 3 + x2 / 3 + x3 / 3 не превосходит максимального значения numeric_limits<uint32_t >::max(), а если быть ещё точнее: это будет максимальное из чисел, делимых на три и не превосходящих при этом numeric_limits<uint32_t >::max(). Пожалуй, я его посчитаю:

 cout << dec << numeric_limits<uint32_t>::max() - numeric_limits<uint32_t>::max() % 3 << ' ' << numeric_limits<uint32_t>::max() << endl;

Вывод:

4294967295 4294967295

То есть numeric_limits<uint32_t>::max() делится на 3 без остатка (если не верите — сосчитайте сумму его цифр и убедитесь, что она делится на 3), x1 / 3 + x2 / 3 + x3 / 3 будет равно исходному числу, а вторая часть с остатками будет нулевой. Случай, если один или более элементов меньше максимального значения, — довольно технический и проверяемый такими же рассуждениями — не будем тратить на это время.

N значений. Первая попытка

Обобщая код из нашего последнего примера, получаем:

uint32_t average(const vector<uint32_t>& nums)
{
    uint32_t  res = 0;
    uint32_t  rem = 0;
    uint32_t  size = nums.size();

    if (!size)
    {
        return 0;
    }

    for (const auto& x: nums)
    {
        res += x / size;
    }

    for (const auto& x: nums)
    {
        rem += x % size;
    }

    return res + rem / size;
}

Отлично, код выглядит, прямо скажем, незамысловато! А вот работает ли?

 cout  << average({1, 2, 3, 4, 5, 6, 7, 8, 9}) << endl;

Эта строчка предсказуемо выведет 5. Попробуем другой пример:

 cout  << hex <<  average(vector<uint32_t>(0x10001, 0x10000)) << endl;

Напомню, что первый аргумент в конструкторе вектора — это количество элементов, и второй — значение этих элементов, вывод будет:

0

Похоже, мы опять столкнулись с переполнением, но ведь выше мы разбирали, почему переполнения возникнуть не может? Или может?

Переполнение не может возникнуть здесь:

    for (const auto& x: nums)
    {
        res += x / size;
    }

И переполнение не может возникнуть здесь:

 return res + rem / size;

При условии, что rem сам по себе посчитан корректно.

Но оно может возникнуть и возникает здесь:

    for (const auto& x: nums)
    {
        rem += x % size;
    }

Проанализируем работу функции с данным вектором:

x / size равно 0, x % size равно самому x т.е. 0x10000, в то же время 0x10000 * 0x10001 = (0x10000 * 0x10000) + 0x10000 — первое слагаемое не влезает в 4 байта и даёт 0 — значит итоговый rem будет 0x10000 и при делении на 0x10001 мы получим 0.

N значений. Чиним алгоритм

В предыдущем абзаце была выявлена проблема с циклом накопления остатка rem — рассмотрим его подробнее. Например, не утруждаясь оперированием большими числами, возьмём вектор из трёх элементов — пусть каждый элемент будет равен 8, тогда x % size равно 2. А есть ли вообще смысл в накоплении суммарного остатка? Нет. Достаточно факта, что остаток превышает или равен модулю (модуль — он же размер вектора в данном контексте) и накопления остатка уже взятого по модулю. Так, на первом шаге мы получим rem = 2, на втором шаге мы не получим 4, вместо этого rem = 1, и инкрементация частного (quotient++).

Воплотим идею в код:

uint32_t average(const vector<uint32_t>& nums)
{
    if (!nums.size())
    {
        return 0;
    }
    
    uint32_t size = nums.size();
    uint32_t quotient = 0;
    uint32_t reminder = 0;
    
    for (const auto& x: nums)
    {
        quotient += x / size;
        uint32_t rem = x % size;
        
        if ((size - reminder) > rem)
        {
            reminder += rem;
        }
        else
        {
            quotient++;
            reminder += rem - size;
        }
    }
    return quotient;
}

Здесь в if мы контролируем переполнение, но не переполнение uint32_t, а переполнение нашего модуля (size). Если его нет, то просто накапливаем остаток — иначе пересчитываем остаток уже по модулю (size) и инкрементируем частное (quotient). Тут нужно понимать, что за одну итерацию частное (quotient) может увеличиться не больше чем на 1, потому что и reminder, и rem оба меньше модуля (size) — следовательно, их сумма меньше 2 * size.

Можно убедиться, что предыдущий пример отработает теперь верно, выведя 0x10000:

 cout  << hex <<  average(vector<uint32_t>(0x10001, 0x10000)) << endl;

Анализ алгоритма и замеры

В терминах О-большого полученный алгоритм ничем не отличается от своей наивной и страдающей от переполнений версии:

uint32_t naive_average(const vector<uint32_t>& nums)
{
    uint32_t  sum = 0;
    
    for (const auto& x: nums)
    {
        sum += x;
    }
    
    return sum / nums.size();
}

В обоих случаях количество итераций равно размеру вектора, и в каждой итерации выполняется конечное, не зависящее от размера вектора, число операций — то есть это линейная сложность O(n). Однако в реальном мире ~~чистогана и наживы~~, где процессорные такты превращаются в секунды выполнения — те секунды, которые пользователи проведут в ожидании результата, — в этом мире эти два алгоритма существенно отличаются. Почему? Потому что деление довольно медленная операция, и, если в naive_average на весь алгоритм у нас одно деление, то в алгоритме без переполнений — столько, сколько элементов в векторе.

Однако на данный момент это всего лишь гипотеза — нужно подкрепить её результатами измерений.

#include <iostream>
#include <vector>
#include <chrono>
#include <cstdint>
#include <random>

using namespace std;
using namespace std::chrono;

uint32_t average(const vector<uint32_t>& nums) {
    if (nums.empty()) {
        return 0;
    }
    
    uint32_t size = nums.size();
    uint32_t quotient = 0;
    uint32_t reminder = 0;
    
    for (const auto& x : nums) {
        quotient += x / size;
        uint32_t rem = x % size;
        
        if ((size - reminder) > rem) {
            reminder += rem;
        } else {
            quotient++;
            reminder += rem - size;
        }
    }
    return quotient;
}

uint32_t naive_average(const vector<uint32_t>& nums) {
    uint32_t sum = 0;
    for (const auto& x : nums) {
        sum += x;
    }
    return sum / nums.size();
}

vector<uint32_t> generate_test_data(size_t size) {
    vector<uint32_t> data;
    data.reserve(size);
    random_device rd;
    mt19937 gen(rd());
    uniform_int_distribution<uint32_t> dist(0, UINT32_MAX);
    for (size_t i = 0; i < size; ++i) {
        data.push_back(dist(gen));
    }
    return data;
}

template<typename Func>
long long measure_time(Func func, const vector<uint32_t>& data, uint32_t& result) {
    auto start = high_resolution_clock::now();
    result = func(data);
    auto end = high_resolution_clock::now();
    return duration_cast<milliseconds>(end - start).count();
}

int main() {
    vector<size_t> sizes = {100'000, 200'000, 400'000, 800'000, 1'600'000};
    
    cout << "Testing different vector sizes:\n";
    for (size_t size : sizes) {
        auto data = generate_test_data(size);
        uint32_t res_opt, res_naive;
        
        auto time_opt = measure_time(average, data, res_opt);
        auto time_naive = measure_time(naive_average, data, res_naive);
        
        cout << "Size: " << size << endl;
        cout << "  Optimized average:  " << time_opt << " ms\t Result: " << res_opt << endl;
        cout << "  Naive average:      " << time_naive << " ms\t Result: " << res_naive << endl;
        cout << "----------------------------------------\n";
    }
    
    return 0;
}

Собрав этот код с флагами -O3 в релизе, получим следующие результаты:

Testing different vector sizes:
Size: 10000000
Optimized Average: 14 ms Result: 2147686718
Naive average: 2 ms Result: 348
— Size: 20000000
Optimized Average: 25 ms Result: 2147422219
Naive average: 3 ms Result: 203
— Size: 40000000
Optimized Average: 49 ms Result: 2147451718
Naive average: 6 ms Result: 67
— Size: 80000000
Optimized Average: 103 ms Result: 2147364889
Naive average: 13 ms Result: 50
— Size: 160000000
Optimized Average: 200 ms Result: 2147464194
Naive average: 31 ms Result: 8
—

Видно, что версия алгоритма, учитывающая переполнение, медленнее в 7-8 раз наивного алгоритма. Может показаться, что разница весьма внушительна, но на самом деле такая разница вызывает вопросы, потому что деление обычно более чем в 10 раз медленнее сложения.
Железо, на котором запускали тесты: 13th Gen Intel® Core(TM) i7-13700H, 32 Гб ОЗУ.

Если посмотреть, как компилятор оптимизировал naive_average, например, в godbolt, то можно заметить, что компилятор использовал в ней SIMD-инструкции, в то время как в average — нет. В свете этого не столь большая разница в их времени выполнения становится ещё более загадочной, тем не менее — оставим исследование этого феномена до следующего раза (или до ответа в комментариях).

Выводы

Даже простые формулы и привычные математические операции часто утрачивают свои аксиоматические свойства, будучи реализованными на реальном «железе».

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

Ссылки

Telegram-канал со скидками, розыгрышами призов и новостями IT 💻

Хабы: