Все мы в школе проходили деление «столбиком» — простой алгоритм, который несложно реализовать, вот только не очень быстрый.
В прошлый раз мы рассматривали, как компилятор оптимизирует деление в случаях, когда делитель известен во время компиляции, но применение его напрямую, чтоб оптимизировать деление для делителей, определямых в run-time, невозможно: вычисление констант сдвига и умножения само по себе требует деления.
В этот раз хочется поговорить о другом методе, сводящем деление к умножениям и битовым сдвигам, основанном на методе поиска корней функции.
Метод касательных и как это связано с делением
Метод касательных — он же метод Ньютона-Рафсона — итеративный численный метод нахождения корня функции, основанный на замене функции её касательной и нахождения уже её корня. Важно понимать: для нашей задачи это всего лишь инструмент, поэтому всеобъемлющего описания приведено не будет — сосредоточимся на применении этого метода для деления чисел.
Итак, у нас есть функция
, чей корень нужно найти, и у нас есть начальное приближение
к искомому корню. Шаги алгоритма получаются заменой функции её касательной, начиная с точки
, и решением уже линейного относительно x уравнения:
. Конечно, наивно полагать, что любая
с наугад выбранным
сойдётся к корню… а нам это и не очень важно — самое время вооружиться конкретной функцией и понять её свойства.
В качестве такой функции возьмём
, где D — наш делитель. Почему такая функция? — Во-первых, очевидно, что она имеет ноль при
, во-вторых, если её продифференцировать и подставить производную в выражение
и упростить, получится формула без единого деления:
. Но на этом не всё: рассмотрим, так называемую, невязку функции:
, которая показывает, насколько близко наш
к истинному значению
, и подставим туда значение для
, получим
Проговорим словами: невязка следующего шага равна квадрату предыдущего шага. Тут важно заметить, что при получении этой формулы
никто не постардал были использованы только алгебраические средства: не звучали слова про «малые окрестности», «разложение в ряд Тейлора» и т. д. Это значит, что равенство носит глобальный (для данной функции характер), и может служить очевидным критерием, когда последовательность сходится, а именно, если
, иначе говоря,
. Другая важная вещь: после первого же шага, невязка, как квадрат числа, станет положительной — наши приближения всегда будут меньше истинного значения
— к этому мы ещё вернёмся.
▍ Демонстрация работы метода
Посмотрим, как работает метод касательной на функции
при разных начальных приближениях
, заодно потестируем наши теоретические выводы.
Возьмём
и
, что даёт
Зеленая точка — корень функции
, видно, что пересечения касательных с осью X стягиваются в эту точку — метод сходится, как и было расчитано.
D возьмём тоже самое, а
, что даёт
Видно, что метод быстро расходится.
Реализация метода
Прежде чем писать код, уточним постановку задачи:
- Делим целые числа без знака, а именно uint16_t.
- Результатом считаем целое частное того же типа.
- Остаток не возвращаем (хотя это возможно).
- Пытаемся написать быстрый код, чтоб посоревноваться со скоростью встроенного оператора деления.
- Алгоритм не будет отдельно обрабатывать единицу и другие степени двойки.
- Алгоритм не будет отдельно обрабатывать ноль в числителе или знаменателе.
Ну вот, казалось бы, всё просто и понятно — есть формулы для шагов, как-нибудь угадаем начальные приближения и готово — ничего сложного. Небольшой нюанс: все числа
, кроме
— дробные. Таким образом, либо мы в реализацию тянем double, либо fixed point numbers. Второе предпочтительнее: в этом случае мы формально не выходим в реализации за integral types. Если вы никогда не слышали про числа с фиксированной точкой — советую почитать
Википедию или одну из
моих предыдущих статей.
▍ Трюки, трюки… очень много трюков
Взглянем опять на формулу «шага» алгоритма:
, если D «пробегает» весь uint16_t от 1го 65535, то наши fixed point numbers должны хранить 16 разрядов до запятой (точки) и столько же после… при этом чаще всего старшая часть будет нулевой. Что, если «нормализовать» D, как это сделано в double/float числах, а именно самый старший ненулевой бит сдвинуть максимально налево и запомнить сдвиг, иначе говоря, хранить мантиссу и степень? В этом случае мантисса
и соответствующая мантисса искомого
. Это преобразование также позволяет искать начальные приближения только в диапазоне
.
Чувствую, пора добавить кода в статью, а именно взглянем на формирование массива приближений
:
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 учтём, что x>=1, y>1, получим, что разница «проваливается» в 4-й разряд.
- Для D = 1 мантисса специально округлена вниз до 0xFF00, чтоб влезть в этот же формат.
- Оценим изначальную невязку, для этого поделим её на D: , с учётом вышеполученного мы отбрасываем 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го бита. Заодно, таким образом, мы избегаем цикла, в котором бы контролировалось значение
. Также, учитывая точность на уровне 16го бита и отбрасывание более младших разрядов, нужно сделать 2 «шага коррекции»:
больше нуля, значит, приближение меньше истинного значения, то есть остаток может уменьшиться и частное увеличиться (а вот в этой
реализации это не так).
Результаты измерений и оптимизаций
Определимся с методикой измерений:
- Время каждой отдельной операции слишком мало — измеряем продолжительность всего цикла по всем делимым и делителям.
- Используем 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 полная таблица дробей будет заниматься неоправданно много памяти, а значит, реализация и применения метода касательных будет иметь смысл.
Источники
© 2024 ООО «МТ ФИНАНС»
Telegram-канал со скидками, розыгрышами призов и новостями IT 💻