Eigen/src/Core/SafeScaling.h

// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0

#ifndef EIGEN_SAFE_SCALING_H
#define EIGEN_SAFE_SCALING_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {
namespace internal {

// Select factors for a positive finite scale. Supported binary floating-point scalars use normal powers of two;
// other scalar types use division for matrix scaling and can request a clamped reciprocal for stableNorm(). General
// scaling rounds down so it does not discard a representable tail that division by the original value would preserve.
// Stable reductions separately round up to keep scaled magnitudes at most one. At the upper exponent boundary, the
// scale is clamped to keep its reciprocal normal; for the standard binary formats this leaves magnitudes below four.
template <typename Scalar>
struct supports_power_of_two_scaling
    : bool_constant<(std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value) &&
                    std::numeric_limits<Scalar>::is_iec559 && std::numeric_limits<Scalar>::radix == 2 &&
                    (sizeof(Scalar) == sizeof(numext::uint32_t) || sizeof(Scalar) == sizeof(numext::uint64_t))> {};

#if !defined(EIGEN_GPU_COMPILE_PHASE)
template <>
struct supports_power_of_two_scaling<long double> : bool_constant<std::numeric_limits<long double>::radix == 2> {};
#endif

template <>
struct supports_power_of_two_scaling<half> : true_type {};

template <>
struct supports_power_of_two_scaling<bfloat16> : true_type {};

template <typename Scalar, bool = supports_power_of_two_scaling<Scalar>::value>
struct safe_scaling;

template <typename Scalar>
struct safe_scaling_factors {
  Scalar scale = Scalar(1);
  // Division-only arithmetic scaling leaves this unused reciprocal at one.
  Scalar invScale = Scalar(1);
};

// value * 2^exponent through the integer significand of the value with bit pattern valueBits, so FTZ/DAZ hardware and
// ARMv7 NEON can flush neither a subnormal input nor a subnormal result. Exact wherever the result is representable;
// a subnormal result rounds to nearest, ties to even, as IEEE 754 multiplication does; below denorm_min / 2 the result
// is a signed zero and above the largest finite value a signed infinity. Zeros, infinities and NaNs pass through. Keep
// the integer-only ABI so compilers cannot replace the reconstruction with FTZ/DAZ-sensitive floating-point arithmetic.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE typename binary_floating_point_traits<Scalar>::Bits scale_binary_bits_by_exponent(
    const typename binary_floating_point_traits<Scalar>::Bits valueBits, int exponent) {
  using Binary = binary_floating_point_traits<Scalar>;
  using Bits = typename Binary::Bits;
  constexpr int kInfinityExponent = int(Binary::kExponentMask >> Binary::kFractionBits);

  const Bits valueExponentBits = valueBits & Binary::kExponentMask;
  const Bits fraction = valueBits & Binary::kFractionMask;
  if (valueExponentBits == Binary::kExponentMask || (valueExponentBits == 0 && fraction == 0)) return valueBits;
  const Bits sign = valueBits & Binary::kSignBit;

  Bits significand = fraction;
  int resultExponent = int(valueExponentBits >> Binary::kFractionBits);
  if (resultExponent == 0) {
    resultExponent = 1;
    while (significand < Binary::kExponentUnit) {
      significand <<= 1;
      --resultExponent;
    }
  }
  // Every finite value saturates beyond this magnitude of exponent, so the clamp changes no result.
  constexpr int kSaturatingExponent = kInfinityExponent + Binary::kFractionBits + 2;
  resultExponent += numext::mini(numext::maxi(exponent, -kSaturatingExponent), kSaturatingExponent);
  if (resultExponent >= kInfinityExponent) return sign | Binary::kExponentMask;
  if (resultExponent > 0)
    return sign | (Bits(resultExponent) << Binary::kFractionBits) | (significand & Binary::kFractionMask);

  // Subnormal result m * 2^(resultExponent - 1) denorm_min, m = significand | kExponentUnit < 2 * kExponentUnit: round
  // to nearest, ties to even, as IEEE 754 multiplication does. A carry into kExponentUnit encodes the smallest normal;
  // for shift > kFractionBits + 1 the result is below denorm_min / 2 and rounds to a signed zero.
  const int shift = 1 - resultExponent;
  if (shift > Binary::kFractionBits + 1) return sign;
  const Bits mantissa = significand | Binary::kExponentUnit;
  const Bits halfway = Bits(1) << (shift - 1);
  const Bits remainder = mantissa & ((halfway << 1) - 1);
  Bits rounded = mantissa >> shift;
  if (remainder > halfway || (remainder == halfway && (rounded & Bits(1)) != 0)) ++rounded;
  return sign | rounded;
}

