Eigen/src/Core/arch/Default/Half.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
//
// The conversion routines are Copyright (c) Fabian Giesen, 2016.
// The original license follows:
//
// Copyright (c) Fabian Giesen, 2016
// All rights reserved.
// Redistribution and use in source and binary forms, with or without
// modification, are permitted.
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
// SPDX-License-Identifier: MPL-2.0

// Standard 16-bit float type, mostly useful for GPUs. Defines a new
// type Eigen::half (inheriting either from CUDA's or HIP's __half struct) with
// operator overloads such that it behaves basically as an arithmetic
// type. It will be quite slow on CPUs (so it is recommended to stay
// in fp32 for CPUs, except for simple parameter conversions, I/O
// to disk and the likes), but fast on GPUs.

#ifndef EIGEN_HALF_H
#define EIGEN_HALF_H

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

// When compiling with GPU support, the "__half_raw" base class as well as
// some other routines are defined in the GPU compiler header files
// (cuda_fp16.h, hip_fp16.h), and they are not tagged constexpr
// As a consequence, we get compile failures when compiling Eigen with
// GPU support. Hence the need to disable EIGEN_CONSTEXPR when building
// Eigen with GPU support.
// Any functions that require `numext::bit_cast` may also not be constexpr,
// including any native types when setting via raw bit values.
#if defined(EIGEN_HAS_GPU_FP16) || EIGEN_HAS_ARM64_FP16 || defined(EIGEN_HAS_BUILTIN_FLOAT16)
#define _EIGEN_MAYBE_CONSTEXPR
#else
#define _EIGEN_MAYBE_CONSTEXPR constexpr
#endif

#define F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, METHOD)                                                  \
  template <>                                                                                              \
  EIGEN_UNUSED EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC PACKET_F16 METHOD<PACKET_F16>(const PACKET_F16& _x) { \
    return float2half(METHOD<PACKET_F>(half2float(_x)));                                                   \
  }

#define EIGEN_INSTANTIATE_GENERIC_MATH_FUNCS_F16(PACKET_F, PACKET_F16) \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pcos)                      \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, psin)                      \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, psinh)                     \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pcosh)                     \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pasinh)                    \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pacosh)                    \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pexp)                      \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pexp2)                     \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pexpm1)                    \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, plog)                      \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, plog1p)                    \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, plog2)                     \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, plog10)                    \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, preciprocal)               \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, prsqrt)                    \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pcbrt)                     \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, psqrt)                     \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, ptanh)

// F16 wrappers for contrib/SpecialFunctions.
#define EIGEN_INSTANTIATE_SPECIAL_FUNCS_F16(PACKET_F, PACKET_F16) \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, perf)                 \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pndtri)

#define EIGEN_INSTANTIATE_BESSEL_FUNCS_F16(PACKET_F, PACKET_F16) \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_i0)          \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_i0e)         \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_i1)          \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_i1e)         \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_j0)          \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_j1)          \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_k0)          \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_k0e)         \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_k1)          \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_k1e)         \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_y0)          \
  F16_PACKET_FUNCTION(PACKET_F, PACKET_F16, pbessel_y1)

