使用固定运行时除数进行快速整数除法和取模。

Question

int n_attrs = some_input_from_other_function() // [2..5000]
vector<int> corr_indexes; // size = n_attrs * n_attrs
vector<char> selected; // szie = n_attrs
vector<pair<int,int>> selectedPairs; // size = n_attrs / 2
// vector::reserve everything here
...
// optimize the code below
const int npairs = n_attrs * n_attrs;
selectedPairs.clear();
for (int i = 0; i < npairs; i++) {
    const int x = corr_indexes[i] / n_attrs;
    const int y = corr_indexes[i] % n_attrs;
    if (selected[x] || selected[y]) continue; // fit inside L1 cache
    
    // below lines are called max 2500 times, so they're insignificant
    selected[x] = true;
    selected[y] = true;
    selectedPairs.emplace_back(x, y);
    if (selectedPairs.size() == n_attrs / 2) break;
}

我有一个像这样的函数。瓶颈在于

    const int x = corr_indexes[i] / n_attrs;
    const int y = corr_indexes[i] % n_attrs;

n_attrs 在循环期间是 const 的，因此我希望找到一种加速此循环的方法。corr_indexes[i]，n_attrs > 0，< max_int32。编辑：请注意，n_attrs 不是编译时常量。

如何优化此循环？不允许使用额外的库。此外，是否有办法使此循环并行化（CPU 或 GPU 都可以，所有内容在此循环之前已经存储在 GPU 内存中）。

- Huy Le

---

评论不适合进行长时间的讨论；此对话已被移至聊天室。 - Samuel Liew

类似相关的：如果你只需要将一个值（如散列）映射到一个范围内，但是它不需要有序（因此实际上你不需要模运算，只需要任何一对多半均匀映射），请参见有没有一种方法将整数值包装到整数范围[min，max]中，而不使用除法或模运算？以获取更快的技巧，只需使用扩展值的高半部分value*(uint64_t)range。因此，value必须在完整的32位范围内相对均匀地分布。 - Peter Cordes

@PeterCordes 哦，我的最终解决方案不使用任何 div/mod。我不是将 (i,j) -> i * N + j 映射，而是将其映射为 (i,j) -> i * 2^K + j，其中 2^K >= N。然后就是通过 2 的幂进行除法/模运算，这比任何整数除法算法都要快得多。因此，我的问题的解决方案与这个问题完全不同 :D - Huy Le

Answer 1

我将限制我的评论在整数除法上，因为从一级角度来看，C++中的模运算可以被视为整数除法加上反向乘法和减法实现，尽管在某些情况下，有更便宜的直接计算模数的方法，例如计算模2ⁿ时。

在大多数平台上，基于软件仿真或迭代硬件实现，整数除法非常缓慢。但是去年根据对苹果M1的微基准测试广泛报道称，它具有极快的整数除法，可能是通过使用专用电路实现的。

自从Torbjörn Granlund和Peter Montgomery在近30年前发表了一篇开创性论文以来，人们就已经广泛知道如何使用整数乘法加上可能的移位和/或其他校正步骤来替换具有恒定除数的整数除法。这个算法通常被称为魔术乘法器技术。它需要预先计算一些与整数除数相关的参数，以用于基于乘法的仿真序列中。

Torbjörn Granlund和Peter L. Montgomery，“Division by invariant integers using multiplication”，ACM SIGPLAN Notices，Vol. 29，June 1994，pp. 61-72 (online)。

目前，所有主要的工具链在处理整数除数为编译时常量时都采用Granlund-Montgomery算法的变体。预计算发生在编译器内部的编译时，然后编译器使用计算出的参数生成代码。一些工具链也可能使用这个算法来处理重复使用的运行时常量除数的除法。对于循环中的运行时常量除数，这可能涉及在循环之前发出一个预计算块来计算必要的参数，然后在循环内使用这些参数进行除法模拟代码。

