c - 将 4 个 uint16_t 封装在一个 uint64_t 中的快速模 12 算法

Question

考虑以下 union :

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

编译器会将除以常数优化为 multiplication by the reciprocal or multiplicative inverse 。例如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进行比较。您可以使用 vector extension 以便 gcc 自动为您生成代码

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

下面是output from 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 不支持 vector 类型上的 !!，因此您需要更改为 (x - 1) % 12 == 0。不幸的是，编译器似乎不支持 __attribute__((vector_size(8)) 发出 MMX 指令。但现在无论如何你都应该使用 SSE 或 AVX

x % 12 == 1 的输出较短，如您在上面的同一个 Godbolt 链接中所见，但您需要一个包含 1 的表进行比较，这可能更好也可能不好。编译器也可能无法像手写代码那样完全优化，因此您可以尝试使用内部函数手动对代码进行矢量化。检查哪个在您的情况下工作得更快

更好的方法是 ((x * 43691) & 0x7ffff) < 43691 或 nwellnhof's answer 中提到的 x * 357913942 < 357913942，它也应该易于矢量化

---

或者，对于像这样的小输入范围，您可以使用查找表。基础版需要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的整除率可以通过多种方式计算

• 一个正在计算奇数和与偶数和之间的差，例如 גלעד ברקן's answer 和检查它是否可以被 3 整除

• 或者您可以检查以 2 为基数的数字总和^2k(例如以 4、16、64 为基数...)是否正确是否能被 3 整除。

  这是有效的，因为在基数 b 中检查 b - 1 的任何除数 n 的整除性，只需检查数字的总和是否可以被 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;
        }
    }
  

被 3 整除到 Roland Illig 的积分

由于自动向量化汇编输出太长，可以在Godbolt link上查看

另见

• How to know if a binary number divides by 3?
• Determine whether or not a binary number is divisible by 3
• Bit representation and divisibility by 3
• building circuit for divisibility by 3
• Check if a number is divisible by 3
• Logic to check the number is divisible by 3 or not?