Foundation/Math/Common.inl

/*
 * Copyright(c) Sophist Solutions, Inc. 1990-2026.  All rights reserved
 */
#include <cmath>
#include <cstdlib>
 
#include "Stroika/Foundation/Common/Concepts.h"
#include "Stroika/Foundation/Common/StdCompat.h"
#include "Stroika/Foundation/Debug/Assertions.h"
 
namespace Stroika::Foundation::Math {
 
    /*
     ********************************************************************************
     ************************************ Math::nan *********************************
     ********************************************************************************
     */
    template <floating_point T>
    constexpr T nan ()
    {
        return numeric_limits<T>::quiet_NaN ();
    }
 
    /*
     ********************************************************************************
     ******************************* Math::infinity *********************************
     ********************************************************************************
     */
    template <floating_point T>
    constexpr T infinity ()
    {
        return numeric_limits<T>::infinity ();
    }
 
    /*
     ********************************************************************************
     ****************************** Math::RoundUpTo *********************************
     ********************************************************************************
     */
    namespace Private {
        template <unsigned_integral T>
        constexpr T RoundUpTo_UnSignedHelper_ (T x, T toNearest)
        {
            return (((x + toNearest - 1u) / toNearest) * toNearest);
        }
        template <unsigned_integral T>
        constexpr T RoundDownTo_UnSignedHelper_ (T x, T toNearest)
        {
            return ((x / toNearest) * toNearest);
        }
        template <signed_integral T>
        constexpr T RoundUpTo_SignedHelper_ (T x, T toNearest)
        {
            using UNSIGNED_T = make_unsigned_t<T>;
            Require (toNearest > 0);
            if (x < 0) {
                return (-static_cast<T> (RoundDownTo_UnSignedHelper_ (static_cast<UNSIGNED_T> (-x), static_cast<UNSIGNED_T> (toNearest))));
            }
            else {
                return static_cast<T> (RoundUpTo_UnSignedHelper_<UNSIGNED_T> (x, toNearest));
            }
        }
        template <signed_integral T>
        constexpr T RoundDownTo_SignedHelper_ (T x, T toNearest)
        {
            using UNSIGNED_T = make_unsigned_t<T>;
            Require (toNearest > 0);
            if (x < 0) {
                return (-static_cast<T> (RoundUpTo_UnSignedHelper_ (static_cast<UNSIGNED_T> (-x), static_cast<UNSIGNED_T> (toNearest))));
            }
            else {
                return (RoundDownTo_UnSignedHelper_ (static_cast<UNSIGNED_T> (x), static_cast<UNSIGNED_T> (toNearest)));
            }
        }
        template <signed_integral T>
        constexpr T RoundUpTo_ (T x, T toNearest)
        {
            return Private::RoundUpTo_SignedHelper_<T> (x, toNearest);
        }
        template <unsigned_integral T>
        constexpr T RoundUpTo_ (T x, T toNearest)
        {
            return Private::RoundUpTo_UnSignedHelper_<T> (x, toNearest);
        }
        template <signed_integral T>
        constexpr T RoundDownTo_ (T x, T toNearest)
        {
            return Private::RoundDownTo_SignedHelper_<T> (x, toNearest);
        }
        template <unsigned_integral T>
        constexpr T RoundDownTo_ (T x, T toNearest)
        {
            return Private::RoundDownTo_UnSignedHelper_<T> (x, toNearest);
        }
    }
    template <Common::IBuiltinArithmetic T>
    constexpr T RoundUpTo (T x, T toNearest)
    {
        // @todo could simplify with if constexpr in this function, and lose the helpers above
        return Private::RoundUpTo_ (x, toNearest);
    }
 
    /*
     ********************************************************************************
     **************************** Math::RoundDownTo *********************************
     ********************************************************************************
     */
    template <Common::IBuiltinArithmetic T>
    constexpr T RoundDownTo (T x, T toNearest)
    {
        // @todo could simplify with if constexpr in this function, and lose the helpers above
        return Private::RoundDownTo_ (x, toNearest);
    }
 
