sqrt の大幅に高速なバージョンをロールすることは可能ですか?

Question

私がプロファイリングしているアプリでは、一部のシナリオでは、この関数が合計実行時間の 10% 以上を占める可能性があることがわかりました。

私は何年にもわたって、卑劣な浮動小数点トリックを使用したより高速な sqrt 実装に関する議論を見てきましたが、そのようなものが最新の CPU で時代遅れかどうかはわかりません。

参考までに、MSVC++ 2008 コンパイラが使用されています...ただし、sqrt はそれほどオーバーヘッドを追加しないと思います。

に関する同様の議論については、ここも参照してください。 モフ 関数。

編集：参考のため、 これ は広く使用されている方法の 1 つですが、実際にははるかに高速なのでしょうか?それにしても最近のSQRTって何サイクルなんですか？

Answer 1

はい、トリックなしでも可能です。

1) 速度のために精度を犠牲にする:sqrt アルゴリズムは反復的であり、より少ない反復で再実装されます。

2) ルックアップテーブル:反復の開始点だけを指定することも、補間と組み合わせてそこまで到達することもできます。

3) キャッシュ:いつも同じ限られた値のセットを二乗しているのでしょうか?そうであれば、キャッシュはうまく機能する可能性があります。これは、同じサイズの多数の形状に対して同じことが計算されるグラフィックス アプリケーションで便利であることがわかりました。そのため、結果を便利にキャッシュできます。

Answer 2

ここに素晴らしい比較表があります。http://assemblyrequired.crashworks.org/timing-square-root/

簡単に言うと、SSE2 の ssqrts は FPU fsqrt より約 2 倍高速で、近似 + 反復はそれより約 4 倍 (全体で 8 倍) 高速です。

また、単精度 sqrt を取得しようとしている場合は、それが実際に得られるものであることを確認してください。float 引数を double に変換し、倍精度 sqrt を呼び出してから float に戻すコンパイラが少なくとも 1 つあると聞いたことがあります。

Answer 3

変更することで速度がさらに向上する可能性が非常に高くなります。 アルゴリズム 変更するよりも 実装:電話してみる sqrt() 通話を速くする代わりに、通話を少なくします。(そして、これは不可能だと思われる場合は、 sqrt() あなたが言及しているのは、まさにそのことです：の改善 アルゴリズム 平方根を計算するために使用されます。)

非常に頻繁に使用されるため、標準ライブラリの実装はおそらく sqrt() 一般的なケースではほぼ最適です。アルゴリズムがいくつかのショートカットを使用できる制限されたドメイン (たとえば、それほど精度が必要ない場合) がない限り、誰かがより高速な実装を思いつく可能性はほとんどありません。

この関数は実行時間の 10% を使用するため、たとえ実行時間の 75% しかかからない実装を思いついたとしても、注意してください。 std::sqrt(), 、これでも実行時間はわずかに短縮されます 2,5%. 。ほとんどのアプリケーションでは、時計を使用して測定する場合を除いて、ユーザーはこれに気づきません。

Answer 4

どれくらい正確な情報が必要ですか sqrt することが？妥当な近似値を非常に迅速に得ることができます。Quake3 の優れた点を参照してください 逆平方根 インスピレーションのための関数 (コードは GPL 化されているため、直接統合したくない場合があることに注意してください)。

Answer 5

これを修正したかどうかはわかりませんが、以前にこの問題について読んだことがあります。一番早いのは、 sqrt インラインアセンブリバージョンの関数。

多数の代替案の説明が表示されます ここ.

最も優れているのは、この魔法の断片です。

double inline __declspec (naked) __fastcall sqrt(double n)
{
    _asm fld qword ptr [esp+4]
    _asm fsqrt
    _asm ret 8
} 

標準よりも約 4.7 倍高速です sqrt 同じ精度で呼び出します。

Answer 6

ここでは、わずか 8KB のルックアップ テーブルを使用する高速な方法を示します。間違いは結果の約 0.5% です。表を簡単に拡大できるので、間違いが減ります。通常の sqrt() よりも約 5 倍高速に実行されます。

// LUT for fast sqrt of floats. Table will be consist of 2 parts, half for sqrt(X) and half for sqrt(2X).
const int nBitsForSQRTprecision = 11;                       // Use only 11 most sagnificant bits from the 23 of float. We can use 15 bits instead. It will produce less error but take more place in a memory. 
const int nUnusedBits   = 23 - nBitsForSQRTprecision;       // Amount of bits we will disregard
const int tableSize     = (1 << (nBitsForSQRTprecision+1)); // 2^nBits*2 because we have 2 halves of the table.
static short sqrtTab[tableSize]; 
static unsigned char is_sqrttab_initialized = FALSE;        // Once initialized will be true

// Table of precalculated sqrt() for future fast calculation. Approximates the exact with an error of about 0.5%
// Note: To access the bits of a float in C quickly we must misuse pointers.
// More info in: http://en.wikipedia.org/wiki/Single_precision
void build_fsqrt_table(void){
    unsigned short i;
    float f;
    UINT32 *fi = (UINT32*)&f;

    if (is_sqrttab_initialized)
        return;

    const int halfTableSize = (tableSize>>1);
    for (i=0; i < halfTableSize; i++){
         *fi = 0;
         *fi = (i << nUnusedBits) | (127 << 23); // Build a float with the bit pattern i as mantissa, and an exponent of 0, stored as 127

         // Take the square root then strip the first 'nBitsForSQRTprecision' bits of the mantissa into the table
         f = sqrtf(f);
         sqrtTab[i] = (short)((*fi & 0x7fffff) >> nUnusedBits);

         // Repeat the process, this time with an exponent of 1, stored as 128
         *fi = 0;
         *fi = (i << nUnusedBits) | (128 << 23);
         f = sqrtf(f);
         sqrtTab[i+halfTableSize] = (short)((*fi & 0x7fffff) >> nUnusedBits);
    }
    is_sqrttab_initialized = TRUE;
}

// Calculation of a square root. Divide the exponent of float by 2 and sqrt() its mantissa using the precalculated table.
float fast_float_sqrt(float n){
    if (n <= 0.f) 
        return 0.f;                           // On 0 or negative return 0.
    UINT32 *num = (UINT32*)&n;
    short e;                                  // Exponent
    e = (*num >> 23) - 127;                   // In 'float' the exponent is stored with 127 added.
    *num &= 0x7fffff;                         // leave only the mantissa 

    // If the exponent is odd so we have to look it up in the second half of the lookup table, so we set the high bit.
    const int halfTableSize = (tableSize>>1);
    const int secondHalphTableIdBit = halfTableSize << nUnusedBits;
    if (e & 0x01) 
        *num |= secondHalphTableIdBit;  
    e >>= 1;                                  // Divide the exponent by two (note that in C the shift operators are sign preserving for signed operands

    // Do the table lookup, based on the quaternary mantissa, then reconstruct the result back into a float
    *num = ((sqrtTab[*num >> nUnusedBits]) << nUnusedBits) | ((e + 127) << 23);
    return n;
}