news.shamcode.ru | Поделить нельзя — умножить или алгоритм быстрого деления по методу Ньютона-Рафсона

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

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

Метод касательных и как это связано с делением

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

Итак, у нас есть функция $inline$ , чей корень нужно найти, и у нас есть начальное приближение $inline$ к искомому корню. Шаги алгоритма получаются заменой функции её касательной, начиная с точки $inline$ , и решением уже линейного относительно x уравнения: $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$ . Конечно, наивно полагать, что любая $inline$ с наугад выбранным $inline$ сойдётся к корню… а нам это и не очень важно — самое время вооружиться конкретной функцией и понять её свойства.

В качестве такой функции возьмём $inline$ , где D — наш делитель. Почему такая функция? — Во-первых, очевидно, что она имеет ноль при $inline$ , во-вторых, если её продифференцировать и подставить производную в выражение $x_{n+1}$ и упростить, получится формула без единого деления: $X_{i+1} = X_i(2-DX_i)$ . Но на этом не всё: рассмотрим, так называемую, невязку функции: $\epsilon_{i+1} = 1 - DX_{i+1}$ , которая показывает, насколько близко наш $X_{i+1}$ к истинному значению $inline$ , и подставим туда значение для $X_{i+1}$ , получим $\epsilon_{i+1} = 1 - D X_{i+1} = 1 - D (X_i(2-DX_i)) = 1 - 2DX_i + D^2X_i^2 = (1 - DX_i)^2 = {\epsilon_i}^2.$ Проговорим словами: невязка следующего шага равна квадрату предыдущего шага. Тут важно заметить, что при получении этой формулы ~~никто не постардал~~ были использованы только алгебраические средства: не звучали слова про «малые окрестности», «разложение в ряд Тейлора» и т. д. Это значит, что равенство носит глобальный (для данной функции характер), и может служить очевидным критерием, когда последовательность сходится, а именно, если $|\epsilon_{0}| = |1 - D X_{0}| < 1.0$ , иначе говоря, $X_0 \in (0; 2/D)$ . Другая важная вещь: после первого же шага, невязка, как квадрат числа, станет положительной — наши приближения всегда будут меньше истинного значения $inline$ — к этому мы ещё вернёмся.

▍ Демонстрация работы метода

Посмотрим, как работает метод касательной на функции $inline$ при разных начальных приближениях $inline$ , заодно потестируем наши теоретические выводы.
Возьмём $inline$ и $inline$ , что даёт $\epsilon_0 = 0.985$

Зеленая точка — корень функции $inline$ , видно, что пересечения касательных с осью X стягиваются в эту точку — метод сходится, как и было расчитано.
D возьмём тоже самое, а $inline$ , что даёт $\epsilon_0 = -1.0999$

Видно, что метод быстро расходится.

Реализация метода

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

Делим целые числа без знака, а именно uint16_t.
Результатом считаем целое частное $\lfloor n/m \rfloor$ того же типа.
Остаток не возвращаем (хотя это возможно).
Пытаемся написать быстрый код, чтоб посоревноваться со скоростью встроенного оператора деления.
Алгоритм не будет отдельно обрабатывать единицу и другие степени двойки.
Алгоритм не будет отдельно обрабатывать ноль в числителе или знаменателе.

Ну вот, казалось бы, всё просто и понятно — есть формулы для шагов, как-нибудь угадаем начальные приближения и готово — ничего сложного. Небольшой нюанс: все числа $inline$ , кроме $inline$ — дробные. Таким образом, либо мы в реализацию тянем double, либо fixed point numbers. Второе предпочтительнее: в этом случае мы формально не выходим в реализации за integral types. Если вы никогда не слышали про числа с фиксированной точкой — советую почитать Википедию или одну из моих предыдущих статей.

▍ Трюки, трюки… очень много трюков

Взглянем опять на формулу «шага» алгоритма: $X_{i+1} = X_i(2-DX_i)$ , если D «пробегает» весь uint16_t от 1го 65535, то наши fixed point numbers должны хранить 16 разрядов до запятой (точки) и столько же после… при этом чаще всего старшая часть будет нулевой. Что, если «нормализовать» D, как это сделано в double/float числах, а именно самый старший ненулевой бит сдвинуть максимально налево и запомнить сдвиг, иначе говоря, хранить мантиссу и степень? В этом случае мантисса $D \in [1.0; 2)$ и соответствующая мантисса искомого $1/D \in (0.5, 1]$ . Это преобразование также позволяет искать начальные приближения только в диапазоне $inline$ .

Чувствую, пора добавить кода в статью, а именно взглянем на формирование массива приближений $inline$ :