如果一个工具链不能优化具有运行时常量除数的除法，则可以手动使用下面代码所示的相同方法。然而，这不太可能达到与基于编译器的解决方案相同的效率，因为在C++级别上不能以便携方式高效地表达所需仿真序列中使用的所有机器操作。这尤其适用于算术右移和带进位的加法。

下面的代码演示了通过乘法进行参数预计算和整数除法仿真的原理。通过投入更多的时间进行设计，很可能可以找到更有效的参数预计算和仿真实现。

#include <cstdio>
#include <cstdlib>
#include <cstdint>

#define PORTABLE  (1)

uint32_t ilog2 (uint32_t i)
{
    uint32_t t = 0;
    i = i >> 1;
    while (i) {
        i = i >> 1;
        t++;
    }
    return (t);
}

/* Based on: Granlund, T.; Montgomery, P.L.: "Division by Invariant Integers 
   using Multiplication". SIGPLAN Notices, Vol. 29, June 1994, pp. 61-72
*/
void prepare_magic (int32_t divisor, int32_t &multiplier, int32_t &add_mask, int32_t &sign_shift)
{
    uint32_t divisoru, d, n, i, j, two_to_31 = uint32_t (1) << 31;
    uint64_t m_lower, m_upper, k, msb, two_to_32 = uint64_t (1) << 32;

    divisoru = uint32_t (divisor);
    d = (divisor < 0) ? (0 - divisoru) : divisoru;
    i = ilog2 (d);
    j = two_to_31 % d;
    msb = two_to_32 << i;
    k = msb / (two_to_31 - j);
    m_lower = msb / d;
    m_upper = (msb + k) / d;
    n = ilog2 (uint32_t (m_lower ^ m_upper));
    n = (n > i) ? i : n;
    m_upper = m_upper >> n;
    i = i - n;
    multiplier = int32_t (uint32_t (m_upper));
    add_mask = (m_upper >> 31) ? (-1) : 0;
    sign_shift = int32_t ((divisoru & two_to_31) | i);
}

int32_t arithmetic_right_shift (int32_t a, int32_t s)
{
    uint32_t msb = uint32_t (1) << 31;
    uint32_t ua = uint32_t (a);
    ua = ua >> s;
    msb = msb >> s;
    return int32_t ((ua ^ msb) - msb);
}

int32_t magic_division (int32_t dividend, int32_t multiplier, int32_t add_mask, int32_t sign_shift)
{
    int64_t prod = int64_t (dividend) * multiplier;
    int32_t quot = (int32_t)(uint64_t (prod) >> 32);
    quot = int32_t (uint32_t (quot) + (uint32_t (dividend) & uint32_t (add_mask)));
#if PORTABLE
    const int32_t byte_mask = 0xff;
    quot = arithmetic_right_shift (quot, sign_shift & byte_mask);
#else // PORTABLE
    quot = quot >> sign_shift; // must mask shift count & use arithmetic right shift
#endif // PORTABLE
    quot = int32_t (uint32_t (quot) + (uint32_t (dividend) >> 31));
    if (sign_shift < 0) quot = -quot;
    return quot;
}

int main (void)
{
    int32_t multiplier;
    int32_t add_mask;
    int32_t sign_shift;
    int32_t divisor;
    
    for (divisor = -20; divisor <= 20; divisor++) {
        /* avoid division by zero */
        if (divisor == 0) {
            divisor++;
            continue;
        }
        printf ("divisor=%d\n", divisor);
        prepare_magic (divisor, multiplier, add_mask, sign_shift);
        printf ("multiplier=%d add_mask=%d sign_shift=%d\n", 
                multiplier, add_mask, sign_shift);
        printf ("exhaustive test of dividends ... ");
        uint32_t dividendu = 0;
        do {
            int32_t dividend = (int32_t)dividendu;
            /* avoid overflow in signed integer division */
            if ((divisor == (-1)) && (dividend == ((-2147483647)-1))) {
                dividendu++;
                continue;
            }
            int32_t res = magic_division (dividend, multiplier, add_mask, sign_shift);
            int32_t ref = dividend / divisor;
            if (res != ref) {
                printf ("\nERR dividend=%d (%08x) divisor=%d  res=%d  ref=%d\n",
                        dividend, (uint32_t)dividend, divisor, res, ref);
                return EXIT_FAILURE;
            }
            dividendu++;
        } while (dividendu);
        printf ("PASSED\n");
    }
    return EXIT_SUCCESS;
}