namespace Eigen {

struct half;

namespace half_impl {

// We want to use the __half_raw struct from the HIP header file only during the device compile phase.
// This is required because of a quirk in the way TensorFlow GPU builds are done.
// When compiling TensorFlow source code with GPU support, files that
//  * contain GPU kernels (i.e. *.cu.cc files) are compiled via hipcc
//  * do not contain GPU kernels ( i.e. *.cc files) are compiled via gcc (typically)
//
// Tensorflow uses the Eigen::half type as its FP16 type, and there are functions that
//  * are defined in a file that gets compiled via hipcc AND
//  * have Eigen::half as a pass-by-value argument AND
//  * are called in a file that gets compiled via gcc
//
// In the scenario described above the caller and callee will see different versions
// of the Eigen::half base class __half_raw, and they will be compiled by different compilers
//
// There appears to be an ABI mismatch between gcc and clang (which is called by hipcc) that results in
// the callee getting corrupted values for the Eigen::half argument.
//
// Making the host side compile phase of hipcc use the same Eigen::half impl, as the gcc compile, resolves
// this error, and hence the following convoluted #if condition
#if !defined(EIGEN_GPUCC) || !defined(EIGEN_GPU_COMPILE_PHASE)

// Make our own __half_raw definition that is similar to CUDA's.
struct __half_raw {
  struct construct_from_rep_tag {};
#if (defined(EIGEN_GPUCC) && !defined(EIGEN_GPU_COMPILE_PHASE))
  // Eigen::half can be used as the datatype for shared memory declarations (in Eigen and TF)
  // The element type for shared memory cannot have non-trivial constructors
  // and hence the following special casing (which skips the zero-initialization).
  // Note that this check gets done even in the host compilation phase, and
  // hence the need for this
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR __half_raw() {}
#else
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR __half_raw() : x(0) {}
#endif

#if EIGEN_HAS_ARM64_FP16
  explicit EIGEN_DEVICE_FUNC __half_raw(numext::uint16_t raw) : x(numext::bit_cast<__fp16>(raw)) {}
  EIGEN_DEVICE_FUNC constexpr __half_raw(construct_from_rep_tag, __fp16 rep) : x{rep} {}
  __fp16 x;
#elif defined(EIGEN_HAS_BUILTIN_FLOAT16)
  explicit EIGEN_DEVICE_FUNC __half_raw(numext::uint16_t raw) : x(numext::bit_cast<_Float16>(raw)) {}
  EIGEN_DEVICE_FUNC constexpr __half_raw(construct_from_rep_tag, _Float16 rep) : x{rep} {}
  _Float16 x;
#else
  explicit EIGEN_DEVICE_FUNC constexpr __half_raw(numext::uint16_t raw) : x(raw) {}
  EIGEN_DEVICE_FUNC constexpr __half_raw(construct_from_rep_tag, numext::uint16_t rep) : x{rep} {}
  numext::uint16_t x;
#endif
};

#elif defined(EIGEN_HIPCC)
// HIP GPU compile phase: nothing to do here.
// HIP fp16 header file has a definition for __half_raw
#elif defined(EIGEN_CUDACC)

// CUDA GPU compile phase.

#elif defined(EIGEN_SYCL_DEVICE_ONLY)
using __half_raw = cl::sycl::half;
#endif

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR __half_raw raw_uint16_to_half(numext::uint16_t x);
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC numext::uint16_t raw_half_as_uint16(const __half_raw& h);
template <bool AssumeArgumentIsNormalOrInfinityOrZero = false>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC __half_raw float_to_half_rtne(float ff);
template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC __half_raw float_to_half_rtne<false>(float ff);
template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC __half_raw float_to_half_rtne<true>(float ff);
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC float half_to_float(__half_raw h);

// Narrowing a wider floating-point value or an integer of more than 24 bits through float rounds twice. Rounding to
// odd in the first step keeps the second rounding, to a format at least two bits narrower than float, correct.
// Paths: 1 for wider floating-point types; 2 for integers of at most 32 bits, rounded to odd on the multiples of 2^8
// where |x| >= 2^24, which float holds exactly and which keep a second rounding to p <= 15 bits correct (its midpoints
// there are multiples of 2^(24 - p) >= 2^9), with no double arithmetic in device code; 3 for 64-bit integers, rounded
// to odd on the multiples of 2^11 where |x| >= 2^53, which double holds exactly; 0 casts directly.
template <typename T>
struct narrow_to_float_path
    : std::integral_constant<int, std::is_floating_point<T>::value ? (sizeof(T) > sizeof(float) ? 1 : 0)
                                  : !std::is_integral<T>::value || std::numeric_limits<T>::digits <= 24 ? 0
                                  : std::numeric_limits<T>::digits <= 32                                ? 2
                                  : std::numeric_limits<T>::digits <= 64                                ? 3
                                                                                                        : 0> {};

template <typename T, int Path = narrow_to_float_path<T>::value>
struct narrow_to_float_impl {
  EIGEN_DEVICE_FUNC static float run(const T& val) { return static_cast<float>(val); }
};
template <typename T>
struct narrow_to_float_impl<T, 1> {
  EIGEN_DEVICE_FUNC static float run(const T& val) {
    const float f = static_cast<float>(val);
    if (static_cast<T>(f) == val || (numext::isnan)(val)) return f;
    // Inexact: truncate toward zero, then set the last bit.
    const bool away = numext::abs(static_cast<T>(f)) > numext::abs(val);
    return numext::bit_cast<float>((numext::bit_cast<numext::uint32_t>(f) - (away ? 1u : 0u)) | 1u);
  }
};
// In two's complement, clearing the low bits rounds toward -inf for either sign; setting the next bit then picks the
// odd neighbor.
template <typename T>
struct narrow_to_float_impl<T, 2> {
  EIGEN_DEVICE_FUNC static float run(const T& val) {
    using U = std::make_unsigned_t<T>;
    const U u = static_cast<U>(val);
    const U low = 0xff;
    const U odd = (u & ~low) | ((u & low) != 0 ? U(0x100) : U(0));
    const U magnitude = std::is_signed<T>::value && (u >> (sizeof(U) * CHAR_BIT - 1)) ? U(0) - u : u;
    return static_cast<float>(magnitude >= (U(1) << 24) ? numext::bit_cast<T>(odd) : val);
  }
};
template <typename T>
struct narrow_to_float_impl<T, 3> {
  EIGEN_DEVICE_FUNC static float run(const T& val) {
    using U = std::make_unsigned_t<T>;
    const U u = static_cast<U>(val);
    const U low = 0x7ff;
    const U odd = (u & ~low) | ((u & low) != 0 ? U(0x800) : U(0));
    const U magnitude = std::is_signed<T>::value && (u >> (sizeof(U) * CHAR_BIT - 1)) ? U(0) - u : u;
    return narrow_to_float_impl<double>::run(
        static_cast<double>(magnitude >= (U(1) << 53) ? numext::bit_cast<T>(odd) : val));
  }
};
template <typename T>
EIGEN_DEVICE_FUNC inline float narrow_to_float(const T& val) {
  return narrow_to_float_impl<T>::run(val);
}

struct half_base : public __half_raw {
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half_base() {}
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half_base(const __half_raw& h) : __half_raw(h) {}

#if defined(EIGEN_GPUCC)
#if defined(EIGEN_HIPCC)
  // Delegate to raw_uint16_to_half, which reinterprets the raw bits for every storage type of
  // __half_raw::x. In the host compile phase on platforms with a native fp16 type (e.g. __fp16 on
  // arm64), a direct "x = __half_as_ushort(h)" would perform a numeric integer-to-float conversion.
  // numext::bit_cast is used to extract the bits because hip_fp16.h defines __half_as_ushort for
  // the device compile phase only (host references fail to link against ROCm 6.3).
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half_base(const __half& h)
      : __half_raw(raw_uint16_to_half(numext::bit_cast<numext::uint16_t>(h))) {}
#elif defined(EIGEN_CUDACC)
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half_base(const __half& h) : __half_raw(*(__half_raw*)&h) {}
#endif
#endif
};

}  // namespace half_impl

// Class definition.
struct half : public half_impl::half_base {
  // Writing this out as separate #if-else blocks to make the code easier to follow
  // The same applies to most #if-else blocks in this file
#if !defined(EIGEN_GPUCC) || !defined(EIGEN_GPU_COMPILE_PHASE)
  // Use the same base class for the following two scenarios
  // * when compiling without GPU support enabled
  // * during host compile phase when compiling with GPU support enabled
  using __half_raw = half_impl::__half_raw;
#elif defined(EIGEN_HIPCC)
  // Nothing to do here
  // HIP fp16 header file has a definition for __half_raw
#elif defined(EIGEN_CUDACC)
  // Nothing to do here.
#endif

  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half() {}

  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half(const __half_raw& h) : half_impl::half_base(h) {}

#if defined(EIGEN_GPUCC)
#if defined(EIGEN_HIPCC)
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half(const __half& h) : half_impl::half_base(h) {}
#elif defined(EIGEN_CUDACC)
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half(const __half& h) : half_impl::half_base(h) {}
#endif
#endif

// In the device compile phase __half_raw is the vendor type, which has no construct_from_rep_tag,
// so these constructors are restricted to the host compile phase and non-GPU builds.
#if EIGEN_HAS_ARM64_FP16 && !defined(EIGEN_GPU_COMPILE_PHASE)
  explicit EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half(__fp16 b)
      : half(__half_raw(__half_raw::construct_from_rep_tag(), b)) {}
#elif defined(EIGEN_HAS_BUILTIN_FLOAT16) && !defined(EIGEN_GPU_COMPILE_PHASE)
  explicit EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half(_Float16 b)
      : half(__half_raw(__half_raw::construct_from_rep_tag(), b)) {}
#endif

  explicit EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR half(bool b)
      : half_impl::half_base(half_impl::raw_uint16_to_half(b ? 0x3c00 : 0)) {}
  template <class T, std::enable_if_t<!std::is_integral<T>::value, int> = 0>
  explicit EIGEN_DEVICE_FUNC half(T val)
      : half_impl::half_base(half_impl::float_to_half_rtne<false>(half_impl::narrow_to_float(val))) {}
  // Float rounds only integers beyond 2^24, which half overflows either way.
  template <class T, std::enable_if_t<std::is_integral<T>::value, int> = 0>
  explicit EIGEN_DEVICE_FUNC half(T val)
      : half_impl::half_base(half_impl::float_to_half_rtne<true>(static_cast<float>(val))) {}
  explicit EIGEN_DEVICE_FUNC half(float f) : half_impl::half_base(half_impl::float_to_half_rtne<false>(f)) {}

  // Following the convention of numpy, converting between complex and
  // float will lead to loss of imag value.
  template <typename RealScalar>
  explicit EIGEN_DEVICE_FUNC half(std::complex<RealScalar> c)
      : half_impl::half_base(half_impl::float_to_half_rtne<false>(half_impl::narrow_to_float(c.real()))) {}

