c - C の高速ユークリッド除算

Question

ユークリッド除算の剰余を取得することに興味があります。つまり、整数のペア (i、n) について、次のような r を見つけます。

i = k * n + r, 0 <= r < |k|

簡単な解決策は次のとおりです。

int euc(int i, int n)
{
    int r;

    r = i % n;
    if ( r < 0) {
        r += n;
    }
    return r;
}

しかし、これを何千万回も実行する必要があるので（多次元配列のイテレータ内で使用されます）、できれば分岐は避けたいです。要件：

• 分岐するが高速であることも望ましい。
• 正の n に対してのみ機能するソリューションは許容されます (ただし、負の i に対して機能する必要があります)。
• n は前もってわからず、> 0 かつ < MAX_INT の任意の値にすることができます

編集

実際には間違った結果を得るのは非常に簡単なので、期待される結果の例を次に示します。

• euc(0, 3) = 0
• euc(1, 3) = 1
• euc(2, 3) = 2
• euc(3, 3) = 0
• euc(-1, 3) = 2
• euc(-2, 3) = 1
• euc(-3, 3) = 0

これを最適化しても意味がないのではないかと心配する人もいます。境界外のアイテムが元の配列を繰り返す「仮想配列」内のアイテムに置き換えられる多次元イテレータにこれが必要です。したがって、配列 x が [1, 2, 3, 4] の場合、仮想配列は [...., 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4]、たとえば x[-2] は x 1など...

次元 d の nd 配列の場合、すべての点について d ユークリッド除算が必要です。am^d カーネルを使用して an^d 配列を相関させる必要がある場合は、n^d * m^d * d ユークリッド除算が必要です。100x100x100 ポイントの 3D 画像と 5*5*5 ポイントのカーネルの場合、それはすでに ~ 4 億のユークリッド分割です。

Answer 1

編集：乗算や分岐はありません。

int euc(int i, int n)
{
    int r;

    r = i % n;
    r += n & (-(r < 0));

    return r;
}

生成されたコードは次のとおりです。MSVC++ 計測プロファイラー (私のテスト) と OP のテストによると、それらはほぼ同じように動作します。

; Original post
00401000  cdq              
00401001  idiv        eax,ecx 
00401003  mov         eax,edx 
00401005  test        eax,eax 
00401007  jge         euc+0Bh (40100Bh) 
00401009  add         eax,ecx 
0040100B  ret              

; Mine
00401020  cdq              
00401021  idiv        eax,ecx 
00401023  xor         eax,eax 
00401025  test        edx,edx 
00401027  setl        al   
0040102A  neg         eax  
0040102C  and         eax,ecx 
0040102E  add         eax,edx 
00401030  ret              

Answer 2

280Z28 と Christopher は、アセンブラー ゴルフを私よりもよくカバーしていると思います。それはランダム アクセスを処理します。

ただし、実際に行っていることは、配列全体を処理しているようです。キャッシュ ミスを回避することは、小さな分岐を回避するよりも何倍も優れた最適化であるため、明らかにメモリ キャッシュの理由から、可能であればこれを順番に実行したいと考えています。

その場合、最初に適切な境界チェックを行うことで、私が「ダッシュ」と呼ぶ内部ループを実行できます。次の k インクリメントがいずれかの配列の最小次元でオーバーフローにならないことを確認し、代わりに「物理」インデックスを毎回 1 ずつインクリメントする新しいさらに内側のループを使用して k ステップを「ダッシュ」します。別のidivを行うこと。あなたまたはコンパイラは、このループを展開したり、Duff のデバイスを使用したりできます。

カーネルが小さい場合、特に固定サイズの場合は、それ (または、追加する代わりに時々減算する適切な展開を伴うその倍数) が、おそらく「ダッシュ」の長さに使用する値です。コンパイル時の一定のダッシュ長がおそらく最適です。これにより、ユーザー (またはコンパイラ) がダッシュ ループを完全に展開し、継続条件を除外できるからです。これによりコードが大きくなりすぎて高速にならない限り、本質的に正のモジュロ演算全体が整数のインクリメントに置き換えられます。

カーネルが固定サイズではなく、最後の次元が非常に小さい場合が多い場合は、最も一般的なサイズに対して異なるバージョンの比較関数を使用し、それぞれでダッシュ ループを完全に展開することを検討してください。

もう 1 つの可能性は、(いずれかの配列で) オーバーフローが発生する次のポイントを計算し、その値までダッシュすることです。ダッシュ ループにはまだ継続条件がありますが、インクリメントのみを使用して可能な限り長くなります。

または、実行している操作が数値の等価性またはその他の単純な操作である場合 (「相関」が何であるかはわかりません)、SIMD 命令などを調べることができます。この場合、ダッシュの長さはの倍数にする必要があります。アーキテクチャ上で最も広い単一命令比較 (または適切な SIMD op)。ただし、これは私が経験したことではありません。

Answer 3

ブランチはありませんが、少しいじっています:

int euc2(int i, int n)
{
    int r;
    r = i % n;
    r += (((unsigned int)r) >> 31) * n;
    return r;
}

乗算なし:

