Switches to using precomputed reciprocal math.

See issue #337 #337 for
more info.

Author: louis-langholtz

PlayRho/Collision/AABB.cpp

@@ -45,7 +45,7 @@ AABB ComputeAABB(const DistanceProxy& proxy, const Transformation& xf) noexcept
 }
 
 AABB ComputeAABB(const DistanceProxy& proxy,
-                   const Transformation& xfm0, const Transformation& xfm1) noexcept
+                 const Transformation& xfm0, const Transformation& xfm1) noexcept
 {
     assert(IsValid(xfm0));
     assert(IsValid(xfm1));

PlayRho/Collision/AABB.hpp

@@ -236,7 +236,8 @@ PLAYRHO_CONSTEXPR inline Vector<Length, N> GetDimensions(const AABB<N>& aabb) no
 template <std::size_t N>
 PLAYRHO_CONSTEXPR inline Vector<Length, N> GetExtents(const AABB<N>& aabb) noexcept
 {
-    return GetDimensions(aabb) / 2;
+    PLAYRHO_CONSTEXPR const auto RealInverseOfTwo = Real{1} / Real{2};
+    return GetDimensions(aabb) * RealInverseOfTwo;
 }
 
 /// @brief Checks whether the first AABB fully contains the second AABB.

PlayRho/Collision/MassData.cpp

@@ -145,16 +145,19 @@ MassData GetMassData(Length vertexRadius, NonNegative<AreaDensity> density,
 
         const auto D = Cross(e1, e2);
 
-        const auto triangleArea = D / Real{2};
+        PLAYRHO_CONSTEXPR const auto RealReciprocalOfTwo = Real{1} / Real{2}; // .5
+        const auto triangleArea = D * RealReciprocalOfTwo;
         area += triangleArea;
 
         // Area weighted centroid
-        center += StripUnit(triangleArea) * (e1 + e2) / Real{3};
+        PLAYRHO_CONSTEXPR const auto RealReciprocalOfThree = Real{1} / Real{3}; // .3333333...
+        center += StripUnit(triangleArea) * (e1 + e2) * RealReciprocalOfThree;
 
         const auto intx2 = Square(GetX(e1)) + GetX(e2) * GetX(e1) + Square(GetX(e2));
         const auto inty2 = Square(GetY(e1)) + GetY(e2) * GetY(e1) + Square(GetY(e2));
 
-        const auto triangleI = D * (intx2 + inty2) / Real{3 * 4};
+        PLAYRHO_CONSTEXPR const auto RealReciprocalOfTwelve = Real{1} / Real{3 * 4}; // .083333..
+        const auto triangleI = D * (intx2 + inty2) * RealReciprocalOfTwelve;
         I += triangleI;
     }
 

PlayRho/Collision/RayCastOutput.cpp

@@ -93,8 +93,9 @@ RayCastOutput RayCast(const AABB& aabb, const RayCastInput& input) noexcept
         }
         else
         {
-            auto t1 = Real{(range.GetMin() - p1i) / pdi};
-            auto t2 = Real{(range.GetMax() - p1i) / pdi};
+            const auto reciprocalOfPdi = Real{1} / pdi;
+            auto t1 = Real{(range.GetMin() - p1i) * reciprocalOfPdi};
+            auto t2 = Real{(range.GetMax() - p1i) * reciprocalOfPdi};
             auto s = -1; // Sign of the normal vector.
             if (t1 > t2)
             {
@@ -187,11 +188,13 @@ RayCastOutput RayCast(const DistanceProxy& proxy, const RayCastInput& input,
             const auto v0off = v0 + offset;
             const auto q_sub_p = v0off - ray0;
 
+            const auto reciprocalRayCrossEdge = Real{1} / ray_cross_edge;
+
             // t = ((q - p) x s) / (r x s)
-            const auto t = Cross(q_sub_p, edge) / ray_cross_edge;
+            const auto t = Cross(q_sub_p, edge) * reciprocalRayCrossEdge;
 
             // u = ((q - p) x r) / (r x s)
-            const auto u = Cross(q_sub_p, ray) / ray_cross_edge;
+            const auto u = Cross(q_sub_p, ray) * reciprocalRayCrossEdge;
 
             if ((t >= 0) && (t <= 1) && (u >= 0) && (u <= 1))
             {

PlayRho/Collision/Simplex.cpp

@@ -111,8 +111,8 @@ Simplex Simplex::Get(const SimplexEdge& s0, const SimplexEdge& s1) noexcept
     }
 
     // Must be in e12 region.
-    const auto d12_sum = d12_1 + d12_2;
-    return Simplex{{s0, s1}, {d12_1 / d12_sum, d12_2 / d12_sum}};
+    const auto inv_sum = Real{1} / (d12_1 + d12_2);
+    return Simplex{{s0, s1}, {d12_1 * inv_sum, d12_2 * inv_sum}};
 }
 
 Simplex Simplex::Get(const SimplexEdge& s0, const SimplexEdge& s1, const SimplexEdge& s2) noexcept
@@ -179,31 +179,31 @@ Simplex Simplex::Get(const SimplexEdge& s0, const SimplexEdge& s1, const Simplex
     const auto d123_3 = n123 * cp_w1_w2;
     if ((d12_1 > 0_m2) && (d12_2 > 0_m2) && (d123_3 <= 0 * SquareMeter * SquareMeter))
     {
-        const auto d12_sum = d12_1 + d12_2;
-        return Simplex{{s0, s1}, {d12_1 / d12_sum, d12_2 / d12_sum}};
+        const auto inv_sum = Real{1} / (d12_1 + d12_2);
+        return Simplex{{s0, s1}, {d12_1 * inv_sum, d12_2 * inv_sum}};
     }
 
     // e13
     const auto cp_w3_w1 = Cross(w3, w1);
     const auto d123_2 = n123 * cp_w3_w1;
     if ((d13_1 > 0_m2) && (d13_2 > 0_m2) && (d123_2 <= 0 * SquareMeter * SquareMeter))
     {
-        const auto d13_sum = d13_1 + d13_2;
-        return Simplex{{s0, s2}, {d13_1 / d13_sum, d13_2 / d13_sum}};
+        const auto inv_sum = Real{1} / (d13_1 + d13_2);
+        return Simplex{{s0, s2}, {d13_1 * inv_sum, d13_2 * inv_sum}};
     }
 
     // e23
     const auto cp_w2_w3 = Cross(w2, w3);
     const auto d123_1 = n123 * cp_w2_w3;
     if ((d23_1 > 0_m2) && (d23_2 > 0_m2) && (d123_1 <= 0 * SquareMeter * SquareMeter))
     {
-        const auto d23_sum = d23_1 + d23_2;
-        return Simplex{{s2, s1}, {d23_2 / d23_sum, d23_1 / d23_sum}};
+        const auto inv_sum = Real{1} / (d23_1 + d23_2);
+        return Simplex{{s2, s1}, {d23_2 * inv_sum, d23_1 * inv_sum}};
     }
 
     // Must be in triangle123
-    const auto d123_sum = d123_1 + d123_2 + d123_3;
-    return Simplex{{s0, s1, s2}, {d123_1 / d123_sum, d123_2 / d123_sum, d123_3 / d123_sum}};
+    const auto inv_sum = Real{1} / (d123_1 + d123_2 + d123_3);
+    return Simplex{{s0, s1, s2}, {d123_1 * inv_sum, d123_2 * inv_sum, d123_3 * inv_sum}};
 }
 
 Simplex Simplex::Get(const SimplexEdges& edges) noexcept

PlayRho/Common/Acceleration.hpp

@@ -132,17 +132,8 @@ namespace d2 {
     /// @relatedalso Acceleration
     PLAYRHO_CONSTEXPR inline Acceleration operator/ (const Acceleration& lhs, const Real rhs)
     {
-        /*
-         * While it can be argued that division operations shouldn't be supported due to
-         * hardware support for division typically being significantly slower than hardware
-         * support for multiplication, it can also be argued that it shouldn't be the
-         * software developer's role to attempt to optimize what the compiler should be
-         * much better at knowing how to optimize. So here the code chooses the latter
-         * strategy which allows the intention to be clearer, and just passes the division
-         * on down to the Vec2 and Angle types (rather than manually rewriting the divisions
-         * as multiplications).
-         */
-        return Acceleration{lhs.linear / rhs, lhs.angular / rhs};
+        const auto inverseRhs = Real{1} / rhs;
+        return Acceleration{lhs.linear * inverseRhs, lhs.angular * inverseRhs};
     }
 
 } // namespace d2

PlayRho/Common/Math.cpp

@@ -62,11 +62,13 @@ Length2 ComputeCentroid(const Span<const Length2>& vertices)
         const auto e1 = p2 - p1;
         const auto e2 = p3 - p1;
 
-        const auto triangleArea = Area{Cross(e1, e2) / Real(2)};
+        PLAYRHO_CONSTEXPR const auto RealInverseOfTwo = Real{1} / Real{2};
+        const auto triangleArea = Area{Cross(e1, e2) * RealInverseOfTwo};
         area += triangleArea;
 
         // Area weighted centroid
-        const auto aveP = (p1 + p2 + p3) / Real{3};
+        PLAYRHO_CONSTEXPR const auto RealInverseOfThree = Real{1} / Real{3};
+        const auto aveP = (p1 + p2 + p3) * RealInverseOfThree;
         c += triangleArea * aveP;
     }
 
@@ -129,7 +131,8 @@ NonNegative<Area> GetAreaOfPolygon(Span<const Length2> vertices)
 
     // Note that using the absolute value isn't necessary for vertices in counter-clockwise
     // ordering; only needed for clockwise ordering.
-    return abs(sum) / Real{2};
+    PLAYRHO_CONSTEXPR const auto RealInverseOfTwo = Real{1} / Real{2};
+    return abs(sum) * RealInverseOfTwo;
 }
 
 SecondMomentOfArea GetPolarMoment(Span<const Length2> vertices)