  EIGEN_DEVICE_FUNC operator float() const {  // NOLINT: Allow implicit conversion to float, because it is lossless.
    return half_impl::half_to_float(*this);
  }

#if defined(EIGEN_HAS_GPU_FP16) && !defined(EIGEN_GPU_COMPILE_PHASE)
  EIGEN_DEVICE_FUNC operator __half() const {
    ::__half_raw hr;
    // raw_half_as_uint16 reinterprets the raw bits for every storage type of __half_raw::x.
    // A direct "hr.x = x" would perform a numeric float-to-integer conversion when x has a
    // native fp16 type (e.g. __fp16 on arm64), since the vendor ::__half_raw::x is an integer.
    hr.x = half_impl::raw_half_as_uint16(*this);
    return __half(hr);
  }
#endif
};

// TODO(majnemer): Get rid of this once we can rely on C++17 inline variables do
// solve the ODR issue.
namespace half_impl {
template <typename = void>
struct numeric_limits_half_impl {
  static constexpr const bool is_specialized = true;
  static constexpr const bool is_signed = true;
  static constexpr const bool is_integer = false;
  static constexpr const bool is_exact = false;
  static constexpr const bool has_infinity = true;
  static constexpr const bool has_quiet_NaN = true;
  static constexpr const bool has_signaling_NaN = true;
  EIGEN_DIAGNOSTICS(push)
  EIGEN_DISABLE_DEPRECATED_WARNING
  static constexpr const std::float_denorm_style has_denorm = std::denorm_present;
  static constexpr const bool has_denorm_loss = false;
  EIGEN_DIAGNOSTICS(pop)
  static constexpr const std::float_round_style round_style = std::round_to_nearest;
  static constexpr const bool is_iec559 = true;
  // The C++ standard defines this as "true if the set of values representable
  // by the type is finite." Half has finite precision.
  static constexpr const bool is_bounded = true;
  static constexpr const bool is_modulo = false;
  static constexpr const int digits = 11;
  static constexpr const int digits10 =
      3;  // according to http://half.sourceforge.net/structstd_1_1numeric__limits_3_01half__float_1_1half_01_4.html
  static constexpr const int max_digits10 =
      5;  // according to http://half.sourceforge.net/structstd_1_1numeric__limits_3_01half__float_1_1half_01_4.html
  static constexpr const int radix = std::numeric_limits<float>::radix;
  static constexpr const int min_exponent = -13;
  static constexpr const int min_exponent10 = -4;
  static constexpr const int max_exponent = 16;
  static constexpr const int max_exponent10 = 4;
  static constexpr const bool traps = std::numeric_limits<float>::traps;
  // IEEE754: "The implementer shall choose how tininess is detected, but shall
  // detect tininess in the same way for all operations in radix two"
  static constexpr const bool tinyness_before = std::numeric_limits<float>::tinyness_before;

  static _EIGEN_MAYBE_CONSTEXPR Eigen::half(min)() { return Eigen::half_impl::raw_uint16_to_half(0x0400); }
  static _EIGEN_MAYBE_CONSTEXPR Eigen::half lowest() { return Eigen::half_impl::raw_uint16_to_half(0xfbff); }
  static _EIGEN_MAYBE_CONSTEXPR Eigen::half(max)() { return Eigen::half_impl::raw_uint16_to_half(0x7bff); }
  static _EIGEN_MAYBE_CONSTEXPR Eigen::half epsilon() { return Eigen::half_impl::raw_uint16_to_half(0x1400); }
  static _EIGEN_MAYBE_CONSTEXPR Eigen::half round_error() { return Eigen::half_impl::raw_uint16_to_half(0x3800); }
  static _EIGEN_MAYBE_CONSTEXPR Eigen::half infinity() { return Eigen::half_impl::raw_uint16_to_half(0x7c00); }
  static _EIGEN_MAYBE_CONSTEXPR Eigen::half quiet_NaN() { return Eigen::half_impl::raw_uint16_to_half(0x7e00); }
  static _EIGEN_MAYBE_CONSTEXPR Eigen::half signaling_NaN() { return Eigen::half_impl::raw_uint16_to_half(0x7d00); }
  static _EIGEN_MAYBE_CONSTEXPR Eigen::half denorm_min() { return Eigen::half_impl::raw_uint16_to_half(0x0001); }
};

// Redundant out-of-class definitions are required pre-C++17 but deprecated since.
#if EIGEN_COMP_CXXVER < 17
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::is_specialized;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::is_signed;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::is_integer;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::is_exact;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::has_infinity;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::has_quiet_NaN;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::has_signaling_NaN;
EIGEN_DIAGNOSTICS(push)
EIGEN_DISABLE_DEPRECATED_WARNING
template <typename T>
constexpr const std::float_denorm_style numeric_limits_half_impl<T>::has_denorm;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::has_denorm_loss;
EIGEN_DIAGNOSTICS(pop)
template <typename T>
constexpr const std::float_round_style numeric_limits_half_impl<T>::round_style;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::is_iec559;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::is_bounded;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::is_modulo;
template <typename T>
constexpr const int numeric_limits_half_impl<T>::digits;
template <typename T>
constexpr const int numeric_limits_half_impl<T>::digits10;
template <typename T>
constexpr const int numeric_limits_half_impl<T>::max_digits10;
template <typename T>
constexpr const int numeric_limits_half_impl<T>::radix;
template <typename T>
constexpr const int numeric_limits_half_impl<T>::min_exponent;
template <typename T>
constexpr const int numeric_limits_half_impl<T>::min_exponent10;
template <typename T>
constexpr const int numeric_limits_half_impl<T>::max_exponent;
template <typename T>
constexpr const int numeric_limits_half_impl<T>::max_exponent10;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::traps;
template <typename T>
constexpr const bool numeric_limits_half_impl<T>::tinyness_before;
#endif
}  // end namespace half_impl
}  // end namespace Eigen

namespace std {
// If std::numeric_limits<T> is specialized, should also specialize
// std::numeric_limits<const T>, std::numeric_limits<volatile T>, and
// std::numeric_limits<const volatile T>
// https://stackoverflow.com/a/16519653/
template <>
class numeric_limits<Eigen::half> : public Eigen::half_impl::numeric_limits_half_impl<> {};
template <>
class numeric_limits<const Eigen::half> : public numeric_limits<Eigen::half> {};
template <>
class numeric_limits<volatile Eigen::half> : public numeric_limits<Eigen::half> {};
template <>
class numeric_limits<const volatile Eigen::half> : public numeric_limits<Eigen::half> {};
}  // end namespace std

