phynexis-deps/eigen
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Core: Implement numext::nextafter for all floating point types

Browse Source 
libeigen/eigen!2834

Closes #2937

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-08-17 21:01:28 -07:00
co-authored by Rasmus Munk Larsen
parent e857597ca0
commit 5593b98453
6 changed files with 133 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Eigen/src/Core/MathFunctions.h
+22
View File
@@ -1928,6 +1928,28 @@ EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double fmod(const double& a, const double&
}
#endif
template <typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T nextafter(const T& from, const T& to) {
EIGEN_USING_STD(nextafter);
return nextafter(from, to);
}
#if defined(SYCL_DEVICE_ONLY)
SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(nextafter, nextafter)
#endif
#if defined(EIGEN_GPUCC)
template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float nextafter(const float& from, const float& to) {
return ::nextafterf(from, to);
}
template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double nextafter(const double& from, const double& to) {
return ::nextafter(from, to);
}
#endif
#if defined(SYCL_DEVICE_ONLY)
#undef SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_BINARY
#undef SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_UNARY
Eigen/src/Core/arch/Default/BFloat16.h
+6 -6
View File
@@ -826,13 +826,13 @@ EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 nextafter(const bfloat16& from, c
}
uint16_t from_bits = numext::bit_cast<uint16_t>(from);
bool from_sign = from_bits >> 15;
// Whether we are adjusting toward the infinity with the same sign as from.
bool toward_inf = (to > from) == !from_sign;
if (toward_inf) {
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 if ((from_bits & 0x7fff) == 0) {
// Adjusting away from inf, but from is zero, so just toggle the sign.
from_bits ^= 0x8000;
} else {
--from_bits;
}
Eigen/src/Core/arch/Default/Half.h
+25
View File
@@ -968,6 +968,31 @@ EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC uint16_t bit_cast<uint16_t, Eigen::half>(c
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,
test/bfloat16_float.cpp
+4 -2
View File
@@ -384,9 +384,11 @@ void test_nextafter() {
std::numeric_limits<bfloat16>::infinity());
VERIFY(numext::nextafter(-(std::numeric_limits<bfloat16>::max)(), -std::numeric_limits<bfloat16>::infinity()) ==
-std::numeric_limits<bfloat16>::infinity());
// The neighbors of ±0 are the smallest subnormals with the sign of the
// direction, matching std::nextafter (IEEE-754 nextUp/nextDown of zero).
VERIFY_BFLOAT16_BITS_EQUAL(numext::nextafter(bfloat16(0.0f), bfloat16(1.0f)), 0x0001);
VERIFY_BFLOAT16_BITS_EQUAL(numext::nextafter(bfloat16(-0.0f), bfloat16(1.0f)), 0x0000);
VERIFY_BFLOAT16_BITS_EQUAL(numext::nextafter(bfloat16(0.0f), bfloat16(-1.0f)), 0x8000);
VERIFY_BFLOAT16_BITS_EQUAL(numext::nextafter(bfloat16(-0.0f), bfloat16(1.0f)), 0x0001);
VERIFY_BFLOAT16_BITS_EQUAL(numext::nextafter(bfloat16(0.0f), bfloat16(-1.0f)), 0x8001);
VERIFY_BFLOAT16_BITS_EQUAL(numext::nextafter(bfloat16(-0.0f), bfloat16(-1.0f)), 0x8001);
VERIFY_BFLOAT16_BITS_EQUAL(numext::nextafter(bfloat16(0.0f), bfloat16(-0.0f)), 0x8000);
VERIFY_BFLOAT16_BITS_EQUAL(numext::nextafter(bfloat16(-0.0f), bfloat16(0.0f)), 0x0000);
test/half_float.cpp
+36
View File
@@ -342,6 +342,41 @@ void test_product() {
VERIFY_IS_APPROX(Ch.noalias() += Ah * Bh, (Cf.noalias() += Af * Bf).cast<half>());
}
void test_nextafter() {
VERIFY((numext::isnan)(numext::nextafter(std::numeric_limits<half>::quiet_NaN(), half(1.0f))));
VERIFY((numext::isnan)(numext::nextafter(half(1.0f), std::numeric_limits<half>::quiet_NaN())));
VERIFY(numext::nextafter(half(0.0f), half(0.0f)) == half(0.0f));
VERIFY(numext::nextafter(half(1.0f), half(1.0f)) == half(1.0f));
VERIFY(numext::nextafter(half(-1.0f), half(-1.0f)) == half(-1.0f));
VERIFY(numext::nextafter(std::numeric_limits<half>::infinity(), std::numeric_limits<half>::infinity()) ==
std::numeric_limits<half>::infinity());
VERIFY(numext::nextafter(std::numeric_limits<half>::infinity(), half(0.0f)) == (std::numeric_limits<half>::max)());
VERIFY(numext::nextafter(-std::numeric_limits<half>::infinity(), half(0.0f)) == -(std::numeric_limits<half>::max)());
VERIFY(numext::nextafter(half(1.0f), std::numeric_limits<half>::infinity()) ==
half(1.0f) + std::numeric_limits<half>::epsilon());
VERIFY(numext::nextafter(half(1.0f), -std::numeric_limits<half>::infinity()) ==
half(1.0f) - std::numeric_limits<half>::epsilon() / half(2.0f));
VERIFY(numext::nextafter(half(-1.0f), -std::numeric_limits<half>::infinity()) ==
half(-1.0f) - std::numeric_limits<half>::epsilon());
VERIFY(numext::nextafter(half(-1.0f), std::numeric_limits<half>::infinity()) ==
half(-1.0f) + std::numeric_limits<half>::epsilon() / half(2.0f));
VERIFY(numext::nextafter((std::numeric_limits<half>::max)(), std::numeric_limits<half>::infinity()) ==
std::numeric_limits<half>::infinity());
VERIFY(numext::nextafter(-(std::numeric_limits<half>::max)(), -std::numeric_limits<half>::infinity()) ==
-std::numeric_limits<half>::infinity());
// The neighbors of ±0 are the smallest subnormals with the sign of the
// direction, matching std::nextafter (IEEE-754 nextUp/nextDown of zero).
VERIFY_HALF_BITS_EQUAL(numext::nextafter(half(0.0f), half(1.0f)), 0x0001);
VERIFY_HALF_BITS_EQUAL(numext::nextafter(half(-0.0f), half(1.0f)), 0x0001);
VERIFY_HALF_BITS_EQUAL(numext::nextafter(half(0.0f), half(-1.0f)), 0x8001);
VERIFY_HALF_BITS_EQUAL(numext::nextafter(half(-0.0f), half(-1.0f)), 0x8001);
// from == to returns to, preserving the sign of zero.
VERIFY_HALF_BITS_EQUAL(numext::nextafter(half(0.0f), half(-0.0f)), 0x8000);
VERIFY_HALF_BITS_EQUAL(numext::nextafter(half(-0.0f), half(0.0f)), 0x0000);
VERIFY_HALF_BITS_EQUAL(numext::nextafter(half(0.0f), half(0.0f)), 0x0000);
VERIFY_HALF_BITS_EQUAL(numext::nextafter(half(-0.0f), half(-0.0f)), 0x8000);
}
EIGEN_DECLARE_TEST(half_float) {
CALL_SUBTEST(test_numtraits());
for (int i = 0; i < g_repeat; i++) {
@@ -352,5 +387,6 @@ EIGEN_DECLARE_TEST(half_float) {
CALL_SUBTEST(test_trigonometric_functions());
CALL_SUBTEST(test_array());
CALL_SUBTEST(test_product());
CALL_SUBTEST(test_nextafter());
}
}
test/numext.cpp
+40
View File
@@ -422,6 +422,40 @@ void check_signbit() {
check_signbit_impl<T>::run();
}
template <typename T>
void check_nextafter() {
const T zero(0);
const T one(1);
const T two(2);
const T eps = std::numeric_limits<T>::epsilon();
const T denorm_min = std::numeric_limits<T>::denorm_min();
const T inf = std::numeric_limits<T>::infinity();
const T nan = std::numeric_limits<T>::quiet_NaN();
const T max = (std::numeric_limits<T>::max)();
// from == to returns to.
VERIFY(numext::equal_strict(numext::nextafter(one, one), one));
// One-ulp steps around 1.
VERIFY(numext::equal_strict(numext::nextafter(one, two), one + eps));
VERIFY(numext::equal_strict(numext::nextafter(one + eps, zero), one));
// The neighbors of ±0 are the smallest subnormals, with the sign of the direction.
VERIFY(numext::equal_strict(numext::nextafter(zero, one), denorm_min));
VERIFY(numext::equal_strict(numext::nextafter(zero, -one), -denorm_min));
VERIFY(numext::equal_strict(numext::nextafter(-zero, one), denorm_min));
VERIFY(numext::equal_strict(numext::copysign(one, numext::nextafter(zero, -one)), -one));
// Stepping the smallest subnormals toward the other sign lands on the zero
// of the starting sign (IEEE-754 nextUp/nextDown).
VERIFY(numext::equal_strict(numext::nextafter(denorm_min, -one), zero));
VERIFY(numext::equal_strict(numext::copysign(one, numext::nextafter(denorm_min, -one)), one));
VERIFY(numext::equal_strict(numext::copysign(one, numext::nextafter(-denorm_min, one)), -one));
// Infinities saturate and unsaturate by one step.
VERIFY(numext::equal_strict(numext::nextafter(max, inf), inf));
VERIFY(numext::equal_strict(numext::nextafter(inf, zero), max));
// NaNs propagate.
VERIFY((numext::isnan)(numext::nextafter(nan, one)));
VERIFY((numext::isnan)(numext::nextafter(one, nan)));
}
template <typename T>
void check_shift() {
using SignedT = typename numext::get_integer_by_size<sizeof(T)>::signed_type;
@@ -530,6 +564,12 @@ EIGEN_DECLARE_TEST(numext) {
CALL_SUBTEST(check_signbit<int32_t>());
CALL_SUBTEST(check_signbit<int64_t>());
CALL_SUBTEST(check_nextafter<half>());
CALL_SUBTEST(check_nextafter<bfloat16>());
CALL_SUBTEST(check_nextafter<float>());
CALL_SUBTEST(check_nextafter<double>());
CALL_SUBTEST(check_nextafter<long double>());
CALL_SUBTEST(check_shift<int8_t>());
CALL_SUBTEST(check_shift<int16_t>());
CALL_SUBTEST(check_shift<int32_t>());