Further implementation of core math functions.

include/J3ML/J3ML.h

@@ -17,6 +17,10 @@
 #include <cassert>
 #include <vector>
 
+
+#include <J3ML/Algorithm/Reinterpret.hpp>
+
+
 /// Swaps two elements in-place without copying their data.
 template <typename T>
 void Swap(T &a, T &b)
@@ -42,6 +46,8 @@ namespace J3ML::SizedIntegralTypes
 
 namespace J3ML::SizedFloatTypes
 {
+    // TODO: Use C++23 <stdfloat>
+    using f16  = float;
     using f32 = float;
     using f64 = double;
     using f128 = long double;
@@ -85,7 +91,8 @@ namespace J3ML::Math::Constants {
     /// ln 10
     constexpr float NaturalLog10 = 2.3025850929940456840179914546843642076011014886287729760333279009675726096773524802359972050895982983;
     constexpr float Infinity = INFINITY;
-    constexpr float NegativeInfinity = INFINITY;
+    constexpr float NegativeInfinity = -INFINITY;
+    constexpr float NotANumber = NAN;
 }
 
 /// This set of functions may be set to use lookup tables or SIMD operations.
@@ -97,6 +104,112 @@ namespace J3ML::Math::Constants {
 
 namespace J3ML::Math::Functions {
 
+    /// Clamps the given input value to the range [min, max].
+    /** @see Clamp01(), Min(), Max(). */
+    template<typename T>
+    T Clamp(const T &val, const T &floor, const T &ceil)
+    {
+        assert(floor <= ceil);
+        return val <= ceil ? (val >= floor ? val : floor) : ceil;
+    }
+    /// Clamps the given input value to the range [0, 1].
+    /** @see Clamp(), Min(), Max(). */
+    template<typename T>
+    T Clamp01(const T &val) { return Clamp(val, T(0), T(1)); }
+
+    /// Computes the smaller of the two values.
+    /** @see Clamp(), Clamp01(), Max() */
+    template <typename T>
+    T Min(const T& a, const T& b) {
+        return a <= b ? a : b;
+    }
+
+    /// Computes the larger of two values.
+    /** @see Clamp(), Clamp01(), Max() */
+    template <typename T>
+    T Max(const T& a, const T& b) {
+        return a >= b ? a : b;
+    }
+
+    /// Computes the smallest in an arbitrary list of values.
+    /** @see Clamp(), Clamp01(), Max() */
+    template <typename T>
+    T Min(const std::initializer_list<T>& list) {
+        T minimum = list[0];
+
+        for (T entry : list) {
+            if (entry <= minimum)
+                minimum = entry;
+        }
+        return minimum;
+    }
+
+    /// Computes the largest in an arbitrary list of values.
+    /** @see Clamp(), Clamp01(), Max() */
+    template <typename T>
+    T Max(const std::initializer_list<T>& list) {
+        T maximum = list[0];
+        for (T entry : list) {
+            if (entry >= maximum)
+                maximum = entry;
+        }
+        return maximum;
+    }
+
+    /** @return True if a > b. */
+    template <typename T>
+    bool GreaterThan(const T& a, const T& b) {
+        return a > b;
+    }
+
+    /** @return True if a < b. */
+    template <typename T>
+    bool LessThan(const T& a, const T& b) {
+        return a < b;
+    }
+
+    /** @return The absolute value of a. */
+    template <typename T>
+    T Abs(const T& a) {
+        return a >= 0 ? a : -a;
+    }
+
+    template<> inline float Abs(const float& a) {
+#ifdef USE_SSE
+#else
+        return a >= 0 ? a : -a;
+#endif
+    }
+
+    /// @return True if a and b are equal, using operator ==().
+    template <typename T>
+    bool EqualExact(const T& a, const T& b) {
+        return a == b;
+    }
+
+    /** Compares the two values for equality up to a small epsilon. */
+    bool Equal(float a, float b, float epsilon = 1e-3f);
+
+    /** Compares the two values for equality up to a small epsilon. */
+    bool Equal(double a, double b, float epsilon = 1e-3f);
+
+    /** Compares the two values for equality, allowing the given amount of absolute error. */
+    bool EqualAbs(float a, float b, float epsilon = 1e-3f);
+
+    /// Computes the relative error of the two variables.
+    float RelativeError(float a, float b);
+
+    inline bool IsFinite(float f) { return (ReinterpretAs<u32>(f) << 1) < 0xFF000000u; }
+    inline bool IsFinite(double d) { return (ReinterpretAs<u64>(d) << 1) < 0xFFE0000000000000ULL; }
+
+    inline bool IsNotANumber(float f) { return (ReinterpretAs<u32>(f) << 1) > 0xFF000000u; }
+    inline bool IsNotANumber(double d) { return (ReinterpretAs<u64>(d) << 1) > 0xFFE0000000000000ULL; }
+
+    inline bool IsInfinite(float f) { return (ReinterpretAs<u32>(f) << 1) == 0xFF000000u; }
+    inline bool IsInfinite(double d) { return (ReinterpretAs<u64>(d) << 1) == 0xFFE0000000000000ULL; }
+
+
+
     float Radians(float deg); /// Converts the given amount of degrees into radians.
     float Degrees(float rad); /// Converts the given amount of radians into degrees.
 
@@ -128,34 +241,29 @@ namespace J3ML::Math::Functions {
     float Ln(float value); /// Computes a logarithm in the natural base (using e as the base).
     float Log10(float value); /// Computes a logarithm in base-10;
 
-    float Ceil(float f);
-    int CeilInt(float f);
-    float Floor(float f);
-    int FloorInt(float f);
-    float Round(float f);
-    int RoundInt(float f);
+    float Ceil(float f);    /// Returns f rounded up to the next integer, as float.
+    int CeilInt(float f);   /// Returns f rounded up to the next integer, as integer.
+    float Floor(float f);   /// Returns f rounded down to the previous integer, as float.
+    int FloorInt(float f);  /// Returns f rounded down to the previous integer, as integer.
+    float Round(float f);   /// Returns f rounded to the nearest integer, as float.
+    int RoundInt(float f);  /// Returns f rounded to the nearest integer, as int.
 
-    float Round(float f, float decimalPlaces);
 
-    float Sign(float f);
+    float Round(float f, float decimalPlaces); /// Returns f rounded to the given decimal places.
+
+
+    float Sign(float f); /// Returns -1 or 1 depending on the sign of f.
+    ///
+    float SignOrZero(float f, float epsilon = 1e-8f);
 
     /// Formats larger numbers into shortened 'Truncated' string representations.
     /// 2241 -> 2.2k, 55421 -> 55.4k, 1000000 -> 1.0M
     std::string Truncate(float input);
 
-    /// Clamps the given input value to the range [min, max].
-    /** @see Clamp01(), Min(), Max(). */
-    template<typename T>
-    T Clamp(const T &val, const T &floor, const T &ceil)
-    {
-        assert(floor <= ceil);
-        return val <= ceil ? (val >= floor ? val : floor) : ceil;
-    }
-    /// Clamps the given input value to the range [0, 1].
-    /** @see Clamp(), Min(), Max(). */
-    template<typename T>
-    T Clamp01(const T &val) { return Clamp(val, T(0), T(1)); }
-    bool EqualAbs(float a, float b, float epsilon = 1e-3f);
+
+
+
+
     float RecipFast(float x);
 
 
@@ -170,29 +278,53 @@ namespace J3ML::Math::Functions {
     float NormalizeToRange(float input, const NumberRange& from, const NumberRange& to);
     // auto rotation_normalized = NormalizeToRange(inp, {0, 360}, {-1, 1});
 
+    /// Linearly interpolates between a and b.
+    /** @param t A value between [0,1].
+        @param a The first endpoint to lerp between.
+        @param b The second endpoint to lerp between.
+        @return This function computes a + t*(b-a). That is, if t==0, this function returns a. If t==1, this function returns b.
+            Otherwise, the returned value linearly moves from a to b as t ranges from 0 to 1.
+        @see LerpMod(), InvLerp(), Step(), SmoothStep(), PingPongMod(), Mod(), ModPos(), Frac(). */
     float Lerp(float a, float b, float t);
+
     /// Linearly interpolates from a to b, under the modulus mod.
-    /// This function takes evaluates a and b in the range [0, mod] and takes the shorter path to reach from a to b.
+    /** This function takes evaluates a and b in the range [0, mod] and takes the shorter path to reach from a to b.
+        @see Lerp(), InvLerp(), Step(), SmoothStep(), PingPongMod(), Mod(), ModPos(), Frac(). */
     float LerpMod(float a, float b, float mod, float t);
+
     /// Computes the lerp factor a and b have to be Lerp()ed to get x.
+    /// /** @see Lerp(), LerpMod(), Step(), SmoothStep(), PingPongMod(), Mod(), ModPos(), Frac(). */
     float InverseLerp(float a, float b, float x);
     /// See http://msdn.microsoft.com/en-us/library/bb509665(v=VS.85).aspx
     float Step(float y, float x);
     /// See http://msdn.microsoft.com/en-us/library/bb509658(v=vs.85).aspx
     float Ramp(float min, float max, float x);
-
+    /// Limits x to the range [0, mod], but instead of wrapping around from mod to 0, the result will move back
+    /// from mod to 0 as x goes from mod to 2*mod.
+    /** @see Lerp(), LerpMod(), InvLerp(), Step(), SmoothStep(), Mod(), ModPos(), Frac(). */
     float PingPongMod(float x, float mod);
 
-    float Sqrt(float x);
-
-    float FastSqrt(float x);
-
-    /// Returns 1/Sqrt(x). (The reciprocal of the square root of x)
-    float RSqrt(float x);
-
-    float FastRSqrt(float x);
-
+    /// Computes a floating-point modulus.
+    /// This function returns a value in the range ]-mod, mod[.
+    /** @see Lerp(), LerpMod(), InvLerp(), Step(), SmoothStep(), PingPongMod(), ModPos(), Frac(). */
+    float Mod(float x, float mod);
+    /// Computes a floating-point modulus using an integer as the modulus.
+    float Mod(float x, int mod);
 
+    /// Computes a floating-point modulus, but restricts the output to the range [0, mod[.
+    /** @see Lerp(), LerpMod(), InvLerp(), Step(), SmoothStep(), PingPongMod(), Mod(), Frac(). */
+    float ModPos(float x, float mod);
+    /// Computes a floating-point modulus, but restricts the output to the range [0, mod[.
+    float ModPos(float x, int mod);
+    /// 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).
 
 
     namespace Interp

src/J3ML/J3ML.cpp

@@ -57,6 +57,11 @@ namespace J3ML
         return 1.f / FastSqrt(x);
     }
 
+    float Math::Functions::PingPongMod(float x, float mod) {
+        x = Mod(x, mod * 2.f);
+        return x >= mod ? (2.f * mod - x) : x;
+    }
+
     float Math::Functions::Sqrt(float x) {
         return std::sqrt(x);
     }
@@ -119,8 +124,22 @@ namespace J3ML
 
     }
 
+    bool Math::Functions::Equal(float a, float b, float epsilon) {
+        return Abs(a-b) <= epsilon;
+    }
+
+    bool Math::Functions::Equal(double a, double b, float epsilon) {
+        return Abs(a - b) <= epsilon;
+    }
+
+    float Math::Functions::RelativeError(float a, float b) {
+        if (a == b) return 0.f; // Handles the special case where approximation and real are both zero.
+        return Abs((a-b) / Max(Abs(a), Abs(b)));
+    }
+
     bool Math::Functions::EqualAbs(float a, float b, float epsilon) {
-        return std::abs(a - b) < epsilon;
+        // TODO: No different from Equal?
+        return Abs(a - b) < epsilon;
     }
 
     float Math::Functions::RecipFast(float x) {
@@ -128,6 +147,64 @@ namespace J3ML
         return 1.f / x;
     }
 
+    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) {
+        a = ModPos(a, mod);
+        b = ModPos(b, mod);
+        if (Abs(b - a) * 2.f <= mod)
+            return Lerp(a, b, t);
+        else {
+            if (a < b)
+                return ModPos(Lerp(a + mod, b, t), mod);
+            else
+                return ModPos(Lerp(a, b + mod, t), mod);
+        }
+    }
+
+    float Math::Functions::InverseLerp(float a, float b, float x) {
+        assert(Abs(b - a) > 1e-5f);
+        return (x - a) / (b - a);
+    }
+
+
+    float Math::Functions::Step(float y, float x) {
+        return (x >= y) ? 1.f : 0.f;
+    }
+
+    float Math::Functions::Ramp(float min, float max, float x) {
+
+        return x <= min ? 0.f : (x >= max ? 1.f : (x - min) / (max - min));
+    }
+
+    float Math::Functions::Mod(float x, float mod) {
+        return std::fmod(x, mod);
+    }
+
+    float Math::Functions::Mod(float x, int mod) {
+        // TODO: Optimize?
+        return std::fmod(x,(float)mod);
+    }
+
+    float Math::Functions::ModPos(float x, float mod) {
+        float m = fmod(x, mod);
+        return m >= 0.f ? m : (m + mod);
+    }
+
+    float Math::Functions::ModPos(float x, int mod) {
+        return ModPos(x, (float)mod);
+    }
+
+    float Math::Functions::Frac(float x) {
+        return x - Floor(x);
+    }
+
+    float Math::Functions::Recip(float x) {
+        return 1.f / x;
+    }
+
+
+
     Math::Rotation::Rotation() : valueInRadians(0) {}
 
     Math::Rotation::Rotation(float value) : valueInRadians(value) {}