Comments 28
Когда-то давно написал программу на ассемблере по извлечению квадратного корня, побитовыми операциями (сдвиги, сравнения и т.п.).
То, что деление по возможности заменяеться умножением это давно не секрет. В приведенном сниппете меня же интересует почему компилятор вставил две инструкции сдвига shr ax, 8 и shr al, 3 вместо одной shr ax, 11 ?
Сдвиг на 8 — это побайтовый сдвиг, может так работает быстрее? (Только предположение.)
Интересный вопрос, я это тоже заметил, пока не могу ответить, скажу только, что их может быть даже три сдвига, посмотрите деление на 7, например, но там понятно почему их нельзя заменить общим сдвигом
Предположу, что связано с установками флагов операций для целых разной длины, например флаг SF должен скинуться для результата сдвига именно 8 битного числа. Так как устанавливается в значение старшего бита именно 8 битного числа в al, а не изначально 16 битного ax.
Вопрос можно считать закрытым: при флаге -O2 компилятор генерирует как раз shr ax, 11
Автор, приводи ассемблерный листинг с -O2. Так гораздо более понятно - всё лишнее удалено.
_Z5div10h:
mov edx, -51
mov eax, edi
mul dl
shr ax, 11
ret
Для тех, кто дочитал статью до конца, но не нашел абзаца про то, зачем вообще весь этот геморрой) Дело в том, что в универсальных процессорах операция ‘div’ в несколько раз медленнее, чем ‘mul’. Почему так, это тема отдельной статьи)
Циферки по операциям можно посмотреть тут: https://www.agner.org/optimize/instruction_tables.pdf
Можно еще упомянуть, что во-первых знаковое деление обрабатывается немного по-другому, нежели беззнаковое. Во-вторых, похожий трюк с умножением на мультипликативное обратное используется и при взятии остатка от деления.
p.s. умножение в свою очередь можно заменить на сложения (вычитания) и сдвиги - даже на x86 компилятор так иногда делает, ну и вообще бывает полезно на контроллерах, где ни команды деления, ни команды умножения нет.
Извиняюсь за оффтоп. Подскажите, существует ли быстрый целочисленный (fixed point) алгоритм вычисления обратного квадратного корня (для нормализации 8- и 16-битных векторов)? Не могу найти ничего специально для фиксированной арифметики, везде предлагают вычислять через float.
Судя по всему существуют, вот что я нашёл: https://github.com/chmike/fpsqrt
Но вот мое (спонтанное) мнение: сдвинуть влево на чётное число разрядов, что б получить целое - потом взять примерный корень, и сдвинуть обратно вправо, на половину первоначальных разрядов
Спасибо, теперь понятно, в каком направлении искать! Stackoverflow знает всё (почти)
Я так понял, они реализуют алгоритм вычисления квадратного корня «столбиком» (я о нём узнал ещё в школе где-то в году 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
Я ерунду сказал насчёт сдвигов: там и так целое, сдвигать надо только для повышения точности когда число не квадрат (видимо) а потом надо учесть все же сдвинуть потому что если исходно число например было 2 знака после точки (т.е. количество 0.25 по сути) после корня оно будет число 0.5 - надо додвигать. А так взятие квадратного корня в бинарной системе вроде не оч сложно - тоже почитал вчера
Да, в десятичной системе перебор из 10 вариантов. В двоичной выбор из двух, что реализуется простым сравнением. Я собственно, этот алгоритм для целого числа и реализовал на ассемблере в студенческие годы (писал о нём в первом комментарии к этой статье).
Я ерунду сказал насчёт сдвигов:
Я сразу понял, что ерунда, но не стал комментировать, с этим разобраться не так уж и сложно. Если извлекать таким алгоритмом корень из вещественного числа, то нужно извлечь корень из мантиссы, а порядок разделить пополам (сдвиг на один бит). Понятно, что там есть определённые детали, но это уже нужно смотреть на конкретный формат записи вещественного числа (там есть неявная единица в начале мантиссы, также мантиссу нужно сдвинуть на 1 бит, если порядок нечётный, и т.п.), нужно выдать ошибку, если знак числа отрицательный и т.п. Но это уже совсем технические детали...
Какие то сложные вычисления. На сколько я могу судить обратное значение делителя и сдвиг вычисляются достаточно просто. Для деления 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
Подумаю над этим на днях
Ну вообще я не прав. Для 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-битной машине за одно умножение результат уже не получишь.
А вы можете дать какие-то ссылки где подробнее можно почитать? Я вот так сходу не могу врубиться откуда формула и тд
Не знаю, сам вывел. Написал тестовую программку для проверки, формула работает. Рассуждал примерно так:
Перевод дробного числа в представление с фиксированной точкой производится путем его умножения на 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
ret
Как поделить не деля или оптимизация деления компиляторам(и)