4 uint16_t的快速模12算法打包在uint64_t中

Question

考虑以下联合:

union Uint16Vect {
    uint16_t _comps[4];
    uint64_t _all;
};

是否有快速算法来确定每个组件是否等于1模12？

一个天真的代码序列是:

Uint16Vect F(const Uint16Vect a) {
    Uint16Vect r;
    for (int8_t k = 0; k < 4; k++) {
        r._comps[k] = (a._comps[k] % 12 == 1) ? 1 : 0;
    }
    return r;
}

Answer 1

这样的常数除法应该乘以乘法逆.如您所见,编译器优化x/12到x*43691 >> 19

bool h(uint16_t x)
{
    return x % 12 == 1;
}
h(unsigned short):
        movzx   eax, di
        imul    eax, eax, 43691 ; = 0xFFFF*8/12 + 1
        shr     eax, 19
        lea     eax, [rax+rax*2]
        sal     eax, 2
        sub     edi, eax
        cmp     di, 1
        sete    al
        ret

因为SSE ​​/ AVX中有乘法指令,所以可以很容易地进行矢量化.此外,x = (x % 12 == 1) ? 1 : 0;可以简化为x = (x % 12 == 1)然后转换为x = (x - 1) % 12 == 0避免值1的常量表进行比较.您可以使用向量扩展名,以便gcc自动为您生成代码

typedef uint16_t ymm __attribute__((vector_size(32)));
ymm mod12(ymm x)
{
    return !!((x - 1) % 12);
}

以下是gcc的输出

mod12(unsigned short __vector(16)):
        vpcmpeqd    ymm3, ymm3, ymm3  ; ymm3 = -1
        vpaddw      ymm0, ymm0, ymm3
        vpmulhuw    ymm1, ymm0, YMMWORD PTR .LC0[rip] ; multiply with 43691
        vpsrlw      ymm2, ymm1, 3
        vpsllw      ymm1, ymm2, 1
        vpaddw      ymm1, ymm1, ymm2
        vpsllw      ymm1, ymm1, 2
        vpcmpeqw    ymm0, ymm0, ymm1
        vpandn      ymm0, ymm0, ymm3
        ret

Clang和ICC不支持!!矢量类型,因此您需要更改为(x - 1) % 12 == 0.不幸的是,似乎编译器不支持__attribute__((vector_size(8))发出MMX指令.但是现在你应该使用SSE或AVX

x % 12 == 1你可以在上面的相同Godbolt链接中看到输出更短,但你需要一个包含1s的表来进行比较,这可能更好或没有.在您的情况下检查哪一个更快

或者对于像这样的小输入范围,您可以使用查找表.基本版本需要65536个元素的数组

#define S1(x) ((x) + 0) % 12 == 1, ((x) + 1) % 12 == 1, ((x) + 2) % 12 == 1, ((x) + 3) % 12 == 1, \
              ((x) + 4) % 12 == 1, ((x) + 4) % 12 == 1, ((x) + 6) % 12 == 1, ((x) + 7) % 12 == 1
#define S2(x) S1((x + 0)*8), S1((x + 1)*8), S1((x + 2)*8), S1((x + 3)*8), \
              S1((x + 4)*8), S1((x + 4)*8), S1((x + 6)*8), S1((x + 7)*8)
#define S3(x) S2((x + 0)*8), S2((x + 1)*8), S2((x + 2)*8), S2((x + 3)*8), \
              S2((x + 4)*8), S2((x + 4)*8), S2((x + 6)*8), S2((x + 7)*8)
#define S4(x) S3((x + 0)*8), S3((x + 1)*8), S3((x + 2)*8), S3((x + 3)*8), \
              S3((x + 4)*8), S3((x + 4)*8), S3((x + 6)*8), S3((x + 7)*8)

bool mod12e1[65536] = {
    S4(0U), S4(8U), S4(16U), S4(24U), S4(32U), S4(40U), S4(48U), S4(56U)
}

在使用时只需更换x % 12 == 1用mod12e1[x].这当然可以是矢量化的

但由于结果只有1或0,因此您还可以使用65536位数组将大小减小到仅8KB

您还可以通过4和3将可分性检验为12除以4的可分性显然是微不足道的.可以通过多种方式计算3的可分性

• 一个是计算差异之间的奇数位数字之和与偶数数字的和像גלעדברקן的回答,并检查它是否是由3或不能整除

• 或者您可以检查基数2 ^2k(如基数4,16,64 ......)中的数字之和是否可以被3整除.

  这是有效的,因为在b检查任何除数n的可除性的基础上b - 1,只需检查数字的总和是否可被n整除.这是它的一个实现

  void modulo12equals1(uint16_t d[], uint32_t size) {
      for (uint32_t i = 0; i < size; i++)
      {
          uint16_t x = d[i] - 1;
          bool divisibleBy4 = x % 4 == 0;
          x = (x >> 8) + (x & 0x00ff); // max 1FE
          x = (x >> 4) + (x & 0x000f); // max 2D
          bool divisibleBy3 = !!((01111111111111111111111ULL >> x) & 1);
          d[i] = divisibleBy3 && divisibleBy4;
      }
  }
  

Roland Illig将 3 分为 3分的可信度

由于自动向量化程序集输出太长,您可以在Godbolt链接上进行检查

也可以看看

• 如何知道二进制数是否除以3？
• 确定二进制数是否可以被3整除
• 比特表示和可分性3
• 建筑电路的可分性3
• 检查一个数字是否可以被3整除
• 检查号码的逻辑可以被3整除吗？

Answer 2

请注意x mod 12 == 1意味着x mod 4 == 1，而后者非常便宜。

所以：

is_mod_12 = ((input & 3) == 1) && ((input % 12) == 1);

input mod 4在经常不 的情况下1，将节省您大量的模运算。然而，如果input mod 4通常是这样1，这将使性能稍微变差。