void print_reciprocal(uint8_t bit_number)
{
    if (!bit_number || bit_number > CHAR_BIT - 1)
    {
        throw std::invalid_argument("It's expected to fit ony byte");
    }

    const uint8_t count = 1 << bit_number;

    std::cout << "uint16_t reciprocals_" << (uint16_t)count << "[] = {";

    for(uint8_t i = 0; i < count; ++i)
    {
        uint8_t denominator = count;
        denominator |= i;

        double reciprocal = static_cast<double>(1.) / static_cast<double>(denominator);

        uint8_t first_byte = std::scalbln(reciprocal, CHAR_BIT + bit_number);

        if (!first_byte)
        {
            first_byte = 0xFF;
        }

        std::cout  << ' ' << std::uppercase << std::hex << "0x" << ((uint16_t)first_byte << 8) << ", ";

        if (!(count % 15))
        {
            std::cout << std::endl;
        }
    }

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

Вывод:

uint16_t reciprocals_8[] = { 0xFF00,  0xE300,  0xCC00,  0xBA00,  0xAA00,  0x9D00,  0x9200,  0x8800, };

Очевидно, что данная функция подготавливает искомый массив, формируя его сразу синтаксически корректной конструкцией языка C++ :), но что такое параметр bit_number и почему второй байт всегда нулевой? Чтобы ответить на эти вопросы, нужно понять, как дальше этот массив будет использован. Кстати, далее по тексту для краткости я буду называть этот массив LUT (Look-Up-Table), следуя нотациям источников.

Итак, у нас есть «нормализованный» делитель D, то есть его старший бит всего равен 1, возьмём его следующие 3 (bit_number) бит и используем их как индекс в LUT. Заметим, что в шестнадцатеричной системе счисления старшая цифра никогда не отличается больше чем на 1, иначе говоря, отличается в 4-ом разряде. В этом и причина, почему мы не храним второй байт: по сути из первого байта мы используем только старшую тетраду.

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

Почему элементы массива LUT отличаются не более чем на 1/16 и почему метод сходится?

Числа в LUT, конечно, это дроби, иначе говоря, количества 1/256.
Рассмотрим разницу двух соседних элементов в LUT $1/x - 1/y = \frac{y-x}{yx} = \frac{1}{8xy}$ учтём, что x>=1, y>1, получим, что разница «проваливается» в 4-й разряд.
Для D = 1 мантисса специально округлена вниз до 0xFF00, чтоб влезть в этот же формат.
Оценим изначальную невязку, для этого поделим её на D: $\frac{1}{D} - X_0$ , с учётом вышеполученного мы отбрасываем 4ый разряд и всё за ним — в любом случае не больше, чем 1/8 (удвоенный 4-й разряд), умножим обратно на D, всё равно получим значение меньше 1го.

Уже можно привести кусок кода, в котором мы «нормализуем» делитель и извлекаем начальное приближение из массива

    
// __builtin_clz returns the number of the first not zero bit counting from the left, and the argument is widened to 4 bytes int
    int shift_to_left =  __builtin_clz(v) - 16;

    // the first not zero bit should be the most left in uint16_t
    v <<= shift_to_left;

    // to look it up in the table we should first move the significant for us part to the right and then to zero the most significant (the fourth) bit
    uint16_t x = reciprocals_8[(v >> (8 + 3 + 1)) - 8];

Прежде чем показать реализацию «шага» обновления x, нужно понять, что именно хранится в числах x и v. Для этого будем использовать следующую нотацию: Q X.Y означает число с фиксированной точкой, где X — число бит для целой части, Y — число бит для дробной части, если целая часть совсем отсутствует пишем Q .Y, очевидно, для нашего примера X+Y=16.
Получается, v будет в формате Q 1.15, x — Q .16, произведение vx — нужное нам по формуле шага — снова Q 1.15 (16 младших разрядов мы отбрасываем), далее, (2 — vx) — опять Q 1.15.

Последнее выражение (в рамках Q 1.15) упрощается до -(vx) — 2 для нас тут то же самое, что 256 для однобайтовых целых. Осталось ещё пара шагов… но проще уже показать код.


// simple helper function
uint16_t multHigherHalf(uint16_t a, uint16_t b)
{
    uint32_t res = a * b;
    res >>= 16;
    return res;
}

x = multHigherHalf(static_cast<uint16_t>((-(v * x)) >> 16), x) << 1;

Сдвиг налево нужен, потому что формально результат получается Q 1.15, но первый бит всегда нулевой, да и для следующего шага нужно иметь x в том же самом формате.

Соберём всё вместе:

