Fast Inverse Square Root on x64

Viewed 9005
3 Answers

Yes, it works if using the correct magic number and corresponding integer type. In addition to the answers above, here's a C++11 implementation that works for both double and float. Conditionals should optimise out at compile time.

template <typename T, char iterations = 2> inline T inv_sqrt(T x) {
    static_assert(std::is_floating_point<T>::value, "T must be floating point");
    static_assert(iterations == 1 or iterations == 2, "itarations must equal 1 or 2");
    typedef typename std::conditional<sizeof(T) == 8, std::int64_t, std::int32_t>::type Tint;
    T y = x;
    T x2 = y * 0.5;
    Tint i = *(Tint *)&y;
    i = (sizeof(T) == 8 ? 0x5fe6eb50c7b537a9 : 0x5f3759df) - (i >> 1);
    y = *(T *)&i;
    y = y * (1.5 - (x2 * y * y));
    if (iterations == 2)
        y = y * (1.5 - (x2 * y * y));
    return y;
}

As for testing, I use the following doctest in my project:

#ifdef DOCTEST_LIBRARY_INCLUDED
    TEST_CASE_TEMPLATE("inv_sqrt", T, double, float) {
        std::vector<T> vals = {0.23, 3.3, 10.2, 100.45, 512.06};
        for (auto x : vals)
            CHECK(inv_sqrt<T>(x) == doctest::Approx(1.0 / std::sqrt(x)));
    }
#endif
Related