Come i bit de-interleavi in modo efficiente (morton inverso)
-
29-10-2019 - |
Domanda
Questa domanda: Come de-interleave Bits (smorzando?) Ha una buona risposta per estrarre una delle due metà di un numero di Morton (solo i bit dispari), ma ho bisogno di una soluzione che estrae entrambe le parti (i bit dispari e i bit uniformi) in poche operazioni possibili. Per il mio uso avrei bisogno di prendere un INT e estrarre un 20 bit ed estrarre due intercetti a 16 bit, dove si sono i bit uniformi e l'altro è i bit dispari spostati direttamente da 1 bit, ad esempio.
input, z: 11101101 01010111 11011011 01101110
output, x: 11100001 10110111 // odd bits shifted right by 1
y: 10111111 11011010 // even bits
.
Sembra che ci siano molte soluzioni utilizzando turni e maschere con numeri magici per generare numeri di morton (cioè bit interleaving), ad es. interleave bit di numeri magici binari , ma non ho ancora trovato nulla per fare il contrario (cioè de-interleaving).
Aggiornamento
Dopo aver riletturato la sezione dalla delizia del hacker su forze perfette / non trasparenti ho trovato alcuni esempi utili che ho adattato come segue:
.
// 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);
}
Penso che questo possa ancora essere migliorato, sia nella sua attuale forma scalare e anche sfruttando il vantaggio del Simd, quindi sono ancora interessato a soluzioni migliori (Scalare o SIMD).
Soluzione
Se il tuo processore gestisce in modo efficace 64 bit, è possibile combinare le operazioni ...
int64 w = (z &0xAAAAAAAA)<<31 | (z &0x55555555 )
w = (w | (w >> 1)) & 0x3333333333333333;
w = (w | (w >> 2)) & 0x0F0F0F0F0F0F0F0F;
...
. Altri suggerimenti
Codice per Intel Haswell e CPU successive.È possibile utilizzare il set di istruzioni BMI2 che contiene le istruzioni di PEXT e PDEP.Questi possono essere utilizzati (tra le altre grandi cose) per costruire le tue funzioni.
#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);
}
. Nel caso in cui qualcuno stia utilizzando i codici Morton in 3D, quindi ha bisogno di leggere un bit ogni 3, e 64 bit qui è la funzione che ho usato:
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)
. Se è necessario utilizzare la velocità di quella che è possibile utilizzare la ricerca della tabella per una conversione di un byte contemporaneamente (la tabella di due byte è più veloce ma a grande).La procedura è realizzata sotto Delphi IDE ma l'assembler / algorithem è lo stesso.
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;
. Non volevo essere limitato a un numero intero di dimensioni fisse e ad esempio elenchi di comandi simili con costanti hardcoded, quindi ho sviluppato una soluzione C ++ 11 che fa uso del metaprogrammazione del modello per generare le funzioni e le costanti. Il codice di assemblaggio generato con -O3
sembra stretto come può ottenere senza utilizzare BMI:
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 sorgente repo e demo live .
.
Implementazione
Fondamentalmente ogni passo nella funzione morton1
funziona spostando e aggiungendo a una sequenza di costanti che assomigliano a questo:
- .
-
0b0101010101010101
(alternativo 1 e 0) -
0b0011001100110011
(alternativo 2x 1 e 0) -
0b0000111100001111
(alternativo 4x 1 e 0) -
0b0000000011111111
(alternativo 8x 1 e 0)Se dovessimo usare dimensioni
D
, avremmo un modello con zeriD-1
e1
. Quindi per generare questi è abbastanza per generare quelli consecutivi e applicare un po 'di 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))> { };
Ora che possiamo generare le costanti al momento della compilazione per le dimensioni arbitrarie con quanto segue:
.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 lo stesso tipo di ricorsione, possiamo generare funzioni per ciascuna delle fasi dell'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
Resta per rispondere alla domanda "Quanti passaggi abbiamo bisogno?". Questo dipende anche dal numero di dimensioni. In generale,
k
Passaggi Compute2^k - 1
Bit di uscita; Il numero massimo di bit significativi per ciascuna dimensione è dato daz = sizeof(T) * 8 / dimensions
, quindi è sufficiente prendere gradini1 + log_2 z
. Il problema è ora che abbiamo bisogno di questo comeconstexpr
per usarlo come parametro modello. Il modo migliore che ho trovato per aggirare questo è per definirelog2
tramite metaprogrammazione:
.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>;
E infine, possiamo eseguire una singola chiamata:
./// @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); }