| // SPDX-FileCopyrightText: The Eigen Authors |
| // SPDX-License-Identifier: MPL-2.0 |
| |
| #include "main.h" |
| #include "fp_control.h" |
| #define EIGEN_TEST_ANNOYING_SCALAR_DONT_THROW |
| #include "AnnoyingScalar.h" |
| |
| template <typename T> |
| void check_power_of_two_scaling_factor() { |
| static_assert(internal::supports_power_of_two_scaling<T>::value, "power-of-two scaling must be enabled"); |
| using Factors = internal::safe_scaling_factors<T>; |
| |
| const auto check_factor = [](const T& value) { |
| const auto factors = internal::safe_scaling<T>::compute_floor_factors(value); |
| VERIFY(factors.scale >= (std::numeric_limits<T>::min)()); |
| VERIFY(factors.invScale >= (std::numeric_limits<T>::min)()); |
| VERIFY_IS_EQUAL(factors.scale * factors.invScale, T(1)); |
| return factors; |
| }; |
| |
| const auto check_round_trip = [&](const T& value) { |
| const auto factors = check_factor(value); |
| Matrix<T, 1, 1> input; |
| input(0) = value; |
| Matrix<T, 1, 1> scaled; |
| const auto scaleToFactors = internal::safe_scaling<T>::scale_to(scaled, input, value); |
| VERIFY_IS_EQUAL(scaleToFactors.scale, factors.scale); |
| VERIFY_IS_EQUAL(scaleToFactors.invScale, factors.invScale); |
| Matrix<T, 1, 1> scaledInPlace; |
| scaledInPlace(0) = value; |
| internal::safe_scaling<T>::scale_in_place(scaledInPlace, value, factors); |
| VERIFY_IS_EQUAL(scaledInPlace, scaled); |
| scaledInPlace(0) = value; |
| const auto inPlaceFactors = internal::safe_scaling<T>::scale_in_place(scaledInPlace, value); |
| VERIFY_IS_EQUAL(inPlaceFactors.scale, factors.scale); |
| VERIFY_IS_EQUAL(inPlaceFactors.invScale, factors.invScale); |
| VERIFY_IS_EQUAL(scaledInPlace, scaled); |
| |
| T restored; |
| internal::safe_scaling<T>::unscale_to(restored, scaled(0), factors); |
| VERIFY_IS_EQUAL(restored, value); |
| Matrix<T, 1, 1> restoredMatrix; |
| internal::safe_scaling<T>::unscale_to(restoredMatrix, scaled, value, factors); |
| VERIFY_IS_EQUAL(restoredMatrix, input); |
| internal::safe_scaling<T>::unscale_in_place(scaled, factors); |
| VERIFY_IS_EQUAL(scaled(0), value); |
| }; |
| |
| const auto denormFactors = check_factor(std::numeric_limits<T>::denorm_min()); |
| VERIFY_IS_EQUAL(denormFactors.scale, (std::numeric_limits<T>::min)()); |
| Factors normalReciprocalFactors; |
| VERIFY(internal::safe_scaling<T>::try_compute_ceiling_factors_with_normal_reciprocal( |
| std::numeric_limits<T>::denorm_min(), (std::numeric_limits<T>::min)(), normalReciprocalFactors)); |
| VERIFY_IS_EQUAL(normalReciprocalFactors.scale, denormFactors.scale); |
| VERIFY_IS_EQUAL(normalReciprocalFactors.invScale, denormFactors.invScale); |
| |
| // bfloat16 arithmetic widens to float, whose subnormal inputs may be flushed before scaling can recover them. |
| if (!std::is_same<T, bfloat16>::value || !ScopedFlushToZero::hardwareFlushesSubnormalInputs()) { |
| using InputScalar = std::conditional_t<std::is_floating_point<T>::value, T, float>; |
| volatile InputScalar denormInput = static_cast<InputScalar>(std::numeric_limits<T>::denorm_min()); |
| const T denormMin = static_cast<T>(denormInput); |
| Matrix<T, 1, 1> denormInputMatrix; |
| denormInputMatrix(0) = denormMin; |
| Matrix<T, 1, 1> scaledDenorm; |
| const auto factors = internal::safe_scaling<T>::scale_to(scaledDenorm, denormInputMatrix, denormMin); |
| Matrix<T, 1, 1> scaledDenormInPlace; |
| scaledDenormInPlace(0) = denormMin; |
| internal::safe_scaling<T>::scale_in_place(scaledDenormInPlace, denormMin, factors); |
| VERIFY_IS_EQUAL(scaledDenormInPlace, scaledDenorm); |
| VERIFY(scaledDenorm(0) > T(0)); |
| } |
| |
| check_round_trip((std::numeric_limits<T>::min)()); |
| check_round_trip(T(0.75)); |
| check_round_trip(T(1.5)); |
| check_round_trip(T(3)); |
| check_round_trip((std::numeric_limits<T>::max)()); |
| |
| const auto floorFactors = internal::safe_scaling<T>::compute_floor_factors(T(1.5)); |
| const auto ceilingFactors = internal::safe_scaling<T>::compute_ceiling_factors_with_normal_reciprocal(T(1.5)); |
| VERIFY_IS_EQUAL(floorFactors.scale, T(1)); |
| VERIFY_IS_EQUAL(ceilingFactors.scale, T(2)); |
| |
| const T highest = (std::numeric_limits<T>::max)(); |
| const auto highestCeilingFactors = internal::safe_scaling<T>::compute_ceiling_factors_with_normal_reciprocal(highest); |
| VERIFY_IS_EQUAL(highestCeilingFactors.scale, check_factor(highest).scale); |
| VERIFY_IS_EQUAL(highestCeilingFactors.scale * highestCeilingFactors.invScale, T(1)); |
| const T scaledHighest = highest * highestCeilingFactors.invScale; |
| VERIFY(scaledHighest >= T(1)); |
| if (std::numeric_limits<T>::max_exponent + std::numeric_limits<T>::min_exponent <= 3) { |
| VERIFY(scaledHighest < T(4)); |
| } |
| } |
| |
| template <typename T> |
| void check_safe_scaling_special_values() { |
| using Scaling = internal::safe_scaling<T>; |
| using Factors = internal::safe_scaling_factors<T>; |
| |
| const Factors identity; |
| VERIFY_IS_EQUAL(identity.scale, T(1)); |
| VERIFY_IS_EQUAL(identity.invScale, T(1)); |
| |
| Matrix<T, 1, 1> zero = Matrix<T, 1, 1>::Zero(); |
| Matrix<T, 1, 1> scaledZero; |
| const auto zeroFactors = Scaling::scale_to(scaledZero, zero, T(0)); |
| VERIFY_IS_EQUAL(zeroFactors.scale, T(1)); |
| VERIFY_IS_EQUAL(zeroFactors.invScale, T(1)); |
| VERIFY_IS_EQUAL(scaledZero(0), T(0)); |
| |
| for (const T special : {std::numeric_limits<T>::infinity(), std::numeric_limits<T>::quiet_NaN()}) { |
| Matrix<T, 2, 1> input; |
| input << special, T(2); |
| Matrix<T, 2, 1> scaled; |
| const auto factors = Scaling::scale_to(scaled, input, special); |
| VERIFY_IS_EQUAL(factors.scale, T(1)); |
| VERIFY_IS_EQUAL(factors.invScale, T(1)); |
| VERIFY_IS_EQUAL(scaled(1), T(2)); |
| if ((numext::isnan)(special)) { |
| VERIFY((numext::isnan)(scaled(0))); |
| } else { |
| VERIFY((numext::isinf)(scaled(0))); |
| } |
| } |
| } |
| |
| void check_safe_scaling_special_value_frontends() { |
| using Scaling = internal::safe_scaling<float>; |
| using Factors = internal::safe_scaling_factors<float>; |
| |
| for (const float special : {std::numeric_limits<float>::infinity(), std::numeric_limits<float>::quiet_NaN()}) { |
| Vector2f input(special, 2.0f); |
| Vector2f scaledTo; |
| const auto scaleToFactors = Scaling::scale_to(scaledTo, input, special); |
| VERIFY_IS_EQUAL(scaleToFactors.scale, 1.0f); |
| VERIFY_IS_EQUAL(scaleToFactors.invScale, 1.0f); |
| |
| Vector2f scaledExpression; |
| const Factors expressionFactors = |
| Scaling::with_scaled(input, special, [&](const auto& expression) { scaledExpression = expression; }); |
| VERIFY_IS_EQUAL(expressionFactors.scale, 1.0f); |
| VERIFY_IS_EQUAL(expressionFactors.invScale, 1.0f); |
| VERIFY_IS_EQUAL(scaledTo(1), 2.0f); |
| VERIFY_IS_EQUAL(scaledExpression(1), 2.0f); |
| if ((numext::isnan)(special)) { |
| VERIFY((numext::isnan)(scaledTo(0))); |
| VERIFY((numext::isnan)(scaledExpression(0))); |
| } else { |
| VERIFY((numext::isinf)(scaledTo(0))); |
| VERIFY((numext::isinf)(scaledExpression(0))); |
| } |
| } |
| } |
| |
| void check_arithmetic_safe_scaling_fallback() { |
| using Scaling = internal::safe_scaling<double, false>; |
| using Factors = internal::safe_scaling_factors<double>; |
| |
| Matrix<double, 2, 1> input; |
| input << 3.0, 6.0; |
| Matrix<double, 2, 1> scaled; |
| const auto factors = Scaling::scale_to(scaled, input, 3.0); |
| VERIFY_IS_EQUAL(factors.scale, 3.0); |
| VERIFY_IS_EQUAL(factors.invScale, 1.0); |
| VERIFY_IS_EQUAL(scaled(0), 1.0); |
| VERIFY_IS_EQUAL(scaled(1), 2.0); |
| Matrix<double, 2, 1> restored; |
| Scaling::unscale_to(restored, scaled, 3.0, factors); |
| VERIFY_IS_EQUAL(restored, input); |
| Scaling::unscale_in_place(scaled, factors); |
| VERIFY_IS_EQUAL(scaled, input); |
| |
| const auto check_extreme_unscale = [](const double value) { |
| const Factors extremeFactors = Scaling::compute_floor_factors(value); |
| double restored; |
| Scaling::unscale_to(restored, 1.0, extremeFactors); |
| VERIFY_IS_EQUAL(restored, value); |
| restored = 1.0; |
| Scaling::unscale_in_place(restored, extremeFactors); |
| VERIFY_IS_EQUAL(restored, value); |
| }; |
| check_extreme_unscale((std::numeric_limits<double>::max)()); |
| volatile double denormMinInput = std::numeric_limits<double>::denorm_min(); |
| const double denormMin = denormMinInput; |
| if (denormMin > 0.0) check_extreme_unscale(denormMin); |
| |
| Factors normalReciprocalFactors; |
| const double normalMin = (std::numeric_limits<double>::min)(); |
| VERIFY(Scaling::try_compute_ceiling_factors_with_normal_reciprocal(3.0, normalMin, normalReciprocalFactors)); |
| VERIFY_IS_EQUAL(normalReciprocalFactors.scale, 3.0); |
| VERIFY_IS_EQUAL(normalReciprocalFactors.invScale, 1.0 / 3.0); |
| VERIFY(!Scaling::try_compute_ceiling_factors_with_normal_reciprocal(std::numeric_limits<double>::denorm_min(), |
| normalMin, normalReciprocalFactors)); |
| VERIFY(!Scaling::try_compute_ceiling_factors_with_normal_reciprocal((std::numeric_limits<double>::max)(), normalMin, |
| normalReciprocalFactors)); |
| } |
| |
| void check_custom_scalar_scaling_exceptions() { |
| using Scalar = AnnoyingScalar; |
| using Scaling = internal::safe_scaling<Scalar>; |
| static_assert(!internal::supports_power_of_two_scaling<Scalar>::value, "exercise arithmetic scaling"); |
| volatile float denormInput = std::numeric_limits<float>::denorm_min(); |
| const float denorm = denormInput; |
| const float twiceDenorm = denorm + denorm; |
| if (!(denorm > 0.0f && twiceDenorm > 0.0f)) return; // No custom-scalar subnormal recovery is promised. |
| |
| Matrix<Scalar, 2, 1> input; |
| input << Scalar(denorm), Scalar(twiceDenorm); |
| Matrix<Scalar, 2, 1> scaled; |
| std::fenv_t savedEnvironment; |
| if (std::feholdexcept(&savedEnvironment) != 0) return; |
| const auto factors = Scaling::scale_to(scaled, input, Scalar(twiceDenorm)); |
| const int overflow = std::fetestexcept(FE_OVERFLOW); |
| std::fesetenv(&savedEnvironment); |
| VERIFY_IS_EQUAL(overflow, 0); |
| VERIFY_IS_EQUAL(scaled(0), Scalar(0.5)); |
| VERIFY_IS_EQUAL(scaled(1), Scalar(1)); |
| Matrix<Scalar, 2, 1> restored; |
| Scaling::unscale_to(restored, scaled, Scalar(twiceDenorm), factors); |
| VERIFY_IS_EQUAL(restored, input); |
| Scaling::unscale_in_place(scaled, factors); |
| VERIFY_IS_EQUAL(scaled, input); |
| } |
| |
| template <typename Scalar> |
| void check_scale_binary_by_power_of_two() { |
| using Binary = internal::binary_floating_point_traits<Scalar>; |
| using Bits = typename Binary::Bits; |
| constexpr int kMinFactorExponent = 1 - std::numeric_limits<Scalar>::min_exponent - Binary::kFractionBits; |
| constexpr int kMaxFactorExponent = 1 - std::numeric_limits<Scalar>::min_exponent; |
| const Bits magnitudes[] = {Bits(0), |
| Bits(1), |
| Bits(3), |
| Binary::kFractionMask, |
| Binary::kExponentUnit, |
| Binary::kExponentUnit + (Binary::kExponentUnit >> 1)}; |
| |
| ScopedFlushToZero flushToZero; |
| // Cover the entire factor range selected by tiny-input recovery. Form the reference from integer significands, |
| // not subnormal floating-point inputs, so it remains valid under FTZ/DAZ. |
| for (int exponent = kMinFactorExponent; exponent <= kMaxFactorExponent; ++exponent) { |
| const Scalar factor = numext::ldexp(Scalar(1), exponent); |
| for (Bits magnitude : magnitudes) { |
| const Scalar expected = numext::ldexp(Scalar(magnitude), std::numeric_limits<Scalar>::min_exponent - |
| std::numeric_limits<Scalar>::digits + exponent); |
| for (Bits sign : {Bits(0), Bits(Binary::kSignBit)}) { |
| const Scalar value = numext::bit_cast<Scalar>(sign | magnitude); |
| const Scalar actual = internal::scale_binary_by_power_of_two(value, factor); |
| VERIFY_IS_EQUAL(Binary::bits(actual), sign | Binary::bits(expected)); |
| } |
| } |
| } |
| } |
| |
| // Rounds a subnormal result exactly where hardware arithmetic may flush it: the exact product m * 2^-k denorm_min is |
| // normal in double, and its nearest integer, ties to even, is the expected magnitude's bit pattern (a carry into the |
| // exponent field encodes the smallest normal). |
| template <typename Scalar> |
| void check_scale_binary_by_power_of_two_rounding() { |
| using Binary = internal::binary_floating_point_traits<Scalar>; |
| using Bits = typename Binary::Bits; |
| // Leading-bit significands in [kExponentUnit, 2 * kExponentUnit): ties and their neighbours at small shifts, and an |
| // all-ones tail that carries into the smallest normal. |
| const Bits significands[] = {Binary::kExponentUnit, |
| Binary::kExponentUnit | Bits(1), |
| Binary::kExponentUnit | Bits(2), |
| Binary::kExponentUnit | Bits(3), |
| Binary::kExponentUnit | Bits(5), |
| Binary::kExponentUnit + (Binary::kExponentUnit >> 1), |
| (Binary::kExponentUnit << 1) - Bits(1)}; |
| // A biased input exponent e scales m by 2^(e - 1); the factor 2^-(e - 1 + k) cancels it and shifts by k. |
| const int exponentFields[] = {1, 2, Binary::kExponentBias - Binary::kFractionBits - 3}; |
| |
| ScopedFlushToZero flushToZero; |
| for (int field : exponentFields) { |
| for (int k = 0; k <= Binary::kFractionBits + 3; ++k) { |
| const Scalar factor = numext::ldexp(Scalar(1), -(field - 1 + k)); |
| for (Bits m : significands) { |
| const Bits expected = Bits(std::nearbyint(std::ldexp(double(m), -k))); |
| for (Bits sign : {Bits(0), Bits(Binary::kSignBit)}) { |
| const Scalar value = |
| numext::bit_cast<Scalar>(sign | (Bits(field) << Binary::kFractionBits) | (m & Binary::kFractionMask)); |
| VERIFY_IS_EQUAL(Binary::bits(internal::scale_binary_by_power_of_two(value, factor)), sign | expected); |
| } |
| } |
| } |
| } |
| // Subnormal inputs normalize before they round. |
| for (Bits magnitude : {Bits(1), Bits(3), Binary::kExponentUnit >> 1, Binary::kFractionMask}) { |
| for (int k = 1; k <= Binary::kFractionBits + 2; ++k) { |
| const Bits expected = Bits(std::nearbyint(std::ldexp(double(magnitude), -k))); |
| for (Bits sign : {Bits(0), Bits(Binary::kSignBit)}) { |
| const Scalar actual = internal::scale_binary_by_power_of_two(numext::bit_cast<Scalar>(sign | magnitude), |
| numext::ldexp(Scalar(1), -k)); |
| VERIFY_IS_EQUAL(Binary::bits(actual), sign | expected); |
| } |
| } |
| } |
| } |
| |
| // Under DAZ a comparison reads a subnormal as zero, so compare representations. |
| template <typename RealScalar> |
| bool same_bits(const RealScalar& a, const RealScalar& b) { |
| using Binary = internal::binary_floating_point_traits<RealScalar>; |
| return Binary::bits(a) == Binary::bits(b); |
| } |
| |
| template <typename RealScalar> |
| bool same_bits(const std::complex<RealScalar>& a, const std::complex<RealScalar>& b) { |
| return same_bits(a.real(), b.real()) && same_bits(a.imag(), b.imag()); |
| } |
| |
| template <typename Scalar> |
| struct scaling_test_value { |
| using RealScalar = typename NumTraits<Scalar>::Real; |
| static Scalar run(RealScalar real, RealScalar) { return real; } |
| }; |
| |
| template <typename RealScalar> |
| struct scaling_test_value<std::complex<RealScalar>> { |
| static std::complex<RealScalar> run(RealScalar real, RealScalar imag) { return std::complex<RealScalar>(real, imag); } |
| }; |
| |
| template <typename Scalar> |
| void check_arithmetic_scaling_expression() { |
| using RealScalar = typename NumTraits<Scalar>::Real; |
| using Scaling = internal::safe_scaling<RealScalar, false>; |
| using Vector2 = Matrix<Scalar, 2, 1>; |
| using Value = scaling_test_value<Scalar>; |
| |
| // Lazy scaling and its adjoint must retain the division policy even below the binary-recovery threshold. |
| for (const RealScalar maxCoeff : |
| {RealScalar(3), numext::ldexp(RealScalar(1.5), std::numeric_limits<RealScalar>::min_exponent + 7)}) { |
| Vector2 input; |
| input << Value::run(maxCoeff, RealScalar(0)), Value::run(maxCoeff / RealScalar(2), -maxCoeff / RealScalar(4)); |
| Vector2 expected; |
| expected << Value::run(RealScalar(1), RealScalar(0)), Value::run(RealScalar(0.5), RealScalar(-0.25)); |
| Vector2 materialized; |
| const auto factors = Scaling::scale_to(materialized, input, maxCoeff); |
| Vector2 lazy; |
| Matrix<Scalar, 1, 2> adjoint; |
| const auto expressionFactors = Scaling::with_scaled(input, maxCoeff, [&](const auto& expression) { |
| lazy = expression; |
| adjoint = expression.adjoint(); |
| }); |
| VERIFY_IS_EQUAL(materialized, expected); |
| VERIFY_IS_EQUAL(lazy, expected); |
| VERIFY_IS_EQUAL(adjoint, expected.adjoint()); |
| VERIFY_IS_EQUAL(expressionFactors.scale, factors.scale); |
| VERIFY_IS_EQUAL(expressionFactors.invScale, factors.invScale); |
| Vector2 inPlace = input; |
| const auto inPlaceFactors = Scaling::scale_in_place(inPlace, maxCoeff); |
| VERIFY_IS_EQUAL(inPlace, expected); |
| VERIFY_IS_EQUAL(inPlaceFactors.scale, factors.scale); |
| VERIFY_IS_EQUAL(inPlaceFactors.invScale, factors.invScale); |
| } |
| } |
| |
| template <typename Scalar> |
| void check_subnormal_preserving_scaling() { |
| using RealScalar = typename NumTraits<Scalar>::Real; |
| using Binary = internal::binary_floating_point_traits<RealScalar>; |
| using Bits = typename Binary::Bits; |
| const RealScalar maxCoeff = numext::ldexp((std::numeric_limits<RealScalar>::min)(), 8); |
| // Arithmetic or narrowing conversions could flush these inputs before they reach the scaling helper. |
| const RealScalar subnormalMaxCoeff = numext::bit_cast<RealScalar>(Bits(64)); |
| const auto factors = internal::safe_scaling<RealScalar>::compute_floor_factors(maxCoeff); |
| |
| Matrix<Scalar, 2, 1> input; |
| input(0) = scaling_test_value<Scalar>::run(maxCoeff, -maxCoeff); |
| input(1) = scaling_test_value<Scalar>::run(numext::bit_cast<RealScalar>(Binary::kExponentUnit >> 1), |
| numext::bit_cast<RealScalar>(Binary::kExponentUnit >> 2)); |
| Matrix<Scalar, 1, 1> subnormalInput; |
| subnormalInput(0) = scaling_test_value<Scalar>::run(subnormalMaxCoeff, RealScalar(0)); |
| ScopedFlushToZero flushToZero; |
| Matrix<Scalar, 2, 1> scaled; |
| internal::safe_scaling<RealScalar>::scale_to(scaled, input, maxCoeff, factors); |
| |
| Matrix<Scalar, 2, 1> scaledExpression; |
| const auto expressionFactors = internal::safe_scaling<RealScalar>::with_scaled( |
| input, maxCoeff, [&](const auto& expression) { scaledExpression = expression; }); |
| |
| VERIFY_IS_EQUAL(scaled(0), scaling_test_value<Scalar>::run(RealScalar(1), RealScalar(-1))); |
| VERIFY_IS_EQUAL(scaled(1), |
| scaling_test_value<Scalar>::run(RealScalar(1) / RealScalar(512), RealScalar(1) / RealScalar(1024))); |
| VERIFY_IS_EQUAL(expressionFactors.scale, factors.scale); |
| VERIFY_IS_EQUAL(expressionFactors.invScale, factors.invScale); |
| VERIFY_IS_EQUAL(scaledExpression, scaled); |
| |
| Matrix<Scalar, 2, 1> scaledInPlace = input; |
| const auto inPlaceFactors = internal::safe_scaling<RealScalar>::scale_in_place(scaledInPlace, maxCoeff); |
| VERIFY_IS_EQUAL(inPlaceFactors.scale, factors.scale); |
| VERIFY_IS_EQUAL(inPlaceFactors.invScale, factors.invScale); |
| VERIFY_IS_EQUAL(scaledInPlace, scaled); |
| |
| Matrix<Scalar, 1, 1> scaledSubnormal; |
| const auto subnormalFactors = |
| internal::safe_scaling<RealScalar>::scale_to(scaledSubnormal, subnormalInput, subnormalMaxCoeff); |
| VERIFY(subnormalFactors.scale != RealScalar(1) || subnormalFactors.invScale != RealScalar(1)); |
| VERIFY(numext::abs(scaledSubnormal(0)) > RealScalar(0)); |
| |
| // Unscaling restores the subnormal components exactly, still under FTZ. |
| Matrix<Scalar, 2, 1> restored; |
| internal::safe_scaling<RealScalar>::unscale_to(restored, scaled, maxCoeff, factors); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(restored(i), input(i))); |
| internal::safe_scaling<RealScalar>::unscale_in_place(scaled, maxCoeff, factors); |
| internal::safe_scaling<RealScalar>::unscale_in_place(scaledSubnormal, subnormalMaxCoeff, subnormalFactors); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(scaled(i), input(i))); |
| VERIFY(same_bits(scaledSubnormal(0), subnormalInput(0))); |
| |
| // A NaN, infinite or zero maxCoeff selects identity factors, and unscaling leaves the coefficients alone. |
| for (const RealScalar special : |
| {std::numeric_limits<RealScalar>::quiet_NaN(), std::numeric_limits<RealScalar>::infinity(), RealScalar(0)}) { |
| Matrix<Scalar, 2, 1> unchanged; |
| const auto specialFactors = internal::safe_scaling<RealScalar>::scale_to(unchanged, input, special); |
| internal::safe_scaling<RealScalar>::unscale_in_place(unchanged, special, specialFactors); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(unchanged(i), input(i))); |
| const auto inPlaceSpecialFactors = internal::safe_scaling<RealScalar>::scale_in_place(unchanged, special); |
| VERIFY_IS_EQUAL(inPlaceSpecialFactors.scale, RealScalar(1)); |
| VERIFY_IS_EQUAL(inPlaceSpecialFactors.invScale, RealScalar(1)); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(unchanged(i), input(i))); |
| } |
| } |
| |
| // Unscaling rounds through integer significands exactly when 0 < maxCoeff < min / eps. With floor factors, eps / 2 |
| // unscales to min / 2 at the threshold and to min / 4 one ulp below it; FTZ flushes only the former. |
| template <typename Scalar> |
| void check_unscale_recovery_threshold() { |
| using Binary = internal::binary_floating_point_traits<Scalar>; |
| using Scaling = internal::safe_scaling<Scalar>; |
| const Scalar threshold = (std::numeric_limits<Scalar>::min)() / NumTraits<Scalar>::epsilon(); |
| const Scalar belowThreshold = numext::nextafter(threshold, Scalar(0)); |
| // Read through volatile so the products are formed at run time, under FTZ. |
| volatile Scalar halfEpsilon = NumTraits<Scalar>::epsilon() / Scalar(2); |
| Matrix<Scalar, 1, 1> atThreshold, below; |
| atThreshold(0) = halfEpsilon; |
| below(0) = halfEpsilon; |
| |
| ScopedFlushToZero flushToZero; |
| Scaling::unscale_in_place(below, belowThreshold, Scaling::compute_floor_factors(belowThreshold)); |
| VERIFY(same_bits(below(0), numext::bit_cast<Scalar>(Binary::kExponentUnit >> 2))); |
| Scaling::unscale_in_place(atThreshold, threshold, Scaling::compute_floor_factors(threshold)); |
| if (flushToZero.isSupported()) VERIFY(same_bits(atThreshold(0), Scalar(0))); |
| } |
| |
| template <typename Scalar> |
| void check_identity_scaling() { |
| using RealScalar = typename NumTraits<Scalar>::Real; |
| using Binary = internal::binary_floating_point_traits<RealScalar>; |
| using Scaling = internal::safe_scaling<RealScalar>; |
| Matrix<Scalar, 2, 1> input; |
| input(0) = scaling_test_value<Scalar>::run(RealScalar(1), RealScalar(-1)); |
| input(1) = scaling_test_value<Scalar>::run(numext::bit_cast<RealScalar>(Binary::kExponentUnit >> 1), |
| numext::bit_cast<RealScalar>(+Binary::kSignBit)); |
| const ScopedFlushToZero flush; |
| for (RealScalar maxCoeff : {RealScalar(0), RealScalar(1), RealScalar(1.5), NumTraits<RealScalar>::infinity(), |
| NumTraits<RealScalar>::quiet_NaN()}) { |
| // Keep the factor selection at run time, where a compiler can merge a copy with multiplication by one. |
| volatile RealScalar runtimeMaxCoeff = maxCoeff; |
| Matrix<Scalar, 2, 1> scaled, restored; |
| const auto factors = Scaling::scale_to(scaled, input, RealScalar(runtimeMaxCoeff)); |
| VERIFY_IS_EQUAL(factors.scale, RealScalar(1)); |
| VERIFY_IS_EQUAL(factors.invScale, RealScalar(1)); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(scaled(i), input(i))); |
| Scaling::unscale_to(restored, input, factors); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(restored(i), input(i))); |
| Scaling::with_scaled(input, RealScalar(runtimeMaxCoeff), [&](const auto& expression) { scaled = expression; }); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(scaled(i), input(i))); |
| Scaling::scale_in_place(scaled, RealScalar(runtimeMaxCoeff)); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(scaled(i), input(i))); |
| Scaling::unscale_in_place(scaled, factors); |
| for (Index i = 0; i < input.size(); ++i) VERIFY(same_bits(scaled(i), input(i))); |
| } |
| } |
| |
| // A reduction that flushes subnormal inputs returns zero; the rescan reads the largest real or imaginary magnitude from |
| // the representation instead, and leaves a nonzero maximum alone. |
| template <typename Scalar> |
| void check_recover_flushed_max_coeff() { |
| using RealScalar = typename NumTraits<Scalar>::Real; |
| using Binary = internal::binary_floating_point_traits<RealScalar>; |
| using Bits = typename Binary::Bits; |
| using Scaling = internal::safe_scaling<RealScalar>; |
| using Value = scaling_test_value<Scalar>; |
| const RealScalar small = numext::bit_cast<RealScalar>(Bits(3)); |
| const RealScalar large = numext::bit_cast<RealScalar>(Bits(Binary::kFractionMask)); |
| Matrix<Scalar, 2, 2> m; |
| m << Value::run(small, RealScalar(0)), Value::run(-small, -large), Value::run(RealScalar(0), small), |
| Value::run(-small, RealScalar(0)); |
| const RealScalar expected = NumTraits<Scalar>::IsComplex ? large : small; |
| VERIFY(same_bits(Scaling::recover_flushed_max_coeff(m, RealScalar(0)), expected)); |
| VERIFY(same_bits(Scaling::recover_flushed_max_coeff(m, RealScalar(1)), RealScalar(1))); |
| VERIFY(same_bits(Scaling::recover_flushed_max_coeff(Matrix<Scalar, 2, 2>::Zero(), RealScalar(0)), RealScalar(0))); |
| } |
| |
| // A strided view reduces coefficient by coefficient, and under Arm FZ a scalar max(a, b) is a compare and a select |
| // that keeps the first subnormal operand rather than producing a zero: the rescan must not trust a subnormal maximum. |
| // Every flush-to-zero mode, with the largest magnitude last. |
| template <typename Scalar> |
| void check_recover_flushed_scalar_max_coeff() { |
| using Binary = internal::binary_floating_point_traits<Scalar>; |
| using Bits = typename Binary::Bits; |
| using Scaling = internal::safe_scaling<Scalar>; |
| using VectorType = Matrix<Scalar, Dynamic, 1>; |
| VectorType storage = VectorType::Zero(6); |
| const Bits significands[3] = {Bits(5), Bits(8), Bits(23)}; |
| for (Index i = 0; i < 3; ++i) storage(2 * i) = numext::bit_cast<Scalar>(significands[i]); |
| const Map<const VectorType, 0, InnerStride<2>> strided(storage.data(), 3); |
| forEachFlushToZeroMode([&](FlushToZeroMode) { |
| const Scalar maxCoeff = strided.cwiseAbs().maxCoeff(); |
| VERIFY(same_bits(Scaling::recover_flushed_max_coeff(strided, maxCoeff), numext::bit_cast<Scalar>(Bits(23)))); |
| VERIFY_IS_EQUAL(internal::binary_frexp_exponent(Scaling::recover_flushed_max_coeff(strided, maxCoeff)), |
| 5 + std::numeric_limits<Scalar>::min_exponent - std::numeric_limits<Scalar>::digits); |
| }); |
| } |
| |
| // Each flush-to-zero mode the host selects does what its name says, read from the results of two probes through the |
| // representation, since a comparison under DAZ reads any subnormal as zero; the modes restore the ambient state. |
| template <typename Scalar> |
| void check_flush_to_zero_modes() { |
| using Binary = internal::binary_floating_point_traits<Scalar>; |
| const auto resultsFlushed = []() { return Binary::bits(underflowProbe<Scalar>()) == 0; }; |
| const auto inputsFlushed = []() { return Binary::bits(subnormalInputProbe<Scalar>()) == 0; }; |
| const bool resultsBefore = resultsFlushed(), inputsBefore = inputsFlushed(); |
| int modes = 0; |
| forEachFlushToZeroMode([&](FlushToZeroMode mode) { |
| ++modes; |
| VERIFY_IS_EQUAL(resultsFlushed(), flushToZeroModeIncludes(mode, FlushToZeroMode::Results)); |
| VERIFY_IS_EQUAL(inputsFlushed(), flushToZeroModeIncludes(mode, FlushToZeroMode::Inputs)); |
| VERIFY_IS_EQUAL(ScopedFlushToZero::hardwareFlushesSubnormalInputs(), |
| flushToZeroModeIncludes(mode, FlushToZeroMode::Inputs)); |
| }); |
| VERIFY(modes >= 1); |
| VERIFY_IS_EQUAL(resultsFlushed(), resultsBefore); |
| VERIFY_IS_EQUAL(inputsFlushed(), inputsBefore); |
| // The default request flushes at least results, as it always did. |
| { |
| const ScopedFlushToZero flushToZero; |
| if (flushToZero.isSupported()) { |
| VERIFY(flushToZeroModeIncludes(flushToZero.mode(), FlushToZeroMode::Results)); |
| VERIFY(resultsFlushed()); |
| } |
| } |
| VERIFY_IS_EQUAL(resultsFlushed(), resultsBefore); |
| } |
| |
| // value * 2^exponent through integer significands agrees bit for bit with the correctly rounded ldexp for every input |
| // class and every exponent that changes the result, and keeps doing so in every flush-to-zero mode, where the |
| // references are taken beforehand. Beyond the exponent range it saturates to the signed zero or infinity ldexp returns. |
| template <typename Scalar> |
| void check_scale_binary_by_exponent() { |
| using Binary = internal::binary_floating_point_traits<Scalar>; |
| using Bits = typename Binary::Bits; |
| using Scaling = internal::safe_scaling<Scalar>; |
| // 2^recovery is the min / eps threshold check_unscale_recovery_threshold() exercises. |
| VERIFY_IS_EQUAL(numext::ldexp(Scalar(1), Scaling::subnormal_recovery_exponent()), |
| (std::numeric_limits<Scalar>::min)() / NumTraits<Scalar>::epsilon()); |
| |
| // Subnormals down to denorm_min, the normal boundary, full fractions, and the largest finite value. |
| const Bits magnitudes[] = {Bits(1), |
| Bits(3), |
| Binary::kExponentUnit >> 1, |
| Binary::kFractionMask, |
| Binary::kExponentUnit, |
| Binary::kExponentUnit | Bits(1), |
| Binary::kExponentUnit + (Binary::kExponentUnit >> 1), |
| (Binary::kExponentUnit << 1) - Bits(1), |
| Binary::bits(Scalar(0.75)), |
| Binary::bits(Scalar(1)), |
| Binary::bits(Scalar(1.5)), |
| Binary::bits((std::numeric_limits<Scalar>::max)())}; |
| // Every exponent from full underflow to full overflow, then values that only the clamp can handle. |
| const int span = std::numeric_limits<Scalar>::max_exponent - std::numeric_limits<Scalar>::min_exponent + |
| std::numeric_limits<Scalar>::digits + 3; |
| std::vector<int> exponents; |
| for (int e = -span; e <= span; ++e) exponents.push_back(e); |
| for (int e : {1 << 24, 1 << 30, (std::numeric_limits<int>::max)()}) { |
| exponents.push_back(e); |
| exponents.push_back(-e); |
| } |
| exponents.push_back((std::numeric_limits<int>::min)()); |
| |
| struct Case { |
| Bits value; |
| int exponent; |
| Bits expected; |
| Bits expectedUnit; // ldexp(1, exponent) |
| int frexpExponent; // of value |
| Bits frexpSignificand; |
| }; |
| std::vector<Case> cases; |
| for (Bits magnitude : magnitudes) { |
| for (int e : exponents) { |
| for (Bits sign : {Bits(0), Bits(Binary::kSignBit)}) { |
| const Scalar value = numext::bit_cast<Scalar>(sign | magnitude); |
| int frexpExponent = 0; |
| EIGEN_USING_STD(frexp); |
| const Scalar frexpSignificand = frexp(value, &frexpExponent); |
| cases.push_back(Case{sign | magnitude, e, Binary::bits(numext::ldexp(value, e)), |
| Binary::bits(numext::ldexp(Scalar(1), e)), frexpExponent, Binary::bits(frexpSignificand)}); |
| } |
| } |
| } |
| for (const Scalar special : {Scalar(0), -Scalar(0), std::numeric_limits<Scalar>::infinity(), |
| -std::numeric_limits<Scalar>::infinity(), std::numeric_limits<Scalar>::quiet_NaN()}) { |
| for (int e : {-1000, -1, 0, 1, 1000}) { |
| cases.push_back(Case{Binary::bits(special), e, Binary::bits(special), Binary::bits(numext::ldexp(Scalar(1), e)), |
| 0, Binary::bits(special)}); |
| } |
| } |
| |
| const auto check = [&]() { |
| for (const Case& c : cases) { |
| const Scalar actual = internal::scale_binary_by_exponent(numext::bit_cast<Scalar>(c.value), c.exponent); |
| VERIFY_IS_EQUAL(Binary::bits(actual), c.expected); |
| const std::complex<Scalar> complexActual = internal::scale_binary_by_exponent( |
| std::complex<Scalar>(numext::bit_cast<Scalar>(c.value), Scalar(1)), c.exponent); |
| VERIFY_IS_EQUAL(Binary::bits(complexActual.real()), c.expected); |
| VERIFY_IS_EQUAL(Binary::bits(complexActual.imag()), c.expectedUnit); |
| // The representation-based classifiers agree with the IEEE ones on every value. |
| const Scalar value = numext::bit_cast<Scalar>(c.value); |
| VERIFY_IS_EQUAL(Binary::bits(internal::ldexp_preserving_subnormals(value, c.exponent)), c.expected); |
| const Bits magnitude = c.value & ~Binary::kSignBit; |
| VERIFY_IS_EQUAL(numext::is_exactly_zero_no_flush(value), magnitude == 0); |
| VERIFY_IS_EQUAL(internal::is_subnormal_magnitude(value), magnitude != 0 && magnitude < Binary::kExponentUnit); |
| if (magnitude == c.value && |
| (numext::isfinite)(value)) { // non-negative: the maximum against every larger power of two |
| VERIFY_IS_EQUAL(Binary::bits(internal::max_preserving_subnormals(value, numext::bit_cast<Scalar>(Bits(1)))), |
| magnitude == 0 ? Bits(1) : magnitude); |
| } |
| // frexp from the representation agrees with frexp on every finite value; zeros give exponent 0. |
| if ((numext::isfinite)(value)) { |
| int frexpExponent = 0; |
| const Scalar significand = internal::binary_frexp(value, frexpExponent); |
| VERIFY_IS_EQUAL(frexpExponent, c.frexpExponent); |
| VERIFY_IS_EQUAL(Binary::bits(significand), c.frexpSignificand); |
| VERIFY_IS_EQUAL(internal::binary_frexp_exponent(value), c.frexpExponent); |
| } |
| } |
| // A normal power-of-two factor is the exponent form of the same scaling. |
| for (int e = std::numeric_limits<Scalar>::min_exponent - 1; e < std::numeric_limits<Scalar>::max_exponent; ++e) { |
| const Scalar factor = numext::bit_cast<Scalar>(Bits(e + Binary::kExponentBias) << Binary::kFractionBits); |
| for (Bits magnitude : magnitudes) { |
| const Scalar value = numext::bit_cast<Scalar>(magnitude); |
| VERIFY_IS_EQUAL(Binary::bits(internal::scale_binary_by_power_of_two(value, factor)), |
| Binary::bits(internal::scale_binary_by_exponent(value, e))); |
| } |
| } |
| // The functor is the coefficient-wise scalar form, complex components included. |
| Matrix<std::complex<Scalar>, 3, 1> input; |
| input << std::complex<Scalar>(numext::bit_cast<Scalar>(Bits(3)), numext::bit_cast<Scalar>(Binary::kFractionMask)), |
| std::complex<Scalar>(Scalar(1.5), -Scalar(0.75)), |
| std::complex<Scalar>((std::numeric_limits<Scalar>::max)(), numext::bit_cast<Scalar>(Bits(1))); |
| for (int e : {-(std::numeric_limits<Scalar>::digits + 1), -3, 5, std::numeric_limits<Scalar>::max_exponent}) { |
| const Matrix<std::complex<Scalar>, 3, 1> scaled = input.unaryExpr(internal::scale_by_exponent_op<Scalar>(e)); |
| for (Index i = 0; i < input.size(); ++i) { |
| VERIFY(same_bits(scaled(i), internal::scale_binary_by_exponent(input(i), e))); |
| } |
| } |
| }; |
| forEachFlushToZeroMode([&](FlushToZeroMode) { check(); }); |
| } |
| |
| // The scalars without a representation path take the C library. |
| void check_preserving_subnormals_fallbacks() { |
| VERIFY_IS_EQUAL(internal::ldexp_preserving_subnormals(1.5L, 3), 12.0L); |
| VERIFY_IS_EQUAL(internal::frexp_exponent_preserving_subnormals(12.0L), 4); |
| int exponent = 0; |
| VERIFY_IS_EQUAL(internal::frexp_preserving_subnormals(-12.0L, exponent), -0.75L); |
| VERIFY_IS_EQUAL(exponent, 4); |
| VERIFY_IS_EQUAL(internal::frexp_exponent_preserving_subnormals(0.0L), 0); |
| VERIFY(numext::is_exactly_zero_no_flush(0.0L)); |
| VERIFY(!internal::is_subnormal_magnitude(1.0L)); |
| VERIFY(internal::is_subnormal_magnitude(std::numeric_limits<long double>::denorm_min())); |
| VERIFY_IS_EQUAL(internal::max_preserving_subnormals(1.5L, 2.5L), 2.5L); |
| } |
| |
| EIGEN_DECLARE_TEST(safe_scaling) { |
| CALL_SUBTEST(check_flush_to_zero_modes<float>()); |
| CALL_SUBTEST(check_flush_to_zero_modes<double>()); |
| CALL_SUBTEST(check_recover_flushed_scalar_max_coeff<float>()); |
| CALL_SUBTEST(check_recover_flushed_scalar_max_coeff<double>()); |
| CALL_SUBTEST(check_scale_binary_by_exponent<float>()); |
| CALL_SUBTEST(check_scale_binary_by_exponent<double>()); |
| CALL_SUBTEST(check_preserving_subnormals_fallbacks()); |
| CALL_SUBTEST(check_identity_scaling<float>()); |
| CALL_SUBTEST(check_identity_scaling<double>()); |
| CALL_SUBTEST(check_identity_scaling<std::complex<float>>()); |
| CALL_SUBTEST(check_identity_scaling<std::complex<double>>()); |
| CALL_SUBTEST(check_power_of_two_scaling_factor<half>()); |
| CALL_SUBTEST(check_power_of_two_scaling_factor<bfloat16>()); |
| CALL_SUBTEST(check_power_of_two_scaling_factor<float>()); |
| CALL_SUBTEST(check_power_of_two_scaling_factor<double>()); |
| CALL_SUBTEST(check_power_of_two_scaling_factor<long double>()); |
| CALL_SUBTEST(check_safe_scaling_special_values<half>()); |
| CALL_SUBTEST(check_safe_scaling_special_values<bfloat16>()); |
| CALL_SUBTEST(check_safe_scaling_special_values<float>()); |
| CALL_SUBTEST(check_safe_scaling_special_values<double>()); |
| CALL_SUBTEST(check_safe_scaling_special_value_frontends()); |
| CALL_SUBTEST(check_arithmetic_safe_scaling_fallback()); |
| CALL_SUBTEST(check_custom_scalar_scaling_exceptions()); |
| CALL_SUBTEST(check_arithmetic_scaling_expression<float>()); |
| CALL_SUBTEST(check_arithmetic_scaling_expression<double>()); |
| CALL_SUBTEST(check_arithmetic_scaling_expression<std::complex<float>>()); |
| CALL_SUBTEST(check_arithmetic_scaling_expression<std::complex<double>>()); |
| CALL_SUBTEST(check_scale_binary_by_power_of_two<float>()); |
| CALL_SUBTEST(check_scale_binary_by_power_of_two<double>()); |
| CALL_SUBTEST(check_scale_binary_by_power_of_two_rounding<float>()); |
| CALL_SUBTEST(check_scale_binary_by_power_of_two_rounding<double>()); |
| CALL_SUBTEST(check_subnormal_preserving_scaling<float>()); |
| CALL_SUBTEST(check_subnormal_preserving_scaling<std::complex<float>>()); |
| CALL_SUBTEST(check_subnormal_preserving_scaling<double>()); |
| CALL_SUBTEST(check_subnormal_preserving_scaling<std::complex<double>>()); |
| CALL_SUBTEST(check_unscale_recovery_threshold<float>()); |
| CALL_SUBTEST(check_unscale_recovery_threshold<double>()); |
| CALL_SUBTEST(check_recover_flushed_max_coeff<float>()); |
| CALL_SUBTEST(check_recover_flushed_max_coeff<std::complex<float>>()); |
| CALL_SUBTEST(check_recover_flushed_max_coeff<double>()); |
| CALL_SUBTEST(check_recover_flushed_max_coeff<std::complex<double>>()); |
| } |
