Implement fmod and tweak fmodf

Author: P1n3appl3

library/compiler-builtins/libm/src/math/fmod.rs

@@ -0,0 +1,80 @@
+use core::u64;
+
+#[inline]
+pub fn fmod(x: f64, y: f64) -> f64 {
+    let mut uxi = x.to_bits();
+    let mut uyi = y.to_bits();
+    let mut ex = (uxi >> 52 & 0x7ff) as i64;
+    let mut ey = (uyi >> 52 & 0x7ff) as i64;
+    let sx = uxi >> 63;
+    let mut i;
+
+    if uyi << 1 == 0 || y.is_nan() || ex == 0x7ff {
+        return (x * y) / (x * y);
+    }
+    if uxi << 1 <= uyi << 1 {
+        if uxi << 1 == uyi << 1 {
+            return 0.0 * x;
+        }
+        return x;
+    }
+
+    /* normalize x and y */
+    if ex == 0 {
+        i = uxi << 12;
+        while i >> 63 == 0 {
+            ex -= 1;
+            i <<= 1;
+        }
+        uxi <<= -ex + 1;
+    } else {
+        uxi &= u64::MAX >> 12;
+        uxi |= 1 << 52;
+    }
+    if ey == 0 {
+        i = uyi << 12;
+        while i >> 63 == 0 {
+            ey -= 1;
+            i <<= 1;
+        }
+        uyi <<= -ey + 1;
+    } else {
+        uyi &= u64::MAX >> 12;
+        uyi |= 1 << 52;
+    }
+
+    /* x mod y */
+    while ex > ey {
+        i = uxi - uyi;
+        if i >> 63 == 0 {
+            if i == 0 {
+                return 0.0 * x;
+            }
+            uxi = i;
+        }
+        uxi <<= 1;
+        ex -= 1;
+    }
+    i = uxi - uyi;
+    if i >> 63 == 0 {
+        if i == 0 {
+            return 0.0 * x;
+        }
+        uxi = i;
+    }
+    while uxi >> 52 == 0 {
+        uxi <<= 1;
+        ex -= 1;
+    }
+
+    /* scale result */
+    if ex > 0 {
+        uxi -= 1 << 52;
+        uxi |= (ex as u64) << 52;
+    } else {
+        uxi >>= -ex + 1;
+    }
+    uxi |= (sx as u64) << 63;
+
+    f64::from_bits(uxi)
+}

library/compiler-builtins/libm/src/math/fmodf.rs

@@ -1,7 +1,5 @@
 use core::u32;
 
-use super::isnanf;
-
 #[inline]
 pub fn fmodf(x: f32, y: f32) -> f32 {
     let mut uxi = x.to_bits();
@@ -11,7 +9,7 @@ pub fn fmodf(x: f32, y: f32) -> f32 {
     let sx = uxi & 0x80000000;
     let mut i;
 
-    if uyi << 1 == 0 || isnanf(y) || ex == 0xff {
+    if uyi << 1 == 0 || y.is_nan() || ex == 0xff {
         return (x * y) / (x * y);
     }
 

library/compiler-builtins/libm/src/math/mod.rs

@@ -7,49 +7,31 @@ macro_rules! force_eval {
 }
 
 mod ceilf;
+mod expf;
 mod fabs;
 mod fabsf;
+mod floor;
 mod floorf;
+mod fmod;
 mod fmodf;
+mod hypot;
+mod hypotf;
+mod logf;
 mod powf;
 mod round;
 mod roundf;
 mod scalbn;
 mod scalbnf;
 mod sqrt;
 mod sqrtf;
-mod logf;
-mod expf;
-mod floor;
 mod trunc;
 mod truncf;
-mod hypot;
-mod hypotf;
 
 //mod service;
 
 pub use self::{
-    ceilf::ceilf,
-    fabs::fabs,
-    fabsf::fabsf,
-    floorf::floorf,
-    fmodf::fmodf,
-    powf::powf,
-    round::round,
-    roundf::roundf,
-    scalbn::scalbn,
-    scalbnf::scalbnf,
-    sqrt::sqrt,
-    sqrtf::sqrtf,
-    logf::logf,
-    expf::expf,
-    floor::floor,
-    trunc::trunc,
+    ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmod::fmod,
+    fmodf::fmodf, hypot::hypot, hypotf::hypotf, logf::logf, powf::powf, round::round,
+    roundf::roundf, scalbn::scalbn, scalbnf::scalbnf, sqrt::sqrt, sqrtf::sqrtf, trunc::trunc,
     truncf::truncf,
-    hypot::hypot,
-    hypotf::hypotf,
 };
-
-fn isnanf(x: f32) -> bool {
-    x.to_bits() & 0x7fffffff > 0x7f800000
-}

library/compiler-builtins/libm/test-generator/src/main.rs

@@ -724,7 +724,7 @@ f64_f64! {
 f64f64_f64! {
     // atan2,
     // fdim,
-    // fmod,
+    fmod,
     hypot,
     // pow,
 }