- njuffa

---

谢谢，这就是我要寻找的算法和关键字。在有人能够提供更快的实现之前，我会将其标记为已解决。 - Huy Le

Answer 2

如何优化这个循环？

这是使用libdivide的完美案例。该库旨在通过使用编译器在编译时使用的策略，在运行时加速常数除法。该库是仅头文件，因此不会创建任何运行时依赖关系。它还支持除法的向量化(即使用SIMD指令)，这绝对是在这种情况下需要使用的，可以大大加快计算速度，而编译器无法在不显着更改循环的情况下完成(最终由于运行时定义的除数，效率不如使用该库)。请注意，libdivide的许可证非常宽松(zlib)，因此您可以轻松地将其包含在项目中，没有强制性的限制(如果您更改了它，您基本上只需要标记它已修改)。

如果仅头文件库不适用，则需要重新实现轮子。想法是将常数除法转换为移位和乘法序列。@njuffa的非常好的答案说明了如何做到这一点。您还可以查看高度优化的libdivide代码。

对于小正除数和小正被除数，没有必要进行长序列的操作。您可以使用基本序列进行欺骗。

uint64_t dividend = corr_indexes[i]; // Must not be too big
uint64_t divider = n_attrs;
uint64_t magic_factor = 4294967296 / n_attrs + 1; // Must be precomputed once
uint32_t result = (dividend * magic_factor) >> 32;

这种方法适用于uint16_t的被除数/除数，但对于更大的值则不安全。实际上，如果dividend的值超过了~800,000，则会失败。更大的被除数需要更复杂的序列，通常也更慢。

有什么方法可以并行化此循环吗？

只有除法/取模可以安全地并行化。在其余循环中存在“循环依赖”，阻止任何并行化（除非做出额外的假设）。因此，循环可以分成两个部分：一个计算除法并将uint16_t结果放入临时数组的部分，稍后以串行方式计算。数组不需要太大，否则计算会受到内存限制，并且产生的并行代码可能比当前的代码更慢。因此，您需要操作适合至少L3高速缓存的小块。如果块太小，则线程同步也可能是问题。最好的解决方案肯定是使用滚动窗口的块。所有这些都肯定有点繁琐/棘手。

请注意，可以使用SIMD指令进行除法部分的计算（使用libdivide很容易）。您还需要分割循环并使用块，但是由于没有同步开销，因此块不需要很大。类似64个整数的东西就足够了。

---

请注意，最近的处理器可以高效地计算此类除法，尤其是针对32位整数（64位的则往往更昂贵）。这是Alder lake、Zen3和M1处理器（P-cores）的情况。请注意，x86/x86-64处理器可以在一个指令中计算余数和商。还要注意，虽然除法具有相当大的延迟，但许多处理器可以将多个除法流水线化，以获得合理的吞吐量。例如，32位div指令在Skylake上的延迟为23~28个周期，但倒数吞吐量为4~6。Zen1/Zen2上似乎不是这种情况。

- Jérôme Richard

---

谢谢！我会考虑使用 libdivide，因为它只有一个头文件。 - Huy Le

“操作小块”: 是的，这也是一个非常好的观点。通过进行 5000^2 次迭代并仅使用 2500 个 emplace_back，可以在临时数组上使用 SIMD 分割，然后检查 selected[x]，selected[y]。谢谢。 - Huy Le

Answer 3

我会通过以下方式优化// optimize the code below后面的部分：

• 使用n_attrs
• 生成如下的函数字符串：

