Если вы никогда не пробовали смотреть как код на C++ разворачивается компилятором в код Assembly – вас ждёт много сюрпризов, причём, не нужно смотреть какой-то замудренный исходный код полный templates или других сложных конструкций: рассмотрите следующий snippet:
uint8_t div10(uint8_t x)
{
return x/10;
}Конечно, я это уже сделал, и приведу результаты прямо здесь, хотя, советую и самим сходить на замечательный ресурс https://godbolt.org/ – выставить там, например, x86-64 gcc 14.1, добавить опцию -O2 и убедиться в том, что результаты крайне интересные:
div10(unsigned char):
mov edx, -51
mov eax, edi
mul dl
shr ax, 11
retДействительно, чего совсем не видно в этом куске кода так это инструкции div, которая несомненно существует и собственно осуществляет деление на x86 архитектуре. Зато – откуда-то взялась магическая константа, ещё и отрицательная! Ну ладно, на самом деле она положительная и равна 205 (т.к. 205+51=256 = mod 256), так что этот вопрос закрыт, но как же всё таки это всё работает?
Базовый принцип
Работает это следующим образом после умножения на 205 мы делаем суммарный сдвиг вправо на 11 разрядов, что эквивалентно делению на 2^11 = 2048 с отбрасыванием остатка. Последнее – очень важно. Заметим, что 205/2048 = 0.10009765625, иначе говоря чуть-чуть больше чем 0.1, если вы умножите последнее число на калькуляторе (или в Python) на 255 вы получите 25.524, иначе говоря, после отбрасывания дробной части – это правильный ответ.
А можно было взять число чуть-чуть меньше чем 0.1? Нет – тривиально умножив на 10 мы бы получили "чуть-чуть меньше чем 1", и просто отбрасывая дробную часть получили бы 0 – немного не тот результат которого ожидаешь деля 10 на 10. А насколько чуть-чуть больше должно быть число, что б трюк сработал? Для uint8_t максимальное число – 255, соответственно, число должно отличаться от 1/10 не больше чем на 1/256. А с любым ли числом (хотя бы из базового типа uint8_t это возможно) – да, с любым.
Тут, я напомню про такие математические функции как floor, ceil и trunc: мы работаем только с положительными числами, поэтому, без лишних сложностей floor=trunc и просто отбрасывает дробную часть, a ceil всегда округляет вверх. Пусть d – наш делитель, N – наша разрядность (у нас 8), мы хотим получить выражение вида: m / 2 ^ (N + k) такое что оно чуть больше 1/d (тут обо всём думаем в вещественных числах), а если точнее то:
m / (2 ^ (N + k)) - (1 / d) >= 1 / 2^N (просто обобщение предыдущего абзаца).
Как найти подходящие числа
Утверждение: m = ceil(2^(ceil(log(d)) + N) / d), k = ceil(log(d)) – как это всё понять, что тут такое написано? ... Это можно понять, например, таким образом: число внизу по условию это степень 2ки, и эта степень точно не меньше чем 2^N, далее предположим, я как-то нашёл k – как теперь по заданному k как подобрать m ? Я уже знаю знаменатель – 2 ^ (N + k), я хочу подобрать целое число, что б оно было чуть больше чем 1/d, что если я возьму просто trunc(2 ^ (N + k) / d) ?
Давайте в числах из примера: я хочу делить uint8_t на 10: т.е. N=8, d = 10, а k пока возьмем равным 2, например, тогда trunc(2 ^ (N + k) / d) = trunk(2^10 / 10) = trunk (1024 /10) = 102, а всё выражение m / (2 ^ (N + k)) = 102/ 1024 = 0,099609375, в общем, близко, но чуть меньше чем нам надо. А вот ceil - всегда будет больше потому что: ceil(x) >= x для положительных чисел. Я ещё не сделал эту оговорку, но сделаю, что d > 1 и d не является точной степенью двойки, то есть точных делений у нас тут не будет, вторая оговорка trunc (x/y) в C++ это просто обычное целочисленное деление.
Итак, я надеюсь, к этому моменту стало понятно, что m в том виде как я ищу действительно аппроксимирует 1/d сверху. Теперь посмотрим почему я выбрал такое k: ceil(2^(ceil(log(d)) + N) / d, вот здесь делая внешний ceil я прибавляю к числителю число не более чем d – потому что остаток не может быть ни больше ни равен d, и понятно что 2^(ceil(log(d)) > d.
Думаю, время показать код:
template<typename InputInteger, typename OutputInteger>
std::pair<OutputInteger, uint8_t> getDivisionMultiplier(InputInteger divisor)
{
if (!divisor)
{
throw std::invalid_argument("Division by zero is impossible");
}
if (divisor == 1)
{
return {1,0};
}
constexpr uint8_t n = sizeof(InputInteger) * CHAR_BIT;
const double log_d_temp = std::log2(static_cast<double>(divisor));
const uint8_t log_d = std::ceil(log_d_temp);
if (log_d == std::floor(log_d_temp))
{
return {1, log_d};
}
OutputInteger res = std::ceil(static_cast<double>(static_cast<OutputInteger>(1) << (log_d + n)) / double(divisor));
return {res, n + log_d};
}
// somewhere in the main function
for(uint8_t divisor = 1; divisor > 0; divisor++)
{
auto [multiplier, shift] = getDivisionMultiplier<uint8_t, uint16_t>(divisor);
for(uint8_t numenator = 1; numenator > 0; numenator++)
{
uint32_t res = static_cast<uint32_t>(numenator * multiplier) >> shift;
if (res != numenator / divisor)
{
std::cout << "panic: did something went wrong?" << std::endl;
}
}
}Наверное, очевидно, что фразу про панику мы никогда не увидим – значит всё? Работает и статью пора заканчивать? – Нет.
Что делает компилятор
Компилятор, точно делает по-другому, действительно, выведем, что код, приведенный выше дает для d=10:
auto p = getDivisionMultiplier<uint8_t, uint16_t>(static_cast<uint8_t>(10));
std::cout << p.first << " " << (uint16_t)(p.second) << std::endl;410 12
А из куска Assembly из начала статьи понятно, что должно было быть 205 и 11... Дело в том, что gcc хочет получить константу m того же размера, что и d и использует для этого более хитрый алгоритм (если вы приглядитесь к моему коду выше – я предусмотрительно использовал тип uint16_t).
Алгоритм gcc основан на подсчёте пары констант m_low = trunc(2^(ceil(log(d)) + N) / d) и m_high = trunc(2^(ceil(log(d)) + N) + 2^(ceil(log(d)) / d), тут важно понимать что m_low – не может быть настоящим m (показывал это выше), а m_high – может (извините, я не буду это тоже пытаться тут расписать как и все оставшиеся выкладки), но m_high, вообще говоря, даже больше m из моего наивного алгоритма. Да, но, и m_low и m_high – точно меньше чем 2^(N + 1), то есть для нашего случая это было бы 9 бит, и ещё точно в целых числах m_high > m_low (без равенства). Что же делает этот алгоритм дальше, чтобы сделать из 9-ти битного числа 8-ми битное число? Правильно, он просто сдвинет m_high на разряд вправо: тут важно понимать, что таким образом он делает trunc(m_high/2), и конечно параллельно, он уменьшит k (число сдвигов направо в конечном коде), но ...если число нечетное, это же не тоже самое? Да, в этом случае есть риск что мы начнём аппроксимировать снизу...поэтому, компилятор так делает и trunc(m_low / 2) и сравнивает эти два числа, потому что m_low уже изначально слишком мала – если m_high деградировала до неё – то нельзя дальше уменьшать m_high и соответствующий сдвиг. А вот и код делающий это:
template<typename InputInteger>
std::tuple<InputInteger, uint8_t, bool> getDivisionMultiplier(InputInteger divisor)
{
if (!divisor)
{
throw std::invalid_argument("Division by zero is impossible");
}
if (divisor == 1)
{
return {1, 0, false};
}
constexpr uint8_t n = sizeof(InputInteger) * CHAR_BIT;
const double log_d_temp = std::log2(static_cast<double>(divisor));
const uint8_t log_d = std::ceil(log_d_temp);
if (log_d == std::floor(log_d_temp))
{
return {1, log_d, false};
}
uint64_t temp_low = (1UL << (log_d + n));
uint64_t temp_hight = (1UL << log_d) | (1UL << (log_d + n));
temp_hight /= divisor;
temp_low /= divisor;
uint8_t additionla_shift = log_d;
while (additionla_shift)
{
if (temp_low /2 >= temp_hight/2)
{
break;
}
temp_low /= 2;
temp_hight /= 2;
--additionla_shift;
}
return {temp_hight, n + additionla_shift, temp_hight > std::numeric_limits<uint8_t>::max()};
}
// somewhere in the main function
auto [coeff, shift, _] = getDivisionMultiplier(static_cast<uint8_t>(10));
std::cout << (uint16_t)coeff << " " << (uint16_t)(shift) << std::endl;205 11
Успех! Коэффициенты совпали! Всё ли на этом? Увы: нет гарантии, что цикл, который редуцирует temp_hight, не выйдет сразу после первой итерации. То есть нет гарантии получить 8-ми битное число, но тогда у нас возможен срез по модулю в return. Gcc умеет успешно использовать это срезанное значение – но это явно отдельная тема.
А если очень интересно или просто не терпится, то я оставлю ссылки, на которые я опирался
Комментарии (27)
checkpoint
02.08.2024 18:47+4То, что деление по возможности заменяеться умножением это давно не секрет. В приведенном сниппете меня же интересует почему компилятор вставил две инструкции сдвига
shr ax, 8иshr al, 3вместо однойshr ax, 11?
StarPilgrim Автор
02.08.2024 18:47Интересный вопрос, я это тоже заметил, пока не могу ответить, скажу только, что их может быть даже три сдвига, посмотрите деление на 7, например, но там понятно почему их нельзя заменить общим сдвигом
lamerok
02.08.2024 18:47+1Предположу, что связано с установками флагов операций для целых разной длины, например флаг SF должен скинуться для результата сдвига именно 8 битного числа. Так как устанавливается в значение старшего бита именно 8 битного числа в al, а не изначально 16 битного ax.
StarPilgrim Автор
02.08.2024 18:47+1Вопрос можно считать закрытым: при флаге -O2 компилятор генерирует как раз shr ax, 11
Panzerschrek
02.08.2024 18:47+12Автор, приводи ассемблерный листинг с -O2. Так гораздо более понятно - всё лишнее удалено.
_Z5div10h: mov edx, -51 mov eax, edi mul dl shr ax, 11 ret
santanico
02.08.2024 18:47+12Для тех, кто дочитал статью до конца, но не нашел абзаца про то, зачем вообще весь этот геморрой) Дело в том, что в универсальных процессорах операция ‘div’ в несколько раз медленнее, чем ‘mul’. Почему так, это тема отдельной статьи)
Циферки по операциям можно посмотреть тут: https://www.agner.org/optimize/instruction_tables.pdf
titbit
02.08.2024 18:47Можно еще упомянуть, что во-первых знаковое деление обрабатывается немного по-другому, нежели беззнаковое. Во-вторых, похожий трюк с умножением на мультипликативное обратное используется и при взятии остатка от деления.
p.s. умножение в свою очередь можно заменить на сложения (вычитания) и сдвиги - даже на x86 компилятор так иногда делает, ну и вообще бывает полезно на контроллерах, где ни команды деления, ни команды умножения нет.
StarPilgrim Автор
02.08.2024 18:47Да, в целом согласен, просто для одной статьи перебор, а в микроконтроллерах я не очень силён
iShrimp
02.08.2024 18:47+1Извиняюсь за оффтоп. Подскажите, существует ли быстрый целочисленный (fixed point) алгоритм вычисления обратного квадратного корня (для нормализации 8- и 16-битных векторов)? Не могу найти ничего специально для фиксированной арифметики, везде предлагают вычислять через float.
StarPilgrim Автор
02.08.2024 18:47+1Судя по всему существуют, вот что я нашёл: https://github.com/chmike/fpsqrt
Но вот мое (спонтанное) мнение: сдвинуть влево на чётное число разрядов, что б получить целое - потом взять примерный корень, и сдвинуть обратно вправо, на половину первоначальных разрядов
iShrimp
02.08.2024 18:47Спасибо, теперь понятно, в каком направлении искать! Stackoverflow знает всё (почти)
askv
02.08.2024 18:47Я так понял, они реализуют алгоритм вычисления квадратного корня «столбиком» (я о нём узнал ещё в школе где-то в году 1990-м от учительницы математики), только в двоичной системе он сильно проще, чем в десятичной.
https://ru.wikipedia.org/wiki/%D0%9A%D0%B2%D0%B0%D0%B4%D1%80%D0%B0%D1%82%D0%BD%D1%8B%D0%B9_%D0%BA%D0%BE%D1%80%D0%B5%D0%BD%D1%8C#%D0%A1%D1%82%D0%BE%D0%BB%D0%B1%D0%B8%D0%BA%D0%BE%D0%BC
StarPilgrim Автор
02.08.2024 18:47Я ерунду сказал насчёт сдвигов: там и так целое, сдвигать надо только для повышения точности когда число не квадрат (видимо) а потом надо учесть все же сдвинуть потому что если исходно число например было 2 знака после точки (т.е. количество 0.25 по сути) после корня оно будет число 0.5 - надо додвигать. А так взятие квадратного корня в бинарной системе вроде не оч сложно - тоже почитал вчера
askv
02.08.2024 18:47+1Да, в десятичной системе перебор из 10 вариантов. В двоичной выбор из двух, что реализуется простым сравнением. Я собственно, этот алгоритм для целого числа и реализовал на ассемблере в студенческие годы (писал о нём в первом комментарии к этой статье).
Я ерунду сказал насчёт сдвигов:
Я сразу понял, что ерунда, но не стал комментировать, с этим разобраться не так уж и сложно. Если извлекать таким алгоритмом корень из вещественного числа, то нужно извлечь корень из мантиссы, а порядок разделить пополам (сдвиг на один бит). Понятно, что там есть определённые детали, но это уже нужно смотреть на конкретный формат записи вещественного числа (там есть неявная единица в начале мантиссы, также мантиссу нужно сдвинуть на 1 бит, если порядок нечётный, и т.п.), нужно выдать ошибку, если знак числа отрицательный и т.п. Но это уже совсем технические детали...
StarPilgrim Автор
02.08.2024 18:47Есть формальное описание как брать этот корень в двоичной системе, но мне оно не очень нравится, лучше покрутить в голове примеры
avdx
02.08.2024 18:47+3Какие то сложные вычисления. На сколько я могу судить обратное значение делителя и сдвиг вычисляются достаточно просто. Для деления N-битного числа X на число D (при условии, что D не является степенью двойки), нам нужно F бит дробной части fixed point представления обратного значения D:
F = N + Nd - 1, где Nd - число бит необходимое для представления числа D
fixed point представление для обратного значения делителя D
R = [2^F / D] + 1, где [] - взятие целой части
соответственно для примера из статьи, N = 8, D =10
Nd = 4
F = 8 + 4 - 1 = 11
R = [2048/10] + 1 = [204.8] + 1 = 205
StarPilgrim Автор
02.08.2024 18:47Подумаю над этим на днях
avdx
02.08.2024 18:47+2Ну вообще я не прав. Для D = 11 приведенный алгоритм не работает. Точное условие, при котором деление будет точным для всех X, это выполнение неравенства:
2^F > Z/(1 + (1/D - (R - 1)/2^F) * Z), где Z = D * (2^N - 1) (вывод приводить здесь не буду)
для N = 8, D =11, это неравенство выполняется для F=12, а для F=11 уже не выполняется. Но тут получается интересная ситуация, в этом случае нам потребуется больше чем 2*N бит для представления R * X, а это может уже приводить к тому, что такую оптимизацию делением для некоторых чисел выполнить просто невозможно, потому что, если например для 32 битных чисел потребуется больше 64 бит для промежуточного результата, то на 64-битной машине за одно умножение результат уже не получишь.
StarPilgrim Автор
02.08.2024 18:47А вы можете дать какие-то ссылки где подробнее можно почитать? Я вот так сходу не могу врубиться откуда формула и тд
avdx
02.08.2024 18:47+1Не знаю, сам вывел. Написал тестовую программку для проверки, формула работает. Рассуждал примерно так:
Перевод дробного числа в представление с фиксированной точкой производится путем его умножения на 2^F (обозначим M = 2^F) и отбрасывания дробной части. В нашем случае мы умножаем на обратное значение для нашего делителя (D) равное 1/D. Пусть его представление с фиксированной точкой будет Rfixed = [1/D * M] ([] - взятие целой части).
Реальное дробное значение этого представления будет R = Rfixed / M
Таким образом 1/D = R + e, где (e) ошибка округления, причем 0 <= e < 1/M
Соответественно R = 1/D - e, а e = 1/D - R
Когда мы умножаем любое число X на R то получаем:
X * R = X/D - e * X
Таким образом вычисленное значение будет всегда меньше реального и при округлении у нас может получаться значение меньше на 1 чем должно быть.
Поэтому прибавим к Rfixed единицу:
Rcfixed = Rfixed + 1
Его дробное значение будет соответственно
Rc = Rcfixed / M = R + 1/M = X / D + X * (1/M - e)
обозначим
ec = 1/M - e, при этом 0 < ec <= 1/M
Соответственно
X * Rс = X/D + eс * X
Нам нужно обеспечить выполнение условия:
X * Rс < [X/D] + 1
иначе при умножении X на Rc результат будет больше на 1, чем должно быть. Отсюда:
X/D + ec * X < [X/D] + 1
Заметим, что максимальным значением для X/D - [X/D] будет (D-1)/D
Таким образом:
(D-1)/D + eс * X < 1, откуда:
eс * X < 1/D
Это условие будет выполняться для любого X, если оно выполняется для максимального значения X = Xmax, т.е.
eс * Xmax < 1/D, откуда:
eс < 1/(D * Xmax), обозначив Z = D * Xmax
eс < 1/Z
Здесь можно поступить просто, используя факт, что ec <= 1/M, обеспечить выполнение условия
1/M < 1/Z, откуда:
M > Z, или
2^F > Z
Для примера из статьи, где D = 10, а Xmax = 255, Z = D * Xmax = 2550, получаем F=12:
2^12 = 4096 > 2550
Но, в этом случае будет нужно больше 2*N бит для промежуточного результата, это в общем случае проблема, поэтому нужно проверять более строгое ограничение:
eс < 1/Z
1/M - e < 1/Z
M * (1/Z + e) > 1
M > Z/(1 + e * Z)
M > Z/(1 + (1/D - R) * Z) (это неравенство можно преобразовать в (M + (M - Rfixed * D) * Xmax) > Z и в коде лучше использовать его, т.к. тут все вычисления целочисленные)
то видим что, для F=11
M = 2048
R = 204
Z/(1 + (1/D - R) * Z) = 2550/(1 + (1/10 - 204/2048)*2550) = 1277.4951076320849
это условие выполняется, а для F = 10
M = 1024
R = 102
Z/(1 + (1/D - R) * Z) = 2550/(1 + (1/10 - 102/1024)*2550) = 1277.4951076320849
уже не выполняется, т.ч. минимальное значение F=11Как я уже писал в предыдущем комментарии, есть числа, при оптимизации деления на которые указанным способом, нужно больше чем 2 * N бит (для них формула дает F = N + Nd). Для таких чисел компилятор использует другую оптимизацию. Вот пример для деления 32-битных чисел на 25 и 27:
#include <cstdint> uint32_t div25(uint32_t x) { return x / 25; } uint32_t div27(uint32_t x) { return x / 27; }clang генерит:
div25(unsigned int): mov eax, edi imul rax, rax, 1374389535 shr rax, 35 ret div27(unsigned int): mov eax, edi imul rax, rax, 795364315 shr rax, 32 sub edi, eax shr edi add eax, edi shr eax, 4 retStarPilgrim Автор
02.08.2024 18:47Да, очень круто для самостоятельного вывода...Эти констатнты они влезают в 32 бита, насколько я вижу. Там я в конце писал, что по алгоритму из статьи может M не влезть в базовый тип, но тогда он срежется по модулю...пример
div27очень похож на те случаиavdx
02.08.2024 18:47Для 27 не влезает. Для 27
F = 37
R = 5090331611А вот если взять F = 34, то R = 636291452, а 795364315 = 636291452 * 5/4. Суть метода для 27 я до конца не понял, конкретно что дает умножение на этот коэффициент и как потом от него избавляются.
askv
Когда-то давно написал программу на ассемблере по извлечению квадратного корня, побитовыми операциями (сдвиги, сравнения и т.п.).