int euc2(int i, int n)
{
    int r;
    r = i % n;
    r += (r >> 31) & n;
    return r;
}

これは与える ：

; _i$ = eax
; _n$ = ecx

cdq
idiv   ecx
mov eax, edx
sar eax, 31
and eax, ecx
add eax, edx

Answer 4

整数の乗算は、除算よりもはるかに高速です。N が既知の多数の呼び出しの場合、N による除算を N の疑似逆数による乗算に置き換えることができます。

これを例に説明します。N=29 を取ります。次に、疑似逆数 2^16/N を 1 回計算します: K=2259 (2259.86 から切り捨てられます...)。I は正で、I*K は 32 ビットに収まると仮定します。

Quo = (I*K)>>16;   // replaces the division, Quo <= I/N
Mod = I - Quo*N;   // Mod >= I%N
while (Mod >= N) Mod -= N;  // compensate for the approximation

私の例では、I=753 を考えてみましょう。Quo=25 と Mod=28 が得られます。（補償は必要ありません）

編集。

3D畳み込みの例では、 i%n への呼び出しのほとんどは 0..n-1 の i で行われるため、ほとんどの場合、次のような最初の行

if (i>=0 && i<n) return i;

コストがかかり、ここでは役に立たない idiv をバイパスします。

また、十分な RAM がある場合は、すべての次元を 2 の累乗に揃え、除算の代わりにビット操作 (シフト、および) を使用します。

編集2。

実際に 10^9 の呼び出しで試してみました。i%n: 2.93 秒、私のコード: 1.38 秒。これは I の境界を意味することに注意してください (I*K は 32 ビットに収まる必要があります)。

別の考え: 値が x+dx で、x が 0..n-1 で dx が小さい場合、以下はすべてのケースをカバーします。

if (i<0) return i+n; else if (i>=n) return i-n;
return i;

Answer 5

iが-nを下回らないことも保証できる場合は、モジュロの前にオプションの加算を置くだけです。そうすれば、ブランチは必要なく、モジュロは、必要がなければ追加したものを切り取ります。

int euc(int i, int n)
{
    return (i + n) % n;
}

iが-n未満の場合でも、この方法を使用できます。このような場合、値がどの範囲にあるかを正確に知っている可能性があります。したがって、nをiに追加する代わりに、x * nをiに追加できます。ここで、xは十分な範囲を与える任意の整数です。速度を上げるために（単一サイクルの乗算がないプロセッサの場合）、乗算の代わりに左シフトすることができます。

Answer 6

TSC を使用して gcc -O3 で全員の提案の時間を測定しました (定数 N の提案を除く)。すべて同じ時間 (1% 以内) でした。

私の考えでは、 ((i%n)+n)%n (分岐なし) または (i+(n<<16))%n (大きな n または極端に負の i の場合は明らかに失敗) の方が高速ですが、すべて同じ時間がかかりました。

Answer 7

int euc(int i, int n)
{
    return (i % n) + (((i % n) < 0) * n);
}

Answer 8

私は次の表現がとても好きです。

r = ((i%n)+n)%n; 

分解は非常に短いです。

r = ((i%n)+n)%n;

004135AC  mov         eax,dword ptr [i] 
004135AF  cdq              
004135B0  idiv        eax,dword ptr [n] 
004135B3  add         edx,dword ptr [n] 
004135B6  mov         eax,edx 
004135B8  cdq              
004135B9  idiv        eax,dword ptr [n] 
004135BC  mov         dword ptr [r],edx 

ジャンプがなく (コストがかかる可能性のある 2 つの idiv)、完全にインライン化できるため、関数呼び出しのオーバーヘッドが回避されます。

どう思いますか？

Answer 9

右シフトが算術でない場合、 JasonのバージョンにフォールバックするChristopher のバージョンを次に示します。

#include <limits.h>
static inline int euc(int i, int n)
{
    // check for arithmetic shift
    #if (-1 >> 1) == -1
        #define OFFSET ((i % n >> (sizeof(int) * CHAR_BIT - 1)) & n)
    #else
        #define OFFSET ((i % n < 0) * n)
    #endif

    return i % n + OFFSET;
}

imulの代わりに使用するため、フォールバック バージョンは遅くなるはずですand。

Answer 10

配列の次元が常に 2 のべき乗であることを保証できる場合は、次のようにすることができます。

r = (i & (n - 1));

ディメンションが特定のサブセットからのものであることをさらに保証できる場合は、次のことができます。

template<int n>
int euc(int i) {
    return (i & (n - 1));
}

int euc(int i, int n) {
    switch (n) {
        case 2: return euc<2>(i);
        case 4: return euc<4>(i);
    }
}

Answer 11

あなたは Eric Bainville への回答で、ほとんどの0 <= i < n場合、

if (i>=0 && i<n) return i;

euc()とにかくあなたの最初の行として。

とにかく比較を行っているので、それらを使用することもできます。

int euc(int i, int n)
{
    if (n <= i)            return i % n;
    else if (i < 0)        return ((i + 1) % n) + n - 1;
    else /* 0 <= i < n */  return i;  // fastest possible response for common case
}