基于RISC-V CLMUL指令优化3D Morton编码的性能问询

本研究针对RISC-V平台展开探索，该平台可选择性支持无进位乘法指令CLMUL。此类指令也适用于其他提供通用寄存器无进位乘法的架构，但目前暂未找到其他相关架构的应用案例。

3D Morton编码的基本需求 #

3D Morton编码需要将三个整数x、y、z的位进行交错排列，满足规则：r₃ᵢ₊₀ = xᵢ，r₃ᵢ₊₁ = yᵢ，r₃ᵢ₊₂ = zᵢ。如果处理器架构支持位存入指令（比如支持BMI2扩展的x86-64平台），实现会非常简洁：

/* 针对范围[0, 1023]的x、y、z，通过位交错计算30位3D Morton编码。
   其中i∈[0,9]：r<3*i+0> = x<i>, r<3*i+1> = y<i>, r<3*i+2> = z<i>。
*/
uint32_t encode_morton3d_10_bits (uint32_t x, uint32_t y, uint32_t z)
{
    return (_pdep_u32 (x, 01111111111) |  // 注意是八进制常量！
            _pdep_u32 (y, 02222222222) |
            _pdep_u32 (z, 04444444444));
}

无位存入指令架构的模拟实现（以RISC-V为例） #

对于不支持位存入操作的架构（如RISC-V），需要模拟将源位存入结果中每第三位的功能。相比常用方法，下述基于乘法扩展位的算法缩短了依赖链长度，提升了指令级并行度：低端CPU上性能不逊于常用方法，高端CPU上性能提升明显。所用乘法因子足够简单，编译器可针对特定微架构替换为更高效的指令：

/* 将输入x（范围[0,1023]）的最低10位分散到结果的每第三位。
   其中i∈[0,9]：r<3*i+0> = x<i>, r<3*i+1> = 0, r<3*i+2> = 0。
*/
uint32_t deposit_every_third (uint32_t x)
{
    const uint32_t sm1  = 0x333ul;      // 源掩码：选中位0,1,4,5,8,9
    const uint32_t sm2  = 0x0ccul;      // 源掩码：选中位2,3,6,7
    const uint32_t ml1a = 0x50005ul;    // 乘法因子：将位分散到0-11,16-27位
    const uint32_t ml1b = 0x00500ul;    // 乘法因子：将位分散到8-19位
    const uint32_t ml2  = 0x05050ul;    // 乘法因子：将位分散到6-13,14-21位
    const uint32_t rm1  = 0x09009009ul; // 结果掩码：保留位0,3,12,15,24,27
    const uint32_t rm2  = 0x00240240ul; // 结果掩码：保留位6,9,18,21
    uint32_t t = x & sm1;
    return ((((t * ml1a) | (t * ml1b)) & rm1) | (((x & sm2) * ml2) & rm2));
}

为计算3D Morton编码，可将移位操作整合到位存入功能中，示例实现如下：

uint32_t deposit_every_third_shifted (uint32_t x, uint32_t s)
{
    const uint32_t sm1  = 0x333ul;           // 源掩码：选中位0,1,4,5,8,9
    const uint32_t sm2  = 0x0ccul;           // 源掩码：选中位2,3,6,7
    const uint32_t ml1a = 0x50005ul << s;    // 乘法因子（带移位）
    const uint32_t ml1b = 0x00500ul << s;
    const uint32_t ml2  = 0x05050ul << s;
    const uint32_t rm1  = 0x09009009ul << s; // 结果掩码（带移位）
    const uint32_t rm2  = 0x00240240ul << s;
    uint32_t t = x & sm1;
    return ((((t * ml1a) | (t * ml1b)) & rm1) | (((x & sm2) * ml2) & rm2));
}

/* 针对范围[0, 1023]的x、y、z，通过位交错计算30位3D Morton编码。
   其中i∈[0,9]：r<3*i+0> = x<i>, r<3*i+1> = y<i>, r<3*i+2> = z<i>。
*/
uint32_t encode_morton3d_10_bits (uint32_t x, uint32_t y, uint32_t z)
{
    return (deposit_every_third_shifted (x, 0) |
            deposit_every_third_shifted (y, 1) |
            deposit_every_third_shifted (z, 2));
}

在32位RISC-V平台上以-O3级别编译上述代码，会生成高效的指令序列：clang 20.1与gcc 15.1均生成58条线性指令，但指令序列差异较大。值得注意的是两款编译器对乘法的优化不同：源码中有9次乘法，clang生成7条MUL指令，而gcc仅生成5条。

利用RISC-V Zbc扩展的CLMUL指令优化 #

在上述deposit_every_third_shifted()的实现中，使用普通整数乘法时需要将第一个乘法因子拆分为ml1a和ml1b两部分，避免进位传播导致结果错误。而支持标准化ISA扩展Zbc（提供无进位乘法指令CLMUL）的RISC-V处理器可以轻松解决这个问题。目前未找到对应的编译器内置函数，因此采用内联汇编实现：

#if defined (PORTABLE)
uint32_t clmul (uint32_t rs1, uint32_t rs2)
{
    uint32_t x = 0;
    for (int i = 0; i < 32; i++)
        if ((rs2 >> i) & 1)
            x ^= rs1 << i;
    return x;
}
#else // !PORTABLE
long clmul (long a, long b) // 注意：假设long的大小与RISC-V寄存器大小匹配
{
  long r;
  __asm__("clmul\t%0,%1,%2" : "=r"(r) : "r"(a), "r"(b));
  return r;
}
#endif // PORTABLE

uint32_t deposit_every_third_clmul_shifted (uint32_t x, uint32_t s)
{
    const uint32_t sm1 = 0x333ul;           // 源掩码：选中位0,1,4,5,8,9
    const uint32_t sm2 = 0x0ccul;           // 源掩码：选中位2,3,6,7
    const uint32_t ml1 = 0x50505ul << s;    // 乘法因子（带移位）
    const uint32_t ml2 = 0x05050ul << s;    
    const uint32_t rm1 = 0x09009009ul << s; // 结果掩码（带移位）
    const uint32_t rm2 = 0x00240240ul << s; 
    return (clmul (x & sm1, ml1) & rm1) | (clmul (x & sm2, ml2) & rm2);
}

uint32_t encode_morton3d_10_bits_clmul (uint32_t x, uint32_t y, uint32_t z)
{
    return (deposit_every_third_clmul_shifted (x, 0) |
            deposit_every_third_clmul_shifted (y, 1) |
            deposit_every_third_clmul_shifted (z, 2));
}

以-O3级别编译时，clang 20.1与gcc 15.1均生成47条指令，指令序列差异显著。由于CLMUL的执行速度不慢于MUL，甚至更快（例如CLMUL延迟为2周期，MUL为4周期），指令数量的减少可显著提升Morton编码的生成速度。

但如果Zbc扩展无其他重要用途，仅靠这一性能提升可能不足以证明在RISC-V处理器中集成该扩展的合理性。

核心疑问 #

在计算3D Morton编码时，是否存在进一步利用CLMUL指令获取更高性能收益的方法？已知clmul(x,x)可在连续源位之间插入一个零位，适用于2D Morton编码，但尚未找到将此思路扩展到3D Morton编码所需的、在连续源位之间插入两个零位的方法。

内容的提问来源于stack exchange，提问作者njuffa