// The same scaling by a positive normal power of two given as the bit pattern of its exponent field.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename binary_floating_point_traits<Scalar>::Bits
scale_binary_bits_by_power_of_two(const typename binary_floating_point_traits<Scalar>::Bits valueBits,
                                  const typename binary_floating_point_traits<Scalar>::Bits factorExponentBits) {
  using Binary = binary_floating_point_traits<Scalar>;
  return scale_binary_bits_by_exponent<Scalar>(
      valueBits, int(factorExponentBits >> Binary::kFractionBits) - Binary::kExponentBias);
}

// value * 2^exponent for the binary floating-point scalars, with the rounding and saturation of
// scale_binary_bits_by_exponent(): the ldexp that FTZ/DAZ hardware cannot flush, for any int exponent.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar scale_binary_by_exponent(const Scalar& value, int exponent) {
  using Binary = binary_floating_point_traits<Scalar>;
  return numext::bit_cast<Scalar>(scale_binary_bits_by_exponent<Scalar>(Binary::bits(value), exponent));
}

template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<Scalar> scale_binary_by_exponent(const std::complex<Scalar>& value,
                                                                                    int exponent) {
  return std::complex<Scalar>(scale_binary_by_exponent(value.real(), exponent),
                              scale_binary_by_exponent(value.imag(), exponent));
}

// frexp's exponent e, 2^(e - 1) <= |value| < 2^e, read from the representation of a finite value so that FTZ/DAZ
// hardware cannot flush a subnormal; 0 for a zero. Infinities and NaNs are not classified: they come back one above
// the largest finite exponent. The classification takes the magnitude bits out of line, like
// scale_binary_bits_by_exponent(): inlined, a compiler folds an integer test on the bit pattern of a float back into
// a floating-point comparison (`magnitude == 0` into `value == 0`), which FTZ/DAZ hardware flushes exactly like the
// arithmetic the bits were meant to bypass.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE int binary_frexp_exponent_of_magnitude(
    const typename binary_floating_point_traits<Scalar>::Bits magnitude) {
  using Binary = binary_floating_point_traits<Scalar>;
  using Bits = typename Binary::Bits;
  if (magnitude == 0) return 0;
  const Bits exponentBits = magnitude & Binary::kExponentMask;
  if (exponentBits != 0) return int(exponentBits >> Binary::kFractionBits) - Binary::kExponentBias + 1;
  // A subnormal is its fraction times 2^(min_exponent - digits).
  return log_2_impl<Bits>::run_floor(magnitude) + 1 + std::numeric_limits<Scalar>::min_exponent -
         std::numeric_limits<Scalar>::digits;
}

template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE int binary_frexp_exponent(const Scalar& value) {
  return binary_frexp_exponent_of_magnitude<Scalar>(binary_floating_point_traits<Scalar>::magnitude(value));
}

// The larger of two magnitude bit patterns, out of line for the same reason.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE typename binary_floating_point_traits<Scalar>::Bits larger_magnitude_bits(
    const typename binary_floating_point_traits<Scalar>::Bits a,
    const typename binary_floating_point_traits<Scalar>::Bits b) {
  return a < b ? b : a;
}

// |value| < the smallest normal from the magnitude bits, out of line for the same reason.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE bool is_below_normal_magnitude_bits(
    const typename binary_floating_point_traits<Scalar>::Bits magnitude) {
  return magnitude < binary_floating_point_traits<Scalar>::kExponentUnit;
}