    /*
     ********************************************************************************
     ********************************* Math::Round **********************************
     ********************************************************************************
     */
    template <integral INT_TYPE, floating_point FLOAT_TYPE>
    constexpr INT_TYPE Round (FLOAT_TYPE n)
    {
        FLOAT_TYPE tmp = ::round (n);
        if (tmp > 0) {
            DISABLE_COMPILER_CLANG_WARNING_START ("clang diagnostic ignored \"-Wimplicit-int-float-conversion\""); // warning: implicit conversion from 'std::__1::numeric_limits<long>::type' (aka 'long') to 'double' changes value from 9223372036854775807 to 9223372036854775808
            return tmp >= numeric_limits<INT_TYPE>::max () ? numeric_limits<INT_TYPE>::max () : static_cast<INT_TYPE> (tmp);
            DISABLE_COMPILER_CLANG_WARNING_END ("clang diagnostic ignored \"-Wimplicit-int-float-conversion\"");
        }
        else {
            return tmp <= numeric_limits<INT_TYPE>::min () ? numeric_limits<INT_TYPE>::min () : static_cast<INT_TYPE> (tmp);
        }
    }
    template <floating_point FLOAT_TYPE>
    inline FLOAT_TYPE Round (FLOAT_TYPE n, unsigned int nDigitsOfPrecision)
    {
        using Common::StdCompat::isinf;
        using Common::StdCompat::isnan;
        Require (nDigitsOfPrecision >= 1);
        if (isnan (n) or isinf (n)) [[unlikely]] {
            return n;
        }
        auto         absN                = fabs (n);
        unsigned int digitsBeforeDecimal = 0;
        if (absN >= 1) {
            digitsBeforeDecimal = Trunc<unsigned int> (log10 (absN)) + 1;
        }
        FLOAT_TYPE pow10Shifter = pow (10, static_cast<int> (nDigitsOfPrecision) - static_cast<int> (digitsBeforeDecimal));
        return round (n * pow10Shifter) / pow10Shifter;
    }
 
    /*
     ********************************************************************************
     ******************************** Math::Trunc ***********************************
     ********************************************************************************
     */
    namespace Private_ {
        void ThrowTruncOfNAN_ ();
    }
    template <integral INT_TYPE, floating_point FLOAT_TYPE>
    constexpr INT_TYPE Trunc (FLOAT_TYPE n)
    {
        using Common::StdCompat::isnan;
        if (isnan (n)) [[unlikely]] {
            Private_::ThrowTruncOfNAN_ ();
        }
        FLOAT_TYPE tmp = ::trunc (n);
        if (tmp > 0) {
            DISABLE_COMPILER_CLANG_WARNING_START ("clang diagnostic ignored \"-Wimplicit-int-float-conversion\""); // warning: implicit conversion from 'std::__1::numeric_limits<long>::type' (aka 'long') to 'double' changes value from 9223372036854775807 to 9223372036854775808
            return tmp >= numeric_limits<INT_TYPE>::max () ? numeric_limits<INT_TYPE>::max () : static_cast<INT_TYPE> (tmp);
            DISABLE_COMPILER_CLANG_WARNING_END ("clang diagnostic ignored \"-Wimplicit-int-float-conversion\"");
        }
        else {
            return tmp <= numeric_limits<INT_TYPE>::min () ? numeric_limits<INT_TYPE>::min () : static_cast<INT_TYPE> (tmp);
        }
    }
 
    /*
     ********************************************************************************
     **************************** Math::NearlyEquals ********************************
     ********************************************************************************
     */
    template <Common::IBuiltinArithmetic T1, Common::IBuiltinArithmetic T2, typename EPSILON_TYPE>
    constexpr bool NearlyEquals (T1 l, T2 r, EPSILON_TYPE epsilon)
    {
        using Common::StdCompat::isinf;
        using Common::StdCompat::isnan;
        if (isnan (l) or isnan (r)) [[unlikely]] {
            return isnan (l) and isnan (r);
        }
        else if (isinf (l) or isinf (r)) [[unlikely]] {
            // only 'equal' if inf and of same sign
            // https://stackoverflow.com/questions/41834621/c-ieee-floats-inf-equal-inf
            return l == r;
        }
        auto diff = l - r;
        Require (epsilon >= 0); // other cases we ignore epsilon
        Assert (not isnan (l) and not isnan (r) and not isinf (l) and not isinf (r));
        return Abs (diff) <= epsilon;
    }
    template <Common::IBuiltinArithmetic T1, Common::IBuiltinArithmetic T2>
    constexpr bool NearlyEquals (T1 l, T2 r)
    {
        using TC = common_type_t<T1, T2>;
        if constexpr (floating_point<TC>) {
            constexpr TC kEpsilon_ = 10000 * numeric_limits<TC>::epsilon (); // pick more than epsilon cuz some math functions have more error than a single bit... - even 1000x not enuf sometimes
            // \see https://realtimecollisiondetection.net/blog/?p=89
            // using relTol = absTol
            TC useEpsilon = kEpsilon_ * std::max<TC> ({static_cast<TC> (1.0), static_cast<TC> (Abs (l)), static_cast<TC> (Abs (r))});
            return NearlyEquals (l, r, useEpsilon);
        }
        else {
            return l == r;
        }
    }
 
