Fast inverse square root (slower by 10%)

Author: knapeczadam

source/Macros.h

@@ -0,0 +1,3 @@
+#pragma once
+
+#define QRSQRT 0

source/MathHelpers.cpp

@@ -0,0 +1,20 @@
+#include "MathHelpers.h"
+
+/**
+ * \brief https://www.youtube.com/watch?v=p8u_k2LIZyo
+ */
+float dae::Q_rsqrt(float number)
+{
+    long i;
+    float x2, y;
+    const float threehalfs = 1.5F;
+
+    x2 = number * 0.5F;
+    y = number;
+    i = *reinterpret_cast<long*>(&y); // evil floating point bit level hacking
+    i = 0x5f3759df - (i >> 1);
+    y = *reinterpret_cast<float*>(&i);
+    y = y * (threehalfs - (x2 * y * y)); // 1st iteration
+    y = y * (threehalfs - (x2 * y * y)); // 2nd iteration, this can be removed
+    return y;
+}

source/MathHelpers.h

@@ -19,6 +19,8 @@ namespace dae
         return a * a;
     }
 
+    float Q_rsqrt(float number);
+
     inline float Lerpf(float a, float b, float factor)
     {
         return ((1 - factor) * a) + (factor * b);

source/RayTracer.vcxproj

@@ -97,6 +97,7 @@ xcopy "$(SolutionDir)..\lib\vld\x64\Microsoft.DTfW.DHL.manifest" "$(OutDir)" /y
     <ClInclude Include="Camera.h" />
     <ClInclude Include="ColorRGB.h" />
     <ClInclude Include="DataTypes.h" />
+    <ClInclude Include="Macros.h" />
     <ClInclude Include="Material.h" />
     <ClInclude Include="MathHelpers.h" />
     <ClInclude Include="Matrix.h" />
@@ -110,6 +111,7 @@ xcopy "$(SolutionDir)..\lib\vld\x64\Microsoft.DTfW.DHL.manifest" "$(OutDir)" /y
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="Camera.cpp" />
+    <ClCompile Include="MathHelpers.cpp" />
     <ClCompile Include="Matrix.cpp" />
     <ClCompile Include="Renderer.cpp" />
     <ClCompile Include="Scene.cpp" />

source/RayTracer.vcxproj.filters

@@ -53,6 +53,7 @@
     <ClInclude Include="DataTypes.h">
       <Filter>Misc</Filter>
     </ClInclude>
+    <ClInclude Include="Macros.h" />
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="main.cpp" />
@@ -75,5 +76,8 @@
     <ClCompile Include="Camera.cpp">
       <Filter>Misc</Filter>
     </ClCompile>
+    <ClCompile Include="MathHelpers.cpp">
+      <Filter>Math</Filter>
+    </ClCompile>
   </ItemGroup>
 </Project>

source/Vector3.cpp

@@ -1,9 +1,11 @@
 #include "Vector3.h"
 
-#include <cassert>
-
 #include "Vector4.h"
+#include "MathHelpers.h"
+#include "Macros.h"
+
 #include <cmath>
+#include <cassert>
 
 namespace dae
 {
@@ -26,7 +28,12 @@ namespace dae
 
     float Vector3::Magnitude() const
     {
-        return sqrtf(x * x + y * y + z * z);
+#if QRSQRT
+        const float magnitudeSq {x * x + y * y + z * z};
+        return 1.0f / Q_rsqrt(magnitudeSq);
+#else
+        return std::sqrt(x * x + y * y + z * z);
+#endif
     }
 
     float Vector3::SqrMagnitude() const
@@ -36,18 +43,30 @@ namespace dae
 
     float Vector3::Normalize()
     {
-        const float m = Magnitude();
+#if QRSQRT
+        const float m {Q_rsqrt(x * x + y * y + z * z)};
+        x *= m;
+        y *= m;
+        z *= m;
+        return m;
+#else
+        const float m{Magnitude()};
         x /= m;
         y /= m;
         z /= m;
-
         return m;
+#endif
     }
 
     Vector3 Vector3::Normalized() const
     {
-        const float m = Magnitude();
+#if QRSQRT
+        const float m {Q_rsqrt(x * x + y * y + z * z)};
+        return {x * m, y * m, z * m};
+#else
+        const float m{Magnitude()};
         return {x / m, y / m, z / m};
+#endif
     }
 
     float Vector3::Dot(const Vector3& v1, const Vector3& v2)

source/Vector4.cpp

@@ -1,9 +1,11 @@
 #include "Vector4.h"
 
-#include <cassert>
-
 #include "Vector3.h"
+#include "MathHelpers.h"
+#include "Macros.h"
+
 #include <cmath>
+#include <cassert>
 
 namespace dae
 {
@@ -17,7 +19,12 @@ namespace dae
 
     float Vector4::Magnitude() const
     {
-        return sqrtf(x * x + y * y + z * z + w * w);
+#if QRSQRT
+        const float magnitudeSq {x * x + y * y + z * z + w * w};
+        return 1.0f / Q_rsqrt(magnitudeSq);
+#else
+        return std::sqrt(x * x + y * y + z * z + w * w);
+#endif
     }
 
     float Vector4::SqrMagnitude() const
@@ -27,19 +34,32 @@ namespace dae
 
     float Vector4::Normalize()
     {
-        const float m = Magnitude();
+#if QRSQRT
+        const float m {Q_rsqrt(x * x + y * y + z * z + w * w)};
+        x *= m;
+        y *= m;
+        z *= m;
+        w *= m;
+        return m;
+#else
+        const float m{Magnitude()};
         x /= m;
         y /= m;
         z /= m;
         w /= m;
-
         return m;
+#endif
     }
 
     Vector4 Vector4::Normalized() const
     {
-        const float m = Magnitude();
+#if QRSQRT
+        const float m {Q_rsqrt(x * x + y * y + z * z + w * w)};
+        return {x * m, y * m, z * m, w * m};
+#else
+        const float m{Magnitude()};
         return {x / m, y / m, z / m, w / m};
+#endif
     }
 
     float Vector4::Dot(const Vector4& v1, const Vector4& v2)