c++ - 在 C++(如 SQRT)中实现自定义数学函数有什么意义？

Question

在过去的几周里，我想知道人们试图重新发明轮子并花费数小时编写自己的 sqrt 函数的意义何在。内置版本优化良好，足够精确和稳定。

例如，我说的是卡马克式平方根。重点是什么？它会在近似过程中失去精度，并使用转换。

英特尔风格的 SSE 平方根给出了精确的结果，但我的计算速度比标准 SQRT 慢。

平均而言，上述所有技巧均被标准 SQRT 击败。所以我的问题是，这有什么意义？

我的电脑有以下 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;
}

请注意，我不是来这里批评任何人的，我来这里只是为了学习、试验并试图找出我自己的方法和最佳工具集。

我正在将自己的游戏引擎写入我的作品集之一。感谢您的友好回答，我愿意接受任何建议。

祝你有美好的一天。

Answer 1

快速平方根倒数的技巧基本上已经过时了。 SSE 内置了自 Pentium 3 在 PC 平台上完全取代它以来一直存在的近似平方根倒数。其他平台通常有自己的平方根倒数，例如 ARM 有 VRSQRTE 和一个执行牛顿步的便捷指令。

顺便说一句，将结果转换为非倒数平方根通常会降低它的用处:主要用例是规范化 vector ，其中“直”平方根很烦人(必须除以)而倒数平方根恰好符合(那么它是乘法)。

您的基准通常不太准确。前段时间刚好做了一些相关的测试，相关的部分是这样的:

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);
}

( quick-bench )

如您所见，我分别测量了吞吐量和延迟，两者之间的区别非常重要。使用牛顿步的平方根倒数需要很长时间，大约与使用普通平方根一样长，但可以以更高的吞吐量进行处理。没有牛顿步骤，单个 vector 归一化操作从开始到结束花费的时间也更少，并且吞吐量变得比以前更好。无论如何，这应该表明对您的平方根做一些事情是有意义的。

顺便说一下，上面的代码并不是好的做法，它会同时规范化 4 个 vector ，以免在计算 单个 (倒数)时浪费 4 宽 SIMD 操作平方根。不过，这并不是真正的问题。