@@ -161,7 +164,8 @@ SecondMomentOfArea GetPolarMoment(Span<const Length2> vertices)
     }
     const auto secondMomentOfAreaX = SecondMomentOfArea{sum_x};
     const auto secondMomentOfAreaY = SecondMomentOfArea{sum_y};
-    return (secondMomentOfAreaX + secondMomentOfAreaY) / Real{12};
+    PLAYRHO_CONSTEXPR const auto RealInverseOfTwelve = Real{1} / Real{12};
+    return (secondMomentOfAreaX + secondMomentOfAreaY) * RealInverseOfTwelve;
 }
 
 namespace d2 {

PlayRho/Common/Math.hpp

@@ -451,11 +451,12 @@ template <typename T, typename U>
 PLAYRHO_CONSTEXPR inline auto Solve(const Matrix22<U> mat, const Vector2<T> b) noexcept
 {
     const auto cp = Cross(get<0>(mat), get<1>(mat));
+    const auto inverseCp = Real{1} / cp;
     using OutType = decltype((U{} * T{}) / cp);
     return (!AlmostZero(StripUnit(cp)))?
         Vector2<OutType>{
-            (get<1>(mat)[1] * b[0] - get<1>(mat)[0] * b[1]) / cp,
-            (get<0>(mat)[0] * b[1] - get<0>(mat)[1] * b[0]) / cp
+            (get<1>(mat)[1] * b[0] - get<1>(mat)[0] * b[1]) * inverseCp,
+            (get<0>(mat)[0] * b[1] - get<0>(mat)[1] * b[0]) * inverseCp
         }: Vector2<OutType>{};
 }
 
@@ -465,10 +466,11 @@ PLAYRHO_CONSTEXPR inline auto Invert(const Matrix22<IN_TYPE> value) noexcept
 {
     const auto cp = Cross(get<0>(value), get<1>(value));
     using OutType = decltype(get<0>(value)[0] / cp);
+    const auto inverseCp = Real{1} / cp;
     return (!AlmostZero(StripUnit(cp)))?
         Matrix22<OutType>{
-            Vector2<OutType>{ get<1>(get<1>(value)) / cp, -get<1>(get<0>(value)) / cp},
-            Vector2<OutType>{-get<0>(get<1>(value)) / cp,  get<0>(get<0>(value)) / cp}
+            Vector2<OutType>{ get<1>(get<1>(value)) * inverseCp, -get<1>(get<0>(value)) * inverseCp},
+            Vector2<OutType>{-get<0>(get<1>(value)) * inverseCp,  get<0>(get<0>(value)) * inverseCp}
         }:
         Matrix22<OutType>{};
 }
@@ -478,7 +480,7 @@ PLAYRHO_CONSTEXPR inline auto Invert(const Matrix22<IN_TYPE> value) noexcept
 PLAYRHO_CONSTEXPR inline Vec3 Solve33(const Mat33& mat, const Vec3 b) noexcept
 {
     const auto dp = Dot(GetX(mat), Cross(GetY(mat), GetZ(mat)));
-    const auto det = (dp != 0)? 1 / dp: dp;
+    const auto det = (dp != 0)? Real{1} / dp: dp;
     const auto x = det * Dot(b, Cross(GetY(mat), GetZ(mat)));
     const auto y = det * Dot(GetX(mat), Cross(b, GetZ(mat)));
     const auto z = det * Dot(GetX(mat), Cross(GetY(mat), b));
@@ -492,7 +494,7 @@ template <typename T>
 PLAYRHO_CONSTEXPR inline T Solve22(const Mat33& mat, const T b) noexcept
 {
     const auto cp = GetX(GetX(mat)) * GetY(GetY(mat)) - GetX(GetY(mat)) * GetY(GetX(mat));
-    const auto det = (cp != 0)? 1 / cp: cp;
+    const auto det = (cp != 0)? Real{1} / cp: cp;
     const auto x = det * (GetY(GetY(mat)) * GetX(b) - GetX(GetY(mat)) * GetY(b));
     const auto y = det * (GetX(GetX(mat)) * GetY(b) - GetY(GetX(mat)) * GetX(b));
     return T{x, y};
@@ -623,7 +625,7 @@ inline Real Normalize(Vec2& vector)
     const auto length = GetMagnitude(vector);
     if (!AlmostZero(length))
     {
-        const auto invLength = 1 / length;
+        const auto invLength = Real{1} / length;
         vector[0] *= invLength;
         vector[1] *= invLength;
         return length;
@@ -734,7 +736,8 @@ PLAYRHO_CONSTEXPR inline Vector2<T> operator* (const UnitVec u, const T s) noexc
 /// @brief Division operator.
 PLAYRHO_CONSTEXPR inline Vec2 operator/ (const UnitVec u, const UnitVec::value_type s) noexcept
 {
-    return Vec2{GetX(u) / s, GetY(u) / s};
+    const auto inverseS = Real{1} / s;
+    return Vec2{GetX(u) * inverseS, GetY(u) * inverseS};
 }
 
 /// @brief Rotates a vector by a given angle.