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

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

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

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


Два значе��ия


Сразу оговорюсь: всюду под средним будет пониматься результат деления суммы аргументов нацело на их количество, иначе говоря, average(1, 1) это 1, и average(1, 2) это тоже 1, потому что$ 3 / 2$ будет 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, что потребует кода сравнения и перестановки двух чисел: $a + (b - a) / 2$.

А ещё есть способ избавиться от обоих делений, используя немного битовой магии: 
(a & b) + ((a ^ b) >> 1).

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

Три значения


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

$a = qb + r$

Где

$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 — нет. В свете этого не столь большая разница в их времени выполнения становится ещё более загадочной, тем не менее — оставим исследование этого феномена до следующего раза (или до ответа в комментариях).





Выводы


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

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

Ссылки



© 2025 ООО «МТ ФИНАНС»

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