// frexp through the representation: the significand in [0.5, 1), with the sign of value, and its exponent.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar binary_frexp(const Scalar& value, int& exponent) {
  exponent = binary_frexp_exponent(value);
  return scale_binary_by_exponent(value, -exponent);
}

// Multiplication by a positive normal power of two through integer significands, so FTZ/DAZ cannot flush a subnormal
// input or result; a subnormal result rounds as IEEE 754 multiplication does, and a result above the largest finite
// value saturates to infinity.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar scale_binary_by_power_of_two(const Scalar& value, const Scalar& factor) {
  using Binary = binary_floating_point_traits<Scalar>;
  using Bits = typename Binary::Bits;
  const Bits factorExponentBits = Binary::bits(factor);
  eigen_internal_assert(factorExponentBits != 0 && factorExponentBits < Binary::kExponentMask &&
                        (factorExponentBits & ~Binary::kExponentMask) == 0);
  return numext::bit_cast<Scalar>(scale_binary_bits_by_power_of_two<Scalar>(Binary::bits(value), factorExponentBits));
}

template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<Scalar> scale_binary_by_power_of_two(
    const std::complex<Scalar>& value, const Scalar& factor) {
  return std::complex<Scalar>(scale_binary_by_power_of_two(value.real(), factor),
                              scale_binary_by_power_of_two(value.imag(), factor));
}

// Coefficient-wise value * 2^exponent through scale_binary_by_exponent(). Scalar only: the integer path is the point.
template <typename Scalar>
struct scale_by_exponent_op {
  EIGEN_DEVICE_FUNC explicit scale_by_exponent_op(int exponent) : m_exponent(exponent) {}

  template <typename CoeffScalar>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffScalar operator()(const CoeffScalar& value) const {
    return scale_binary_by_exponent(value, m_exponent);
  }

  int m_exponent;
};

template <typename Scalar>
struct functor_traits<scale_by_exponent_op<Scalar>> {
  static constexpr int Cost = 10 * NumTraits<Scalar>::MulCost;
  static constexpr bool PacketAccess = false;
  static constexpr bool IsRepeatable = true;
};

template <typename FactorScalar>
struct scale_by_power_of_two_op {
  EIGEN_DEVICE_FUNC explicit scale_by_power_of_two_op(const FactorScalar& factor) : m_factor(factor) {}

  template <typename CoeffScalar>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffScalar operator()(const CoeffScalar& value) const {
    return scale_binary_by_power_of_two(value, m_factor);
  }

  FactorScalar m_factor;
};

template <typename FactorScalar>
struct functor_traits<scale_by_power_of_two_op<FactorScalar>> {
  // Integer significand scaling is scalar and calls an out-of-line helper.
  static constexpr int Cost = 10 * NumTraits<FactorScalar>::MulCost;
  static constexpr bool PacketAccess = false;
  static constexpr bool IsRepeatable = true;
};

template <typename FactorScalar, typename CoeffScalar>
struct use_subnormal_preserving_scaling
    : bool_constant<
          (std::is_same<FactorScalar, float>::value &&
           (std::is_same<CoeffScalar, float>::value || std::is_same<CoeffScalar, std::complex<float>>::value)) ||
          (std::is_same<FactorScalar, double>::value &&
           (std::is_same<CoeffScalar, double>::value || std::is_same<CoeffScalar, std::complex<double>>::value))> {};

// Subnormal-magnitude tests read the representation for float and double, where a comparison under DAZ
// reads a subnormal as zero as well; the other scalars compare arithmetically.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_subnormal_magnitude_impl(const Scalar& value, true_type) {
  // The smallest normal, 2^(min_exponent - 1), has frexp exponent min_exponent.
  const int exponent = binary_frexp_exponent(value);
  return exponent != 0 && exponent < std::numeric_limits<Scalar>::min_exponent;
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_subnormal_magnitude_impl(const Scalar& value, false_type) {
  const Scalar magnitude = numext::abs(value);
  return magnitude > Scalar(0) && magnitude < (std::numeric_limits<Scalar>::min)();
}
// 0 < |value| < the smallest normal.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_subnormal_magnitude(const Scalar& value) {
  return is_subnormal_magnitude_impl(value, bool_constant<use_subnormal_preserving_scaling<Scalar, Scalar>::value>());
}

