emit_x64_floating_point: Optimize 32-bit EmitFPRSqrtEstimate

CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.8)
-project(dynarmic C CXX)
+project(dynarmic C CXX ASM)
 
 # Determine if we're built as a subproject (using add_subdirectory)
 # or if this is the master project.

src/backend/x64/emit_x64_floating_point.cpp

@@ -995,6 +995,111 @@ static void EmitFPRSqrtEstimate(BlockOfCode& code, EmitContext& ctx, IR::Inst* i
             ctx.reg_alloc.DefineValue(inst, result);
             return;
         }
+
+        // TODO: VRSQRT14SS implementation (AVX512F)
+
+        auto args = ctx.reg_alloc.GetArgumentInfo(inst);
+
+        const Xbyak::Xmm operand = ctx.reg_alloc.UseXmm(args[0]);
+        const Xbyak::Xmm result = ctx.reg_alloc.ScratchXmm();
+        const Xbyak::Xmm value = ctx.reg_alloc.ScratchXmm();
+        [[maybe_unused]] const Xbyak::Reg32 tmp = ctx.reg_alloc.ScratchGpr().cvt32();
+
+        Xbyak::Label fallback, bad_values, end;
+
+        if constexpr (fsize == 64) {
+            code.cvtsd2ss(value, operand);
+
+            if (ctx.FPCR().RMode() == FP::RoundingMode::TowardsMinusInfinity || ctx.FPCR().RMode() == FP::RoundingMode::TowardsZero) {
+                code.ucomiss(value, code.MConst(xword, FP::FPInfo<u32>::MaxNormal(false)));
+                code.je(bad_values, code.T_NEAR);
+            }
+        } else {
+            code.movaps(value, operand);
+        }
+
+        code.movaps(xmm0, code.MConst(xword, 0xFFFF8000));
+        code.pand(value, xmm0);
+        code.por(value, code.MConst(xword, 0x00008000));
+
+        // Detect NaNs, negatives, zeros, denormals and infinities
+        code.ucomiss(value, code.MConst(xword, 0x00800000));
+        code.jna(bad_values, code.T_NEAR);
+
+        code.sqrtss(value, value);
+
+        code.movd(result, code.MConst(xword, 0x3F800000));
+        code.divss(result, value);
+
+        code.paddd(result, code.MConst(xword, 0x00004000));
+        code.pand(result, xmm0);
+
+        if constexpr (fsize == 64) {
+            code.cvtss2sd(result, result);
+        }
+        code.L(end);
+
+        code.SwitchToFarCode();
+
+        code.L(bad_values);
+        bool needs_fallback = false;
+        if constexpr (fsize == 32) {
+            Xbyak::Label default_nan;
+
+            code.movd(tmp, operand);
+
+            if (!ctx.FPCR().FZ()) {
+                if (ctx.FPCR().DN()) {
+                    // a > 0x80000000
+                    code.cmp(tmp, 0x80000000);
+                    code.ja(default_nan, code.T_NEAR);
+                }
+
+                // a > 0 && a < 0x00800000;
+                code.sub(tmp, 1);
+                code.cmp(tmp, 0x007FFFFF);
+                code.jb(fallback);
+                needs_fallback = true;
+            }
+
+            code.rsqrtss(result, operand);
+
+            if (ctx.FPCR().DN()) {
+                code.ucomiss(result, result);
+                code.jnp(end, code.T_NEAR);
+            } else {
+                // FZ ? (a >= 0x80800000 && a <= 0xFF800000) : (a >= 0x80000001 && a <= 0xFF800000)
+                // !FZ path takes into account the subtraction by one from the earlier block
+                code.add(tmp, ctx.FPCR().FZ() ? 0x7F800000 : 0x80000000);
+                code.cmp(tmp, ctx.FPCR().FZ() ? 0x7F000001 : 0x7F800000);
+                code.jnb(end, code.T_NEAR);
+            }
+
+            code.L(default_nan);
+            code.movd(result, code.MConst(xword, 0x7FC00000));
+            code.jmp(end, code.T_NEAR);
+        } else {
+            needs_fallback = true;
+        }
+
+        code.L(fallback);
+        if (needs_fallback) {
+            code.sub(rsp, 8);
+            ABI_PushCallerSaveRegistersAndAdjustStackExcept(code, HostLocXmmIdx(result.getIdx()));
+            code.movq(code.ABI_PARAM1, operand);
+            code.mov(code.ABI_PARAM2.cvt32(), ctx.FPCR().Value());
+            code.lea(code.ABI_PARAM3, code.ptr[code.r15 + code.GetJitStateInfo().offsetof_fpsr_exc]);
+            code.CallFunction(&FP::FPRSqrtEstimate<FPT>);
+            code.movq(result, rax);
+            ABI_PopCallerSaveRegistersAndAdjustStackExcept(code, HostLocXmmIdx(result.getIdx()));
+            code.add(rsp, 8);
+            code.jmp(end, code.T_NEAR);
+        }
+
+        code.SwitchToNearCode();
+
+        ctx.reg_alloc.DefineValue(inst, result);
+        return;
     }
 
     auto args = ctx.reg_alloc.GetArgumentInfo(inst);