void dynamicFunction(MyType & selectedPairs, Foo & selected)
{
    const int npairs = @@ * @@;
    selectedPairs.clear();
    for (int i = 0; i < npairs; i++) {
        const int x = corr_indexes[i] / @@;
        const int y = corr_indexes[i] % @@;
        if (selected[x] || selected[y]) continue; // fit inside L1 cache
    
        // below lines are called max 2500 times, so they're insignificant
        selected[x] = true;
        selected[y] = true;
        selectedPairs.emplace_back(x, y);
        if (selectedPairs.size() == @@ / 2) 
            break;
    }
}

• 用n_attrs的值替换所有的 @@
• 编译并生成一个DLL
• 链接并调用函数

这样，对于DLL来说，n_attrs是一个编译时常量，并且编译器可以自动对该值进行大部分优化，例如：

• 当x是2的幂时，对n&(x-1)进行操作而不是n%x
• 使用移位和乘法代替除法
• 也许还有其他优化，例如对x和y进行预计算的展开循环（因为已知x）

在紧密循环中进行的一些整数运算更容易被编译器SIMD化/向量化，当大部分部分在编译时已知。

如果你的CPU是AMD，甚至可以尝试使用神奇的浮点运算代替未知/未知除法以获得向量化。

通过缓存所有（或大百分比的）n_attrs值，您可以消除以下延迟：

• 字符串生成
• 编译
• 文件（DLL）读取（假设某个面向对象的DLL包装）

---

如果要优化的部分将在GPU中运行，则很有可能CUDA/OpenCL实现已经以浮点数的方式执行整数除法（以保持SIMD路径被占用而不是在整数除法上串行化），或者直接能够作为SIMD整数操作直接使用。因此，您可以将其直接用于GPU中，除了std::vector之外，它不受所有C++ CUDA编译器（以及不在OpenCL内核中）的支持。这些与主机环境相关的部分可以在内核执行后计算（与不包括emplace_back或用在GPU中的struct交换的部分）。

- huseyin tugrul buyukisik

---

哇，这是一个非常有创意的解决方案。不幸的是，n_attrs 可以在函数调用之间改变，它只在我展示的循环内是常量。更不用说编译的成本太大了（整个循环成本 < 100ms）。 - Huy Le

你说n_attrs只有5000个不同的值。缓存在这里有帮助吗？一个小函数的5000个不同版本不应该占用太多RAM，但如果它太随机，缓存内容可能会被抛弃。 - huseyin tugrul buyukisik

如果 n_attrs = 2^k，我已经使用 num >> k 和 num & (n_attrs - 1) 来替换 / %，这样速度更快。但是我正在尝试解决通用情况。 - Huy Le

关于整数除法的函数式编程黑科技：https://dev59.com/57r4oIgBc1ULPQZFYRUK 可能会对一些 AMD 处理器有所帮助。但在 Intel 上速度要慢得多，并且某些编译器标志（如舍入模式等）可能会完全破坏该算法。请自行决定是否使用。 - huseyin tugrul buyukisik

你的意思是我需要生成5000个不同的.DLL文件吗？ - Huy Le

如果使用仅缓冲区作为后端的虚拟文件系统，则可以防止其暴露5000个动态链接库，并且它会像RAMDISK一样具有快速访问速度。我在CUDA项目中使用了类似的方式，只是内核位于可执行文件的子文件夹中。如果仅用于内核而不是主机，则甚至可以将所有二进制代码合并到单个文本文件中，作为PTX格式或二进制可执行格式，并且操作系统已经可以将文件缓存在RAM中。但我也会使用非常简单的直接映射缓存来摆脱过多的文件API调用。 - huseyin tugrul buyukisik

Answer 4

所以在我的情况下，实际上最好的解决方案是：

不要使用 index = row * n_cols + col 的表示方法，而是使用 index = (row << 16) | col（32 位）或者 index = (row << 32) | col（64 位）。然后，row = index >> 32，col = index & (32 - 1)。或者更好的方法是，只需使用 uint16_t* pairs = reinterpret_cast<uint16_t*>(index_array);，然后对于每个 i % 2 == 0，pair[i], pair[i+1] 就是一对。

这是假设行数/列数小于 2^16（或 2^32）的情况。

我仍然保留最佳答案，因为它仍然回答了必须使用除法的情况。

- Huy Le