Эффективное разделение беззнакового значения на две силы, округление
Я хочу реализовать целочисленное деление без знака a на произвольную степень двух, округляя вверх, эффективно. Итак, математически я хочу ceiling(p/q)
0. В C реализация соломона, которая не использует ограниченный домен q
, может быть похожа на следующую функцию 1:
/** q must be a power of 2, although this version works for any q */
uint64_t divide(uint64_t p, uint64_t q) {
uint64_t res = p / q;
return p % q == 0 ? res : res + 1;
}
... конечно, я действительно не хочу использовать разделение или мод на уровне машины, так как это занимает много циклов даже на современном оборудовании. Я ищу сокращение силы, которое использует сдвиги и/или некоторые другие дешевые операции (операции) - воспользовавшись тем, что q
является степенью 2.
Вы можете предположить, что у нас есть эффективная функция lg(unsigned int x)
, которая возвращает журнал base-2 x
, если x
является степенью двух.
Undefined поведение прекрасное, если q
равно нулю.
Обратите внимание, что простое решение: (p+q-1) >> lg(q)
не работает вообще - попробуйте его с помощью p == 2^64-100 and q == 256
2 например.
Сведения о платформе
Я интересуюсь решениями на C, которые, вероятно, будут хорошо работать на разных платформах, но ради конкретности, присуждая награду, и поскольку любое окончательное обсуждение производительности должно включать целевую архитектуру, Будем говорить о том, как я их протежу:
- Сервер Skylake
-
gcc 5.4.0
с флагами компиляции -O3 -march=haswell
Использование встроенных gcc (таких как встроенные функции bitcan/leading zero) отлично, и, в частности, я реализовал функцию lg()
, которую я сказал, был доступен следующим образом:
inline uint64_t lg(uint64_t x) {
return 63U - (uint64_t)__builtin_clzl(x);
}
inline uint32_t lg32(uint32_t x) {
return 31U - (uint32_t)__builtin_clz(x);
}
Я проверил, что они скомпилированы до одной инструкции bsr
, по крайней мере с -march=haswell
, несмотря на кажущееся участие вычитания. Вы, конечно, можете игнорировать их и использовать любые другие встроенные функции, которые вы хотите в своем решении.
Benchmark
Я написал критерий для существующих ответов и поделился и обновил результаты по мере внесения изменений.
Написание хорошего ориентира для небольшой, потенциально связанной операции довольно сложно. Когда код встроен в сайт вызова, большая часть работы функции может исчезнуть, особенно когда она находится в цикле 3.
Вы могли бы просто избежать всей проблемы вложения, гарантируя, что ваш код не вставлен: объявите его в другой блок компиляции. Я попытался сделать это с помощью bench
двоичного кода, но на самом деле результаты довольно бессмысленны. Почти все реализации привязаны к 4 или 5 циклам за вызов, но даже фиктивный метод, который ничего не делает, кроме return 0
, принимает одно и то же время. Таким образом, вы в основном просто измеряете накладные расходы call + ret
. Кроме того, вы почти никогда не будете использовать такие функции, как это, если вы не испортились, они будут доступны для вложения и что все это изменит.
Итак, два теста, на которых я буду больше всего фокусироваться, многократно вызывают тестируемый метод в цикле, позволяя встраивание, кросс-функциональную оптизацию, подъем петли и даже векторизации.
Существует два основных типа тестов: латентность и пропускная способность. Ключевое различие заключается в том, что в тесте латентности каждый вызов divide
зависит от предыдущего вызова, поэтому общие вызовы не могут быть легко перекрыты 4:
uint32_t bench_divide_latency(uint32_t p, uint32_t q) {
uint32_t total = p;
for (unsigned i=0; i < ITERS; i++) {
total += divide_algo(total, q);
q = rotl1(q);
}
return total;
}
Обратите внимание, что работающий total
зависит от вывода каждого деления, а также является вызовом вызова divide
.
Вариант пропускной способности, с другой стороны, не подает вывод одного деления на последующий. Это позволяет работать на одном вызове с последующим дублированием (как компилятором, так и главным образом процессором) и даже допускает векторизация:
uint32_t bench_divide_throughput(uint32_t p, uint32_t q) {
uint32_t total = p;
for (unsigned i=0; i < ITERS; i++) {
total += fname(i, q);
q = rotl1(q);
}
return total;
}
Обратите внимание, что здесь мы корнем в цикле счетчика как дивиденд - это переменная, но она не зависит от предыдущего деления.
Кроме того, каждый тест имеет три варианта поведения для делителя, q
:
- Постоянный делитель времени компиляции. Например, вызов
divide(p, 8)
. Это обычно на практике, и код может быть намного проще, когда делитель известен во время компиляции.
- Инвариантный делитель. Здесь делитель не знает во время компиляции, но является постоянным для всего цикла бенчмаркинга. Это позволяет подмножество оптимизаций, которые выполняет константа времени компиляции.
- Переменный делитель. На каждой итерации цикла изменяется делитель. Контрольные функции выше показывают этот вариант, используя команду "rotate left 1" для изменения делителя.
Объединяя все, что вы получаете в общей сложности 6 различных тестов.
Результаты
В целом
В целях выбора общего наилучшего алгоритма я просмотрел каждый из 12 подмножеств для предложенных алгоритмов: (latency, throughput) x (constant a, invariant q, variable q) x (32-bit, 64-bit)
и присвоил оценку 2, 1 или 0 на каждый подтест следующим образом:
- Лучший алгоритм (с точностью до 5%) получает оценку 2.
- "Близкие" алгоритмы (не более 50% медленнее, чем лучшие) получают оценку 1.
- Остальные алгоритмы оценивают нуль.
Следовательно, максимальный общий балл равен 24, но алгоритм не достиг этого. Вот общие итоговые результаты:
╔═══════════════════════╦═══════╗
║ Algorithm ║ Score ║
╠═══════════════════════╬═══════╣
║ divide_user23_variant ║ 20 ║
║ divide_chux ║ 20 ║
║ divide_user23 ║ 15 ║
║ divide_peter ║ 14 ║
║ divide_chrisdodd ║ 12 ║
║ stoke32 ║ 11 ║
║ divide_chris ║ 0 ║
║ divide_weather ║ 0 ║
╚═══════════════════════╩═══════╝
Таким образом, для целей этого конкретного тестового кода, с этим конкретным компилятором и на этой платформе, user2357112 "вариант" (с ... + (p & mask) != 0
) лучше всего работает, связанный с предложением chux (что на самом деле является идентичным кодом).
Вот все подзадачи, которые суммируются с приведенным выше:
╔══════════════════════════╦═══════╦════╦════╦════╦════╦════╦════╗
║ ║ Total ║ LC ║ LI ║ LV ║ TC ║ TI ║ TV ║
╠══════════════════════════╬═══════╬════╬════╬════╬════╬════╬════╣
║ divide_peter ║ 6 ║ 1 ║ 1 ║ 1 ║ 1 ║ 1 ║ 1 ║
║ stoke32 ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║
║ divide_chux ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║
║ divide_user23 ║ 8 ║ 1 ║ 1 ║ 2 ║ 2 ║ 1 ║ 1 ║
║ divide_user23_variant ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║
║ divide_chrisdodd ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║
║ divide_chris ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║
║ divide_weather ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║
║ ║ ║ ║ ║ ║ ║ ║ ║
║ 64-bit Algorithm ║ ║ ║ ║ ║ ║ ║ ║
║ divide_peter_64 ║ 8 ║ 1 ║ 1 ║ 1 ║ 2 ║ 2 ║ 1 ║
║ div_stoke_64 ║ 5 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 1 ║
║ divide_chux_64 ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║
║ divide_user23_64 ║ 7 ║ 1 ║ 1 ║ 2 ║ 1 ║ 1 ║ 1 ║
║ divide_user23_variant_64 ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║
║ divide_chrisdodd_64 ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║
║ divide_chris_64 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║
║ divide_weather_64 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║
╚══════════════════════════╩═══════╩════╩════╩════╩════╩════╩════╝
Здесь каждый тест называется XY, с X в {Latency, Throughput} и Y в {Constant Q, Invariant Q, Variable Q}. Так, например, LC - это "Тест латентности с константой q".
Анализ
На самом высоком уровне решения можно условно разделить на две категории: быстрые (верхние 6 финишеров) и медленные (нижние два). Разница больше: все быстрые алгоритмы были самыми быстрыми, по крайней мере, на двух подтестах, и вообще, когда они не закончили сначала, они попали в категорию "достаточно близко" (только исключения были неудачными векторизацией в случае stoke
и chrisdodd
). Однако медленные алгоритмы оценивали 0 (даже не закрывали) при каждом тесте. Таким образом, вы можете в основном устранить медленные алгоритмы из дальнейшего рассмотрения.
Автоматическая векторизация
Среди быстрых алгоритмов большой дифференциатор - это возможность авто-векторизация.
Ни один из алгоритмов не смог авто-векторизовать в тестах на задержку, что имеет смысл, поскольку тесты задержки Vec?, показывающим, включена ли функция или нет:
Size Vec? Name
045 N bench_c_div_stoke_64
049 N bench_c_divide_chrisdodd_64
059 N bench_c_stoke32_64
212 Y bench_c_divide_chux_64
227 Y bench_c_divide_peter_64
220 Y bench_c_divide_user23_64
212 Y bench_c_divide_user23_variant_64
Тенденция ясна - векторизованные функции занимают около 4x размер не-вексельных. Это связано с тем, что сами петли ядра больше (векторные инструкции имеют тенденцию быть больше и их больше), а потому, что настройка цикла и особенно код после цикла намного больше: например, для векторизованной версии требуется сокращение до суммируем все частичные суммы в векторе. Счетчик циклов фиксирован и кратен 8, поэтому код хвоста не генерируется - но если бы они были переменными, сгенерированный код был бы еще больше.
Кроме того, несмотря на значительное улучшение во время выполнения, векторизация gcc
фактически невелика. Здесь выдержка из векторизованной версии программы Питера:
on entry: ymm4 == all zeros
on entry: ymm5 == 0x00000001 0x00000001 0x00000001 ...
4007a4: c5 ed 76 c4 vpcmpeqd ymm0,ymm2,ymm4
4007ad: c5 fd df c5 vpandn ymm0,ymm0,ymm5
4007b1: c5 dd fa c0 vpsubd ymm0,ymm4,ymm0
4007b5: c5 f5 db c0 vpand ymm0,ymm1,ymm0
Этот фрагмент работает независимо от 8 DWORD
элементов, происходящих из ymm2
. Если мы возьмем x
как единственный элемент DWORD
ymm2
, а y
входящее значение ymm1
, эти команды foud соответствуют:
x == 0 x != 0
x = x ? 0 : -1; // -1 0
x = x & 1; // 1 0
x = 0 - x; // -1 0
x = y1 & x; // y1 0
Таким образом, первые три команды могут быть просто заменены первым, так как состояния в обоих случаях одинаковы. Так что два цикла добавлены к этой цепочке зависимостей (которая не переносится циклом) и два дополнительных uops. Очевидно, что фазы оптимизации gcc
как-то плохо взаимодействуют с кодом векторизации здесь, так как такие тривиальные оптимизации редко пропускаются в скалярном коде. Рассмотрение других векторизованных версий аналогичным образом показывает большую производительность, упавшую на пол.
Филиалы vs Без отраслевых
Почти все решения, скомпилированные для кода без ветвей, даже если код C имел условные выражения или явные ветки. Условные части были достаточно малы, чтобы компилятор вообще решил использовать условное перемещение или какой-то вариант. Одним из исключений является chrisdodd
, который скомпилирован с веткой (проверка if p == 0
) во всех тестах пропускной способности, но ни одна из латентных. Здесь типичный пример из теста пропускной способности константы q
:
0000000000400e60 <bench_c_divide_chrisdodd_32>:
400e60: 89 f8 mov eax,edi
400e62: ba 01 00 00 00 mov edx,0x1
400e67: eb 0a jmp 400e73 <bench_c_divide_chrisdodd_32+0x13>
400e69: 0f 1f 80 00 00 00 00 nop DWORD PTR [rax+0x0]
400e70: 83 c2 01 add edx,0x1
400e73: 83 fa 01 cmp edx,0x1
400e76: 74 f8 je 400e70 <bench_c_divide_chrisdodd_32+0x10>
400e78: 8d 4a fe lea ecx,[rdx-0x2]
400e7b: c1 e9 03 shr ecx,0x3
400e7e: 8d 44 08 01 lea eax,[rax+rcx*1+0x1]
400e82: 81 fa 00 ca 9a 3b cmp edx,0x3b9aca00
400e88: 75 e6 jne 400e70 <bench_c_divide_chrisdodd_32+0x10>
400e8a: c3 ret
400e8b: 0f 1f 44 00 00 nop DWORD PTR [rax+rax*1+0x0]
В ветки 400e76
пропущен случай, когда p == 0
. Фактически, компилятор мог просто очистить первую итерацию (вычислять ее результат явно), а затем полностью избежать прыжка, так как после этого он может доказать, что p != 0
. В этих тестах ветвь совершенно предсказуема, что может дать преимущество коду, который фактически компилируется с использованием ветки (поскольку код сравнения и ветвления существенно не соответствует линии и близок к свободному), и это большая часть того, почему chrisdodd
выигрывает пропускная способность, переменный случай q.
Подробные результаты испытаний
Здесь вы можете найти некоторые подробные результаты тестов и некоторые детали самих тестов.
Задержка
Результаты ниже проверяют каждый алгоритм на протяжении 1е9 итераций. Циклы вычисляются просто путем умножения времени/вызова на тактовую частоту. Обычно можно предположить, что что-то вроде 4.01
совпадает с 4.00
, но большие отклонения, такие как 5.11
, кажутся реальными и воспроизводимыми.
Результаты для divide_plusq_32
используют (p + q - 1) >> lg(q)
, но показаны только для справки, так как эта функция не работает при больших p + q
. Результаты для dummy
являются очень простой функцией: return p + q
и позволяют оценить эксплуатационные издержки контрольного пакета 5 (само добавление должно пройти максимум не более).
==============================
Bench: Compile-time constant Q
==============================
Function ns/call cycles
divide_peter_32 2.19 5.67
divide_peter_64 2.18 5.64
stoke32_32 1.93 5.00
stoke32_64 1.97 5.09
stoke_mul_32 2.75 7.13
stoke_mul_64 2.34 6.06
div_stoke_32 1.94 5.03
div_stoke_64 1.94 5.03
divide_chux_32 1.55 4.01
divide_chux_64 1.55 4.01
divide_user23_32 1.97 5.11
divide_user23_64 1.93 5.00
divide_user23_variant_32 1.55 4.01
divide_user23_variant_64 1.55 4.01
divide_chrisdodd_32 1.95 5.04
divide_chrisdodd_64 1.93 5.00
divide_chris_32 4.63 11.99
divide_chris_64 4.52 11.72
divide_weather_32 2.72 7.04
divide_weather_64 2.78 7.20
divide_plusq_32 1.16 3.00
divide_plusq_64 1.16 3.00
divide_dummy_32 1.16 3.00
divide_dummy_64 1.16 3.00
==============================
Bench: Invariant Q
==============================
Function ns/call cycles
divide_peter_32 2.19 5.67
divide_peter_64 2.18 5.65
stoke32_32 1.93 5.00
stoke32_64 1.93 5.00
stoke_mul_32 2.73 7.08
stoke_mul_64 2.34 6.06
div_stoke_32 1.93 5.00
div_stoke_64 1.93 5.00
divide_chux_32 1.55 4.02
divide_chux_64 1.55 4.02
divide_user23_32 1.95 5.05
divide_user23_64 2.00 5.17
divide_user23_variant_32 1.55 4.02
divide_user23_variant_64 1.55 4.02
divide_chrisdodd_32 1.95 5.04
divide_chrisdodd_64 1.93 4.99
divide_chris_32 4.60 11.91
divide_chris_64 4.58 11.85
divide_weather_32 12.54 32.49
divide_weather_64 17.51 45.35
divide_plusq_32 1.16 3.00
divide_plusq_64 1.16 3.00
divide_dummy_32 0.39 1.00
divide_dummy_64 0.39 1.00
==============================
Bench: Variable Q
==============================
Function ns/call cycles
divide_peter_32 2.31 5.98
divide_peter_64 2.26 5.86
stoke32_32 2.06 5.33
stoke32_64 1.99 5.16
stoke_mul_32 2.73 7.06
stoke_mul_64 2.32 6.00
div_stoke_32 2.00 5.19
div_stoke_64 2.00 5.19
divide_chux_32 2.04 5.28
divide_chux_64 2.05 5.30
divide_user23_32 2.05 5.30
divide_user23_64 2.06 5.33
divide_user23_variant_32 2.04 5.29
divide_user23_variant_64 2.05 5.30
divide_chrisdodd_32 2.04 5.30
divide_chrisdodd_64 2.05 5.31
divide_chris_32 4.65 12.04
divide_chris_64 4.64 12.01
divide_weather_32 12.46 32.28
divide_weather_64 19.46 50.40
divide_plusq_32 1.93 5.00
divide_plusq_64 1.99 5.16
divide_dummy_32 0.40 1.05
divide_dummy_64 0.40 1.04
Пропускная способность
Вот результаты тестов пропускной способности. Обратите внимание, что многие из алгоритмов здесь были авто-векторизованы, поэтому производительность относительно хороша для тех, кто во многих случаях является частью цикла. Одним из результатов является то, что в отличие от большинства результатов латентности 64-битные функции значительно медленнее, поскольку векторизация более эффективна при меньших размерах элементов (хотя разрыв больше, чем я ожидал).
==============================
Bench: Compile-time constant Q
==============================
Function ns/call cycles
stoke32_32 0.39 1.00
divide_chux_32 0.15 0.39
divide_chux_64 0.53 1.37
divide_user23_32 0.14 0.36
divide_user23_64 0.53 1.37
divide_user23_variant_32 0.15 0.39
divide_user23_variant_64 0.53 1.37
divide_chrisdodd_32 1.16 3.00
divide_chrisdodd_64 1.16 3.00
divide_chris_32 4.34 11.23
divide_chris_64 4.34 11.24
divide_weather_32 1.35 3.50
divide_weather_64 1.35 3.50
divide_plusq_32 0.10 0.26
divide_plusq_64 0.39 1.00
divide_dummy_32 0.08 0.20
divide_dummy_64 0.39 1.00
==============================
Bench: Invariant Q
==============================
Function ns/call cycles
stoke32_32 0.48 1.25
divide_chux_32 0.15 0.39
divide_chux_64 0.48 1.25
divide_user23_32 0.17 0.43
divide_user23_64 0.58 1.50
divide_user23_variant_32 0.15 0.38
divide_user23_variant_64 0.48 1.25
divide_chrisdodd_32 1.16 3.00
divide_chrisdodd_64 1.16 3.00
divide_chris_32 4.35 11.26
divide_chris_64 4.36 11.28
divide_weather_32 5.79 14.99
divide_weather_64 17.00 44.02
divide_plusq_32 0.12 0.31
divide_plusq_64 0.48 1.25
divide_dummy_32 0.09 0.23
divide_dummy_64 0.09 0.23
==============================
Bench: Variable Q
==============================
Function ns/call cycles
stoke32_32 1.16 3.00
divide_chux_32 1.36 3.51
divide_chux_64 1.35 3.50
divide_user23_32 1.54 4.00
divide_user23_64 1.54 4.00
divide_user23_variant_32 1.36 3.51
divide_user23_variant_64 1.55 4.01
divide_chrisdodd_32 1.16 3.00
divide_chrisdodd_64 1.16 3.00
divide_chris_32 4.02 10.41
divide_chris_64 3.84 9.95
divide_weather_32 5.40 13.98
divide_weather_64 19.04 49.30
divide_plusq_32 1.03 2.66
divide_plusq_64 1.03 2.68
divide_dummy_32 0.63 1.63
divide_dummy_64 0.66 1.71
a По крайней мере, указав unsigned, мы избегаем целой банки червей, связанной с правым сдвигом поведения целых чисел со знаком в C и С++.
0 Конечно, эта нотация на самом деле не работает на C, где /
обрезает результат, поэтому ceiling
ничего не делает. Поэтому рассмотрим, что псевдо-нотация, а не прямая C.
1 Я также интересуюсь решениями, где все типы uint32_t
, а не uint64_t
.
2 В общем случае любые p
и q
, где p + q >= 2^64
вызывает проблему из-за переполнения.
3 Тем не менее, функция должна быть в цикле, потому что производительность микроскопической функции, которая занимает полдюжины циклов, действительно имеет значение, если она вызывается в довольно жесткой петле.
4 Это немного упрощает - только дивиденд p
зависит от вывода предыдущей итерации, поэтому некоторые работы, связанные с обработкой q
, все еще могут перекрываться.
5 Используйте такие оценки с осторожностью, однако накладные расходы не просто аддитивны. Если накладные расходы отображаются в 4 циклах, а какая-то функция f
занимает 5, вероятно, неточно сказать, что стоимость реальной работы в f
равна 5 - 4 == 1
из-за перекрытия выполнения.
Ответы
Ответ 1
Этот ответ о том, какой идеал в asm; что мы хотели бы убедить компилятор испустить для нас. (Я не предлагаю на самом деле использовать inline asm, за исключением точки сравнения при выводе компилятора для сравнения. https://gcc.gnu.org/wiki/DontUseInlineAsm).
Мне удалось получить довольно хороший выход asm из чистого C для ceil_div_andmask
, см. ниже. (Это хуже, чем CMOV на Broadwell/Skylake, но, вероятно, хорошо на Haswell. версия user23/chux выглядит даже лучше для обоих случаев.) Это в основном просто стоит упомянуть как один из немногих случаев, когда я получил компилятор для извлечения asm, который я хотел.
Похоже, что идея Криса Додда о return ((p-1) >> lg(q)) + 1
с обработкой специального случая для d = 0 является одним из лучших вариантов. То есть оптимальная реализация его в asm трудно превзойти с оптимальной реализацией чего-либо еще. У Chux (p >> lg(q)) + (bool)(p & (q-1))
также есть преимущества (например, более низкая латентность от результата p- > ) и более CSE, когда один и тот же q используется для нескольких разделов, Ниже приведен анализ задержки/пропускной способности того, что gcc делает с ним.
Если один и тот же e = lg(q)
повторно используется для нескольких дивидендов или один и тот же дивиденд используется повторно для нескольких делителей, различные реализации могут CSE больше выражения. Они также могут эффективно векторизовать с помощью AVX2.
Филиалы дешевы и очень эффективны, если они прогнозируют очень хорошо, поэтому разветвление на d==0
будет лучше, если оно почти никогда не будет принято. Если d==0
не является редкими, нераспространяющийся asm будет работать лучше в среднем. В идеале мы можем написать что-то в C, которое позволит gcc сделать правильный выбор во время оптимизации на основе профиля и скомпилировать его с хорошим asm для любого случая.
Так как лучшие безветровые реализации asm не добавляют много латентности против разветвленной реализации, веткистая, вероятно, способ пойти, если ветка не будет идти одинаково, может быть, в 99% случаев. Вероятно, это возможно для ветвления на p==0
, но, вероятно, менее вероятно для ветвления на p & (q-1)
.
Трудно направить gcc5.4 на испускание чего-либо оптимального. Это моя работа над Godbolt).
Я думаю, что оптимальные последовательности для Skylake для этого алгоритма заключаются в следующем. (Представлены как автономные функции для AMD64 SysV ABI, но говорят о пропускной способности/задержке в предположении, что компилятор будет выпустить что-то подобное в виде петли без привязки RET).
Ветвь на переносе от вычисления d-1
для обнаружения d == 0 вместо отдельного теста и ветки. Уменьшает счетчик uop красиво, esp на семействе SnB, где JC может сжимать макрос с SUB.
ceil_div_pjc_branch:
xor eax,eax ; can take this uop off the fast path by adding a separate xor-and-return block, but in reality we want to inline something like this.
sub rdi, 1
jc .d_was_zero ; fuses with the sub on SnB-family
tzcnt rax, rsi ; tzcnt rsi,rsi also avoids any false-dep problems, but this illustrates that the q input can be read-only.
shrx rax, rdi, rax
inc rax
.d_was_zero:
ret
- Fused-domain uops: 5 (не считая ret), а один из них - xor-zero (без исполнительного блока)
- Задержка HSW/SKL с успешным предсказанием ветвления:
- (d == 0): нет зависимости от данных d или q, разрывает цепочку dep. (управляющая зависимость от d для обнаружения ошибочных прогнозов и выхода из ветки).
- (d!= 0): q- > result: tzcnt + shrx + inc = 5c
- (d!= 0): d- > результат: sub + shrx + inc = 3c
- Пропускная способность: возможно, просто узкое место по пропускной способности uop
Я пробовал, но не смог получить gcc для перехода на CF из вычитания, но он всегда хочет провести отдельное сравнение. Я знаю, что gcc можно уговорить на ветвление на CF после вычитания двух переменных, но, возможно, это не удается, если вы являетесь константой времени компиляции. (IIRC, это обычно компилируется в CF-тест с неподписанными vars: foo -= bar; if(foo>bar) carry_detected = 1;
)
Без ветвления с ADC/SBB для обработки события d==0
. Zero-handling добавляет только одну инструкцию к критическому пути (по сравнению с версией без специальной обработки для d == 0), но также преобразует еще одну из INC в sbb rax, -1
, чтобы сделать CF отмененным -= -1
. Использование CMOV дешевле на pre-Broadwell, но требует дополнительных инструкций по настройке.
ceil_div_pjc_asm_adc:
tzcnt rsi, rsi
sub rdi, 1
adc rdi, 0 ; d? d-1 : d. Sets CF=CF
shrx rax, rdi, rsi
sbb rax, -1 ; result++ if d was non-zero
ret
- Fused-domain uops: 5 (не считая ret) на SKL. 7 по HSW
- Задержка SKL:
- q- > результат: tzcnt + shrx + sbb = 5c
- d- > result: sub + adc + shrx (здесь начинается деление на q) + sbb = 4c
- Пропускная способность: TZCNT работает на p1. SBB, ADC и SHRX работают только на p06. Таким образом, я думаю, что мы сглаживаем 3 раза за p06 за итерацию, делая этот запуск в лучшем случае одной итерацией за 1.5c.
Если q и d станут готовыми одновременно, обратите внимание, что эта версия может запускать SUB/ADC параллельно с задержкой 3c для TZCNT. Если оба они поступают из одной и той же кешированной кэш-линии, это, безусловно, возможно. В любом случае преимущество в определении q как можно позже в цепочке зависимостей результатов d- > является преимуществом.
Получение этого из С представляется маловероятным с помощью gcc5.4. Существует встроенная функция add-with-carry, но gcc создает полный беспорядок. Он не использует непосредственные операнды для АЦП или SBB и сохраняет перенос в целочисленную регистрацию между каждой операцией. gcc7, clang3.9 и icc17 делают из этого ужасный код.
#include <x86intrin.h>
// compiles to completely horrible code, putting the flags into integer regs between ops.
T ceil_div_adc(T d, T q) {
T e = lg(q);
unsigned long long dm1; // unsigned __int64
unsigned char CF = _addcarry_u64(0, d, -1, &dm1);
CF = _addcarry_u64(CF, 0, dm1, &dm1);
T shifted = dm1 >> e;
_subborrow_u64(CF, shifted, -1, &dm1);
return dm1;
}
# gcc5.4 -O3 -march=haswell
mov rax, -1
tzcnt rsi, rsi
add rdi, rax
setc cl
xor edx, edx
add cl, -1
adc rdi, rdx
setc dl
shrx rdi, rdi, rsi
add dl, -1
sbb rax, rdi
ret
CMOV, чтобы исправить весь результат: худшая латентность от результата q- > , так как она раньше использовалась в цепочке отрезка результата d- > .
ceil_div_pjc_asm_cmov:
tzcnt rsi, rsi
sub rdi, 1
shrx rax, rdi, rsi
lea rax, [rax+1] ; inc preserving flags
cmovc rax, zeroed_register
ret
- Fused-domain uops: 5 (не считая ret) на SKL. 6 по HSW
- Задержка SKL:
- q- > результат: tzcnt + shrx + lea + cmov = 6c (хуже, чем ADC/SBB на 1c)
- d- > результат: sub + shrx (здесь начинается от q) + lea + cmov = 4c
- Пропускная способность: TZCNT работает на p1. LEA - p15. CMOV и SHRX - p06. SUB - p0156. Теоретически только узкое место в пропускной способности uop с плавным доменом, поэтому одна итерация на 1.25c. При наличии большого количества независимых операций конфликты ресурсов из SUB или LEA, кража p1 или p06, не должны быть проблемой пропускной способности, поскольку при 1 инерции на 1.25c ни один порт не насыщен uops, который может работать только на этом порту.
CMOV, чтобы получить операнд для SUB: я надеялся, что смогу найти способ создания операнда для более поздней инструкции, которая при необходимости создавала бы нуль без зависимости ввода от q, e, или результат SHRX. Это помогло бы, если d
готов до q
или в то же время.
Это не достигает этой цели и нуждается в дополнительном 7-байтовом mov rdx,-1
в цикле.
ceil_div_pjc_asm_cmov:
tzcnt rsi, rsi
mov rdx, -1
sub rdi, 1
shrx rax, rdi, rsi
cmovnc rdx, rax
sub rax, rdx ; res += d ? 1 : -res
ret
версия с низкой задержкой для процессоров с предварительным BDW с дорогостоящим CMOV, используя SETCC для создания маски для AND.
ceil_div_pjc_asm_setcc:
xor edx, edx ; needed every iteration
tzcnt rsi, rsi
sub rdi, 1
setc dl ; d!=0 ? 0 : 1
dec rdx ; d!=0 ? -1 : 0 // AND-mask
shrx rax, rdi, rsi
inc rax
and rax, rdx ; zero the bogus result if d was initially 0
ret
Тем не менее 4c задержка от результата d- > (и 6 из результата q- > ), поскольку SETC/DEC происходит параллельно с SHRX/INC. Total uop count: 8. Большинство из этих insns могут работать на любом порту, поэтому он должен быть 1 итера за 2 такта.
Конечно, для pre-HSW вам также необходимо заменить SHRX.
Мы можем получить gcc5.4, чтобы испустить что-то почти хорошее: (по-прежнему использует отдельный TEST вместо установки маски на основе sub rdi, 1
, но в остальном те же инструкции, что и выше). Смотрите в Godbolt.
T ceil_div_andmask(T p, T q) {
T mask = -(T)(p!=0); // TEST+SETCC+NEG
T e = lg(q);
T nonzero_result = ((p-1) >> e) + 1;
return nonzero_result & mask;
}
Когда компилятор знает, что p
отличен от нуля, он использует преимущество и делает хороший код:
// http://stackoverflow.com/questions/40447195/can-i-hint-the-optimizer-by-giving-the-range-of-an-integer
#if defined(__GNUC__) && !defined(__INTEL_COMPILER)
#define assume(x) do{if(!(x)) __builtin_unreachable();}while(0)
#else
#define assume(x) (void)(x) // still evaluate it once, for side effects in case anyone is insane enough to put any inside an assume()
#endif
T ceil_div_andmask_nonzerop(T p, T q) {
assume(p!=0);
return ceil_div_andmask(p, q);
}
# gcc5.4 -O3 -march=haswell
xor eax, eax # gcc7 does tzcnt in-place instead of wasting an insn on this
sub rdi, 1
tzcnt rax, rsi
shrx rax, rdi, rax
add rax, 1
ret
Chux/user23_variant
только 3c латентность от результата p- > , а константа q может много CSE.
T divide_A_chux(T p, T q) {
bool round_up = p & (q-1); // compiles differently from user23_variant with clang: AND instead of
return (p >> lg(q)) + round_up;
}
xor eax, eax # in-place tzcnt would save this
xor edx, edx # target for setcc
tzcnt rax, rsi
sub rsi, 1
test rsi, rdi
shrx rdi, rdi, rax
setne dl
lea rax, [rdx+rdi]
ret
Выполнение SETCC перед TZCNT позволит использовать TZCNT на месте, сохраняя xor eax,eax
. Я не смотрел, как это происходит в цикле.
- Fused-domain uops: 8 (не считая ret) на HSW/SKL
- Задержка HSW/SKL:
- q- > result: (tzcnt + shrx (p) | sub + test (p) + setne) + lea (или add) = 5c
- d- > результат: test (здесь начинается от q) + setne + lea = 3c. (цепочка shrx- > lea короче и, следовательно, не критический путь)
- Пропускная способность: возможно, только узкое место на интерфейсе, на один раз за 2c. Сохранение
xor eax,eax
должно ускорить это до 1 на 1.75c (но, конечно, любые накладные расходы цикла будут частью узкого места, потому что узкие места в интерфейсе подобны).
Ответ 2
uint64_t exponent = lg(q);
uint64_t mask = q - 1;
// v divide
return (p >> exponent) + (((p & mask) + mask) >> exponent)
// ^ round up
Отдельное вычисление части "round up" позволяет избежать проблем с переполнением (p+q-1) >> lg(q)
. В зависимости от того, насколько интеллектуальным является ваш компилятор, может быть возможно выразить часть "round up" как ((p & mask) != 0)
без разветвления.
Ответ 3
Эффективный способ деления на 2 для беззнакового целого числа в C является правым сдвигом сдвига вправо, разделяемым на два (округление вниз), поэтому смещение справа на n делит на 2 n (округление).
Теперь вы хотите округлить, а не вниз, что вы можете сделать, сначала добавив 2 n -1 или эквивалентно вычитая один перед сдвигом и добавив один после (кроме 0). Это работает примерно так:
unsigned ceil_div(unsigned d, unsigned e) {
/* compute ceil(d/2**e) */
return d ? ((d-1) >> e) + 1 : 0;
}
Условие можно удалить, используя логическое значение d для сложения и вычитания единицы:
unsigned ceil_div(unsigned d, unsigned e) {
/* compute ceil(d/2**e) */
return ((d - !!d) >> e) + !!d;
}
Из-за его размера и требования к скорости, функция должна быть сделана статической встроенной. Вероятно, это не будет отличаться для оптимизатора, но параметры должны быть const. Если он должен быть общим для многих файлов, определите его в заголовке:
static inline unsigned ceil_div(const unsigned d, const unsigned e){...
Ответ 4
Эффективное разделение неподписанного значения на степень двух, округление
[Повторно написать], учитывая OP разъяснение относительно мощности-2.
Завершающая или потолочная часть проста, когда переполнение не вызывает беспокойства. Просто добавьте q-1
, затем сдвиньте.
В противном случае, поскольку возможность округления зависит от всех бит p
меньше, чем q
, обнаружение этих бит необходимо прежде, чем они будут сдвинуты.
uint64_t divide_A(uint64_t p, uint64_t q) {
bool round_up = p & (q-1);
return (p >> lg64(q)) + round_up;
}
Это предполагает, что код имеет эффективную функцию lg64(uint64_t x)
, которая возвращает журнал базового 2 x
, если x
является степенью двух.
Ответ 5
Мой старый ответ не сработал, если p
был чем-то большим, чем силой двух (крики). Итак, мое новое решение, использующее функции __builtin_ctzll()
и __builtin_ffsll()
0 доступные в gcc
(который в качестве бонуса предоставляет быстрый логарифм, который вы упомянули!):
uint64_t divide(uint64_t p,uint64_t q) {
int lp=__builtin_ffsll(p);
int lq=__builtin_ctzll(q);
return (p>>lq)+(lp<(lq+1)&&lp);
}
Обратите внимание, что это предполагает, что a long long
- 64 бит. Это должно быть изменено немного иначе.
Идея здесь в том, что если нам нужно переполнение тогда и только тогда, когда p
имеет меньше конечных нулей, чем q
. Обратите внимание, что для мощности двух число конечных нулей равно логарифму, поэтому мы также можем использовать этот встроенный для журнала.
Часть &&lp
предназначена только для случая угла, где p
равно нулю: в противном случае здесь будет выводиться 1
.
0 Нельзя использовать __builtin_ctzll()
для обоих, потому что это undefined, если p==0
.
Ответ 6
Если дивиденд/делитель может быть гарантированно не превышать 63 (или 31) бита, вы можете использовать следующую версию, указанную в вопросе. Обратите внимание, как p + q может переполняться, если они используют все 64 бит. Это было бы хорошо, если команда SHR сдвигалась в флагом переноса, но AFAIK это не делает.
uint64_t divide(uint64_t p, uint64_t q) {
return (p + q - 1) >> lg(q);
}
Если эти ограничения не могут быть гарантированы, вы можете просто сделать метод пола, а затем добавить 1, если он будет округлен. Это можно определить, проверяя, находятся ли какие-либо биты в дивиденде в пределах диапазона делителя.
Примечание: p & -p извлекает младший бит набора на машинах с дополнением 2s или инструкцию BLSI
uint64_t divide(uint64_t p, uint64_t q) {
return (p >> lg(q)) + ( (p & -p ) < q );
}
Какой clang компилируется в:
bsrq %rax, %rsi
shrxq %rax, %rdi, %rax
blsiq %rdi, %rcx
cmpq %rsi, %rcx
adcq $0, %rax
retq
Это немного словесно и использует некоторые новые инструкции, поэтому, возможно, есть способ использовать флаг переноса в исходной версии. Посмотрим:
RCR инструкция, но похоже, что было бы хуже... возможно SHRD... Это было бы что-то вроде этого (невозможно проверить на данный момент)
xor edx, edx ;edx = 0 (will store the carry flag)
bsr rcx, rsi ;rcx = lg(q) ... could be moved anywhere before shrd
lea rax, [rsi-1] ;rax = q-1 (adding p could carry)
add rax, rdi ;rax += p (handle carry)
setc dl ;rdx = carry flag ... or xor rdx and setc
shrd rax, rdx, cl ;rax = rdx:rax >> cl
ret
Еще одна инструкция, но должна быть совместима со старыми процессорами (если она работает... Я всегда получаю исходный/целевой обмен, не стесняйтесь редактировать)
Приложение:
Я реализовал функцию lg(), которую я сказал, был доступен следующим образом:
inline uint64_t lg(uint64_t x) {
return 63U - (uint64_t)__builtin_clzl(x);
}
inline uint32_t lg32(uint32_t x) {
return 31U - (uint32_t)__builtin_clz(x);
}
Функции быстрого журнала не полностью оптимизируются для bsr на clang и ICC, но вы можете это сделать:
#if defined(__x86_64__) && (defined(__clang__) || defined(__INTEL_COMPILER))
static inline uint64_t lg(uint64_t x){
inline uint64_t ret;
//other compilers may want bsrq here
__asm__("bsr %0, %1":"=r"(ret):"r"(x));
return ret;
}
#endif
#if defined(__i386__) && (defined(__clang__) || defined(__INTEL_COMPILER)) static inline uint32_t lg32(uint32_t x){
inline uint32_t ret;
__asm__("bsr %0, %1":"=r"(ret):"r"(x));
return ret;
}
#endif
Ответ 7
К этой проблеме уже приложено много человеческих мозговых возможностей с несколькими вариантами отличных ответов в C
вместе с ответом Питера Кордеса, который охватывает все, на что можно надеяться в asm, с заметками о попытке сопоставить его вернуться к C
.
Таким образом, пока люди пишут в банке, я подумал, что нужно сказать о какой-то грубой вычислительной мощности! С этой целью я использовал Stanford STOKE superoptimizer, чтобы попытаться найти хорошие решения для 32-разрядных и 64-разрядных версий этой проблемы.
Обычно, супероптимизатор обычно является чем-то вроде поиска грубой силы через все возможные последовательности команд, пока вы не найдете лучший по какой-либо метрике. Конечно, с что-то вроде 1000 инструкций, которое быстро выйдет из-под контроля для нескольких команд 1. STOKE, на руке, использует ориентированный рандомизированный подход: он случайным образом превращает мутации в существующую кандидатную программу, оценивая на каждом шаге функцию стоимости, которая выполняет как эффективность, так и правильность. В любом случае однострочный лайнер - много бумаг, если это вызвало ваше любопытство.
Итак, через несколько минут STOKE нашел несколько довольно интересных решений. Он нашел почти все идеи высокого уровня в существующих решениях, а также несколько уникальных. Например, для 32-битной функции STOKE нашел эту версию:
neg rsi
dec rdi
pext rax, rsi, rdi
inc eax
ret
Он не использует никаких указаний о переходе/переходе на нуль или сдвиге вообще. В значительной степени он использует neg esi
, чтобы превратить делитель в маску с 1s в высоких битах, а затем pext
эффективно выполняет сдвиг с использованием этой маски. Вне этого трюка он использует тот же трюк, что и пользователь QuestionC: декремент p, shift, increment p - но он работает даже для нулевого дивиденда, потому что он использует 64-битные регистры, чтобы отличить нулевой регистр от MSB-набора больших p
case.
Я добавил версию C этого алгоритма к эталону и добавил его к результатам. Он конкурирует с другими хорошими алгоритмами, привязывая их сначала к "переменным Q" . Он выполняет векторизация, но не так же хорошо, как и другие 32-битные алгоритмы, поскольку для этого требуется 64-битная математика, и поэтому векторы могут обрабатывать только половину количества элементов одновременно.
Еще лучше, в 32-битном случае он придумал множество решений, которые используют тот факт, что некоторые интуитивные решения, которые терпят неудачу для краевых случаев, "просто работают", если вы используете 64-битные операционные системы для части из этого. Например:
tzcntl ebp, esi
dec esi
add rdi, rsi
sarx rax, rdi, rbp
ret
Это эквивалент предложения return (p + q - 1) >> lg(q)
, упомянутого в вопросе. Это не работает вообще, так как для больших p + q
он переполняется, но для 32-разрядных p
и q
это решение отлично работает, занимая важные части в 64-битной версии. Преобразуйте это обратно в C с помощью некоторых приведений, и на самом деле выясняется, что использование lea
сделает добавление в одной команде 1
stoke_32(unsigned int, unsigned int):
tzcnt edx, esi
mov edi, edi ; goes away when inlining
mov esi, esi ; goes away when inlining
lea rax, [rsi-1+rdi]
shrx rax, rax, rdx
ret
Таким образом, это решение с тремя командами, когда оно встроено в нечто, которое уже имеет значения, равные нулю, в rdi
и rsi
. Для определения автономной функции требуются инструкции mov
для нулевого расширения, поскольку как работает ABI SysV x64.
Для 64-битной функции она не придумала ничего, что сдуло бы существующие решения, но придумал некоторые аккуратные вещи, например:
tzcnt r13, rsi
tzcnt rcx, rdi
shrx rax, rdi, r13
cmp r13b, cl
adc rax, 0
ret
Этот парень учитывает конечные нули обоих аргументов, а затем добавляет 1 к результату, если q
имеет меньше конечных нулей, чем p
, так как это нужно, чтобы округлить. Умно!
В целом, он понял идею, что вам нужно shl
tzcnt
очень быстро (как и большинство людей), а затем придумал массу других решений проблемы настройки результата для учета округления. Ему даже удалось использовать blsi
и bzhi
в нескольких решениях. Здесь было предложено 5-командное решение:
tzcnt r13, rsi
shrx rax, rdi, r13
imul rsi, rax
cmp rsi, rdi
adc rax, 0
ret
Это в основном подход "умножить и проверить" - возьмите усеченный res = p \ q
, умножьте его обратно и, если он отличается от p
, добавьте один: return res * q == p ? ret : ret + 1
. Круто. Однако не лучше, чем решения Peter. У STOKE, по-видимому, есть некоторые недостатки в расчете латентности - он считает, что у вышеупомянутого есть латентность 5 - но это больше похоже на 8 или 9 в зависимости от архитектуры. Поэтому он иногда сужается в решениях, основанных на его ошибочном вычислении латентности.
1 Интересно, хотя этот подход с грубой силой достигает своей выполнимости около 5-6 команд: если вы предполагаете, что можете урезать количество команд, чтобы сказать 300, исключив инструкции SIMD и x87, тогда вам понадобится ~ 28 дней, чтобы попробовать все 300 ^ 5
5 последовательностей команд со скоростью 1 000 000 кандидатов в секунду. Вы могли бы уменьшить это в 1000 раз с различными оптимизациями, что означает менее часа для 5-командных последовательностей и, возможно, неделю для 6-инструкций. Как это бывает, большинство лучших решений этой проблемы попадают в это окно с инструкциями 5-6...
2 Это будет медленным lea
, поэтому последовательность, найденная STOKE, по-прежнему оптимальна для того, для чего я оптимизирован, что было латентностью.
Ответ 8
Вы можете сделать это так, сравнивая разделение n / d
с делением на (n-1) / d
.
#include <stdio.h>
int main(void) {
unsigned n;
unsigned d;
unsigned q1, q2;
double actual;
for(n = 1; n < 6; n++) {
for(d = 1; d < 6; d++) {
actual = (double)n / d;
q1 = n / d;
if(n) {
q2 = (n - 1) / d;
if(q1 == q2) {
q1++;
}
}
printf("%u / %u = %u (%f)\n", n, d, q1, actual);
}
}
return 0;
}
Выход программы:
1 / 1 = 1 (1.000000)
1 / 2 = 1 (0.500000)
1 / 3 = 1 (0.333333)
1 / 4 = 1 (0.250000)
1 / 5 = 1 (0.200000)
2 / 1 = 2 (2.000000)
2 / 2 = 1 (1.000000)
2 / 3 = 1 (0.666667)
2 / 4 = 1 (0.500000)
2 / 5 = 1 (0.400000)
3 / 1 = 3 (3.000000)
3 / 2 = 2 (1.500000)
3 / 3 = 1 (1.000000)
3 / 4 = 1 (0.750000)
3 / 5 = 1 (0.600000)
4 / 1 = 4 (4.000000)
4 / 2 = 2 (2.000000)
4 / 3 = 2 (1.333333)
4 / 4 = 1 (1.000000)
4 / 5 = 1 (0.800000)
5 / 1 = 5 (5.000000)
5 / 2 = 3 (2.500000)
5 / 3 = 2 (1.666667)
5 / 4 = 2 (1.250000)
5 / 5 = 1 (1.000000)
Обновление
Я опубликовал ранний ответ на исходный вопрос, который работает, но не учитывал эффективность алгоритма или что делитель всегда имеет силу 2. Выполнение двух делений было бесполезно дорого.
Я использую 32-разрядный компилятор MSVC в 64-битной системе, поэтому нет никаких шансов, что я смогу обеспечить наилучшее решение для требуемой цели. Но это интересный вопрос, поэтому я попытался найти, что лучшее решение обнаружит бит 2 ** n делителя. Использование функции библиотеки log2
работало, но было настолько медленным. Выполнение моей собственной смены было намного лучше, но все же мой лучший C-решение -
unsigned roundup(unsigned p, unsigned q)
{
return p / q + ((p & (q-1)) != 0);
}
Мое встроенное 32-битное ассемблерное решение выполняется быстрее, но, конечно, это не ответит на вопрос. Я краду несколько циклов, предположив, что eax
возвращается как значение функции.
unsigned roundup(unsigned p, unsigned q)
{
__asm {
mov eax,p
mov edx,q
bsr ecx,edx ; cl = bit number of q
dec edx ; q-1
and edx,eax ; p & (q-1)
shr eax,cl ; divide p by q, a power of 2
sub edx,1 ; generate a carry when (p & (q-1)) == 0
cmc
adc eax,0 ; add 1 to result when (p & (q-1)) != 0
}
} ; eax returned as function value
Ответ 9
Это кажется эффективным и работает для подписанного, если ваш компилятор использует арифметические сдвиги вправо (обычно это правда).
#include <stdio.h>
int main (void)
{
for (int i = -20; i <= 20; ++i) {
printf ("%5d %5d\n", i, ((i - 1) >> 1) + 1);
}
return 0;
}
Используйте >> 2
для разделения на 4, >> 3
, чтобы разделить на 8, & ct. Эффективный lg
выполняет там работу.
Вы можете даже разделить на 1! >> 0