Add float next_toward() method

Author: danakj

subspace/num/__private/float_methods.inc

@@ -489,6 +489,20 @@ sus_pure inline _self mul_add(_self a, _self b) const& noexcept {
   return static_cast<_primitive>(
       ::fma(primitive_value, a.primitive_value, b.primitive_value));
 }
+
+/// Returns the next representable value of the float type after `self` in the
+/// direction of `toward`. If `self == toward`, `toward` is returned. If either
+/// `self` or `toward` is NAN, NAN is returned.
+///
+/// This is implemented by the
+/// [cmath](https://en.cppreference.com/w/c/numeric/math/nextafter) library, see
+/// the documentation for details on errors.
+//
+// TODO: constexpr in C++23.
+sus_pure inline _self next_toward(_self toward) const& noexcept {
+  return __private::next_toward(primitive_value, toward.primitive_value);
+}
+
 /// Raises a number to a floating point power.
 sus_pure inline _self powf(_self n) const& noexcept {
   // MSVC pow(float) is returning a double for some reason.

subspace/num/__private/intrinsics.h

@@ -14,20 +14,20 @@
 
 #pragma once
 
-#include <math.h>  // TODO: Remove this and define nans and truncf ourselves.
 #include <stddef.h>
 #include <stdint.h>
 
 #include <bit>
+#include <cmath>
 #include <type_traits>
 
 #if _MSC_VER
 #include <intrin.h>
 #endif
 
 #include "subspace/assertions/unreachable.h"
-#include "subspace/macros/inline.h"
 #include "subspace/macros/builtin.h"
+#include "subspace/macros/inline.h"
 #include "subspace/macros/pure.h"
 #include "subspace/marker/unsafe.h"
 #include "subspace/mem/size_of.h"
@@ -1511,4 +1511,14 @@ sus_pure_const constexpr inline T float_clamp(marker::UnsafeFnMarker, T x,
   return x;
 }
 
+// TODO: constexpr in C++23.
+template <class T>
+  requires(std::is_floating_point_v<T> && ::sus::mem::size_of<T>() <= 8)
+sus_pure_const inline T next_toward(T from, T to) {
+  if constexpr (::sus::mem::size_of<T>() == 4u)
+    return std::nexttowardf(from, to);
+  else
+    return std::nexttoward(from, to);
+}
+
 }  // namespace sus::num::__private

subspace/num/f32_unittest.cc

@@ -1106,4 +1106,15 @@ TEST(f32, fmt) {
   EXPECT_EQ(fmt::format("{:.4f}", 1234.567_f32), "1234.5670");
 }
 
+TEST(f32, NextToward) {
+  EXPECT_EQ((0_f32).next_toward(f32::NAN).is_nan(), true);
+  EXPECT_EQ(f32::NAN.next_toward(0_f32).is_nan(), true);
+  EXPECT_EQ((1_f32).next_toward(f32::INFINITY) - 1_f32,
+            f32(std::numeric_limits<float>::epsilon()));
+  EXPECT_GT((0_f32).next_toward(f32::INFINITY), 0_f32);
+  EXPECT_GT((0_f32).next_toward(1_f32), 0_f32);
+  EXPECT_LT((0_f32).next_toward(f32::NEG_INFINITY), 0_f32);
+  EXPECT_LT((0_f32).next_toward(-1_f32), 0_f32);
+}
+
 }  // namespace

subspace/num/f64_unittest.cc

@@ -1082,4 +1082,15 @@ TEST(f64, fmt) {
   EXPECT_EQ(fmt::format("{:.4f}", 1234890.567_f64), "1234890.5670");
 }
 
+TEST(f64, NextToward) {
+  EXPECT_EQ((0_f64).next_toward(f64::NAN).is_nan(), true);
+  EXPECT_EQ(f64::NAN.next_toward(0_f64).is_nan(), true);
+  EXPECT_EQ((1_f64).next_toward(f64::INFINITY) - 1_f64,
+            f64(std::numeric_limits<double>::epsilon()));
+  EXPECT_GT((0_f64).next_toward(f64::INFINITY), 0_f64);
+  EXPECT_GT((0_f64).next_toward(1_f64), 0_f64);
+  EXPECT_LT((0_f64).next_toward(f64::NEG_INFINITY), 0_f64);
+  EXPECT_LT((0_f64).next_toward(-1_f64), 0_f64);
+}
+
 }  // namespace