使用const運行時除數的快速整數除法和取模-有解無憂

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 記憶體上）。

uj5u.com熱心網友回復：

我將我的評論限制在整數除法上，因為在 C 中首先對模運算進行排序可以被視為整數除法加上反乘和減法，盡管在某些情況下，有更便宜的直接計算模數的方法，例如當計算模 2 ⁿ時。

基于軟體仿真或迭代硬體實作，大多數平臺上的整數除法都非常慢。但去年廣泛報道稱，基于蘋果 M1 的微基準測驗，它具有極快的整數除法，大概是通過使用專用電路。

自從將近 30 年前 Torbj?rn Granlund 和 Peter Montgomery 的一篇開創性論文以來，眾所周知如何通過使用整數乘法加上可能的移位和/或其他校正步驟來用常數除數代替整數除法。這種演算法通常被稱為魔術乘法器技術。它需要從整數除數中預先計算一些相關引數，以用于基于乘法的仿真序列。

Torbj?rn Granlund 和 Peter L. Montgomery，“使用乘法除以不變整數”，ACM SIGPLAN 通知，卷。29，1994 年 6 月，第 61-72 頁（在線）。

目前，在處理編譯時常量的整數除數時，所有主要工具鏈都包含 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 (x) 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;
}

uj5u.com熱心網友回復：

我將// 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 的編譯時常量值，編譯器可以自動對值進行大部分優化，例如：

做n&(x-1)而不是n%x當 x 是 2 的冪值時
移位和乘法而不是除法
也許還有其他優化，例如使用預先計算的 x 和 y 索引展開回圈（因為 x 是已知的）

當更多部分在編譯時已知時，緊密回圈中的一些整數數學運算更容易通過編譯器進行 SIMDify/矢量化。

如果您的 CPU 是 AMD，您甚至可以嘗試使用魔術浮點運算來代替未知/未知除法來獲得矢量化。

通過快取的所有（或大部分）值n_attrs，您可以消除以下延遲：

字串生成
編譯
檔案（DLL）讀取（假設一些面向物件的 DLL 包裝）

如果要優化的部分將在 GPU 中運行，那么 CUDA/OpenCL 實作很有可能已經以浮點方式進行整數除法（以保持 SIMD 路徑被占用而不是在整數除法上被序列化）或者只是能夠直接作為 SIMD 整數運算，因此您可以像在 GPU 中一樣使用它，除了std::vector所有 C CUDA 編譯器不支持（并且不在 OpenCL 內核中）。這些與主機環境相關的部分可以在內核（不包括 emplace_back 或與在 GPU 中作業的結構交換的部分）執行后計算。

uj5u.com熱心網友回復：

如何優化這個回圈？

這是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 核）的情況。請注意，模數和除法都是在 x86/x86-64 處理器上的一條指令中計算的。另請注意，雖然除法具有相當大的延遲，但許多處理器可以流水線化多個除法以獲得合理的吞吐量。例如，一條 32 位div指令在 Skylake 上的延遲為 23~28 個周期，但倒數吞吐量為 4~6。Zen1/Zen2 顯然不是這種情況。

轉載請註明出處，本文鏈接：https://www.uj5u.com/qiye/507882.html

標籤：C 数学优化库达整数除法

上一篇：如何正確制作相對于其他vector2的vector2？

下一篇：無法實作Maclaurin的系列