Add nextafter and nextafterf from musl

src/math/mod.rs

@@ -146,6 +146,8 @@ mod log2f;
 mod logf;
 mod modf;
 mod modff;
+mod nextafter;
+mod nextafterf;
 mod pow;
 mod powf;
 mod remainder;
@@ -258,6 +260,8 @@ pub use self::log2f::log2f;
 pub use self::logf::logf;
 pub use self::modf::modf;
 pub use self::modff::modff;
+pub use self::nextafter::nextafter;
+pub use self::nextafterf::nextafterf;
 pub use self::pow::pow;
 pub use self::powf::powf;
 pub use self::remainder::remainder;

src/math/nextafter.rs

@@ -0,0 +1,38 @@
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn nextafter(x: f64, y: f64) -> f64 {
+    if x.is_nan() || y.is_nan() {
+        return x + y;
+    }
+
+    let mut ux_i = x.to_bits();
+    let uy_i = y.to_bits();
+    if ux_i == uy_i {
+        return y;
+    }
+
+    let ax = ux_i & !1_u64 / 2;
+    let ay = uy_i & !1_u64 / 2;
+    if ax == 0 {
+        if ay == 0 {
+            return y;
+        }
+        ux_i = (uy_i & 1_u64 << 63) | 1;
+    } else if ax > ay || ((ux_i ^ uy_i) & 1_u64 << 63) != 0 {
+        ux_i -= 1;
+    } else {
+        ux_i += 1;
+    }
+
+    let e = ux_i.wrapping_shr(52 & 0x7ff);
+    // raise overflow if ux.f is infinite and x is finite
+    if e == 0x7ff {
+        force_eval!(x + x);
+    }
+    let ux_f = f64::from_bits(ux_i);
+    // raise underflow if ux.f is subnormal or zero
+    if e == 0 {
+        force_eval!(x * x + ux_f * ux_f);
+    }
+    ux_f
+}

src/math/nextafterf.rs

@@ -0,0 +1,38 @@
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn nextafterf(x: f32, y: f32) -> f32 {
+    if x.is_nan() || y.is_nan() {
+        return x + y;
+    }
+
+    let mut ux_i = x.to_bits();
+    let uy_i = y.to_bits();
+    if ux_i == uy_i {
+        return y;
+    }
+
+    let ax = ux_i & 0x7fff_ffff_u32;
+    let ay = uy_i & 0x7fff_ffff_u32;
+    if ax == 0 {
+        if ay == 0 {
+            return y;
+        }
+        ux_i = (uy_i & 0x8000_0000_u32) | 1;
+    } else if ax > ay || ((ux_i ^ uy_i) & 0x8000_0000_u32) != 0 {
+        ux_i -= 1;
+    } else {
+        ux_i += 1;
+    }
+
+    let e = ux_i.wrapping_shr(0x7f80_0000_u32);
+    // raise overflow if ux_f is infinite and x is finite
+    if e == 0x7f80_0000_u32 {
+        force_eval!(x + x);
+    }
+    let ux_f = f32::from_bits(ux_i);
+    // raise underflow if ux_f is subnormal or zero
+    if e == 0 {
+        force_eval!(x * x + ux_f * ux_f);
+    }
+    ux_f
+}