src/backend/x64/emit_x64_vector_floating_point.cpp

@@ -406,22 +406,17 @@ void EmitThreeOpVectorOperation(BlockOfCode& code, EmitContext& ctx, IR::Inst* i
     ctx.reg_alloc.DefineValue(inst, result);
 }
 
-template<size_t fpcr_controlled_arg_index = 1, typename Lambda>
+template<typename Lambda>
-void EmitTwoOpFallback(BlockOfCode& code, EmitContext& ctx, IR::Inst* inst, Lambda lambda) {
+void EmitTwoOpFallbackWithoutRegAlloc(BlockOfCode& code, EmitContext& ctx, Xbyak::Xmm result, Xbyak::Xmm arg1, Lambda lambda, bool fpcr_controlled) {
     const auto fn = static_cast<mp::equivalent_function_type<Lambda>*>(lambda);
 
-    auto args = ctx.reg_alloc.GetArgumentInfo(inst);
+    const u32 fpcr = ctx.FPCR(fpcr_controlled).Value();
-    const bool fpcr_controlled = args[fpcr_controlled_arg_index].GetImmediateU1();
-    const Xbyak::Xmm arg1 = ctx.reg_alloc.UseXmm(args[0]);
-    const Xbyak::Xmm result = ctx.reg_alloc.ScratchXmm();
-    ctx.reg_alloc.EndOfAllocScope();
-    ctx.reg_alloc.HostCall(nullptr);
 
     constexpr u32 stack_space = 2 * 16;
     code.sub(rsp, stack_space + ABI_SHADOW_SPACE);
     code.lea(code.ABI_PARAM1, ptr[rsp + ABI_SHADOW_SPACE + 0 * 16]);
     code.lea(code.ABI_PARAM2, ptr[rsp + ABI_SHADOW_SPACE + 1 * 16]);
-    code.mov(code.ABI_PARAM3.cvt32(), ctx.FPCR(fpcr_controlled).Value());
+    code.mov(code.ABI_PARAM3.cvt32(), fpcr);
     code.lea(code.ABI_PARAM4, code.ptr[code.r15 + code.GetJitStateInfo().offsetof_fpsr_exc]);
 
     code.movaps(xword[code.ABI_PARAM2], arg1);
@@ -429,6 +424,19 @@ void EmitTwoOpFallback(BlockOfCode& code, EmitContext& ctx, IR::Inst* inst, Lamb
     code.movaps(result, xword[rsp + ABI_SHADOW_SPACE + 0 * 16]);
 
     code.add(rsp, stack_space + ABI_SHADOW_SPACE);
+}
+
+template<size_t fpcr_controlled_arg_index = 1, typename Lambda>
+void EmitTwoOpFallback(BlockOfCode& code, EmitContext& ctx, IR::Inst* inst, Lambda lambda) {
+    auto args = ctx.reg_alloc.GetArgumentInfo(inst);
+    const Xbyak::Xmm arg1 = ctx.reg_alloc.UseXmm(args[0]);
+    const Xbyak::Xmm result = ctx.reg_alloc.ScratchXmm();
+    ctx.reg_alloc.EndOfAllocScope();
+    ctx.reg_alloc.HostCall(nullptr);
+
+    const bool fpcr_controlled = args[fpcr_controlled_arg_index].GetImmediateU1();
+
+    EmitTwoOpFallbackWithoutRegAlloc(code, ctx, result, arg1, lambda, fpcr_controlled);
 
     ctx.reg_alloc.DefineValue(inst, result);
 }
@@ -1435,6 +1443,12 @@ template<size_t fsize>
 static void EmitRSqrtEstimate(BlockOfCode& code, EmitContext& ctx, IR::Inst* inst) {
     using FPT = mp::unsigned_integer_of_size<fsize>;
 
+    const auto fallback_fn = [](VectorArray<FPT>& result, const VectorArray<FPT>& operand, FP::FPCR fpcr, FP::FPSR& fpsr) {
+        for (size_t i = 0; i < result.size(); i++) {
+            result[i] = FP::FPRSqrtEstimate<FPT>(operand[i], fpcr, fpsr);
+        }
+    };
+
     if constexpr (fsize != 16) {
         if (ctx.HasOptimization(OptimizationFlag::Unsafe_ReducedErrorFP)) {
             auto args = ctx.reg_alloc.GetArgumentInfo(inst);
@@ -1454,11 +1468,7 @@ static void EmitRSqrtEstimate(BlockOfCode& code, EmitContext& ctx, IR::Inst* ins
         }
     }
 
-    EmitTwoOpFallback(code, ctx, inst, [](VectorArray<FPT>& result, const VectorArray<FPT>& operand, FP::FPCR fpcr, FP::FPSR& fpsr) {
+    EmitTwoOpFallback(code, ctx, inst, fallback_fn);
-        for (size_t i = 0; i < result.size(); i++) {
-            result[i] = FP::FPRSqrtEstimate<FPT>(operand[i], fpcr, fpsr);
-        }
-    });
 }
 
 void EmitX64::EmitFPVectorRSqrtEstimate16(EmitContext& ctx, IR::Inst* inst) {

src/common/fp/op/FPRSqrtEstimate.cpp

@@ -35,6 +35,7 @@ FPT FPRSqrtEstimate(FPT op, FPCR fpcr, FPSR& fpsr) {
     }
 
     if (type == FPType::Infinity) {
+        // Note: Just +Inf reaches here, negatives are handled in the case above.
         return FPInfo<FPT>::Zero(false);
     }
 

tests/CMakeLists.txt

@@ -19,6 +19,13 @@ add_executable(dynarmic_tests
     rand_int.h
 )
 
+if (NOT MSVC)
+    target_sources(dynarmic_tests PRIVATE
+        rsqrt_test.cpp
+        rsqrt_test_fn.s
+    )
+endif()
+
 if (DYNARMIC_TESTS_USE_UNICORN)
     target_sources(dynarmic_tests PRIVATE
         A32/fuzz_arm.cpp
@@ -47,7 +54,7 @@ create_target_directory_groups(dynarmic_print_info)
 target_link_libraries(dynarmic_tests PRIVATE dynarmic boost catch fmt mp xbyak)
 target_include_directories(dynarmic_tests PRIVATE . ../src)
 target_compile_options(dynarmic_tests PRIVATE ${DYNARMIC_CXX_FLAGS})
-target_compile_definitions(dynarmic_tests PRIVATE FMT_USE_USER_DEFINED_LITERALS=0)
+target_compile_definitions(dynarmic_tests PRIVATE FMT_USE_USER_DEFINED_LITERALS=0 CATCH_CONFIG_ENABLE_BENCHMARKING=1)
 
 target_link_libraries(dynarmic_print_info PRIVATE dynarmic boost catch fmt mp)
 target_include_directories(dynarmic_print_info PRIVATE . ../src)

tests/rsqrt_test.cpp

@@ -0,0 +1,148 @@
+/* This file is part of the dynarmic project.
+ * Copyright (c) 2021 MerryMage
+ * SPDX-License-Identifier: 0BSD
+ */
+
+#include <catch.hpp>
+#include <fmt/printf.h>
+
+#include "common/common_types.h"
+#include "common/fp/fpcr.h"
+#include "common/fp/fpsr.h"
+#include "common/fp/op/FPRSqrtEstimate.h"
+
+extern "C" u32 rsqrt_inaccurate(u32);
+extern "C" u32 rsqrt_full(u32);
+extern "C" u32 rsqrt_full_gpr(u32);
+extern "C" u32 rsqrt_full_nb(u32);
+extern "C" u32 rsqrt_full_nb2(u32);
+extern "C" u32 rsqrt_full_nb_gpr(u32);
+extern "C" u32 rsqrt_newton(u32);
+extern "C" u32 rsqrt_hack(u32);
+
+using namespace Dynarmic;
+
+extern "C" u32 rsqrt_fallback(u32 value) {
+    FP::FPCR fpcr;
+    FP::FPSR fpsr;
+    return FP::FPRSqrtEstimate(value, fpcr, fpsr);
+}
+
+void Test(u32 value) {
+    FP::FPCR fpcr;
+    FP::FPSR fpsr;
+
+    const u32 expect = FP::FPRSqrtEstimate(value, fpcr, fpsr);
+    const u32 full = rsqrt_full(value);
+    const u32 full_gpr = rsqrt_full_gpr(value);
+    const u32 newton = rsqrt_newton(value);
+    const u32 hack = rsqrt_hack(value);
+
+    if (expect != full || expect != full_gpr || expect != newton || expect != hack) {
+        fmt::print("{:08x} = {:08x} : {:08x} : {:08x} : {:08x} : {:08x}\n", value, expect, full, full_gpr, newton, hack);
+
+        REQUIRE(expect == full);
+        REQUIRE(expect == full_gpr);
+        REQUIRE(expect == newton);
+        REQUIRE(expect == hack);
+    }
+}
+
+TEST_CASE("RSqrt Tests", "[fp][.]") {
+    Test(0x00000000);
+    Test(0x80000000);
+    Test(0x7f8b7201);
+    Test(0x7f800000);
+    Test(0x7fc00000);
+    Test(0xff800000);
+    Test(0xffc00000);
+    Test(0xff800001);
+
+    for (u64 i = 0; i < 0x1'0000'0000; i++) {
+        const u32 value = static_cast<u32>(i);
+        Test(value);
+    }
+}
+
+TEST_CASE("Benchmark RSqrt", "[fp][.]") {
+    BENCHMARK("Inaccurate") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_inaccurate(value);
+        }
+        return total;
+    };
+
+    BENCHMARK("Full divss") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_full(value);
+        }
+        return total;
+    };
+
+    BENCHMARK("Full divss (GPR)") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_full_gpr(value);
+        }
+        return total;
+    };
+
+    BENCHMARK("Full divss (NB)") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_full_nb(value);
+        }
+        return total;
+    };
+
+    BENCHMARK("Full divss (NB2)") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_full_nb2(value);
+        }
+        return total;
+    };
+
+    BENCHMARK("Full divss (NB + GPR)") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_full_nb_gpr(value);
+        }
+        return total;
+    };
+
+    BENCHMARK("One Newton iteration") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_newton(value);
+        }
+        return total;
+    };
+
+    BENCHMARK("Ugly Hack") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_hack(value);
+        }
+        return total;
+    };
+
+    BENCHMARK("Softfloat") {
+        u64 total = 0;
+        for (u64 i = 0; i < 0x1'0000'0000; i += 0x1234) {
+            const u32 value = static_cast<u32>(i);
+            total += rsqrt_fallback(value);
+        }
+        return total;
+    };
+}

