Cómo desentrelazar bits de manera eficiente (Morton inverso)
https://stackoverflow.com/questions/4909263
-
29-10-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPregunta
Esta pregunta: Cómo desentrelazar bits (¿DesMortonizar?) tiene una buena respuesta para extraer una de las dos mitades de un número de Morton (solo los bits impares), pero necesito una solución que extraiga ambas partes (los bits impares y los pares) en el menor número de operaciones posible.
Para mi uso, necesitaría tomar un int de 32 bits y extraer dos int de 16 bits, donde uno son los bits pares y el otro son los bits impares desplazados a la derecha 1 bit, por ejemplo.
input, z: 11101101 01010111 11011011 01101110
output, x: 11100001 10110111 // odd bits shifted right by 1
y: 10111111 11011010 // even bits
Parece haber muchas soluciones que utilizan turnos y máscaras con números mágicos para generar números de Morton (es decir,bits entrelazados), p.e. Intercalar bits mediante números mágicos binarios, pero todavía no he encontrado nada para hacer lo contrario (es decir,desintercalado).
ACTUALIZAR
Después de volver a leer la sección de Hacker's Delight sobre mezclas/desordenes perfectas, encontré algunos ejemplos útiles que adapté de la siguiente manera:
// 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);
}
Creo que esto aún se puede mejorar, tanto en su forma escalar actual como aprovechando SIMD, por lo que todavía estoy interesado en mejores soluciones (ya sean escalares o SIMD).
Solución
Si su procesador maneja un INTS de 64 bits de manera eficiente, podría combinar las operaciones ...
int64 w = (z &0xAAAAAAAA)<<31 | (z &0x55555555 )
w = (w | (w >> 1)) & 0x3333333333333333;
w = (w | (w >> 2)) & 0x0F0F0F0F0F0F0F0F;
...
Otros consejos
Código para Intel Haswell y luego CPUs.Puede usar el conjunto de instrucciones BMI2 que contiene las instrucciones de PEXT y PDEP.Estos pueden (entre otras grandes cosas) ser utilizadas para construir sus funciones.
#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);
}
En caso de que alguien esté usando códigos de Morton en 3D, por lo que necesita leer un bit cada 3 y 64 bits aquí es la función que utilizo:
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)
Si necesita velocidad de lo que puede usar la búsqueda de la tabla para una conversión de bytes a la vez (la tabla de dos bytes es más rápida pero a la grande).El procedimiento está hecho bajo Delphi IDE, pero el ensamblador / algorithem es el mismo.
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;
No quería limitarme a un número entero de tamaño fijo y hacer listas de comandos similares con constantes codificadas, así que desarrollé una solución C++11 que utiliza metaprogramación de plantillas para generar las funciones y las constantes.El código ensamblador generado con -O3 Parece lo más ajustado posible sin usar el IMC:
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
TL;DR repositorio de origen y demo en vivo.
Implementación
Básicamente cada paso en el morton1 La función funciona cambiando y agregando a una secuencia de constantes que se ven así:
0b0101010101010101(alternar 1 y 0)0b0011001100110011(alternar 2x 1 y 0)0b0000111100001111(alternar 4x 1 y 0)0b0000000011111111(alternar 8x 1 y 0)
Si tuviéramos que usar D dimensiones, tendríamos un patrón con D-1 ceros y 1 uno.Entonces para generarlos basta con generar unos consecutivos y aplicar algunos bit a bit o:
/// @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))> {
};
Ahora que podemos generar las constantes en tiempo de compilación para dimensiones arbitrarias con lo siguiente:
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)>;
Con el mismo tipo de recursividad podemos generar funciones para cada uno de los pasos del algoritmo. 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
Queda por responder a la pregunta "¿cuántos pasos necesitamos?".Esto también depende del número de dimensiones.En general, k pasos computar 2^k - 1 bits de salida;el número máximo de bits significativos para cada dimensión viene dado por z = sizeof(T) * 8 / dimensions, por lo tanto basta con tomar 1 + log_2 z pasos.El problema ahora es que necesitamos esto como constexpr para usarlo como parámetro de plantilla.La mejor manera que encontré para solucionar esto es definir log2 a través de metaprogramación:
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>;
Y por último, podemos realizar una única llamada:
/// @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);
}
