structural_real.hpp

Go to the documentation of this file.

#ifndef TDLS_CORE_STRUCTURAL_REAL_HPP
#define TDLS_CORE_STRUCTURAL_REAL_HPP
 
 
 
 
 
 
#include <type_traits>
 
#include <tdls/core/macros.hpp>
 
 
 
namespace tdls {
 
 
 
template<typename T>
struct StructuralReal {
 
    static_assert(std::is_floating_point_v<T>, "tdls::StructuralReal: T must be floating-point");
 
    long long mantissa = 0; 
    int exponent       = 0; 
 
    constexpr StructuralReal() noexcept = default;
 
    TDLS_HOST_DEVICE constexpr StructuralReal(const long long m, const int e) noexcept
        : mantissa(m), exponent(e) {
    }
 
    TDLS_HOST_DEVICE constexpr StructuralReal(const T v) noexcept {
        if (v == T(0)) return;
        // A NaN differs from itself; an infinity is the only non-zero
        // value unchanged by halving (zero was excluded above). Neither
        // test produces a NaN, which constant evaluation would reject,
        // and neither needs a library call, which nvcc would reject in
        // this host-and-device constructor.
        if (v != v || v * T(0.5) == v) {
            // Non-finite sentinel: no finite value yields a zero mantissa
            // with a non-zero exponent.
            mantissa = 0;
            exponent = 1;
            return;
        }
        constexpr T big = T(1LL << 62);
        T a             = v;
        int e           = 0;
        while (a >= big || a <= -big) {
            a /= T(2);
            ++e;
        }
        while (T(static_cast<long long>(a)) != a) {
            // A mantissa wider than 63 bits (long double): the sentinel.
            if (a >= big || a <= -big) {
                mantissa = 0;
                exponent = 1;
                return;
            }
            a *= T(2);
            --e;
        }
        long long m = static_cast<long long>(a);
        while ((m & 1LL) == 0) {
            m /= 2;
            ++e;
        }
        mantissa = m;
        exponent = e;
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE constexpr bool is_finite() const noexcept {
        return !(mantissa == 0 && exponent != 0);
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE constexpr operator T() const noexcept {
        T r = T(mantissa);
        for (int e = exponent; e > 0; --e)
            r *= T(2);
        for (int e = exponent; e < 0; ++e)
            r /= T(2);
        return r;
    }
};
 
 
 
} // namespace tdls
 
 
 
#endif // TDLS_CORE_STRUCTURAL_REAL_HPP