我将我的评论限制在整数除法上，因为对于一阶，C中的模运算可以被看作和实现为整数除法加上反乘和减法，尽管在某些情况下，有直接计算模的更便宜的方法，例如当计算模2n时。
整数除法在大多数平台上都相当慢，要么是基于软件仿真，要么是基于迭代硬件实现。但去年有广泛报道称，基于苹果M1的微基准测试，它的整数除法速度极快，可能是因为使用了专用电路。
自从Torbjörn Granlund和Peter蒙哥马利在大约三十年前发表了一篇开创性的论文以来，人们已经广泛地知道如何通过使用整数乘法加上可能的移位和/或其他校正步骤来用常数除数替换整数除法。这种算法通常被称为魔术乘法器技术。它需要预先计算整数除数的一些相关参数，以便在乘法中使用。基于仿真序列。
Torbjörn Granlund和Peter L.蒙哥马利，“使用乘法除以不变整数”，*ACM SIGPLAN通知 *，第29卷，1994年6月，第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;
}

如何优化这个循环？
这是libdivide的一个完美用例。此库旨在通过使用编译器在编译时使用的策略，在运行时加快常量除法。此库仅包含标头，因此不会创建任何运行时依赖项。它还支持除法的矢量化（即使用SIMD指令），在这种情况下，这无疑是为了大幅加快计算速度，而编译器在不显著改变循环的情况下无法做到这一点（最后，由于运行时定义的除数，它的效率会降低）注意libdivide的许可证是非常宽松的（zlib），因此您可以轻松地将其包含在您的项目中，而不需要严格的约束（如果您更改它，基本上只需要将其标记为已修改）。
如果header only-libraries不合适，那么你需要重新实现wheel。这个想法是把一个常数的除法转换成一个移位和乘法的序列。@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缓存的小块进行操作。如果块太小，那么线程同步也会成为一个问题。2最好的解决方案当然是使用块的滚动窗口。3所有这些实现起来肯定有点繁琐/棘手。
注意，SIMD指令可以用于除法部分（使用libdivide很容易）。您还需要拆分循环并使用块，但块不需要很大，因为没有同步开销。64个整数应该足够了。
请注意，最新的处理器可以像这样高效地计算除法，特别是对于32位整数（64位的往往要贵得多）。桤木湖的情况尤其如此，Zen 3和M1处理器（P核）。注意，在x86/x86-64处理器上，模和除都在一个指令中计算。还注意，虽然除具有相当大的延迟，许多处理器可以流水线化多个除法以获得合理的吞吐量。例如，32位div指令在Skylake上的延迟为23~28个周期，而吞吐量的倒数为4~6。这显然不是Zen 1/Zen 2上的情况。

我会通过以下方式优化// 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
• 链接和调用函数

因此，n_attrs是DLL的编译时常量值，编译器可以自动对该值进行大部分优化，如下所示：

• 当x是2的幂数值时，执行n&(x-1)而不是n%x
• 移位和乘法而不是除法
• 也可以使用其他优化，比如使用预先计算的x和y索引展开循环（因为x是已知的）

当在编译时知道更多的部分时，紧循环中的一些整数数学运算更容易被编译器SIMDify/向量化。
如果您的CPU是AMD的，您甚至可以尝试使用神奇的浮点运算来代替未知/未知除法来实现矢量化。
通过缓存n_attrs的所有（或大部分）值，您可以消除以下延迟：

• 字符串生成
• 编译
• 文件（DLL）阅读（假设DLL有一些面向对象的 Package ）

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

所以在我的情况下实际上最好的解决方案。
不要用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);，那么pair[i], pair[i+1] for each i % 2 == 0就是一对。
这是假设行数/列数小于2^16（或2^32）。
我仍然保留上面的答案，因为它仍然回答了必须使用除法的情况。