|
| 1 | +// Copyright (c) Microsoft Corporation. |
| 2 | +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
| 3 | + |
| 4 | +#pragma once |
| 5 | + |
| 6 | +#include <cassert> |
| 7 | +#include <cfenv> |
| 8 | +#include <cmath> |
| 9 | +#include <float.h> |
| 10 | +#include <type_traits> |
| 11 | +#include <xutility> |
| 12 | + |
| 13 | +namespace fputil { |
| 14 | + template <typename T> |
| 15 | + using float_bits_t = typename _STD _Float_traits<T>::type; |
| 16 | + |
| 17 | + template <typename T> |
| 18 | + _INLINE_VAR constexpr float_bits_t<T> magnitude_mask_v = _STD _Float_traits<T>::_Magnitude_mask; |
| 19 | + |
| 20 | + template <typename T> |
| 21 | + _INLINE_VAR constexpr float_bits_t<T> exponent_mask_v = _STD _Float_traits<T>::_Exponent_mask; |
| 22 | + |
| 23 | + template <typename T> |
| 24 | + _INLINE_VAR constexpr float_bits_t<T> significand_mask_v = magnitude_mask_v<T> & ~exponent_mask_v<T>; |
| 25 | + |
| 26 | + template <typename T> |
| 27 | + _INLINE_VAR constexpr float_bits_t<T> sign_mask_v = _STD _Float_traits<T>::_Sign_mask; |
| 28 | + |
| 29 | + template <typename T> |
| 30 | + _INLINE_VAR constexpr float_bits_t<T> norm_min_bits_v = significand_mask_v<T> + 1U; |
| 31 | + |
| 32 | + template <typename T> |
| 33 | + _INLINE_VAR constexpr float_bits_t<T> norm_max_bits_v = exponent_mask_v<T> - 1U; |
| 34 | + |
| 35 | + template <typename T> |
| 36 | + _INLINE_VAR constexpr float_bits_t<T> infinity_bits_v = exponent_mask_v<T>; |
| 37 | + |
| 38 | + // not affected by abrupt underflow |
| 39 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 40 | + constexpr bool iszero(const T& x) { |
| 41 | + return _STD _Float_abs_bits(x) == 0; |
| 42 | + } |
| 43 | + |
| 44 | + // not affected by /fp:fast |
| 45 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 46 | + constexpr bool signbit(const T& x) { |
| 47 | + const auto bits = std::_Bit_cast<float_bits_t<T>>(x); |
| 48 | + return (bits & sign_mask_v<T>) != 0; |
| 49 | + } |
| 50 | + |
| 51 | + enum class rounding_mode { |
| 52 | + to_nearest_ties_even = FE_TONEAREST, |
| 53 | + toward_zero = FE_TOWARDZERO, |
| 54 | + toward_positive = FE_UPWARD, |
| 55 | + toward_negative = FE_DOWNWARD, |
| 56 | + }; |
| 57 | + |
| 58 | + bool is_directed_rounding_mode(const rounding_mode mode) { |
| 59 | + switch (mode) { |
| 60 | + case rounding_mode::to_nearest_ties_even: |
| 61 | + return false; |
| 62 | + |
| 63 | + case rounding_mode::toward_zero: |
| 64 | + case rounding_mode::toward_positive: |
| 65 | + case rounding_mode::toward_negative: |
| 66 | + return true; |
| 67 | + |
| 68 | + default: |
| 69 | + assert(false); |
| 70 | + return false; |
| 71 | + } |
| 72 | + } |
| 73 | + |
| 74 | +#if TEST_FP_ROUNDING |
| 75 | + |
| 76 | +#ifdef __clang__ |
| 77 | +// TRANSITION, should be #pragma STDC FENV_ACCESS ON |
| 78 | +#else // ^^^ clang / MSVC vvv |
| 79 | +// TRANSITION, VSO-923474 -- should be #pragma STDC FENV_ACCESS ON |
| 80 | +#pragma fenv_access(on) |
| 81 | +#endif // ^^^ MSVC ^^^ |
| 82 | + |
| 83 | + constexpr rounding_mode all_rounding_modes[] = { |
| 84 | + rounding_mode::to_nearest_ties_even, |
| 85 | + rounding_mode::toward_zero, |
| 86 | + rounding_mode::toward_positive, |
| 87 | + rounding_mode::toward_negative, |
| 88 | + }; |
| 89 | + |
| 90 | + class rounding_guard { |
| 91 | + public: |
| 92 | + explicit rounding_guard(const rounding_mode mode) : old_mode{static_cast<rounding_mode>(std::fegetround())} { |
| 93 | + const int result = std::fesetround(static_cast<int>(mode)); |
| 94 | + assert(result == 0); |
| 95 | + } |
| 96 | + |
| 97 | + ~rounding_guard() { |
| 98 | + const int result = std::fesetround(static_cast<int>(old_mode)); |
| 99 | + assert(result == 0); |
| 100 | + } |
| 101 | + |
| 102 | + rounding_guard(const rounding_guard&) = delete; |
| 103 | + rounding_guard& operator=(const rounding_guard&) = delete; |
| 104 | + |
| 105 | + private: |
| 106 | + rounding_mode old_mode; |
| 107 | + }; |
| 108 | + |
| 109 | +#else // ^^^ alternative rounding modes / default rounding mode only vvv |
| 110 | + |
| 111 | + constexpr rounding_mode all_rounding_modes[] = {rounding_mode::to_nearest_ties_even}; |
| 112 | + |
| 113 | + class rounding_guard { |
| 114 | + public: |
| 115 | + explicit rounding_guard(const rounding_mode mode) { |
| 116 | + static_cast<void>(mode); |
| 117 | + } |
| 118 | + |
| 119 | + ~rounding_guard() = default; |
| 120 | + |
| 121 | + rounding_guard(const rounding_guard&) = delete; |
| 122 | + rounding_guard& operator=(const rounding_guard&) = delete; |
| 123 | + }; |
| 124 | + |
| 125 | +#endif // ^^^ default rounding mode only ^^^ |
| 126 | + |
| 127 | + // compares whether two floating point values are equal |
| 128 | + // all NaNs are equal, +0.0 and -0.0 are not equal |
| 129 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 130 | + bool precise_equal(const T& actual, const T& expected) { |
| 131 | + if (_STD _Is_nan(actual) || _STD _Is_nan(expected)) { |
| 132 | + return _STD _Is_nan(actual) == _STD _Is_nan(expected); |
| 133 | + } else { |
| 134 | + return actual == expected && fputil::signbit(actual) == fputil::signbit(expected); |
| 135 | + } |
| 136 | + } |
| 137 | + |
| 138 | + namespace detail { |
| 139 | + // 0x80...00 = zero, 0x80...01 = numeric_limits<T>::denorm_min(), 0x7f...ff = -numeric_limits<T>::denorm_min() |
| 140 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 141 | + float_bits_t<T> offset_representation(const T& x) { |
| 142 | + const float_bits_t<T> abs_bits = _STD _Float_abs_bits(x); |
| 143 | + return fputil::signbit(x) ? sign_mask_v<T> - abs_bits : sign_mask_v<T> + abs_bits; |
| 144 | + } |
| 145 | + |
| 146 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 147 | + float_bits_t<T> is_offset_value_subnormal_or_zero(const float_bits_t<T> offset_value) { |
| 148 | + constexpr float_bits_t<T> positive_norm_min_offset = sign_mask_v<T> + norm_min_bits_v<T>; |
| 149 | + constexpr float_bits_t<T> negative_norm_min_offset = sign_mask_v<T> - norm_min_bits_v<T>; |
| 150 | + |
| 151 | + return negative_norm_min_offset < offset_value && offset_value < positive_norm_min_offset; |
| 152 | + } |
| 153 | + |
| 154 | + // number of ulps above zero, if we count [0, numeric_limits<T>::min()) as 1 ulp |
| 155 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 156 | + double abrupt_underflow_ulp(const float_bits_t<T> offset_value) { |
| 157 | + using bits_type = float_bits_t<T>; |
| 158 | + |
| 159 | + constexpr bits_type offset_positive_norm_min = sign_mask_v<T> + norm_min_bits_v<T>; |
| 160 | + constexpr bits_type offset_negative_norm_min = sign_mask_v<T> - norm_min_bits_v<T>; |
| 161 | + |
| 162 | + if (offset_value >= offset_positive_norm_min) { |
| 163 | + return 1.0 + (offset_value - offset_positive_norm_min); |
| 164 | + } else if (offset_value <= offset_negative_norm_min) { |
| 165 | + return -1.0 - (offset_negative_norm_min - offset_value); |
| 166 | + } else if (offset_value >= sign_mask_v<T>) { |
| 167 | + return static_cast<double>(offset_value - sign_mask_v<T>) / norm_min_bits_v<T>; |
| 168 | + } else { |
| 169 | + return -static_cast<double>(sign_mask_v<T> - offset_value) / norm_min_bits_v<T>; |
| 170 | + } |
| 171 | + } |
| 172 | + |
| 173 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 174 | + bool is_within_ulp_tolerance(const T& actual, const T& expected, const int ulp_tolerance) { |
| 175 | + if (_STD _Is_nan(actual) || _STD _Is_nan(expected)) { |
| 176 | + return _STD _Is_nan(actual) == _STD _Is_nan(expected); |
| 177 | + } |
| 178 | + |
| 179 | + if (_STD _Is_inf(expected)) { |
| 180 | + return actual == expected; |
| 181 | + } |
| 182 | + |
| 183 | + if (fputil::signbit(actual) != fputil::signbit(expected)) { |
| 184 | + return false; |
| 185 | + } |
| 186 | + |
| 187 | + using bits_type = float_bits_t<T>; |
| 188 | + |
| 189 | + // compute ulp difference |
| 190 | + const bits_type actual_offset = detail::offset_representation(actual); |
| 191 | + const bits_type expected_offset = detail::offset_representation(expected); |
| 192 | + const bits_type ulp_diff = |
| 193 | + actual_offset < expected_offset ? expected_offset - actual_offset : actual_offset - expected_offset; |
| 194 | + |
| 195 | + if (ulp_diff <= static_cast<unsigned int>(ulp_tolerance) && ulp_tolerance >= 0) { |
| 196 | + return true; |
| 197 | + } |
| 198 | + |
| 199 | +#if WITH_FP_ABRUPT_UNDERFLOW |
| 200 | + // handle abrupt underflow |
| 201 | + if (detail::is_offset_value_subnormal_or_zero<T>(expected_offset) |
| 202 | + || detail::is_offset_value_subnormal_or_zero<T>(actual_offset)) { |
| 203 | + const double adjusted_actual_ulp = detail::abrupt_underflow_ulp<T>(actual_offset); |
| 204 | + const double adjusted_expected_ulp = detail::abrupt_underflow_ulp<T>(expected_offset); |
| 205 | + const double adjusted_ulp_diff = std::abs(adjusted_actual_ulp - adjusted_expected_ulp); |
| 206 | + |
| 207 | + if (adjusted_ulp_diff <= ulp_tolerance) { |
| 208 | + return true; |
| 209 | + } |
| 210 | + } |
| 211 | +#endif // WITH_FP_ABRUPT_UNDERFLOW |
| 212 | + |
| 213 | + return false; |
| 214 | + } |
| 215 | + |
| 216 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 217 | + bool is_within_absolute_tolerance(const T& actual, const T& expected, const double absolute_tolerance) { |
| 218 | + return _STD _Is_finite(actual) && _STD _Is_finite(expected) |
| 219 | + && std::abs(actual - expected) <= absolute_tolerance; |
| 220 | + } |
| 221 | + } // namespace detail |
| 222 | + |
| 223 | + // returns whether floating point result is nearly equal to the expected value |
| 224 | + template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0> |
| 225 | + bool near_equal( |
| 226 | + const T& actual, const T& expected, const int ulp_tolerance = 1, const double absolute_tolerance = 0) { |
| 227 | + if (precise_equal(actual, expected)) { |
| 228 | + return true; |
| 229 | + } |
| 230 | + |
| 231 | + if (ulp_tolerance > 0 && detail::is_within_ulp_tolerance(actual, expected, ulp_tolerance)) { |
| 232 | + return true; |
| 233 | + } |
| 234 | + |
| 235 | + if (absolute_tolerance > 0 && detail::is_within_absolute_tolerance(actual, expected, absolute_tolerance)) { |
| 236 | + return true; |
| 237 | + } |
| 238 | + |
| 239 | + return false; |
| 240 | + } |
| 241 | +} // namespace fputil |
0 commit comments