Add legacy Carmack implementation of Reciprocal Sqrt

include/J3ML/J3ML.hpp

@@ -441,13 +441,23 @@ namespace J3ML::Math::Functions {
     /// Returns the fractional part of x.
     /** @see Lerp(), LerpMod(), InvLerp(), Step(), SmoothStep(), PingPongMod(), Mod(), ModPos(). */
     float Frac(float x);
-    float Sqrt(float x); /// Returns the square root of x.
-    float FastSqrt(float x); /// Computes a fast approximation of the square root of x.
-    float RSqrt(float x); /// Returns 1/Sqrt(x). (The reciprocal of the square root of x)
-    float FastRSqrt(float x); /// SSE implementation of reciprocal square root.
-    float Recip(float x); /// Returns 1/x, the reciprocal of x.
-    float RecipFast(float x); /// Returns 1/x, the reciprocal of x, using a fast approximation (SSE rcp instruction).
-
+    /// Returns the square root of x.
+    float Sqrt(float x);
+    /// Computes a fast approximation of the square root of x.
+    float FastSqrt(float x);
+    /// Returns 1/Sqrt(x). (The reciprocal of the square root of x)
+    float RSqrt(float x);
+    /// SSE implementation of reciprocal square root.
+    float FastRSqrt(float x);
+    /// Returns 1/x, the reciprocal of x.
+    float Recip(float x);
+    /// Returns 1/x, the reciprocal of x, using a fast approximation (SSE rcp instruction).
+    float RecipFast(float x);
+    /// Carmack (Quake) implementation of inverse (reciprocal) square root.
+    /// Relies on funky type-casting hacks to avoid floating-point division and sqrt, which is very slow on legacy hardware.
+    /// This technique is superseded by modern processors having built-in support, but is included for its historical significance.
+    /// https://en.wikipedia.org/wiki/Fast_inverse_square_root
+    float QRSqrt(float x);
 }
 
 namespace J3ML::Math::Functions::Interpolation

src/J3ML/J3ML.cpp

@@ -235,6 +235,20 @@ namespace J3ML {
         return 1.f / x;
     }
 
+    float Math::Functions::QRSqrt(float x) {
+        long i;
+        float x2, y;
+        const float threehalfs = 1.5f;
+        x2 = x * 0.5f;
+        y = x;
+        i = *(long*) &y; // evil floating point bit level hacking
+        i = 0x5f3759df - (i >> 1); // what the fuck?
+        y = *(float*) &i;
+        y = y * (threehalfs - (x2 * y * y)); // increase precision of approximation via Newton's Method
+
+        return y;
+    }
+
     float Math::Functions::Lerp(float a, float b, float t) { return a + t * (b-a);}
 
     float Math::Functions::LerpMod(float a, float b, float mod, float t) {