x86 - sse 命令を使用した複雑な Mul および Div

Question

複雑な乗算と除算を SSE 命令で実行することは有益ですか? SSE を使用すると、足し算と引き算のパフォーマンスが向上することがわかっています。SSE を使用して複雑な乗算を実行し、パフォーマンスを向上させる方法を誰かに教えてもらえますか?

Answer 1

完全を期すために、ここからダウンロードできるインテル® 64 および IA-32 アーキテクチャー最適化リファレンス・マニュアルには、複素数の乗算 (例 6-9) と複素数の除算 (例 6-10) のアセンブリが含まれています。

たとえば、乗算コードは次のとおりです。

// Multiplication of (ak + i bk ) * (ck + i dk )
// a + i b can be stored as a data structure
movsldup xmm0, src1; load real parts into the destination, a1, a1, a0, a0
movaps xmm1, src2; load the 2nd pair of complex values, i.e. d1, c1, d0, c0
mulps xmm0, xmm1; temporary results, a1d1, a1c1, a0d0, a0c0
shufps xmm1, xmm1, b1; reorder the real and imaginary parts, c1, d1, c0, d0
movshdup xmm2, src1; load imaginary parts into the destination, b1, b1, b0, b0
mulps xmm2, xmm1; temporary results, b1c1, b1d1, b0c0, b0d0
addsubps xmm0, xmm2; b1c1+a1d1, a1c1 -b1d1, b0c0+a0d0, ; a0c0-b0d0

アセンブリは、gccs X86 組み込み関数に直接マップされます(各命令を で述語するだけ__builtin_ia32_です)。

Answer 2

複雑な乗算は次のように定義されます。

((c1a * c2a) - (c1b * c2b)) + ((c1b * c2a) + (c1a * c2b))i

したがって、複素数の2つのコンポーネントは次のようになります

((c1a * c2a) - (c1b * c2b)) and ((c1b * c2a) + (c1a * c2b))i

したがって、次のように定義された 4 つの複素数を表すために 8 つの float を使用していると仮定します。

c1a, c1b, c2a, c2b
c3a, c3b, c4a, c4b

そして、(c1 * c3) と (c2 * c4) を同時に実行したい場合、SSE コードは次のような「何か」になります。

（注：WindowsでMSVCを使用しましたが、原則は同じです）。

__declspec( align( 16 ) ) float c1c2[]        = { 1.0f, 2.0f, 3.0f, 4.0f };
__declspec( align( 16 ) ) float c3c4[]          = { 4.0f, 3.0f, 2.0f, 1.0f };
__declspec( align( 16 ) ) float mulfactors[]    = { -1.0f, 1.0f, -1.0f, 1.0f };
__declspec( align( 16 ) ) float res[]           = { 0.0f, 0.0f, 0.0f, 0.0f };

__asm 
{
    movaps xmm0, xmmword ptr [c1c2]         // Load c1 and c2 into xmm0.
    movaps xmm1, xmmword ptr [c3c4]         // Load c3 and c4 into xmm1.
    movaps xmm4, xmmword ptr [mulfactors]   // load multiplication factors into xmm4

    movaps xmm2, xmm1                       
    movaps xmm3, xmm0                       
    shufps xmm2, xmm1, 0xA0                 // Change order to c3a c3a c4a c4a and store in xmm2
    shufps xmm1, xmm1, 0xF5                 // Change order to c3b c3b c4b c4b and store in xmm1
    shufps xmm3, xmm0, 0xB1                 // change order to c1b c1a c2b c2a abd store in xmm3

    mulps xmm0, xmm2                        
    mulps xmm3, xmm1                    
    mulps xmm3, xmm4                        // Flip the signs of the 'a's so the add works correctly.

    addps xmm0, xmm3                        // Add together

    movaps xmmword ptr [res], xmm0          // Store back out
};

float res1a = (c1c2[0] * c3c4[0]) - (c1c2[1] * c3c4[1]);
float res1b = (c1c2[1] * c3c4[0]) + (c1c2[0] * c3c4[1]);

float res2a = (c1c2[2] * c3c4[2]) - (c1c2[3] * c3c4[3]);
float res2b = (c1c2[3] * c3c4[2]) + (c1c2[2] * c3c4[3]);

if ( res1a != res[0] || 
     res1b != res[1] || 
     res2a != res[2] || 
     res2b != res[3] )
{
    _exit( 1 );
}

上で行ったことは、数学を少し単純化したということです。以下を仮定します。

c1a c1b c2a c2b
c3a c3b c4a c4b

再配置すると、次のベクトルになります

0 => c1a c1b c2a c2b
1 => c3b c3b c4b c4b
2 => c3a c3a c4a c4a
3 => c1b c1a c2b c2a

次に、0 と 2 を乗算して取得します。

0 => c1a * c3a, c1b * c3a, c2a * c4a, c2b * c4a

次に、3 と 1 を掛け合わせて次のようにします。

3 => c1b * c3b, c1a * c3b, c2b * c4b, c2a * c4b

最後に、いくつかのフロートの符号を 3 で反転します。

3 => -(c1b * c3b), c1a * c3b, -(c2b * c4b), c2a * c4b

だから私はそれらを一緒に追加して取得することができます

(c1a * c3a) - (c1b * c3b), (c1b * c3a ) + (c1a * c3b), (c2a * c4a) - (c2b * c4b), (c2b * c4a) + (c2a * c4b)

これが私たちが求めていたものです:)

Answer 3

Intel最適化リファレンスのアルゴリズムはNaN、入力のオーバーフローとsを適切に処理しません。

NaN数値の実数部または虚数部の1つが、他の部分に誤って拡散します。

無限大（たとえば、無限大* 0）を使用するいくつかの操作はで終わるためNaN、オーバーフローによりNaN、他の点では正常に動作するデータにsが表示される可能性があります。

オーバーフローとNaNsがまれな場合、これを回避する簡単な方法はNaN、結果をチェックして、コンパイラのIEEE準拠の実装で再計算することです。

float complex a[2], b[2];
__m128 res = simd_fast_multiply(a, b);

/* store unconditionally, can be executed in parallel with the check
 * making it almost free if there is no NaN in data */
_mm_store_ps(dest, res);

/* check for NaN */
__m128 n = _mm_cmpneq_ps(res, res);
int have_nan = _mm_movemask_ps(n);
if (have_nan != 0) {
    /* do it again unvectorized */
    dest[0] = a[0] * b[0];
    dest[1] = a[1] * b[1];
}