num_util.h

Go to the documentation of this file.

 
#ifndef GUL14_NUM_UTIL_H_
#define GUL14_NUM_UTIL_H_
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <type_traits>
 
#include "gul14/internal.h"
 
namespace gul14 {
 
template<typename ValueT>
constexpr auto abs(ValueT n) noexcept -> std::enable_if_t<std::is_unsigned<ValueT>::value, ValueT>
{
    return n;
}
 
template<typename ValueT>
constexpr auto abs(ValueT n) noexcept -> std::enable_if_t<not std::is_unsigned<ValueT>::value, ValueT>
{
    return std::abs(n);
}
 
template<typename NumT, typename OrderT,
    typename = std::enable_if_t<
        std::is_arithmetic<NumT>::value
        and std::is_arithmetic<OrderT>::value
    >>
bool within_orders(const NumT a, const NumT b, const OrderT orders) noexcept(false) {
    // std::pow() is not noexcept, which might or might not be true
    auto difference = gul14::abs(a - b);
    if (difference == NumT{ 0 })
        return true;
    auto maximum = std::max(gul14::abs(a), gul14::abs(b));
    auto limit = maximum / std::pow(static_cast<std::decay_t<NumT>>(10.0), orders);
    return difference < limit;
}
 
template<typename NumT>
bool within_abs(NumT a, NumT b, NumT tol) noexcept {
    tol = gul14::abs(tol); // Negative diff does not make sense
    bool ret{};
    if (a > b) {
        if (std::is_floating_point<NumT>::value)
            ret = a - tol <= b; // different formula needed because of inf/-inf and subnormal values
        else
            ret = a - b <= tol;
    } else {
        if (std::is_floating_point<NumT>::value)
            ret = b - tol <= a;
        else
            ret = b - a <= tol;
    }
    return ret;
}
 
template<typename NumT,
    typename = std::enable_if_t<std::is_floating_point<NumT>::value>>
bool within_ulp(NumT a, NumT b, unsigned int ulp)
{
    if (a == b)
        return true;
    auto diff = std::abs(a - b);
    auto maxval = std::max(std::abs(a), std::abs(b));
 
    // NORMAL cases
    if (std::isnormal(maxval))
    {
        return diff <= std::numeric_limits<NumT>::epsilon()
                       * static_cast<NumT>(ulp) * maxval;
    }
 
    // SUBNORMAL cases
    if (std::isfinite(maxval))
        return diff <= std::numeric_limits<NumT>::denorm_min() * static_cast<NumT>(ulp);
 
    // NAN and INF cases: Two NANs or (a xor b is (-)INF) or (INF, -INF pair)
    return false;
}
 
template<class NumT>
constexpr const NumT& clamp(const NumT& v, const NumT& lo, const NumT& hi) {
    return std::max(lo, std::min(v, hi)); // usually more optimized than actually using operator<()
}
template<class NumT, class Compare>
constexpr const NumT& clamp(const NumT& v, const NumT& lo, const NumT& hi, Compare comp) {
    if (comp(v, lo))
        return lo;
    if (comp(hi, v))
        return hi;
    return v;
}
 
} // namespace gul14
 
#endif
 
// vi:ts=4:sw=4:et