如何有效地去交错位（逆莫顿）

Posted 2023-02-19

【中文标题】如何有效地去交错位（逆莫顿）

【英文标题】：How to efficiently de-interleave bits (inverse Morton)

【发布时间】：2011-06-22 00:14:09

【问题描述】：

这个问题：How to de-interleave bits (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

似乎有很多解决方案使用带有幻数的移位和掩码来生成莫顿数（即交错位），例如Interleave bits by Binary Magic Numbers，但我还没有找到任何相反的方法（即去交错）。

更新

在重新阅读 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）感兴趣。

【问题讨论】：

您链接到的交织解决方案使用的操作是您链接到的去交织解决方案的两倍。如果这是可以接受的，您可以通过两次应用解交织解决方案来实现相同的性能。我认为没有比这更好的了，因为两种解决方案都使用相同的原理，并且有一半位为 0 的阶段，所以它们一次只能处理一半的信息，所以如果你想要所有信息你需要两个去。当然，如果你有 64 位整数，你可以一次性完成；然后您可以将其中一个奇偶校验位移动到高 32 位。

我又玩了一些——我没有想出更好的解决方案，但我确实做了一个有趣的观察：如果可以的话，您可以有效地将 AaBbCcDd.. 更改为 ABabCDcd..有效地将 0aB00cD0.. 更改为 0Ba00Dc0.. -- 因此您可以减少此步骤以有效交换两个位，这意味着实现映射 0->0、3->3、1->2、2->1。我能想到的对两位（mod 4）的可逆运算是：加 0、1、2 或 3，与 1 或 3 进行异或运算，或乘以 3。但这些只会生成 S_4 的 8 元素子组，它不会'不包括所需的排列。

我假设“交错操作”是指将 32 位字的高 16 位视为奇数位，将低 16 位视为偶数位，并通过交错获得新的 32 位字他们？抽象的答案是，是的，它是循环的，因为它是双射运算，并且存在有限数量的不同 32 位字 :-) 但更实际地说，循环长度是 5：交错操作循环二进制中的数字位索引的表示，最低有效位变为最高有效位，对于 32 位字，有 5 位数字循环。

我的另一个想法，开箱即用：您需要以正确顺序排列的奇数位和偶数位吗？或者您是否可以重组其余代码（例如通过使用不同的查找表），以便您可以以不同的顺序接受它们？因为让它们以不同的顺序排列很容易：odd = x & 0xaaaaaaaa;奇数 = (奇数 | (奇数 >>> 17)) & 0xffff;偶数 = x & 0x55555555;偶数 = (偶数 | (偶数 >>> 15)) & 0xffff;

@joriki：不幸的是，我需要以正确的顺序排列这些位 - 我将使用它们作为数组的索引，我需要按 Morton 顺序遍历该数组。

【参考方案1】：

如果您的处理器有效地处理 64 位整数，您可以组合这些操作...

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

【讨论】：

啊，是的——这是个好主意——我只会使用 64 位 CPU，所以这绝对值得。

【参考方案2】：

适用于 Intel Haswell 和更高版本 CPU 的代码。您可以使用包含 pext 和 pdep 指令的 BMI2 指令集。这些可以（以及其他很棒的东西）用于构建您的功能。

#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);

【参考方案3】：

如果有人在 3d 中使用 morton 代码，所以他需要每 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)

【参考方案4】：

如果您需要速度，则可以使用表查找一次进行一个字节转换（两个字节表更快但很大）。程序是在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;

【讨论】：

谢谢 - 我也在考虑一种表查找方法，但使用 SSSE3 的 PSHUFB 指令并行执行多个 4 位查找 - 除了并行性之外，这样做的好处是它不会t 需要内存访问来进行表查找。

很高兴看到这个解决方案... :)

【参考方案5】：

我不想局限于固定大小的整数并使用硬编码常量列出类似命令，因此我开发了一个 C++11 解决方案，它利用模板元编程来生成函数和常量。使用-O3 生成的汇编代码似乎在不使用 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 source repo 和 live demo。

---

实施

morton1 函数中的每个步骤基本上都是通过移位和添加到如下所示的一系列常量来工作的：

0b0101010101010101（备用 1 和 0） 0b0011001100110011（交替 2x 1 和 0） 0b0000111100001111（交替 4x 1 和 0） 0b0000000011111111（交替 8x 1 和 0）

如果我们使用D 维度，我们将有一个D-1 零和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);

【讨论】：

非常好的通用实现。

【参考方案6】：

您可以通过像这样相乘来提取 8 个交错位：

uint8_t deinterleave_even(uint16_t x) 
    return ((x & 0x5555) * 0xC00030000C0003 & 0x0600180060008001) * 0x0101010101010101 >> 56;

uint8_t deinterleave_odd(uint16_t x) 
    return ((x & 0xAAAA) * 0xC00030000C0003 & 0x03000C003000C000) * 0x0101010101010101 >> 56;

将它们组合为 32 位或更大应该是微不足道的。