    /*
     ********************************************************************************
     ************************* Math::PinToSpecialPoint ******************************
     ********************************************************************************
     */
    template <floating_point T>
    constexpr T PinToSpecialPoint (T p, T special)
    {
        if (Math::NearlyEquals (p, special)) {
            return special;
        }
        return p;
    }
    template <floating_point T>
    constexpr T PinToSpecialPoint (T p, T special, T epsilon)
    {
        if (Math::NearlyEquals (p, special, epsilon)) {
            return special;
        }
        return p;
    }
 
    /*
     ********************************************************************************
     ***************************** Math::PinInRange *********************************
     ********************************************************************************
     */
    template <typename T>
    [[deprecated ("Since Stroika v3.0d12 - use std::clamp")]] constexpr T PinInRange (T initialValue, T lowerBound, T upperBound)
    {
        Require (lowerBound <= upperBound);
        auto r  = max (lowerBound, min (upperBound, initialValue));
        auto r2 = clamp (initialValue, lowerBound, upperBound);
        Assert (r2 == r);
        return r2;
    }
 
    /*
     ********************************************************************************
     ******************************** Math::AtLeast *********************************
     ********************************************************************************
     */
    template <typename T>
    constexpr T AtLeast (T initialValue, T lowerBound)
    {
        return max (initialValue, lowerBound);
    }
 
    /*
     ********************************************************************************
     ******************************** Math::AtMost **********************************
     ********************************************************************************
     */
    template <typename T>
    constexpr T AtMost (T initialValue, T upperBound)
    {
        return min (initialValue, upperBound);
    }
 
    /*
     ********************************************************************************
     ************************** Math::PinToMaxForType *******************************
     ********************************************************************************
     */
    template <typename NEW_T, typename T>
    constexpr NEW_T PinToMaxForType (T initialValue)
    {
        using LargerType = decltype (NEW_T{} + T{}); // maybe should use conditional<> for this?
        return static_cast<NEW_T> (min<LargerType> (initialValue, numeric_limits<NEW_T>::max ()));
    }
 
    /*
     ********************************************************************************
     ********************************** Math::Abs ***********************************
     ********************************************************************************
     */
    template <Common::IBuiltinArithmetic T, typename RESULT_TYPE>
    constexpr RESULT_TYPE Abs (T v)
    {
        if constexpr (is_unsigned_v<T>) {
            return v; // and no unary minus on an unsigned type, which MSVC warns about (C4146) even where never evaluated
        }
        else {
            return v < 0 ? -v : v;
        }
    }
 
    /*
     ********************************************************************************
     ********************************** Math::IsOdd *********************************
     ********************************************************************************
     */
    template <integral T>
    constexpr bool IsOdd (T v)
    {
        return v % 2 == 1;
    }
 
    /*
     ********************************************************************************
     ********************************* Math::IsEven *********************************
     ********************************************************************************
     */
    template <integral T>
    constexpr bool IsEven (T v)
    {
        static_assert (integral<T>);
        return v % 2 == 0;
    }
 
    /*
     ********************************************************************************
     ******************************** Math::IsPrime *********************************
     ********************************************************************************
     */
    template <integral T>
    bool IsPrime (T v)
    {
        Require (v >= 0); // no negative numbers
        // @todo - redo this as http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes but this is simpler and
        // has no memory requirements and is good enuf for me to test...
        if (v == 1) {
            return false;
        }
        if (v == 2) {
            return true; // special case
        }
        T checkUpTo = static_cast<T> (::sqrt (v)) + static_cast<T> (1);
        // Check each number from 3 up to checkUpTo and see if its a divisor
        for (T d = 2; d <= checkUpTo; ++d) {
            if (v % d == 0) {
                return false;
            }
        }
        return true;
    }
 
    /*
     ********************************************************************************
     *************************** Math::PrimeAtLeastThisBig **************************
     ********************************************************************************
     */
    template <integral T>
    T PrimeAtLeastThisBig (T v)
    {
        while (not IsPrime (v)) {
            ++v;
        }
        return v;
    }
 
}