ホーム>source

過去数週間、私は人々が車輪を再発明し、例えば自分のsqrt関数を書くのに何時間も費やそうとしている点は何だろうと思っていました。組み込みバージョンは、最適化され、正確で、十分に安定しています。

たとえば、カーマックスタイルの平方根について話しています。ポイントは?近似中に精度が失われ、キャストが使用されます。

IntelスタイルのSSE Square Rootは正確な結果を出していましたが、私の計算では標準のSQRTよりも遅くなりました。

平均して、上記のトリックはすべて標準のSQRTをはるかに上回っています。だから私の質問は、ポイントは何ですか?

私のPCには以下のCPUがあります:

Intel(R)Core(TM)i7-6700HQ CPU @ 2.60GHz。

各メソッドについて次の結果が得られました(有用なコメントによる以下の提案に従って、パフォーマンステストを修正しました。そのn.mに感謝します)。

(ニュートン法のような近似を使用している場合、精度が失われるため、それに応じて計算を調整する必要があることに注意してください。)

参照用に以下のソースコードを見つけることができます。

#include <chrono>
#include <cmath>
#include <deque>
#include <iomanip>
#include <iostream>
#include <immintrin.h>
#include <random>
using f64 = double;
using s64 = int64_t;
using u64 = uint64_t;
static constexpr u64 cycles = 24;
static constexpr u64 sample_max = 1000000;
f64 sse_sqrt(const f64 x) {
    __m128d root = _mm_sqrt_pd(_mm_load_pd(&x));
    return *(reinterpret_cast<f64*>(&root));
}
constexpr f64 carmack_sqrt(const f64 x) {
    union {
        f64 x;
        s64 i;
    } u = {};
    u.x = x;
    u.i = 0x5fe6eb50c7b537a9 - (u.i >> 1);
    f64 xhalf = 0.5 * x;
    u.x = u.x * (1.5 - xhalf * u.x * u.x);
    # u.x = u.x * (1.5 - xhalf * u.x * u.x);
    # u.x = u.x * (1.5 - xhalf * u.x * u.x);
    # ... so on, if you want more precise result ...
    return u.x * x;
}
int main(int /* argc */, char ** /*argv*/) {
    std::random_device r;
    std::default_random_engine e(r());
    std::uniform_real_distribution<f64> dist(1, sample_max);
    std::deque<f64> samples(sample_max);
    for (auto& sample : samples) {
        sample = dist(e);
    }
    // std sqrt
    {
        std::cout << "> Measuring std sqrt.\r\n> Please wait . . .\r\n";
        f64 result = 0;
        auto t1 = std::chrono::high_resolution_clock::now();
        for (auto cycle = 0; cycle < cycles; ++cycle) {
            for (auto& sample : samples) {
                result += std::sqrt(static_cast<f64>(sample));
            }
        }
        auto t2 = std::chrono::high_resolution_clock::now();
        auto dt = t2 - t1;
        std::cout << "> Accumulated result: " << std::setprecision(19) << result << "\n";
        std::cout << "> Total execution time: " <<
        std::chrono::duration_cast<std::chrono::milliseconds>(dt).count() << " ms.\r\n\r\n";
    }
    // sse sqrt
    {
        std::cout << "> Measuring sse sqrt.\r\n> Please wait . . .\r\n";
        f64 result = 0;
        auto t1 = std::chrono::high_resolution_clock::now();
        for (auto cycle = 0; cycle < cycles; ++cycle) {
            for (auto& sample : samples) {
                result += sse_sqrt(static_cast<f64>(sample));
            }
        }
        auto t2 = std::chrono::high_resolution_clock::now();
        auto dt = t2 - t1;
        std::cout << "> Accumulated result: " << std::setprecision(19) << result << "\n";
        std::cout << "> Total execution time: " <<
        std::chrono::duration_cast<std::chrono::milliseconds>(dt).count() << " ms.\r\n\r\n";
    }
    // carmack sqrt
    {
        std::cout << "> Measuring carmack sqrt.\r\n> Please wait . . .\r\n";
        f64 result = 0;
        auto t1 = std::chrono::high_resolution_clock::now();
        for (auto cycle = 0; cycle < cycles; ++cycle) {
            for (auto& sample : samples) {
                result += carmack_sqrt(static_cast<f64>(sample));
           }
        }
        auto t2 = std::chrono::high_resolution_clock::now();
        auto dt = t2 - t1;
        std::cout << "> Accumulated result: " << std::setprecision(19) << result << "\n";
        std::cout << "> Total execution time: " <<
        std::chrono::duration_cast<std::chrono::milliseconds>(dt).count() << " ms.\r\n\r\n";
    }
    std::cout << "> Press any key to exit . . .\r\n";
    std::getchar();
    return 0;
}