namespace Eigen {

namespace half_impl {

#if defined(EIGEN_GPU_COMPILE_PHASE)
// Note: We deliberately do *not* define this to 1 even if we have Arm's native
// fp16 type since GPU half types are rather different from native CPU half types.
#define EIGEN_HAS_NATIVE_GPU_FP16
#endif

// Intrinsics for native fp16 support. Note that on current hardware,
// these are no faster than fp32 arithmetic (you need to use the half2
// versions to get the ALU speed increased), but you do save the
// conversion steps back and forth.

#if defined(EIGEN_HAS_NATIVE_GPU_FP16)
EIGEN_STRONG_INLINE __device__ half operator+(const half& a, const half& b) { return __hadd(::__half(a), ::__half(b)); }
EIGEN_STRONG_INLINE __device__ half operator*(const half& a, const half& b) { return __hmul(a, b); }
EIGEN_STRONG_INLINE __device__ half operator-(const half& a, const half& b) { return __hsub(a, b); }
EIGEN_STRONG_INLINE __device__ half operator/(const half& a, const half& b) { return __hdiv(a, b); }
EIGEN_STRONG_INLINE __device__ half operator-(const half& a) { return __hneg(a); }
EIGEN_STRONG_INLINE __device__ half& operator+=(half& a, const half& b) {
  a = a + b;
  return a;
}
EIGEN_STRONG_INLINE __device__ half& operator*=(half& a, const half& b) {
  a = a * b;
  return a;
}
EIGEN_STRONG_INLINE __device__ half& operator-=(half& a, const half& b) {
  a = a - b;
  return a;
}
EIGEN_STRONG_INLINE __device__ half& operator/=(half& a, const half& b) {
  a = a / b;
  return a;
}
EIGEN_STRONG_INLINE __device__ bool operator==(const half& a, const half& b) { return __heq(a, b); }
// __hne is the ordered comparison, false for a NaN operand; != must be true there.
EIGEN_STRONG_INLINE __device__ bool operator!=(const half& a, const half& b) { return !__heq(a, b); }
EIGEN_STRONG_INLINE __device__ bool operator<(const half& a, const half& b) { return __hlt(a, b); }
EIGEN_STRONG_INLINE __device__ bool operator<=(const half& a, const half& b) { return __hle(a, b); }
EIGEN_STRONG_INLINE __device__ bool operator>(const half& a, const half& b) { return __hgt(a, b); }
EIGEN_STRONG_INLINE __device__ bool operator>=(const half& a, const half& b) { return __hge(a, b); }

#endif  // EIGEN_HAS_NATIVE_GPU_FP16

#if (EIGEN_HAS_ARM64_FP16 || defined(EIGEN_HAS_BUILTIN_FLOAT16)) && !defined(EIGEN_GPU_COMPILE_PHASE)
// nvcc's EDG front end does not promote __fp16 arithmetic to float, and it ranks the mandatory
// __fp16 -> float promotion as an exact match, so constructing a half from a native fp16 expression
// is an ambiguous tie between half(__fp16) and half(float). Every nvcc arm64 translation unit that
// includes Eigen then fails to compile. Building __half_raw names one constructor, and is what
// half(__fp16) and half(_Float16) do internally, so the stored value is unchanged.
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half half_from_rep(decltype(__half_raw::x) rep) {
  return half(__half_raw(__half_raw::construct_from_rep_tag(), rep));
}
#endif

#if defined(EIGEN_HAS_ARM64_FP16_SCALAR_ARITHMETIC) && !defined(EIGEN_GPU_COMPILE_PHASE)
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator+(const half& a, const half& b) {
  return half_from_rep(vaddh_f16(a.x, b.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator*(const half& a, const half& b) {
  return half_from_rep(vmulh_f16(a.x, b.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator-(const half& a, const half& b) {
  return half_from_rep(vsubh_f16(a.x, b.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator/(const half& a, const half& b) {
  return half_from_rep(vdivh_f16(a.x, b.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator-(const half& a) { return half_from_rep(vnegh_f16(a.x)); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator+=(half& a, const half& b) {
  a = half_from_rep(vaddh_f16(a.x, b.x));
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator*=(half& a, const half& b) {
  a = half_from_rep(vmulh_f16(a.x, b.x));
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator-=(half& a, const half& b) {
  a = half_from_rep(vsubh_f16(a.x, b.x));
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator/=(half& a, const half& b) {
  a = half_from_rep(vdivh_f16(a.x, b.x));
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator==(const half& a, const half& b) { return vceqh_f16(a.x, b.x); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator!=(const half& a, const half& b) { return !vceqh_f16(a.x, b.x); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator<(const half& a, const half& b) { return vclth_f16(a.x, b.x); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator<=(const half& a, const half& b) { return vcleh_f16(a.x, b.x); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator>(const half& a, const half& b) { return vcgth_f16(a.x, b.x); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator>=(const half& a, const half& b) { return vcgeh_f16(a.x, b.x); }

#elif (EIGEN_HAS_ARM64_FP16 || defined(EIGEN_HAS_BUILTIN_FLOAT16)) && !defined(EIGEN_GPU_COMPILE_PHASE)

// NOTE: While GCC appears to only apply `-Wdouble-promotion` to single-to-double
// conversions, Clang applies the warning to half-to-single implicit conversion also.  For
// the operations below, implicit conversion is intentional so that---depending on the
// target features---the compiler can choose between using native half-precision
// instructions or converting to single-precision.
#if EIGEN_COMP_CLANG
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wdouble-promotion"
#endif

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator+(const half& a, const half& b) {
  return half_from_rep(static_cast<decltype(a.x)>(a.x + b.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator*(const half& a, const half& b) {
  return half_from_rep(static_cast<decltype(a.x)>(a.x * b.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator-(const half& a, const half& b) {
  return half_from_rep(static_cast<decltype(a.x)>(a.x - b.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator/(const half& a, const half& b) {
  return half_from_rep(static_cast<decltype(a.x)>(a.x / b.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator-(const half& a) {
  return half_from_rep(static_cast<decltype(a.x)>(-a.x));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator+=(half& a, const half& b) {
  a = a + b;
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator*=(half& a, const half& b) {
  a = a * b;
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator-=(half& a, const half& b) {
  a = a - b;
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator/=(half& a, const half& b) {
  a = a / b;
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator==(const half& a, const half& b) { return a.x == b.x; }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator!=(const half& a, const half& b) { return a.x != b.x; }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator<(const half& a, const half& b) { return a.x < b.x; }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator<=(const half& a, const half& b) { return a.x <= b.x; }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator>(const half& a, const half& b) { return a.x > b.x; }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator>=(const half& a, const half& b) { return a.x >= b.x; }

#if EIGEN_COMP_CLANG
#pragma clang diagnostic pop
#endif

// We need to distinguish ‘clang as the CUDA compiler’ from ‘clang as the host compiler,
// invoked by NVCC’ (e.g. on MacOS). The former needs to see both host and device implementation
// of the functions, while the latter can only deal with one of them.
#elif !defined(EIGEN_HAS_NATIVE_GPU_FP16) || (EIGEN_COMP_CLANG && !EIGEN_COMP_NVCC)  // Emulate support for half floats

#if EIGEN_COMP_CLANG && defined(EIGEN_GPUCC)
// We need to provide emulated *host-side* FP16 operators for clang.
#pragma push_macro("EIGEN_DEVICE_FUNC")
#undef EIGEN_DEVICE_FUNC
#if defined(EIGEN_GPUCC) && defined(EIGEN_HAS_NATIVE_GPU_FP16)
#define EIGEN_DEVICE_FUNC __host__
#else  // both host and device need emulated ops.
#define EIGEN_DEVICE_FUNC __host__ __device__
#endif
#endif

// Definitions for CPUs and older HIP+CUDA, mostly working through conversion
// to/from fp32.
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator+(const half& a, const half& b) { return half(float(a) + float(b)); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator*(const half& a, const half& b) { return half(float(a) * float(b)); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator-(const half& a, const half& b) { return half(float(a) - float(b)); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator/(const half& a, const half& b) { return half(float(a) / float(b)); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator-(const half& a) {
  return raw_uint16_to_half(static_cast<numext::uint16_t>(raw_half_as_uint16(a) ^ 0x8000));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator+=(half& a, const half& b) {
  a = half(float(a) + float(b));
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator*=(half& a, const half& b) {
  a = half(float(a) * float(b));
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator-=(half& a, const half& b) {
  a = half(float(a) - float(b));
  return a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half& operator/=(half& a, const half& b) {
  a = half(float(a) / float(b));
  return a;
}

// Non-negative floating point numbers have a monotonic mapping to non-negative integers.
// This property allows floating point numbers to be reinterpreted as integers for comparisons, which is useful if there
// is no native floating point comparison operator. Floating point signedness is handled by the sign-magnitude
// representation, whereas integers typically use two's complement. Converting the bit pattern from sign-magnitude to
// two's complement allows the transformed bit patterns be compared as signed integers. All edge cases (+/-0 and +/-
// infinity) are handled automatically, except NaN.
//
// fp16 uses 1 sign bit, 5 exponent bits, and 10 mantissa bits. The bit pattern conveys NaN when all the exponent
// bits (5) are set, and at least one mantissa bit is set. The sign bit is irrelevant for determining NaN. To check for
// NaN, clear the sign bit and check if the integral representation is greater than 01111100000000. To test
// for non-NaN, clear the sign bit and check if the integral representation is less than or equal to 01111100000000.

// convert sign-magnitude representation to two's complement
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC int16_t mapToSigned(uint16_t a) {
#if EIGEN_COMP_NVHPC
  // NVHPC through at least 26.5 can miscompile repeated inlined comparisons of the transformed integer representations.
  EIGEN_OPTIMIZATION_BARRIER(a)
#endif
  constexpr uint16_t kAbsMask = (1 << 15) - 1;
  // If the sign bit is set, clear the sign bit and return the (integer) negation. Otherwise, return the input.
  return (a >> 15) ? -(a & kAbsMask) : a;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool isOrdered(const half& a, const half& b) {
  constexpr uint16_t kInf = ((1 << 5) - 1) << 10;
  constexpr uint16_t kAbsMask = (1 << 15) - 1;
  return numext::maxi(a.x & kAbsMask, b.x & kAbsMask) <= kInf;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator==(const half& a, const half& b) {
  bool result = mapToSigned(a.x) == mapToSigned(b.x);
  result &= isOrdered(a, b);
  return result;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator!=(const half& a, const half& b) { return !(a == b); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator<(const half& a, const half& b) {
  bool result = mapToSigned(a.x) < mapToSigned(b.x);
  result &= isOrdered(a, b);
  return result;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator<=(const half& a, const half& b) {
  bool result = mapToSigned(a.x) <= mapToSigned(b.x);
  result &= isOrdered(a, b);
  return result;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator>(const half& a, const half& b) {
  bool result = mapToSigned(a.x) > mapToSigned(b.x);
  result &= isOrdered(a, b);
  return result;
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool operator>=(const half& a, const half& b) {
  bool result = mapToSigned(a.x) >= mapToSigned(b.x);
  result &= isOrdered(a, b);
  return result;
}

#if EIGEN_COMP_CLANG && defined(EIGEN_GPUCC)
#pragma pop_macro("EIGEN_DEVICE_FUNC")
#endif

#endif  // Emulate support for half floats

// Division by an index. Do it in full float precision to avoid accuracy
// issues in converting the denominator to half.
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator/(const half& a, Index b) {
  return half(static_cast<float>(a) / static_cast<float>(b));
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator++(half& a) {
  a += half(1);
  return a;
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator--(half& a) {
  a -= half(1);
  return a;
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator++(half& a, int) {
  half original_value = a;
  ++a;
  return original_value;
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half operator--(half& a, int) {
  half original_value = a;
  --a;
  return original_value;
}

// Conversion routines, including fallbacks for the host or older CUDA.
// Note that newer Intel CPUs (Haswell or newer) have vectorized versions of
// these in hardware. If we need more performance on older/other CPUs, they are
// also possible to vectorize directly.

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR __half_raw raw_uint16_to_half(numext::uint16_t x) {
  // In the device compile phase of a GPU build we cannot simply do a "return __half_raw(x)",
  // because there __half_raw is the vendor type (a union in the hip_fp16 header file) that has
  // no uint16 constructor, and that will trigger a compile error. There "h.x = x" assigns the
  // raw bits, since the vendor member x is an integer.
  // In the host compile phase (and in non-GPU builds) Eigen's own __half_raw is in effect and its
  // member x may be a native fp16 type (__fp16 on arm64, _Float16 with AVX512FP16 or riscv-zfh),
  // so "h.x = x" would perform a numeric integer-to-float conversion that corrupts the raw bits;
  // the explicit uint16 constructor reinterprets the bits for every storage type instead.
  // Fortunately, since we need to disable EIGEN_CONSTEXPR for GPU anyway, we can get out
  // of this catch22 by having separate bodies for the GPU device phase / everything else.
#if defined(EIGEN_GPUCC) && defined(EIGEN_GPU_COMPILE_PHASE)
  __half_raw h;
  h.x = x;
  return h;
#else
  return __half_raw(x);
#endif
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC numext::uint16_t raw_half_as_uint16(const __half_raw& h) {
  // HIP/CUDA/Default have a member 'x' of type uint16_t.
  // For ARM64 native half, the member 'x' is of type __fp16, so we need to bit-cast.
  // For SYCL, cl::sycl::half is _Float16, so cast directly.
#if EIGEN_HAS_ARM64_FP16
  return numext::bit_cast<numext::uint16_t>(h.x);
#elif defined(EIGEN_HAS_BUILTIN_FLOAT16)
  return numext::bit_cast<numext::uint16_t>(h.x);
#elif defined(EIGEN_SYCL_DEVICE_ONLY)
  return numext::bit_cast<numext::uint16_t>(h);
#else
  return h.x;
#endif
}

template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC __half_raw float_to_half_rtne<false>(float ff) {
#if defined(EIGEN_GPU_COMPILE_PHASE)
  __half tmp_ff = __float2half(ff);
  return *(__half_raw*)&tmp_ff;

#elif EIGEN_HAS_ARM64_FP16
  __half_raw h;
  h.x = static_cast<__fp16>(ff);
  return h;

#elif defined(EIGEN_HAS_BUILTIN_FLOAT16)
  __half_raw h;
  h.x = static_cast<_Float16>(ff);
  return h;

#elif defined(EIGEN_HAS_FP16_C)
  __half_raw h;
  // Spell out the vector conversion to avoid Clang warning about the C99 compound literal used by _cvtss_sh.
  h.x = static_cast<numext::uint16_t>(_mm_extract_epi16(_mm_cvtps_ph(_mm_set_ss(ff), 0), 0));
  return h;

#else
  uint32_t f_bits = Eigen::numext::bit_cast<uint32_t>(ff);
  const uint32_t f32infty_bits = {255 << 23};
  const uint32_t f16max_bits = {(127 + 16) << 23};
  const uint32_t denorm_magic_bits = {((127 - 15) + (23 - 10) + 1) << 23};
  const uint32_t sign_mask = 0x80000000u;
  __half_raw o;
  o.x = static_cast<uint16_t>(0x0u);

  const uint32_t sign = f_bits & sign_mask;
  f_bits ^= sign;

  // NOTE all the integer compares in this function can be safely
  // compiled into signed compares since all operands are below
  // 0x80000000. Important if you want fast straight SSE2 code
  // (since there's no unsigned PCMPGTD).

  if (f_bits >= f16max_bits) {                         // result is Inf or NaN (all exponent bits set)
    o.x = (f_bits > f32infty_bits) ? 0x7e00 : 0x7c00;  // NaN->qNaN and Inf->Inf
  } else {                                             // (De)normalized number or zero
    if (f_bits < (113 << 23)) {                        // resulting FP16 is subnormal or zero
      // use a magic value to align our 10 mantissa bits at the bottom of
      // the float. as long as FP addition is round-to-nearest-even this
      // just works.
      f_bits = Eigen::numext::bit_cast<uint32_t>(Eigen::numext::bit_cast<float>(f_bits) +
                                                 Eigen::numext::bit_cast<float>(denorm_magic_bits));

      // and one integer subtract of the bias later, we have our final float!
      o.x = static_cast<numext::uint16_t>(f_bits - denorm_magic_bits);
    } else {
      const uint32_t mant_odd = (f_bits >> 13) & 1;  // resulting mantissa is odd

      // update exponent, rounding bias part 1
      // Equivalent to `f.u += ((unsigned int)(15 - 127) << 23) + 0xfff`, but
      // without arithmetic overflow.
      f_bits += 0xc8000fffU;
      // rounding bias part 2
      f_bits += mant_odd;
      // take the bits!
      o.x = static_cast<numext::uint16_t>(f_bits >> 13);
    }
  }

  o.x |= static_cast<numext::uint16_t>(sign >> 16);
  return o;
#endif
}

// float_to_half_rtne template specialization that assumes that its function
// argument (ff) is either a normal floating-point number that does not
// underflow to a half subnormal, +/-infinity, or zero. Used to improve
// compile-time constant folding and conversion from an integer type to half.
template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC __half_raw float_to_half_rtne<true>(float ff) {
#if defined(EIGEN_GPU_COMPILE_PHASE)
  __half tmp_ff = __float2half(ff);
  return *(__half_raw*)&tmp_ff;

#elif EIGEN_HAS_ARM64_FP16
  __half_raw h;
  h.x = static_cast<__fp16>(ff);
  return h;

#elif defined(EIGEN_HAS_BUILTIN_FLOAT16)
  __half_raw h;
  h.x = static_cast<_Float16>(ff);
  return h;

#else
  uint32_t f_bits = Eigen::numext::bit_cast<uint32_t>(ff);
  const uint32_t f16max_bits = {(127 + 16) << 23};
  const uint32_t sign = f_bits & 0x80000000u;
  f_bits ^= sign;

  numext::uint16_t raw = 0x0000u;
  if (f_bits >= f16max_bits) {
    raw = 0x7c00u;
  } else if (f_bits != 0) {
    const uint32_t mant_odd = (f_bits >> 13) & 1u;
    f_bits += 0xc8000fffU + mant_odd;
    raw = static_cast<numext::uint16_t>(f_bits >> 13);
  }
  return raw_uint16_to_half(raw | static_cast<numext::uint16_t>(sign >> 16));
#endif
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC float half_to_float(__half_raw h) {
#if defined(EIGEN_GPU_COMPILE_PHASE)
  return __half2float(h);
#elif EIGEN_HAS_ARM64_FP16 || defined(EIGEN_HAS_BUILTIN_FLOAT16)
  return static_cast<float>(h.x);
#elif defined(EIGEN_HAS_FP16_C)
#if EIGEN_COMP_MSVC
  // MSVC does not have scalar instructions.
  return _mm_cvtss_f32(_mm_cvtph_ps(_mm_set1_epi16(h.x)));
#else
  return _cvtsh_ss(h.x);
#endif
#else
  const float magic = Eigen::numext::bit_cast<float>(static_cast<uint32_t>(113 << 23));
  const uint32_t shifted_exp = 0x7c00 << 13;  // exponent mask after shift
  uint32_t o_bits = (h.x & 0x7fff) << 13;     // exponent/mantissa bits
  const uint32_t exp = shifted_exp & o_bits;  // just the exponent
  o_bits += (127 - 15) << 23;                 // exponent adjust

  // handle exponent special cases
  if (exp == shifted_exp) {      // Inf/NaN?
    o_bits += (128 - 16) << 23;  // extra exp adjust
  } else if (exp == 0) {         // Zero/Denormal?
    o_bits += 1 << 23;           // extra exp adjust
    // renormalize
    o_bits = Eigen::numext::bit_cast<uint32_t>(Eigen::numext::bit_cast<float>(o_bits) - magic);
  }

  o_bits |= (h.x & 0x8000u) << 16;  // sign bit
  return Eigen::numext::bit_cast<float>(o_bits);
#endif
}

// --- standard functions ---

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool(isinf)(const half& a) { return (raw_half_as_uint16(a) & 0x7fff) == 0x7c00; }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool(isnan)(const half& a) {
#if defined(EIGEN_GPU_COMPILE_PHASE)
  return __hisnan(a);
#else
  return (raw_half_as_uint16(a) & 0x7fff) > 0x7c00;
#endif
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool(isfinite)(const half& a) {
  return (raw_half_as_uint16(a) & 0x7fff) < 0x7c00;
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half abs(const half& a) {
#if defined(EIGEN_HAS_ARM64_FP16_SCALAR_ARITHMETIC)
  return half_from_rep(vabsh_f16(a.x));
#else
  return raw_uint16_to_half(static_cast<numext::uint16_t>(raw_half_as_uint16(a) & 0x7FFF));
#endif
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half exp(const half& a) {
#if defined(EIGEN_CUDA_ARCH) || defined(EIGEN_HIP_DEVICE_COMPILE)
  return half(hexp(::__half(a)));
#else
  return half(numext::exp(float(a)));
#endif
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half exp2(const half& a) {
#if defined(EIGEN_CUDA_ARCH) || defined(EIGEN_HIP_DEVICE_COMPILE)
  return half(hexp2(::__half(a)));
#else
  return half(numext::exp2(float(a)));
#endif
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half expm1(const half& a) { return half(numext::expm1(float(a))); }
// float covers the half range, so the single conversion back rounds correctly.
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half ldexp(const half& a, int exponent) {
  return half(numext::ldexp(float(a), exponent));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half log(const half& a) {
#if defined(EIGEN_GPU_COMPILE_PHASE)
  return half(hlog(::__half(a)));
#else
  return half(numext::log(float(a)));
#endif
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half log1p(const half& a) { return half(numext::log1p(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half log10(const half& a) { return half(numext::log10(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half log2(const half& a) {
  return half(static_cast<float>(EIGEN_LOG2E) * numext::log(float(a)));
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half sqrt(const half& a) {
#if defined(EIGEN_CUDA_ARCH) || defined(EIGEN_HIP_DEVICE_COMPILE)
  return half(hsqrt(::__half(a)));
#else
  // numext::sqrt<float> would be instantiated here, before the SSE specialization is declared.
  EIGEN_USING_STD(sqrt);
  return half(sqrt(float(a)));
#endif
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half pow(const half& a, const half& b) {
  return half(numext::pow(float(a), float(b)));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half atan2(const half& a, const half& b) {
  return half(numext::atan2(float(a), float(b)));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half sin(const half& a) { return half(numext::sin(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half cos(const half& a) { return half(numext::cos(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half tan(const half& a) { return half(numext::tan(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half tanh(const half& a) {
  // numext::tanh(float) can select an approximation that changes half rounding.
  EIGEN_USING_STD(tanh);
  return half(tanh(float(a)));
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half asin(const half& a) { return half(numext::asin(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half acos(const half& a) { return half(numext::acos(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half atan(const half& a) { return half(numext::atan(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half atanh(const half& a) { return half(numext::atanh(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half floor(const half& a) {
#if (defined(EIGEN_CUDA_ARCH)) || defined(EIGEN_HIP_DEVICE_COMPILE)
  return half(hfloor(::__half(a)));
#else
  return half(numext::floor(float(a)));
#endif
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half ceil(const half& a) {
#if (defined(EIGEN_CUDA_ARCH)) || defined(EIGEN_HIP_DEVICE_COMPILE)
  return half(hceil(::__half(a)));
#else
  return half(numext::ceil(float(a)));
#endif
}
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half rint(const half& a) { return half(numext::rint(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half round(const half& a) { return half(numext::round(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half trunc(const half& a) { return half(numext::trunc(float(a))); }
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half fmod(const half& a, const half& b) {
  return half(numext::fmod(float(a), float(b)));
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half(min)(const half& a, const half& b) { return b < a ? b : a; }

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half(max)(const half& a, const half& b) { return a < b ? b : a; }

// a * b + c for half or bfloat16 a, b and c, rounded to odd in double: narrowing the result to half or bfloat16, which
// rounds to odd in float and then to nearest, then rounds a * b + c once. A float fma rounds twice instead: for half,
// fma(33/32, 65/64, 2^-24) rounds to the midpoint 2145/2048 in float, then to even, 1 + 3/64, not up to 1 + 49/1024.
// a * b is exact in double, and TwoSum (Knuth) gives the error e of s = a * b + c exactly. Device code does not use
// it: GPUs run double at a fraction of the float rate, and SYCL devices need not support it.
EIGEN_DEVICE_FUNC inline double fma_round_to_odd(float a, float b, float c) {
  const double p = static_cast<double>(a) * static_cast<double>(b);
  const double c_wide = static_cast<double>(c);
  const double s = p + c_wide;
  // p is exact, so s is nonfinite only for a nonfinite operand; TwoSum would then evaluate inf - inf (FE_INVALID).
  if (!(numext::isfinite)(s)) return s;
  const double p_part = s - c_wide;
  const double e = (p - p_part) + (c_wide - (s - p_part));
  if (e == 0) return s;
  // Inexact: truncate toward zero, then set the last bit.
  const bool away = (s > 0) != (e > 0);
  return numext::bit_cast<double>((numext::bit_cast<numext::uint64_t>(s) - (away ? 1u : 0u)) | 1u);
}

EIGEN_DEVICE_FUNC inline half fma(const half& a, const half& b, const half& c) {
#if defined(EIGEN_HAS_NATIVE_GPU_FP16)
  // Rounded once in hardware.
  return __hfma(::__half(a), ::__half(b), ::__half(c));
#elif defined(__SYCL_DEVICE_ONLY__)
  // Rounds twice; see fma_round_to_odd.
  return half(numext::fma(static_cast<float>(a), static_cast<float>(b), static_cast<float>(c)));
#elif defined(EIGEN_HAS_ARM64_FP16_SCALAR_ARITHMETIC)
  return half_from_rep(vfmah_f16(c.x, a.x, b.x));
#elif defined(EIGEN_VECTORIZE_AVX512FP16)
  // Reduces to vfmadd213sh.
  return half(_mm_cvtsh_h(_mm_fmadd_ph(_mm_set_sh(a.x), _mm_set_sh(b.x), _mm_set_sh(c.x))));
#else
  return half(fma_round_to_odd(static_cast<float>(a), static_cast<float>(b), static_cast<float>(c)));
#endif
}

#ifndef EIGEN_NO_IO
EIGEN_ALWAYS_INLINE std::ostream& operator<<(std::ostream& os, const half& v) {
  os << static_cast<float>(v);
  return os;
}
#endif

}  // end namespace half_impl

// import Eigen::half_impl::half into Eigen namespace
// using half_impl::half;

namespace internal {

template <>
struct is_arithmetic<half> : std::true_type {};

template <>
struct random_impl<half> {
  enum : int { MantissaBits = 10 };
  using Impl = random_impl<float>;
  static EIGEN_DEVICE_FUNC inline half run(const half& x, const half& y) {
    float result = Impl::run(x, y, MantissaBits);
    return half(result);
  }
  static EIGEN_DEVICE_FUNC inline half run() {
    float result = Impl::run(MantissaBits);
    return half(result);
  }
};

}  // end namespace internal

template <>
struct NumTraits<Eigen::half> : GenericNumTraits<Eigen::half> {
  enum { IsSigned = true, IsInteger = false, IsComplex = false, RequireInitialization = false };

  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR static EIGEN_STRONG_INLINE Eigen::half epsilon() {
    // 0x1400 is 2^-10, the fp16 machine epsilon, matching std::numeric_limits<Eigen::half>::epsilon().
    return half_impl::raw_uint16_to_half(0x1400);
  }
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR static EIGEN_STRONG_INLINE Eigen::half dummy_precision() {
    return half_impl::raw_uint16_to_half(0x211f);  //  Eigen::half(1e-2f);
  }
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR static EIGEN_STRONG_INLINE Eigen::half highest() {
    return half_impl::raw_uint16_to_half(0x7bff);
  }
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR static EIGEN_STRONG_INLINE Eigen::half lowest() {
    return half_impl::raw_uint16_to_half(0xfbff);
  }
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR static EIGEN_STRONG_INLINE Eigen::half infinity() {
    return half_impl::raw_uint16_to_half(0x7c00);
  }
  EIGEN_DEVICE_FUNC _EIGEN_MAYBE_CONSTEXPR static EIGEN_STRONG_INLINE Eigen::half quiet_NaN() {
    return half_impl::raw_uint16_to_half(0x7e00);
  }
};

}  // end namespace Eigen

#undef _EIGEN_MAYBE_CONSTEXPR

namespace Eigen {
namespace numext {

#if defined(EIGEN_GPU_COMPILE_PHASE)

template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool(isnan)(const Eigen::half& h) {
  return (half_impl::isnan)(h);
}

template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool(isinf)(const Eigen::half& h) {
  return (half_impl::isinf)(h);
}

template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool(isfinite)(const Eigen::half& h) {
  return (half_impl::isfinite)(h);
}

#endif

template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half bit_cast<Eigen::half, uint16_t>(const uint16_t& src) {
  return Eigen::half(Eigen::half_impl::raw_uint16_to_half(src));
}

template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC uint16_t bit_cast<uint16_t, Eigen::half>(const Eigen::half& src) {
  return Eigen::half_impl::raw_half_as_uint16(src);
}

EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half nextafter(const Eigen::half& from, const Eigen::half& to) {
  if (numext::isnan EIGEN_NOT_A_MACRO(from)) {
    return from;
  }
  if (numext::isnan EIGEN_NOT_A_MACRO(to)) {
    return to;
  }
  if (from == to) {
    return to;
  }
  uint16_t from_bits = numext::bit_cast<uint16_t>(from);
  bool from_sign = from_bits >> 15;
  if ((from_bits & 0x7fff) == 0) {
    // From ±0 toward a nonzero value: the neighbor is the smallest subnormal
    // carrying the sign of the direction (IEEE-754 nextUp/nextDown of zero).
    from_bits = (to > from) ? uint16_t(0x0001) : uint16_t(0x8001);
  } else if ((to > from) != from_sign) {
    // Toward the infinity with the same sign as from: increase the magnitude.
    ++from_bits;
  } else {
    --from_bits;
  }
  return numext::bit_cast<Eigen::half>(from_bits);
}

// Specialize multiply-add to match packet operations and reduce conversions to/from float.
template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half madd<Eigen::half>(const Eigen::half& x, const Eigen::half& y,
                                                                    const Eigen::half& z) {
  return Eigen::half(static_cast<float>(x) * static_cast<float>(y) + static_cast<float>(z));
}

}  // namespace numext
}  // namespace Eigen

// Warp shuffle overloads for Eigen::half.
// CUDA uses __shfl_*_sync (with mask); HIP uses __shfl_* (no mask).
// Note that the following are __device__ - only functions.
#if defined(EIGEN_CUDACC) || defined(EIGEN_HIPCC)

#if defined(EIGEN_CUDACC)

__device__ EIGEN_STRONG_INLINE Eigen::half __shfl_sync(unsigned mask, Eigen::half var, int srcLane,
                                                       int width = warpSize) {
  const __half h = var;
  return static_cast<Eigen::half>(__shfl_sync(mask, h, srcLane, width));
}

__device__ EIGEN_STRONG_INLINE Eigen::half __shfl_up_sync(unsigned mask, Eigen::half var, unsigned int delta,
                                                          int width = warpSize) {
  const __half h = var;
  return static_cast<Eigen::half>(__shfl_up_sync(mask, h, delta, width));
}

__device__ EIGEN_STRONG_INLINE Eigen::half __shfl_down_sync(unsigned mask, Eigen::half var, unsigned int delta,
                                                            int width = warpSize) {
  const __half h = var;
  return static_cast<Eigen::half>(__shfl_down_sync(mask, h, delta, width));
}

__device__ EIGEN_STRONG_INLINE Eigen::half __shfl_xor_sync(unsigned mask, Eigen::half var, int laneMask,
                                                           int width = warpSize) {
  const __half h = var;
  return static_cast<Eigen::half>(__shfl_xor_sync(mask, h, laneMask, width));
}

#else  // HIP

__device__ EIGEN_STRONG_INLINE Eigen::half __shfl(Eigen::half var, int srcLane, int width = warpSize) {
  const int ivar = static_cast<int>(Eigen::numext::bit_cast<Eigen::numext::uint16_t>(var));
  return Eigen::numext::bit_cast<Eigen::half>(static_cast<Eigen::numext::uint16_t>(__shfl(ivar, srcLane, width)));
}

__device__ EIGEN_STRONG_INLINE Eigen::half __shfl_up(Eigen::half var, unsigned int delta, int width = warpSize) {
  const int ivar = static_cast<int>(Eigen::numext::bit_cast<Eigen::numext::uint16_t>(var));
  return Eigen::numext::bit_cast<Eigen::half>(static_cast<Eigen::numext::uint16_t>(__shfl_up(ivar, delta, width)));
}

__device__ EIGEN_STRONG_INLINE Eigen::half __shfl_down(Eigen::half var, unsigned int delta, int width = warpSize) {
  const int ivar = static_cast<int>(Eigen::numext::bit_cast<Eigen::numext::uint16_t>(var));
  return Eigen::numext::bit_cast<Eigen::half>(static_cast<Eigen::numext::uint16_t>(__shfl_down(ivar, delta, width)));
}

__device__ EIGEN_STRONG_INLINE Eigen::half __shfl_xor(Eigen::half var, int laneMask, int width = warpSize) {
  const int ivar = static_cast<int>(Eigen::numext::bit_cast<Eigen::numext::uint16_t>(var));
  return Eigen::numext::bit_cast<Eigen::half>(static_cast<Eigen::numext::uint16_t>(__shfl_xor(ivar, laneMask, width)));
}

#endif  // HIP vs CUDA
#endif  // __shfl*

// ldg() has an overload for __half_raw, but we also need one for Eigen::half.
#if defined(EIGEN_CUDACC) || defined(EIGEN_HIPCC)
EIGEN_STRONG_INLINE __device__ Eigen::half __ldg(const Eigen::half* ptr) {
  return Eigen::half_impl::raw_uint16_to_half(__ldg(reinterpret_cast<const Eigen::numext::uint16_t*>(ptr)));
}
#endif  // __ldg

#if EIGEN_HAS_STD_HASH
namespace std {
template <>
struct hash<Eigen::half> {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::size_t operator()(const Eigen::half& a) const {
    // -0 == +0, so both must hash alike.
    const Eigen::numext::uint16_t bits = Eigen::numext::bit_cast<Eigen::numext::uint16_t>(a);
    return static_cast<std::size_t>((bits & 0x7fff) == 0 ? 0 : bits);
  }
};
}  // end namespace std
#endif

namespace Eigen {
namespace internal {

template <>
struct cast_impl<float, half> {
  EIGEN_DEVICE_FUNC static inline half run(const float& a) {
#if defined(EIGEN_GPU_COMPILE_PHASE)
    return __float2half(a);
#else
    return half(a);
#endif
  }
};

template <>
struct cast_impl<int, half> {
  EIGEN_DEVICE_FUNC static inline half run(const int& a) {
#if defined(EIGEN_GPU_COMPILE_PHASE)
    return __float2half(static_cast<float>(a));
#else
    return half(static_cast<float>(a));
#endif
  }
};

template <>
struct cast_impl<half, float> {
  EIGEN_DEVICE_FUNC static inline float run(const half& a) {
#if defined(EIGEN_GPU_COMPILE_PHASE)
    return __half2float(a);
#else
    return static_cast<float>(a);
#endif
  }
};

}  // namespace internal
}  // namespace Eigen

#endif  // EIGEN_HALF_H