tests/rsqrt_test_fn.s

@@ -0,0 +1,303 @@
+.global _rsqrt_inaccurate
+.global rsqrt_inaccurate
+.global _rsqrt_full
+.global rsqrt_full
+.global _rsqrt_full_gpr
+.global rsqrt_full_gpr
+.global _rsqrt_full_nb
+.global rsqrt_full_nb
+.global _rsqrt_full_nb2
+.global rsqrt_full_nb2
+.global _rsqrt_full_nb_gpr
+.global rsqrt_full_nb_gpr
+.global _rsqrt_newton
+.global rsqrt_newton
+.global _rsqrt_hack
+.global rsqrt_hack
+.global _rsqrt_fallback
+
+.text
+.intel_syntax
+
+.align 16
+min_pos_denorm:
+.long 0x00800000,0,0,0
+penultimate_bit:
+.long 0x00008000,0,0,0
+ultimate_bit:
+.long 0x00004000,0,0,0
+top_mask:
+.long 0xFFFF8000,0,0,0
+one:
+.long 0x3f800000,0,0,0
+half:
+.long 0x3f000000,0,0,0
+one_point_five:
+.long 0x3fc00000,0,0,0
+magic1:
+.long 0x60000000,0,0,0
+magic2:
+.long 0x3c000000,0,0,0
+magic3:
+.long 0x000047ff,0,0,0
+
+_rsqrt_inaccurate:
+rsqrt_inaccurate:
+    movd xmm0, edi
+
+    rsqrtss xmm0, xmm0
+
+    movd eax, xmm0
+    ret
+
+_rsqrt_full:
+rsqrt_full:
+    movd xmm0, edi
+
+    pand xmm0, [rip + top_mask]
+    por xmm0, [rip + penultimate_bit]
+
+    vcmpngt_uqss xmm1, xmm0, [rip + min_pos_denorm]
+    ptest xmm1, xmm1
+    jnz rsqrt_full_bad
+
+    sqrtss xmm0, xmm0
+
+    movd xmm1, [rip + one]
+    divss xmm1, xmm0
+
+    paddd xmm1, [rip + ultimate_bit]
+    pand xmm1, [rip + top_mask]
+
+    movd eax, xmm1
+    ret
+
+_rsqrt_full_gpr:
+rsqrt_full_gpr:
+    movd eax, xmm0 # Emulate regalloc mov
+
+    mov eax, edi
+    and eax, 0xFFFF8000
+    or eax, 0x00008000
+
+    movd xmm0, eax
+    vcmpngt_uqss xmm1, xmm0, [rip + min_pos_denorm]
+    ptest xmm1, xmm1
+    jnz rsqrt_full_bad
+
+    sqrtss xmm0, xmm0
+
+    movd xmm1, [rip + one]
+    divss xmm1, xmm0
+    movd eax, xmm1
+
+    add eax, 0x00004000
+    and eax, 0xffff8000
+
+    movd xmm0, eax # Emulate regalloc mov
+    ret
+
+_rsqrt_full_nb2:
+rsqrt_full_nb2:
+    movd xmm0, edi
+
+    pand xmm0, [rip + top_mask]
+    por xmm0, [rip + penultimate_bit]
+
+    ucomiss xmm0, [rip + min_pos_denorm]
+    jna rsqrt_full_bad_new1
+
+    sqrtss xmm0, xmm0
+
+    movd xmm1, [rip + one]
+    divss xmm1, xmm0
+
+    paddd xmm1, [rip + ultimate_bit]
+    pand xmm1, [rip + top_mask]
+
+    movd eax, xmm1
+    ret
+
+_rsqrt_full_nb:
+rsqrt_full_nb:
+    movd xmm0, edi
+
+    pand xmm0, [rip + top_mask]
+    por xmm0, [rip + penultimate_bit]
+
+    vcmpngt_uqss xmm1, xmm0, [rip + min_pos_denorm]
+    ptest xmm1, xmm1
+    jnz rsqrt_full_bad_new1
+
+    sqrtss xmm0, xmm0
+
+    movd xmm1, [rip + one]
+    divss xmm1, xmm0
+
+    paddd xmm1, [rip + ultimate_bit]
+    pand xmm1, [rip + top_mask]
+
+    movd eax, xmm1
+    ret
+
+rsqrt_full_bad_new1:
+    cmp edi, 0x00800000
+    jb rsqrt_full_bad_new_fallback1
+
+    movd xmm0, edi
+    rsqrtss xmm1, xmm0
+
+    ucomiss xmm1, xmm1
+    jp rsqrt_full_bad_new1_nan
+
+    movd eax, xmm1
+    ret
+
+rsqrt_full_bad_new_fallback1:
+    call _rsqrt_fallback
+    ret
+
+rsqrt_full_bad_new1_nan:
+    ucomiss xmm0, xmm0
+    jp rsqrt_full_bad_new1_nan_ret
+
+    mov eax, 0x7FC00000
+    ret
+
+rsqrt_full_bad_new1_nan_ret:
+    ret
+
+_rsqrt_full_nb_gpr:
+rsqrt_full_nb_gpr:
+    movd eax, xmm0 # Emulate regalloc mov
+
+    mov eax, edi
+    and eax, 0xFFFF8000
+    or eax, 0x00008000
+
+    movd xmm0, eax
+    vcmpngt_uqss xmm1, xmm0, [rip + min_pos_denorm]
+    ptest xmm1, xmm1
+    jnz rsqrt_full_bad_new2
+
+    sqrtss xmm0, xmm0
+
+    movd xmm1, [rip + one]
+    divss xmm1, xmm0
+    movd eax, xmm1
+
+    add eax, 0x00004000
+    and eax, 0xffff8000
+
+    movd xmm0, eax # Emulate regalloc mov
+    ret
+
+rsqrt_full_bad_new2:
+    cmp edi, 0x00800000
+    jb rsqrt_full_bad_new_fallback2
+
+    movd xmm0, edi
+    rsqrtss xmm1, xmm0
+
+    test edi, edi
+    js rsqrt_full_bad_new2_nan
+
+    movd eax, xmm1
+    ret
+
+rsqrt_full_bad_new_fallback2:
+    call _rsqrt_fallback
+    ret
+
+rsqrt_full_bad_new2_nan:
+    mov eax, 0x7FC00000
+    ret
+
+rsqrt_full_bad:
+    xorps xmm1, xmm1
+    movd xmm0, edi
+    ucomiss xmm0, xmm1
+    jp rsqrt_full_nan
+    je rsqrt_full_zero
+    jc rsqrt_full_neg
+
+    cmp edi, 0x7F800000
+    je rsqrt_full_inf
+
+    # TODO: Full Denormal Implementation
+    call _rsqrt_fallback
+    ret
+
+rsqrt_full_neg:
+    mov eax, 0x7FC00000
+    ret
+
+rsqrt_full_inf:
+    xor eax, eax
+    ret
+
+rsqrt_full_nan:
+    mov eax, edi
+    or eax, 0x00400000
+    ret
+
+rsqrt_full_zero:
+    mov eax, edi
+    or eax, 0x7F800000
+    ret
+
+_rsqrt_newton:
+rsqrt_newton:
+    movd xmm0, edi
+
+    pand xmm0, [rip + top_mask]
+    por xmm0, [rip + penultimate_bit]
+
+    vcmpngt_uqss xmm1, xmm0, [rip + min_pos_denorm]
+    ptest xmm1, xmm1
+    jnz rsqrt_full_bad
+
+    rsqrtps xmm1, xmm0
+    mulss xmm0, [rip + half]
+    vmulss xmm2, xmm1, xmm1
+    mulss xmm2, xmm0
+    movaps xmm0, [rip + one_point_five]
+    subss xmm0, xmm2
+    mulss xmm0, xmm1
+
+    paddd xmm0, [rip + ultimate_bit]
+    pand xmm0, [rip + top_mask]
+
+    movd eax, xmm0
+    ret
+
+_rsqrt_hack:
+rsqrt_hack:
+    movd xmm9, edi
+
+    vpand xmm0, xmm9, [rip + top_mask]
+    por xmm0, [rip + penultimate_bit]
+
+    # detect NaNs, negatives, zeros, denormals and infinities
+    vcmpngt_uqss xmm1, xmm0, [rip + min_pos_denorm]
+    ptest xmm1, xmm1
+    jnz rsqrt_full_bad
+
+    # calculate x64 estimate
+    rsqrtps xmm0, xmm0
+
+    # calculate correction factor
+    vpslld xmm1, xmm9, 8
+    vpsrad xmm2, xmm1, 31
+    paddd xmm1, [rip + magic1]
+    pcmpgtd xmm1, [rip + magic2]
+    pxor xmm1, xmm2
+    movaps xmm2, [rip + magic3]
+    psubd xmm2, xmm1
+
+    # correct x64 estimate
+    paddd xmm0, xmm2
+    pand xmm0, [rip + top_mask]
+
+    movd eax, xmm0
+    ret