使用SIMD指令的平行二项式系数

Question

背景

我最近一直在使用一些旧代码（〜1998）并重写其中一些代码以提高性能。以前，在状态的基本数据结构中，我将元素存储在多个数组中，现在，我使用的是原始位（对于需要少于64位的情况）。也就是说，在我拥有b个元素的数组之前，现在我已经在单个64位整数中设置了b位，这些位指示该值是否是我的状态的一部分。

使用_pext_u64和_pdep_u64之类的内在函数，我设法使所有操作的速度提高了5-10倍。我正在做最后一个操作，该操作与计算完美的哈希函数有关。

散列函数的确切细节并不是太重要，但是归结为计算二项式系数（对于各种n choose k和n!/((n-k)!k!)，都是n-k。我当前的代码使用了一个较大的查找表这可能很难单独显着加快速度（除了我尚未测量的表中可能发生的高速缓存未命中）。

但是，我在想，使用SIMD指令，我可以直接并行地针对多个状态计算这些指令，从而提高整体性能。

一些约束：

• 在每个64位状态中总是准确地设置了b位（代表小数字）。
• 二项式系数中的k值与b相关，并且在计算中均匀变化。这些值很小（大多数时间<= 5）。
• 最终的哈希将为<1500万（很容易容纳32位）。

因此，我可以很容易地写出并行执行此数学运算，以及将所有运算保持为整数倍/除法而没有余数，同时保持在32位以内的方式。总体流程为：

1. 将这些位提取为适合SIMD指令的值。
2. 以一种避免溢出的方式执行n choose k计算。
3. 从每个条目中提取最终的哈希值

但是，我之前没有写过SIMD代码，因此我仍在加快所有可用功能以及它们的警告/效率。

示例：

以前，我总是将数据放在一个数组中，假设总是有5个元素：

[3 7 19 31 38]

现在我为此使用单个64位值：

0x880080088

这使许多其他操作非常有效。对于完美的哈希，我需要有效地计算出这样的值（使用c作为选择）：

(50c5)-(38c5) + (37c4)-(31c4) + (30c3)-(19c3) + ...

但是，实际上，我需要计算很多，只是值略有不同：

(50c5)-(Xc5) + ((X-1)c4)-(Yc4) + ((Y-1)c3)-(Zc3) + ...

所有X / Y / Z ...将有所不同，但是每个的计算形式都相同。

问题：

1. 通过转换为SIMD操作获得效率的直觉是否合理？ （Some sources suggest "no"，但这是计算单个系数而不是并行执行多个系数的问题。）

2. 是否有比重复的_tzcnt_u64调用更有效的方法，用于将位提取到数据结构中以进行SIMD操作？ （例如，如果有帮助，我可以将64位状态表示形式临时分成32位块，但这样就不能保证每个元素中设置的位数相同。）

3. 当我知道不会发生溢出时，计算二项式系数的几个顺序乘/除运算的最佳内在函数是什么。 （当我浏览英特尔参考资料时，在浏览所有变体时都难以快速解释其命名-尚不清楚我想要的是可用的。）

4. 如果直接计算系数不太可能有效，是否可以将SIMD指令用于并行查找我之前的系数查找表？

（我为将几个问题放在一起而表示歉意，但鉴于具体情况，我认为将它们作为一个整体放在一起会更好。）

Answer 1

这里是一种可能的解决方案，它一次使用一种状态从查找表进行计算。在多个状态下并行执行此操作（而不是使用单个状态）可能会更有效率。

int64_t GetPerfectHash2(State &s)
{
    // 6 values will be used
    __m256i offsetsm1 = _mm256_setr_epi32(6*boardSize-1,5*boardSize-1,
                                          4*boardSize-1,3*boardSize-1,
                                          2*boardSize-1,1*boardSize-1,0,0);
    __m256i offsetsm2 = _mm256_setr_epi32(6*boardSize-2,5*boardSize-2,
                                          4*boardSize-2,3*boardSize-2,
                                          2*boardSize-2,1*boardSize-2,0,0);
    int32_t index[9];
    uint64_t value = _pext_u64(s.index2, ~s.index1);
    index[0] = boardSize-numItemsSet+1;
    for (int x = 1; x < 7; x++)
    {
        index[x] = boardSize-numItemsSet-_tzcnt_u64(value);
        value = _blsr_u64(value);
    }
    index[8] = index[7] = 0;

    // Load values and get index in table
    __m256i firstLookup = _mm256_add_epi32(_mm256_loadu_si256((const __m256i*)&index[0]), offsetsm2);
    __m256i secondLookup = _mm256_add_epi32(_mm256_loadu_si256((const __m256i*)&index[1]), offsetsm1);
    // Lookup in table
    __m256i values1 = _mm256_i32gather_epi32(combinations, firstLookup, 4);
    __m256i values2 = _mm256_i32gather_epi32(combinations, secondLookup, 4);
    // Subtract the terms
    __m256i finalValues = _mm256_sub_epi32(values1, values2);
    _mm256_storeu_si256((__m256i*)index, finalValues);

    // Extract out final sum
    int64_t result = 0;
    for (int x = 0; x < 6; x++)
    {
        result += index[x];
    }
    return result;  
}

请注意，我实际上有两个类似的案例。在第一种情况下，我不需要_pext_u64，并且此代码比我现有的代码慢约3倍。在第二种情况下，我需要它，它快了25％。