uint16_t divide(uint16_t u, uint16_t v)
{
    // __builtin_clz returns the number of the first not zero bit counting from the left, and the argument is widened to 4 bytes int
    int shift_to_left =  __builtin_clz(v) - 16;

    // the first not zero bit should be the most left in uint16_t (normalization)
    v <<= shift_to_left;

    // to look it up in the table we should first move the significant for us part to the right and then to zero the most significant bit
    uint16_t x = reciprocals_8[(v >> (8 + 3 + 1)) - 8];

    // two steps of Newton
    x = multHigherHalf(static_cast<uint16_t>((-(v * x)) >> 16), x) << 1;
    x = multHigherHalf(static_cast<uint16_t>((-(v * x)) >> 16), x) << 1;

    uint16_t q = multHigherHalf(x, u);
    q >>= 16 - shift_to_left - 1;

    //get normalization back
    v >>= shift_to_left;

    uint32_t reminder = u - q * v;

    // the reminder should be always less than v
    if (reminder >= v)
    {
        reminder -= v;
        ++q;

        if (reminder >= v)
        {
            ++q;

        }
    }

    return q;
}

Тут осталось обсудить всего пару моментов:

Почему именно 2 шага метода касательных.
Зачем считается остаток и увеличивается частное.

Нам достаточно 2 шага, потому что у нас точность на уровне 4го бита и метод имеет квадратичную сходимость — после двух шагов точность будет на уровне 16го бита. Заодно, таким образом, мы избегаем цикла, в котором бы контролировалось значение $\epsilon_i$ . Также, учитывая точность на уровне 16го бита и отбрасывание более младших разрядов, нужно сделать 2 «шага коррекции»: $\epsilon_i$ больше нуля, значит, приближение меньше истинного значения, то есть остаток может уменьшиться и частное увеличиться (а вот в этой реализации это не так).

Результаты измерений и оптимизаций

Определимся с методикой измерений:

Время каждой отдельной операции слишком мало — измеряем продолжительность всего цикла по всем делимым и делителям.
Используем std::chrono::high_resolution_clock::now().
Поскольку результаты измерений будут различаться от разу к разу — запустим 5 раз и приведём среднюю величину, округлённую до 2ух знаков.
Весь код однопоточный.
CPU: Intel® Core(TM) i5-9300H CPU @ 2.40GHz.
Измеряем на релизной сборке, собранной с флагами -Wall -Wextra -Wpedantic -Werror -O3.

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

    auto start = std::chrono::high_resolution_clock::now();

    for(uint16_t divisor = 1; divisor > 0; divisor++)
    {
        for(uint16_t numenator = 1; numenator > 0; numenator++)
        {
            volatile uint16_t res = numenator / divisor;
            (void)(res);
        }
    }

    auto stop = std::chrono::high_resolution_clock::now();

    std::cout << "duration = " << std::chrono::duration_cast<std::chrono::microseconds>(stop - start).count() << std::endl;

Получаем средний результат: 9.4 секунды.

Не привожу исходный код измерений для метода divide, потому что он отличается только одной строкой, средний результат: 5.8 секунды — лучше! Причём — значительно!

А можно ещё лучше? — Да, но это будет хак: тип uint16_t не такой уж большой… можно сформировать массив всех дробей 1/v, и это займёт всего 2 * 64Кб — 128Кб памяти.

void print_all_reciprocals(const char* file_name)
{
    if (!file_name || ! std::strlen(file_name))
    {
        return;
    }

    std::ofstream out(file_name);

    uint16_t i = 0;

    out << "uint16_t all_reciprocals" << "[] = {";

    while (true)
    {
        if (!i)
        {
            out << "0x0000, ";
            i++;
            continue;
        }

        if (i == 1)
        {
            out << "0xFFFF, ";
            i++;
            continue;
        }

        double reciprocal = static_cast<double>(1.) / static_cast<double>(i);

        uint16_t two_bytes = std::scalbln(reciprocal, CHAR_BIT * sizeof(uint16_t));

        out << "0x" << std::uppercase << std::hex << std::setfill('0') << std::setw(4) << two_bytes << ", ";

        if (!(i % 15))
        {
            out << std::endl;
        }

        if (i == std::numeric_limits<uint16_t>::max())
        {
            break;
        }

        i++;
    }

    out << "};" << std::endl;

    out.close();
}

uint16_t divide_wo(uint16_t u, uint16_t v)
{
    uint16_t x = all_reciprocals[v];
    uint16_t q = multHigherHalf(x, u);

    uint32_t reminder = u - q * v;

    if (reminder >= v)
    {
        q++;
    }

    return q;
}

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

Выводы

Мы показали, что метод Ньютона-Рафсона может соревноваться по скорости со встроенным оператором деления языка C++. Хотя, для uint16_t он является избыточным для более широких типов uint32_t, uint64_t полная таблица дробей будет заниматься неоправданно много памяти, а значит, реализация и применения метода касательных будет иметь смысл.

Источники

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