私は誰かを批判するためにここにいるのではないことに注意してください。私はここで自分の方法と選択するための最良のツールセットを学び、実験し、理解しようとしています。

ポートフォリオの1つに自分のゲームエンジンを書いています。私はあなたの親切な答えに感謝しており、私はどんな提案に対してもオープンです。

ごきげんよう。

あなたの答え
  • 解決した方法 # 1

    楽しさと利益のために?

    あなたの質問に基づいて、そうする理由はありませんが、言語を学習したい場合は、整数/浮動小数点((ほとんどの場合)すべての言語のプリミティブ)に依存しているため、数学的な問題を解決することをお勧めしますアルゴリズムは十分に文書化されています。

    「実際の」コードでは、libcが提供するメソッドを使用している限り、それを使用する必要があります。通常、組み込みプラットフォームにはlibcがないか、独自のプラットフォームをロールバックするため、独自のプラットフォームを実装する必要があります。

  • 解決した方法 # 2

    その高速の逆数平方根のトリックは、ほとんど時代遅れです。 SSEは、Pentium 3がPCプラットフォーム上で完全に置き換えたために存在する近似の逆平方根に組み込まれています。通常、他のプラットフォームには独自の逆数平方根があります。たとえば、ARMには VRSQRTE があります  そして、ニュートンのステップも行う便利な指示。

    ちなみに、結果を非可逆平方根にすると、通常はあまり役に立たなくなります。主なユースケースは、ベクトルを正規化することです。ここで、「真っ直ぐな」平方根は、逆数ながら、平方根は正確に適合します(それは乗算です)。

    多くの場合、ベンチマークは非常に正確ではありません。私はしばらく前にいくつかの関連するテストを行ったことがありますが、関連する部分は次のようになります。

    std::sqrt  ベース:

    HMM_INLINE float HMM_LengthVec4(hmm_vec4 A)
    {
        float Result = std::sqrt(HMM_LengthSquaredVec4(A));
        return(Result);
    }
    HMM_INLINE hmm_vec4 HMM_NormalizeVec4(hmm_vec4 A)
    {
        hmm_vec4 Result = {0};
        float VectorLength = HMM_LengthVec4(A);
        /* NOTE(kiljacken): We need a zero check to not divide-by-zero */
        if (VectorLength != 0.0f)
        {
            float Multiplier = 1.0f / VectorLength;
    #ifdef HANDMADE_MATH__USE_SSE
            __m128 SSEMultiplier = _mm_set1_ps(Multiplier);
            Result.InternalElementsSSE = _mm_mul_ps(A.InternalElementsSSE, SSEMultiplier);        
    #else 
            Result.X = A.X * Multiplier;
            Result.Y = A.Y * Multiplier;
            Result.Z = A.Z * Multiplier;
            Result.W = A.W * Multiplier;
    #endif
        }
        return (Result);
    }
    
    

    SSEの逆平方根とニュートンステップ:

    HMM_INLINE hmm_vec4 HMM_NormalizeVec4_new(hmm_vec4 A)
    {
        hmm_vec4 Result;
        // square elements and add them together, result is in every lane
        __m128 t0 = _mm_mul_ps(A.InternalElementsSSE, A.InternalElementsSSE);
        __m128 t1 = _mm_add_ps(t0, _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 3, 0, 1)));
        __m128 sq = _mm_add_ps(t1, _mm_shuffle_ps(t1, t1, _MM_SHUFFLE(0, 1, 2, 3)));
        // compute reciprocal square root with Newton step for ~22bit accuracy
        __m128 rLen = _mm_rsqrt_ps(sq);
        __m128 half = _mm_set1_ps(0.5);
        __m128 threehalf = _mm_set1_ps(1.5);
        __m128 t = _mm_mul_ps(_mm_mul_ps(sq, half), _mm_mul_ps(rLen, rLen));
        rLen = _mm_mul_ps(rLen, _mm_sub_ps(threehalf, t));
        // multiply elements by the reciprocal of the vector length
        __m128 normed = _mm_mul_ps(A.InternalElementsSSE, rLen);
        // normalize zero-vector to zero, not to NaN
        __m128 zero = _mm_setzero_ps();
        Result.InternalElementsSSE = _mm_andnot_ps(_mm_cmpeq_ps(A.InternalElementsSSE, zero), normed);
        return (Result);
    }
    
    

    SSEの逆平方根なしで ニュートンステップ:

    HMM_INLINE hmm_vec4 HMM_NormalizeVec4_lowacc(hmm_vec4 A)
    {
        hmm_vec4 Result;
        // square elements and add them together, result is in every lane
        __m128 t0 = _mm_mul_ps(A.InternalElementsSSE, A.InternalElementsSSE);
        __m128 t1 = _mm_add_ps(t0, _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 3, 0, 1)));
        __m128 sq = _mm_add_ps(t1, _mm_shuffle_ps(t1, t1, _MM_SHUFFLE(0, 1, 2, 3)));
        // compute reciprocal square root without Newton step for ~12bit accuracy
        __m128 rLen = _mm_rsqrt_ps(sq);
        // multiply elements by the reciprocal of the vector length
        __m128 normed = _mm_mul_ps(A.InternalElementsSSE, rLen);
        // normalize zero-vector to zero, not to NaN
        __m128 zero = _mm_setzero_ps();
        Result.InternalElementsSSE = _mm_andnot_ps(_mm_cmpeq_ps(A.InternalElementsSSE, zero), normed);
        return (Result);
    }
    
    

    (クイックベンチ)

    ご覧のとおり、スループットとレイテンシを個別に測定しましたが、その違いは非常に重要でした。ニュートンステップの逆平方根は、通常の平方根を使用するのと同じくらい長い時間がかかりますが、より高いスループットで処理できます。ニュートンステップを使用しない場合、1回のベクトル正規化操作でも開始から終了までの時間が短くなり、スループットは以前よりもさらに向上します。とにかく、これはあなたの平方根について何かをするためのポイントがあることを実証する必要があります。

    ちなみに、上記のコードはグッドプラクティスを意味するものではなく、4つのベクトルを同時に正規化して、4つのSIMD操作を無駄に計算しないようにします。シングル (逆数)平方根。ただし、ここでは実際には問題ではありません。

  • 解決した方法 # 3

    Carmackの手法のポイントは、1990年代の浮動小数点演算から得られるよりも優れたパフォーマンスを整数演算から抽出することでした。それ以来、浮動小数点のパフォーマンスは大幅に改善されました!独自のベンチマークでわかるように。ハードウェアで同様の制約に直面した場合を除き、新しいコードでその手法を使用する実用的な理由はありませんが、i7ではそうではありません。

  • 解決した方法 # 4

    プロジェクトの範囲による標準ライブラリの使用制限

    What is the point implementing custom math functions in C++ (like SQRT)?

    他の回答で既に言及されていることに加えて、独自の(カスタム)数学関数を実装するプロジェクトの選択は、次の理由による可能性があります。

    プロジェクトに課された制限、例えば何らかの標準のため、コンパイラに同梱されている数学ライブラリを使用しないようにします。

    一例として、ISO 26262規格に準拠したASIL分類プロジェクトがあります。たとえば、適切な資格w.r.tを提供するコンパイラを使用します。プロジェクトのソースコードを正しくコンパイルする、しかしそれは十分な資格を提供しません付属の標準ライブラリ、ここで数学ライブラリは、ソースコードではなく、オブジェクトによってのみリンクできます(後者の場合、適切なテストとソースコードの資格は、プロジェクト自体によって記述できます)。

関連記事

  • 前へ java - JPAクエリ:サブクエリをグループ化条件に結合する
  • 次へ javascript - ReactボタンのonClickパス引数