Как эффективно устранить чередование битов (обратный Мортон)
https://stackoverflow.com/questions/4909263
-
29-10-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Этот вопрос: Как устранить чередование битов (UnMortonizing?) есть хороший ответ на извлечение одной из двух половин числа Мортона (только нечетные биты), но мне нужно решение, которое извлекает обе части (нечетные и четные биты) за минимально возможное количество операций.
Для моего использования мне нужно было бы взять 32-битное целое число и извлечь два 16-битных целых числа, где один — четные биты, а другой — нечетные биты, смещенные вправо на 1 бит, например
input, z: 11101101 01010111 11011011 01101110
output, x: 11100001 10110111 // odd bits shifted right by 1
y: 10111111 11011010 // even bits
Кажется, существует множество решений, использующих сдвиги и маски с магическими числами для генерации чисел Мортона (т.перемежающиеся биты), например Чередование битов с помощью двоичных магических чисел, но я еще не нашел ничего, что можно было бы сделать наоборот (т.е.деперемежение).
ОБНОВЛЯТЬ
Перечитав раздел Hacker's Delight об идеальных перетасовках/перетасовках, я нашел несколько полезных примеров, которые адаптировал следующим образом:
// morton1 - extract even bits
uint32_t morton1(uint32_t x)
{
x = x & 0x55555555;
x = (x | (x >> 1)) & 0x33333333;
x = (x | (x >> 2)) & 0x0F0F0F0F;
x = (x | (x >> 4)) & 0x00FF00FF;
x = (x | (x >> 8)) & 0x0000FFFF;
return x;
}
// morton2 - extract odd and even bits
void morton2(uint32_t *x, uint32_t *y, uint32_t z)
{
*x = morton1(z);
*y = morton1(z >> 1);
}
Я думаю, что это еще можно улучшить, как в его нынешней скалярной форме, так и за счет использования преимуществ SIMD, поэтому меня все еще интересуют лучшие решения (скалярные или SIMD).
Решение
Если ваш процессор обрабатывает 64 бит эффективно, вы можете объединить операции ...
int64 w = (z &0xAAAAAAAA)<<31 | (z &0x55555555 )
w = (w | (w >> 1)) & 0x3333333333333333;
w = (w | (w >> 2)) & 0x0F0F0F0F0F0F0F0F;
...
Другие советы
код для Intel Haswell и более поздних процессоров.Вы можете использовать набор инструкций BMI2, который содержит инструкции PEXT и PDEP.Это может (среди прочего великих вещей) использоваться для создания ваших функций.
#include <immintrin.h>
#include <stdint.h>
// on GCC, compile with option -mbmi2, requires Haswell or better.
uint64_t xy_to_morton (uint32_t x, uint32_t y)
{
return _pdep_u32(x, 0x55555555) | _pdep_u32(y,0xaaaaaaaa);
}
uint64_t morton_to_xy (uint64_t m, uint32_t *x, uint32_t *y)
{
*x = _pext_u64(m, 0x5555555555555555);
*y = _pext_u64(m, 0xaaaaaaaaaaaaaaaa);
}
В случае, если кто-то использует Morton Codes в 3D, поэтому ему нужно читать один бит каждые 3, а 64 бита здесь - функция, которую я использовал:
uint64_t morton3(uint64_t x) {
x = x & 0x9249249249249249;
x = (x | (x >> 2)) & 0x30c30c30c30c30c3;
x = (x | (x >> 4)) & 0xf00f00f00f00f00f;
x = (x | (x >> 8)) & 0x00ff0000ff0000ff;
x = (x | (x >> 16)) & 0xffff00000000ffff;
x = (x | (x >> 32)) & 0x00000000ffffffff;
return x;
}
uint64_t bits;
uint64_t x = morton3(bits)
uint64_t y = morton3(bits>>1)
uint64_t z = morton3(bits>>2)
Если вам нужна скорость, чем вы можете использовать настольный поиск для одного байтового преобразования сразу (два байта таблицы быстрее, но для большой).Процедура выполнена под Delphi IDE, но ассемблер / алгоритм одинаковы.
const
MortonTableLookup : array[byte] of byte = ($00, $01, $10, $11, $12, ... ;
procedure DeinterleaveBits(Input: cardinal);
//In: eax
//Out: dx = EvenBits; ax = OddBits;
asm
movzx ecx, al //Use 0th byte
mov dl, byte ptr[MortonTableLookup + ecx]
//
shr eax, 8
movzx ecx, ah //Use 2th byte
mov dh, byte ptr[MortonTableLookup + ecx]
//
shl edx, 16
movzx ecx, al //Use 1th byte
mov dl, byte ptr[MortonTableLookup + ecx]
//
shr eax, 8
movzx ecx, ah //Use 3th byte
mov dh, byte ptr[MortonTableLookup + ecx]
//
mov ecx, edx
and ecx, $F0F0F0F0
mov eax, ecx
rol eax, 12
or eax, ecx
rol edx, 4
and edx, $F0F0F0F0
mov ecx, edx
rol ecx, 12
or edx, ecx
end;
Я не хотел ограничиваться целым числом фиксированного размера и составлять списки похожих команд с жестко запрограммированными константами, поэтому я разработал решение на C++11, которое использует метапрограммирование шаблонов для генерации функций и констант.Ассемблерный код, созданный с помощью -O3 кажется настолько точным, насколько это возможно без использования ИМТ:
andl $0x55555555, %eax
movl %eax, %ecx
shrl %ecx
orl %eax, %ecx
andl $0x33333333, %ecx
movl %ecx, %eax
shrl $2, %eax
orl %ecx, %eax
andl $0xF0F0F0F, %eax
movl %eax, %ecx
shrl $4, %ecx
orl %eax, %ecx
movzbl %cl, %esi
shrl $8, %ecx
andl $0xFF00, %ecx
orl %ecx, %esi
ТЛ;ДР исходный репозиторий и живая демонстрация.
Выполнение
По сути, каждый шаг в morton1 Функция работает путем сдвига и добавления к последовательности констант, которая выглядит следующим образом:
0b0101010101010101(чередовать 1 и 0)0b0011001100110011(поочередно 2x 1 и 0)0b0000111100001111(поочередно 4x 1 и 0)0b0000000011111111(поочередно 8x 1 и 0)
Если бы мы использовали D размеров, у нас будет шаблон с D-1 нули и 1 один.Итак, чтобы сгенерировать их, достаточно сгенерировать последовательные и применить некоторые побитовые или:
/// @brief Generates 0b1...1 with @tparam n ones
template <class T, unsigned n>
using n_ones = std::integral_constant<T, (~static_cast<T>(0) >> (sizeof(T) * 8 - n))>;
/// @brief Performs `@tparam input | (@tparam input << @tparam width` @tparam repeat times.
template <class T, T input, unsigned width, unsigned repeat>
struct lshift_add :
public lshift_add<T, lshift_add<T, input, width, 1>::value, width, repeat - 1> {
};
/// @brief Specialization for 1 repetition, just does the shift-and-add operation.
template <class T, T input, unsigned width>
struct lshift_add<T, input, width, 1> : public std::integral_constant<T,
(input & n_ones<T, width>::value) | (input << (width < sizeof(T) * 8 ? width : 0))> {
};
Теперь мы можем генерировать константы во время компиляции для произвольных измерений с помощью следующего:
template <class T, unsigned step, unsigned dimensions = 2u>
using mask = lshift_add<T, n_ones<T, 1 << step>::value, dimensions * (1 << step), sizeof(T) * 8 / (2 << step)>;
Используя тот же тип рекурсии, мы можем генерировать функции для каждого из шагов алгоритма. x = (x | (x >> K)) & M:
template <class T, unsigned step, unsigned dimensions>
struct deinterleave {
static T work(T input) {
input = deinterleave<T, step - 1, dimensions>::work(input);
return (input | (input >> ((dimensions - 1) * (1 << (step - 1))))) & mask<T, step, dimensions>::value;
}
};
// Omitted specialization for step 0, where there is just a bitwise and
Осталось ответить на вопрос «сколько нам нужно шагов?».Это также зависит от количества измерений.В общем, k шаги вычисления 2^k - 1 выходные биты;максимальное количество значимых битов для каждого измерения определяется выражением z = sizeof(T) * 8 / dimensions, поэтому достаточно взять 1 + log_2 z шаги.Проблема в том, что нам нужно это как constexpr чтобы использовать его в качестве параметра шаблона.Лучший способ обойти это — определить log2 с помощью метапрограммирования:
template <unsigned arg>
struct log2 : public std::integral_constant<unsigned, log2<(arg >> 1)>::value + 1> {};
template <>
struct log2<1u> : public std::integral_constant<unsigned, 0u> {};
/// @brief Helper constexpr which returns the number of steps needed to fully interleave a type @tparam T.
template <class T, unsigned dimensions>
using num_steps = std::integral_constant<unsigned, log2<sizeof(T) * 8 / dimensions>::value + 1>;
И, наконец, мы можем выполнить один единственный вызов:
/// @brief Helper function which combines @see deinterleave and @see num_steps into a single call.
template <class T, unsigned dimensions>
T deinterleave_first(T n) {
return deinterleave<T, num_steps<T, dimensions>::value - 1, dimensions>::work(n);
}
