2D morton 代码编码/解码 64 位

Posted 2023-02-19

【中文标题】2D morton 代码编码/解码 64 位

【英文标题】：2D morton code encode/decode 64bits

【发布时间】：2015-08-12 21:31:18

【问题描述】：

如何将给定 [x, y] 的 morton 代码（z 顺序）编码/解码为 32 位无符号整数，生成 64 位 morton 代码，反之亦然？ 我确实有 xy2d 和 d2xy 但仅适用于 16 位宽的坐标，产生 32 位莫顿数。在网上搜了很多，但没有找到。请帮忙。

【问题讨论】：

将32位版本扩展到64位真的不难。将所有蒙版的宽度加倍，并按照与其他蒙版相同的模式添加一个额外的步骤。

【参考方案1】：

如果您可以使用特定于体系结构的指令，那么您可能能够加速操作，超出使用 bit-twiddeling hack 所能达到的速度：

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

这是一个完整的例子（用 GCC 测试）：

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


void morton_to_xy(uint64_t m, uint32_t *x, uint32_t *y)

  *x = _pext_u64(m, 0x5555555555555555);
  *y = _pext_u64(m, 0xaaaaaaaaaaaaaaaa);

如果您必须支持较早的 CPU 或 ARM 平台，则不会丢失所有内容。您至少可以从特定于密码学的说明中获得有关 xy_to_morton 函数的帮助。

如今，许多 CPU 都支持无进位乘法。在 ARM 上，来自 NEON 指令集的 vmul_p8。在 X86 上，您会在 CLMUL 指令集（自 2010 年起提供）中找到 PCLMULQDQ。

这里的技巧是，一个数字与其自身的无进位乘法将返回一个位模式，该模式包含参数的原始位以及交错的零位。所以它与上面显示的 _pdep_u32(x,0x55555555) 相同。例如。它变成以下字节：

 +----+----+----+----+----+----+----+----+
 | b7 | b6 | b5 | b4 | b3 | b2 | b1 | b0 |
 +----+----+----+----+----+----+----+----+

进入：

 +----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+
 | 0  | b7 | 0  | b6 | 0  | b5 | 0  | b4 | 0  | b3 | 0  | b2 | 0  | b1 | 0  | b0 |
 +----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+

现在您可以将 xy_to_morton 函数构建为（此处显示为 CLMUL 指令集）：

#include <wmmintrin.h>
#include <stdint.h>

// on GCC, compile with option -mpclmul

uint64_t carryless_square (uint32_t x)

  uint64_t val[2] = x, 0;
  __m128i *a = (__m128i * )val;
  *a = _mm_clmulepi64_si128 (*a,*a,0);
  return val[0];


uint64_t xy_to_morton (uint32_t x, uint32_t y)

  return carryless_square(x)|(carryless_square(y) <<1);

_mm_clmulepi64_si128 生成一个 128 位的结果，我们只使用低 64 位。因此，您甚至可以改进上述版本并使用单个 _mm_clmulepi64_si128 来完成这项工作。

这是您在主流平台（例如，带有 NEON 和 x86 的现代 ARM）上所能获得的一样好。不幸的是，我不知道使用加密指令加速 morton_to_xy 函数的任何技巧，我努力了几个月。

【讨论】：

真的很棒。满足。

@DawidSzymański 如果您想了解更多信息，我建议您查看此博客：bitmath.blogspot.de 并阅读有关 tesseral 算术的内容（即使用以 morton 顺序存储的数字进行算术运算，而无需对其进行编码/解码）。我很确定你可以用它来填充空间曲线。

@harold，有趣的事实是：我们已经享受了 GF(2'm) 中 x*x 运算的位旋转幂的数学奇点。然而，加密人也喜欢在 GF(2'm) 中有一个快速的 sqrt(x)。他们已经发现这是关于将偶数位与奇数位分开，但他们还不知道位旋转的技巧。我认为每个人都可以从中学习！

@NilsPipenbrinck 在这么久之后才找到这个答案，好奇它们是否存在于 3D 空间中？说将 x,y,z 编码为 Z 顺序，反之亦然。

【参考方案2】：

void xy2d_morton(uint64_t x, uint64_t y, uint64_t *d)

    x = (x | (x << 16)) & 0x0000FFFF0000FFFF;
    x = (x | (x << 8)) & 0x00FF00FF00FF00FF;
    x = (x | (x << 4)) & 0x0F0F0F0F0F0F0F0F;
    x = (x | (x << 2)) & 0x3333333333333333;
    x = (x | (x << 1)) & 0x5555555555555555;

    y = (y | (y << 16)) & 0x0000FFFF0000FFFF;
    y = (y | (y << 8)) & 0x00FF00FF00FF00FF;
    y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0F;
    y = (y | (y << 2)) & 0x3333333333333333;
    y = (y | (y << 1)) & 0x5555555555555555;

    *d = x | (y << 1);


// morton_1 - extract even bits

uint32_t morton_1(uint64_t x)

    x = x & 0x5555555555555555;
    x = (x | (x >> 1))  & 0x3333333333333333;
    x = (x | (x >> 2))  & 0x0F0F0F0F0F0F0F0F;
    x = (x | (x >> 4))  & 0x00FF00FF00FF00FF;
    x = (x | (x >> 8))  & 0x0000FFFF0000FFFF;
    x = (x | (x >> 16)) & 0x00000000FFFFFFFF;
    return (uint32_t)x;


void d2xy_morton(uint64_t d, uint64_t &x, uint64_t &y)

    x = morton_1(d);
    y = morton_1(d >> 1);

【讨论】：

在morton_1 中，最后一个值不应该是0x00000000FFFFFFFF吗？

附注morton_1 可以返回 uint32_t。

【参考方案3】：

无论位数如何，原始代码都是相同的。如果你不需要超快的位旋转版本，这个就可以了

uint32_t x;
uint32_t y;
uint64_t z = 0;

for (int i = 0; i < sizeof(x) * 8; i++)

  z |= (x & (uint64_t)1 << i) << i | (y & (uint64_t)1 << i) << (i + 1);

如果您需要更快的位旋转，那么这个应该可以工作。请注意，x 和 y 必须是 64 位变量。

uint64_t x;
uint64_t y;
uint64_t z = 0;

x = (x | (x << 16)) & 0x0000FFFF0000FFFF;
x = (x | (x << 8)) & 0x00FF00FF00FF00FF;
x = (x | (x << 4)) & 0x0F0F0F0F0F0F0F0F;
x = (x | (x << 2)) & 0x3333333333333333;
x = (x | (x << 1)) & 0x5555555555555555;

y = (y | (y << 16)) & 0x0000FFFF0000FFFF;
y = (y | (y << 8)) & 0x00FF00FF00FF00FF;
y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0F;
y = (y | (y << 2)) & 0x3333333333333333;
y = (y | (y << 1)) & 0x5555555555555555;

z = x | (y << 1);

【讨论】：

对快速方式和反向更感兴趣？