使用 const 运行时除数进行快速整数除法和取模

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 内存上）。

Answer 1

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

\n

基于软件模拟或迭代硬件实现，整数除法在大多数平台上都相当慢。但去年有广泛报道称，基于苹果 M1 的微基准测试，它具有极快的整数除法，大概是通过使用专用电路实现的。

\n

自从近三十年前 Torbj\xc3\xb6rn Granlund 和 Peter Montgomery 发表了一篇开创性论文以来，如何通过使用整数乘法加上可能的移位和/或其他校正步骤来用常除数替换整数除法已经广为人知。该算法通常被称为魔法乘数技术。它需要从整数除数中预先计算一些相关参数，以便在基于乘法的仿真序列中使用。

\n

Torbj\xc3\xb6rn Granlund 和 Peter L. Montgomery，“使用乘法除以不变整数”，ACM SIGPLAN 公告，卷。29，1994 年 6 月，第 61-72 页（在线）。

\n

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

\n

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

\n

下面的代码演示了参数预计算和通过乘法进行整数除法模拟的原理。通过在设计中投入比我愿意为这个答案花费的时间更多的时间，很可能可以识别出参数预计算和仿真的更有效的实现。

\n

#include <cstdio>\n#include <cstdlib>\n#include <cstdint>\n\n#define PORTABLE  (1)\n\nuint32_t ilog2 (uint32_t i)\n{\n    uint32_t t = 0;\n    i = i >> 1;\n    while (i) {\n        i = i >> 1;\n        t++;\n    }\n    return (t);\n}\n\n/* Based on: Granlund, T.; Montgomery, P.L.: "Division by Invariant Integers \n   using Multiplication". SIGPLAN Notices, Vol. 29, June 1994, pp. 61-72\n*/\nvoid prepare_magic (int32_t divisor, int32_t &multiplier, int32_t &add_mask, int32_t &sign_shift)\n{\n    uint32_t divisoru, d, n, i, j, two_to_31 = uint32_t (1) << 31;\n    uint64_t m_lower, m_upper, k, msb, two_to_32 = uint64_t (1) << 32;\n\n    divisoru = uint32_t (divisor);\n    d = (divisor < 0) ? (0 - divisoru) : divisoru;\n    i = ilog2 (d);\n    j = two_to_31 % d;\n    msb = two_to_32 << i;\n    k = msb / (two_to_31 - j);\n    m_lower = msb / d;\n    m_upper = (msb + k) / d;\n    n = ilog2 (uint32_t (m_lower ^ m_upper));\n    n = (n > i) ? i : n;\n    m_upper = m_upper >> n;\n    i = i - n;\n    multiplier = int32_t (uint32_t (m_upper));\n    add_mask = (m_upper >> 31) ? (-1) : 0;\n    sign_shift = int32_t ((divisoru & two_to_31) | i);\n}\n\nint32_t arithmetic_right_shift (int32_t a, int32_t s)\n{\n    uint32_t msb = uint32_t (1) << 31;\n    uint32_t ua = uint32_t (a);\n    ua = ua >> s;\n    msb = msb >> s;\n    return int32_t ((ua ^ msb) - msb);\n}\n\nint32_t magic_division (int32_t dividend, int32_t multiplier, int32_t add_mask, int32_t sign_shift)\n{\n    int64_t prod = int64_t (dividend) * multiplier;\n    int32_t quot = (int32_t)(uint64_t (prod) >> 32);\n    quot = int32_t (uint32_t (quot) + (uint32_t (dividend) & uint32_t (add_mask)));\n#if PORTABLE\n    const int32_t byte_mask = 0xff;\n    quot = arithmetic_right_shift (quot, sign_shift & byte_mask);\n#else // PORTABLE\n    quot = quot >> sign_shift; // must mask shift count & use arithmetic right shift\n#endif // PORTABLE\n    quot = int32_t (uint32_t (quot) + (uint32_t (dividend) >> 31));\n    if (sign_shift < 0) quot = -quot;\n    return quot;\n}\n\nint main (void)\n{\n    int32_t multiplier;\n    int32_t add_mask;\n    int32_t sign_shift;\n    int32_t divisor;\n    \n    for (divisor = -20; divisor <= 20; divisor++) {\n        /* avoid division by zero */\n        if (divisor == 0) {\n            divisor++;\n            continue;\n        }\n        printf ("divisor=%d\\n", divisor);\n        prepare_magic (divisor, multiplier, add_mask, sign_shift);\n        printf ("multiplier=%d add_mask=%d sign_shift=%d\\n", \n                multiplier, add_mask, sign_shift);\n        printf ("exhaustive test of dividends ... ");\n        uint32_t dividendu = 0;\n        do {\n            int32_t dividend = (int32_t)dividendu;\n            /* avoid overflow in signed integer division */\n            if ((divisor == (-1)) && (dividend == ((-2147483647)-1))) {\n                dividendu++;\n                continue;\n            }\n            int32_t res = magic_division (dividend, multiplier, add_mask, sign_shift);\n            int32_t ref = dividend / divisor;\n            if (res != ref) {\n                printf ("\\nERR dividend=%d (%08x) divisor=%d  res=%d  ref=%d\\n",\n                        dividend, (uint32_t)dividend, divisor, res, ref);\n                return EXIT_FAILURE;\n            }\n            dividendu++;\n        } while (dividendu);\n        printf ("PASSED\\n");\n    }\n    return EXIT_SUCCESS;\n}\n

\n