// |value| < the smallest normal, zero included: a maximum that a flushing comparison may have produced. A scalar
// reduction under FTZ/DAZ does not always return a zero: where max(a, b) is a compare and a select (Arm FZ with
// scalar VFP or fcsel, ARMv7's double, which has no packet), every comparison is false and the running maximum keeps
// the bits of the first operand.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_zero_or_subnormal_magnitude_impl(const Scalar& value, true_type) {
  return is_below_normal_magnitude_bits<Scalar>(binary_floating_point_traits<Scalar>::magnitude(value));
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_zero_or_subnormal_magnitude_impl(const Scalar& value, false_type) {
  return numext::abs(value) < (std::numeric_limits<Scalar>::min)();
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_zero_or_subnormal_magnitude(const Scalar& value) {
  return is_zero_or_subnormal_magnitude_impl(value,
                                             bool_constant<use_subnormal_preserving_scaling<Scalar, Scalar>::value>());
}

// The larger of two non-negative values, such as two recovered maxima, read from the representation for float and
// double: an arithmetic maximum (a compiler may lower it to fmaxnm) flushes subnormal operands to zero under FTZ/DAZ.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar max_preserving_subnormals_impl(const Scalar& a, const Scalar& b,
                                                                            true_type) {
  using Binary = binary_floating_point_traits<Scalar>;
  return numext::bit_cast<Scalar>(larger_magnitude_bits<Scalar>(Binary::magnitude(a), Binary::magnitude(b)));
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar max_preserving_subnormals_impl(const Scalar& a, const Scalar& b,
                                                                            false_type) {
  return numext::maxi(a, b);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar max_preserving_subnormals(const Scalar& a, const Scalar& b) {
  return max_preserving_subnormals_impl(a, b, bool_constant<use_subnormal_preserving_scaling<Scalar, Scalar>::value>());
}

// |value| and value < 0 of a finite real scalar, from the representation for float and double: a comparison reads a
// negative subnormal as zero under DAZ, and an abs() that widens (MSVC 19.29 takes abs(float) through double) flushes
// a subnormal result under FTZ when it narrows back.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar abs_preserving_subnormals_impl(const Scalar& value, true_type) {
  return numext::bit_cast<Scalar>(binary_floating_point_traits<Scalar>::magnitude(value));
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar abs_preserving_subnormals_impl(const Scalar& value, false_type) {
  return numext::abs(value);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar abs_preserving_subnormals(const Scalar& value) {
  return abs_preserving_subnormals_impl(value,
                                        bool_constant<use_subnormal_preserving_scaling<Scalar, Scalar>::value>());
}

template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_negative_preserving_subnormals_impl(const Scalar& value, true_type) {
  using Binary = binary_floating_point_traits<Scalar>;
  return (Binary::bits(value) & Binary::kSignBit) != 0 && !is_zero_magnitude_bits<Scalar>(Binary::magnitude(value));
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_negative_preserving_subnormals_impl(const Scalar& value, false_type) {
  return value < Scalar(0);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_negative_preserving_subnormals(const Scalar& value) {
  return is_negative_preserving_subnormals_impl(
      value, bool_constant<use_subnormal_preserving_scaling<Scalar, Scalar>::value>());
}

// The position of the largest magnitude of a real vector with at least one coefficient, read from the representation
// and out of line like the other bit classifications.
template <typename Scalar, typename Derived>
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE Index index_of_largest_magnitude(const DenseBase<Derived>& src) {
  using Binary = binary_floating_point_traits<Scalar>;
  using Bits = typename Binary::Bits;
  const evaluator<Derived> coeffs(src.derived());
  Index position = 0;
  Bits largest = Binary::magnitude(coeffs.coeff(0));
  for (Index i = 1; i < src.size(); ++i) {
    const Bits magnitude = Binary::magnitude(coeffs.coeff(i));
    if (magnitude > largest) {
      largest = magnitude;
      position = i;
    }
  }
  return position;
}

// ldexp and frexp's exponent for a real scalar: float and double through the representation, which FTZ/DAZ cannot
// reach, the other scalars through the C library. frexp's exponent is e with 2^(e - 1) <= |value| < 2^e, 0 for a
// zero.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar ldexp_preserving_subnormals_impl(const Scalar& value, int exponent,
                                                                              true_type) {
  return scale_binary_by_exponent(value, exponent);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar ldexp_preserving_subnormals_impl(const Scalar& value, int exponent,
                                                                              false_type) {
  return numext::ldexp(value, exponent);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar ldexp_preserving_subnormals(const Scalar& value, int exponent) {
  return ldexp_preserving_subnormals_impl(value, exponent,
                                          bool_constant<use_subnormal_preserving_scaling<Scalar, Scalar>::value>());
}

template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE int frexp_exponent_preserving_subnormals_impl(const Scalar& value, true_type) {
  return binary_frexp_exponent(value);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE int frexp_exponent_preserving_subnormals_impl(const Scalar& value, false_type) {
  int exponent = 0;
  EIGEN_USING_STD(frexp);
  frexp(value, &exponent);
  return exponent;
}
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE int frexp_exponent_preserving_subnormals(const Scalar& value) {
  return frexp_exponent_preserving_subnormals_impl(
      value, bool_constant<use_subnormal_preserving_scaling<Scalar, Scalar>::value>());
}
// The significand in [0.5, 1), with the sign of value, and its exponent.
template <typename Scalar>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar frexp_preserving_subnormals(const Scalar& value, int& exponent) {
  exponent = frexp_exponent_preserving_subnormals(value);
  return ldexp_preserving_subnormals(value, -exponent);
}

template <typename Scalar, bool IsPowerOfTwo_>
struct safe_scaling_operations {
 private:
  using Factors = safe_scaling_factors<Scalar>;

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE bool needs_subnormal_recovery(const Scalar& value) {
    using Binary = binary_floating_point_traits<Scalar>;
    using Bits = typename Binary::Bits;
    constexpr Bits kThreshold = Bits(subnormal_recovery_exponent() + Binary::kExponentBias) << Binary::kFractionBits;
    return Binary::magnitude(value) - Bits(1) < kThreshold - Bits(1);
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE bool is_identity(const Factors& factors) {
    return factors.scale == Scalar(1) && factors.invScale == Scalar(1);
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE bool has_normal_reciprocal(const Scalar&, const Scalar&, false_type) {
    return false;
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE bool has_normal_reciprocal(const Scalar& value, const Scalar& normalMin,
                                                                          true_type) {
    return value >= normalMin && value <= Scalar(1) / normalMin;
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE bool is_positive_finite(const Scalar& value, false_type) {
    return value > Scalar(0) && value <= NumTraits<Scalar>::highest();
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE bool is_positive_finite(const Scalar& value, true_type) {
    using Binary = binary_floating_point_traits<Scalar>;
    const typename Binary::Bits bits = Binary::bits(value);
    return bits > 0 && bits < Binary::kExponentMask;
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Factors select_factors(const Scalar& maxCoeff) {
    using IsIeeeBinary = bool_constant<std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value>;
    if (!is_positive_finite(maxCoeff, IsIeeeBinary())) return Factors{};
    return safe_scaling<Scalar, IsPowerOfTwo_>::compute_floor_factors(maxCoeff);
  }

  template <typename Src, typename Func>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void with_scaled_impl(const Src& src, const Scalar&,
                                                                     const Factors& factors, const Func& func,
                                                                     false_type) {
    if (is_identity(factors)) {
      // Keep the copy distinct from multiplication by one, which flushes subnormals under FTZ.
      const Src* unscaled = &src;
      EIGEN_OPTIMIZATION_BARRIER(unscaled);
      func(*unscaled);
    } else EIGEN_IF_CONSTEXPR (IsPowerOfTwo_)
      func(src * factors.invScale);
    else
      func(src / factors.scale);
  }

  template <typename Src, typename Func>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void with_scaled_impl(const Src& src, const Scalar& maxCoeff,
                                                                     const Factors& factors, const Func& func,
                                                                     true_type) {
    if (!is_identity(factors) && needs_subnormal_recovery(maxCoeff))
      func(src.unaryExpr(scale_by_power_of_two_op<Scalar>(factors.invScale)));
    else
      with_scaled_impl(src, maxCoeff, factors, func, false_type());
  }

  template <typename Dest, typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void unscale_to_impl(Dest& dest, const Src& src, const Scalar&,
                                                                    const Factors& factors, false_type) {
    unscale_to(dest, src, factors);
  }

  // Below the recovery threshold, unscaling rounds coefficients that are significant relative to maxCoeff into the
  // subnormal range, where FTZ/DAZ arithmetic would zero them.
  template <typename Dest, typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void unscale_to_impl(Dest& dest, const Src& src, const Scalar& maxCoeff,
                                                                    const Factors& factors, true_type) {
    if (!is_identity(factors) && needs_subnormal_recovery(maxCoeff))
      dest = src.unaryExpr(scale_by_power_of_two_op<Scalar>(factors.scale));
    else
      unscale_to_impl(dest, src, maxCoeff, factors, false_type());
  }

  template <typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar recover_flushed_max_coeff_impl(const Src&, const Scalar& maxCoeff,
                                                                                     false_type) {
    return maxCoeff;
  }

  template <typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar recover_flushed_max_coeff_impl(const Src& src,
                                                                                     const Scalar& maxCoeff,
                                                                                     true_type) {
    using Binary = binary_floating_point_traits<Scalar>;
    using Bits = typename Binary::Bits;
    if (!is_zero_or_subnormal_magnitude(maxCoeff)) return maxCoeff;
    // A product expression has no coefficient access; its evaluator materializes it.
    const evaluator<Src> coeffs(src);
    Bits maxBits = 0;
    for (Index col = 0; col < src.cols(); ++col) {
      for (Index row = 0; row < src.rows(); ++row) {
        const typename Src::Scalar coeff = coeffs.coeff(row, col);
        maxBits = numext::maxi(
            maxBits, numext::maxi(Binary::magnitude(numext::real(coeff)), Binary::magnitude(numext::imag(coeff))));
      }
    }
    return numext::bit_cast<Scalar>(maxBits);
  }

  template <typename Derived>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Index recover_flushed_max_coeff_index_impl(const DenseBase<Derived>&,
                                                                                          Index position, false_type) {
    return position;
  }

  template <typename Derived>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Index recover_flushed_max_coeff_index_impl(const DenseBase<Derived>& src,
                                                                                          Index, true_type) {
    return index_of_largest_magnitude<Scalar>(src);
  }

 public:
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Factors
  compute_ceiling_factors_with_normal_reciprocal(const Scalar& value) {
    Factors factors = safe_scaling<Scalar, IsPowerOfTwo_>::compute_ceiling_factors(value);
    EIGEN_IF_CONSTEXPR (!IsPowerOfTwo_) {
      if (factors.invScale > NumTraits<Scalar>::highest()) {
        factors.invScale = NumTraits<Scalar>::highest();
        factors.scale = Scalar(1) / factors.invScale;
      }
    }
    return factors;
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE bool try_compute_ceiling_factors_with_normal_reciprocal(
      const Scalar& value, const Scalar& normalMin, Factors& factors) {
    if (!IsPowerOfTwo_ &&
        !has_normal_reciprocal(value, normalMin, bool_constant<std::is_floating_point<Scalar>::value>()))
      return false;
    factors = compute_ceiling_factors_with_normal_reciprocal(value);
    return true;
  }

  // Dispatch once per expression; ordinary scaling retains packet access. The callback consumes the lazy
  // expression before its operands expire, without requiring a common C++ type for the different paths.
  template <typename Src, typename Func>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Factors with_scaled(const Src& src, const Scalar& maxCoeff,
                                                                   const Func& func) {
    const Factors factors = select_factors(maxCoeff);
    constexpr bool kPreserveSubnormalInputs =
        IsPowerOfTwo_ && use_subnormal_preserving_scaling<Scalar, typename Src::Scalar>::value;
    with_scaled_impl(src, maxCoeff, factors, func, bool_constant<kPreserveSubnormalInputs>());
    return factors;
  }

  template <typename Dest, typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void unscale_to(Dest&& dest, const Src& src, const Factors& factors) {
    if (factors.scale == Scalar(1)) {
      const Src* unscaled = &src;
      EIGEN_OPTIMIZATION_BARRIER(unscaled);
      dest = *unscaled;
      return;
    }
    dest = src * factors.scale;
  }

  template <typename ValueType>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void unscale_in_place(ValueType& value, const Factors& factors) {
    if (factors.scale == Scalar(1)) return;
    value *= factors.scale;
  }

  template <typename MatrixType>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void unscale_in_place(MatrixType& matrix, const Scalar& maxCoeff,
                                                                     const Factors& factors) {
    if (factors.scale == Scalar(1)) return;
    unscale_to(matrix, matrix, maxCoeff, factors);
  }

  template <typename Dest, typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void unscale_to(Dest&& dest, const Src& src, const Scalar& maxCoeff,
                                                               const Factors& factors) {
    constexpr bool kPreserveSubnormalOutputs =
        IsPowerOfTwo_ && use_subnormal_preserving_scaling<Scalar, typename Src::Scalar>::value;
    unscale_to_impl(dest, src, maxCoeff, factors, bool_constant<kPreserveSubnormalOutputs>());
  }

  // Below 2^subnormal_recovery_exponent(), a coefficient within a factor 2^-kFractionBits of the largest one can be
  // subnormal, so a scaling that must survive FTZ/DAZ goes through integer significands (scale_by_exponent_op,
  // scale_by_power_of_two_op) instead of packet multiplication.
  EIGEN_DEVICE_FUNC static constexpr int subnormal_recovery_exponent() {
    return std::numeric_limits<Scalar>::min_exponent + std::numeric_limits<Scalar>::digits - 2;
  }

  // A SIMD unit that flushes subnormal inputs (ARMv7 NEON, Arm FZ, DAZ) reduces an all-subnormal matrix to a zero
  // maximum, and a flushed scalar reduction to whichever subnormal it held first, so rescan a zero or subnormal
  // maxCoeff from the representation. Any subnormal maximum, including the largest real or imaginary magnitude of a
  // complex matrix, selects the same factors as the true one.
  template <typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar recover_flushed_max_coeff(const Src& src,
                                                                                const Scalar& maxCoeff) {
    constexpr bool kPreserveSubnormalInputs =
        IsPowerOfTwo_ && use_subnormal_preserving_scaling<Scalar, typename Src::Scalar>::value;
    return recover_flushed_max_coeff_impl(src, maxCoeff, bool_constant<kPreserveSubnormalInputs>());
  }

  // The position of the largest magnitude in a real vector whose maxCoeff(&position) read zero or subnormal,
  // rescanned from the representation like recover_flushed_max_coeff(); the position maxCoeff() reported where there
  // is no representation path.
  template <typename Derived>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Index recover_flushed_max_coeff_index(const DenseBase<Derived>& src,
                                                                                     Index position) {
    constexpr bool kPreserveSubnormalInputs =
        IsPowerOfTwo_ && use_subnormal_preserving_scaling<Scalar, typename Derived::Scalar>::value;
    return recover_flushed_max_coeff_index_impl(src, position, bool_constant<kPreserveSubnormalInputs>());
  }

  template <typename MatrixType>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Factors scale_in_place(MatrixType& matrix, const Scalar& maxCoeff) {
    const Factors factors = select_factors(maxCoeff);
    scale_in_place(matrix, maxCoeff, factors);
    return factors;
  }

  template <typename MatrixType>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void scale_in_place(MatrixType& matrix, const Scalar& maxCoeff,
                                                                   const Factors& factors) {
    if (!is_identity(factors)) scale_to(matrix, matrix, maxCoeff, factors);
  }

  template <typename Dest, typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE void scale_to(Dest& dest, const Src& src, const Scalar& maxCoeff,
                                                             const Factors& factors) {
    constexpr bool kPreserveSubnormalInputs =
        IsPowerOfTwo_ && use_subnormal_preserving_scaling<Scalar, typename Src::Scalar>::value;
    with_scaled_impl(
        src, maxCoeff, factors, [&](const auto& scaled) { dest = scaled; }, bool_constant<kPreserveSubnormalInputs>());
  }

  template <typename Dest, typename Src>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Factors scale_to(Dest& dest, const Src& src, const Scalar& maxCoeff) {
    const Factors factors = select_factors(maxCoeff);
    scale_to(dest, src, maxCoeff, factors);
    return factors;
  }
};

template <typename Scalar, bool>
struct safe_scaling : safe_scaling_operations<Scalar, false> {
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE safe_scaling_factors<Scalar> compute_floor_factors(const Scalar& value) {
    return {value, Scalar(1)};
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE safe_scaling_factors<Scalar> compute_ceiling_factors(
      const Scalar& value) {
    return {value, Scalar(1) / value};
  }
};

template <typename Scalar>
struct safe_scaling<Scalar, true> : safe_scaling_operations<Scalar, true> {
 private:
  using Binary = binary_floating_point_traits<Scalar>;
  using Bits = typename Binary::Bits;

  template <bool RoundUp>
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE safe_scaling_factors<Scalar> compute_power_of_two_factors(
      const Scalar& value) {
    constexpr Bits kMaxScale = Binary::kMaxFiniteExponentBits - Binary::kExponentUnit;
    Bits scaleBits = RoundUp ? numext::bit_cast<Bits>(numext::ceil_power_of_two(value))
                             : Binary::bits(value) & Binary::kExponentMask;
    if (scaleBits < Binary::kExponentUnit) scaleBits = Binary::kExponentUnit;
    if (scaleBits > kMaxScale) scaleBits = kMaxScale;
    const Bits invScaleBits = Bits(Binary::kMaxFiniteExponentBits - scaleBits);
    return {numext::bit_cast<Scalar>(scaleBits), numext::bit_cast<Scalar>(invScaleBits)};
  }

 public:
  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE safe_scaling_factors<Scalar> compute_floor_factors(const Scalar& value) {
    return compute_power_of_two_factors<false>(value);
  }

  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE safe_scaling_factors<Scalar> compute_ceiling_factors(
      const Scalar& value) {
    return compute_power_of_two_factors<true>(value);
  }
};

#if !defined(EIGEN_GPU_COMPILE_PHASE)
template <>
struct safe_scaling<long double, true> : safe_scaling_operations<long double, true> {
  static EIGEN_STRONG_INLINE safe_scaling_factors<long double> compute_floor_factors(const long double& value) {
    EIGEN_USING_STD(frexp);
    int exponent = 0;
    frexp(value, &exponent);
    return compute_factors_from_exponent(exponent - 1);
  }

  static EIGEN_STRONG_INLINE safe_scaling_factors<long double> compute_ceiling_factors(const long double& value) {
    EIGEN_USING_STD(frexp);
    int exponent = 0;
    const long double fraction = frexp(value, &exponent);
    if (numext::equal_strict(fraction, 0.5L)) --exponent;
    return compute_factors_from_exponent(exponent);
  }

 private:
  static EIGEN_STRONG_INLINE safe_scaling_factors<long double> compute_factors_from_exponent(int exponent) {
    constexpr int kMinNormalExponent = std::numeric_limits<long double>::min_exponent - 1;
    constexpr int kMaxNormalExponent = std::numeric_limits<long double>::max_exponent - 1;
    constexpr int kMinScaleExponent =
        kMinNormalExponent > -kMaxNormalExponent ? kMinNormalExponent : -kMaxNormalExponent;
    constexpr int kMaxScaleExponent =
        kMaxNormalExponent < -kMinNormalExponent ? kMaxNormalExponent : -kMinNormalExponent;
    if (exponent < kMinScaleExponent) exponent = kMinScaleExponent;
    if (exponent > kMaxScaleExponent) exponent = kMaxScaleExponent;

    EIGEN_USING_STD(ldexp);
    const long double scale = ldexp(1.0L, exponent);
    return {scale, 1.0L / scale};
  }
};
#endif

}  // namespace internal
}  // namespace Eigen

#endif  // EIGEN_SAFE_SCALING_H