bit manipulation - bitwise - Cómo desentrelazar eficientemente los bits(Morton inverso)

bitwise operators c

Question

Esta pregunta: ¿Cómo desentrelazar bits (UnMortonizing?) 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 incluso bits) en tan pocas operaciones como sea posible.

Para mi uso necesitaría tomar un int de 32 bits y extraer dos ints de 16 bits, donde uno es el bit de par y el otro es el bit de impar desplazado a la derecha en 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 que hay muchas soluciones que utilizan turnos y máscaras con números mágicos para generar números Morton (es decir, bits de intercalación), por ejemplo, bits de intercalación por números mágicos binarios , pero todavía no he encontrado nada para hacer lo contrario (es decir, desentrelazado) .

ACTUALIZAR

Después de volver a leer la sección de Hacker''s Delight sobre el orden aleatorio perfecto, 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 puede mejorarse, tanto en su forma escalar actual como también aprovechando SIMD, por lo que todavía estoy interesado en mejores soluciones (ya sea escalar o SIMD).

Answer 1

Código para las CPUs Intel Haswell y posteriores. Puede utilizar el conjunto de instrucciones BMI2 que contiene las instrucciones pext y pdep. Estas pueden (entre otras grandes cosas) ser usadas 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);
}

Answer 2

En caso de que alguien esté usando códigos morton en 3d, entonces necesita leer un bit cada 3, y aquí está la función de 64 bits que utilicé:

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)

Answer 3

No quería limitarme a un entero de tamaño fijo y hacer listas de comandos similares con constantes codificadas, por lo que desarrollé una solución C ++ 11 que utiliza la metaprogramación de plantillas para generar las funciones y las constantes. El código de ensamblaje generado con -O3 parece lo más ajustado posible sin utilizar 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 fuente repo y demostración en vivo .

Implementación

Básicamente, cada paso en la función morton1 funciona cambiando y agregando una secuencia de constantes que se ven así:

1. 0b0101010101010101 (alternativa 1 y 0)
2. 0b0011001100110011 (alterna 2x 1 y 0)
3. 0b0000111100001111 (alternativo 4x 1 y 0)
4. 0b0000000011111111 (alterna 8x 1 y 0)

Si tuviéramos que usar las dimensiones D , tendríamos un patrón con ceros D-1 y 1 uno. Así que para generar estos es suficiente generar unos consecutivos y aplicar algunos bitwise 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 recursión, 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 la pregunta "¿cuántos pasos necesitamos?". Esto depende también del número de dimensiones. En general, los pasos k calculan 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, es suficiente tomar 1 + log_2 z pasos 1 + log_2 z . El problema ahora es que necesitamos esto como constexpr para usarlo como un parámetro de plantilla. La mejor manera de log2 este problema es definir log2 través de la 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 finalmente, podemos realizar una sola 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);
}

Answer 4

Si necesita velocidad, puede utilizar la búsqueda de tablas para la conversión de un byte a la vez (la tabla de dos bytes es más rápida pero más grande). El procedimiento se realiza bajo el IDE de Delphi pero el ensamblador / algoritmo 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;

Answer 5

Si su procesador maneja eficientemente 64 bits, podría combinar las operaciones ...

int64 w = (z &0xAAAAAAAA)<<31 | (z &0x55555555 )
w = (w | (w >> 1)) & 0x3333333333333333;
w = (w | (w >> 2)) & 0x0F0F0F0F0F0F0F0F; 
...