improve impl of __floorh, __floorx, __ceilh and __ceilx

lib/compiler_rt/ceil.zig

@@ -3,6 +3,7 @@
 //!
 //! https://git.musl-libc.org/cgit/musl/tree/src/math/ceilf.c
 //! https://git.musl-libc.org/cgit/musl/tree/src/math/ceil.c
+//! https://git.musl-libc.org/cgit/musl/tree/src/math/ceill.c
 
 const std = @import("std");
 const builtin = @import("builtin");
@@ -27,8 +28,23 @@ comptime {
 }
 
 pub fn __ceilh(x: f16) callconv(.c) f16 {
-    // TODO: more efficient implementation
+    var u: u16 = @bitCast(x);
-    return @floatCast(ceilf(x));
+    const e = @as(i16, @intCast((u >> 10) & 31)) - 15;
+    var m: u16 = undefined;
+
+    if (e >= 10) return x;
+
+    if (e >= 0) {
+        m = @as(u16, 0x03FF) >> @intCast(e);
+        if (u & m == 0) return x;
+        if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1.0p120);
+        if (u >> 15 == 0) u += m;
+        u &= ~m;
+        return @bitCast(u);
+    } else {
+        if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1.0p120);
+        return if (u >> 15 != 0) -0.0 else if (u << 1 != 0) 1.0 else x;
+    }
 }
 
 pub fn ceilf(x: f32) callconv(.c) f32 {
@@ -36,31 +52,18 @@ pub fn ceilf(x: f32) callconv(.c) f32 {
     const e = @as(i32, @intCast((u >> 23) & 0xFF)) - 0x7F;
     var m: u32 = undefined;
 
-    // TODO: Shouldn't need this explicit check.
+    if (e >= 23) return x;
-    if (x == 0.0) {
-        return x;
-    }
 
-    if (e >= 23) {
+    if (e >= 0) {
-        return x;
-    } else if (e >= 0) {
         m = @as(u32, 0x007FFFFF) >> @intCast(e);
-        if (u & m == 0) {
+        if (u & m == 0) return x;
-            return x;
-        }
         if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1.0p120);
-        if (u >> 31 == 0) {
+        if (u >> 31 == 0) u += m;
-            u += m;
-        }
         u &= ~m;
         return @bitCast(u);
     } else {
         if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1.0p120);
-        if (u >> 31 != 0) {
+        return if (u >> 31 != 0) -0.0 else if (u << 1 != 0) 1.0 else x;
-            return -0.0;
-        } else {
-            return 1.0;
-        }
     }
 }
 
@@ -96,8 +99,32 @@ pub fn ceil(x: f64) callconv(.c) f64 {
 }
 
 pub fn __ceilx(x: f80) callconv(.c) f80 {
-    // TODO: more efficient implementation
+    const f80_toint = 1.0 / math.floatEps(f80);
-    return @floatCast(ceilq(x));
+
+    const u: u80 = @bitCast(x);
+    const e = (u >> 64) & 0x7FFF;
+    var y: f80 = undefined;
+
+    if (e >= 0x3FFF + 64 or x == 0) return x;
+
+    if (u >> 79 != 0) {
+        y = x - f80_toint + f80_toint - x;
+    } else {
+        y = x + f80_toint - f80_toint - x;
+    }
+
+    if (e <= 0x3FFF - 1) {
+        if (common.want_float_exceptions) mem.doNotOptimizeAway(y);
+        if (u >> 79 != 0) {
+            return -0.0;
+        } else {
+            return 1.0;
+        }
+    } else if (y < 0) {
+        return x + y + 1;
+    } else {
+        return x + y;
+    }
 }
 
 pub fn ceilq(x: f128) callconv(.c) f128 {
@@ -140,6 +167,12 @@ pub fn ceill(x: c_longdouble) callconv(.c) c_longdouble {
     }
 }
 
+test "ceil16" {
+    try expect(__ceilh(1.3) == 2.0);
+    try expect(__ceilh(-1.3) == -1.0);
+    try expect(__ceilh(0.2) == 1.0);
+}
+
 test "ceil32" {
     try expect(ceilf(1.3) == 2.0);
     try expect(ceilf(-1.3) == -1.0);
@@ -152,12 +185,26 @@ test "ceil64" {
     try expect(ceil(0.2) == 1.0);
 }
 
+test "ceil80" {
+    try expect(__ceilx(1.3) == 2.0);
+    try expect(__ceilx(-1.3) == -1.0);
+    try expect(__ceilx(0.2) == 1.0);
+}
+
 test "ceil128" {
     try expect(ceilq(1.3) == 2.0);
     try expect(ceilq(-1.3) == -1.0);
     try expect(ceilq(0.2) == 1.0);
 }
 
+test "ceil16.special" {
+    try expect(__ceilh(0.0) == 0.0);
+    try expect(__ceilh(-0.0) == -0.0);
+    try expect(math.isPositiveInf(__ceilh(math.inf(f16))));
+    try expect(math.isNegativeInf(__ceilh(-math.inf(f16))));
+    try expect(math.isNan(__ceilh(math.nan(f16))));
+}
+
 test "ceil32.special" {
     try expect(ceilf(0.0) == 0.0);
     try expect(ceilf(-0.0) == -0.0);
@@ -174,6 +221,14 @@ test "ceil64.special" {
     try expect(math.isNan(ceil(math.nan(f64))));
 }
 
+test "ceil80.special" {
+    try expect(__ceilx(0.0) == 0.0);
+    try expect(__ceilx(-0.0) == -0.0);
+    try expect(math.isPositiveInf(__ceilx(math.inf(f80))));
+    try expect(math.isNegativeInf(__ceilx(-math.inf(f80))));
+    try expect(math.isNan(__ceilx(math.nan(f80))));
+}
+
 test "ceil128.special" {
     try expect(ceilq(0.0) == 0.0);
     try expect(ceilq(-0.0) == -0.0);

lib/compiler_rt/floor.zig

@@ -3,6 +3,7 @@
 //!
 //! https://git.musl-libc.org/cgit/musl/tree/src/math/floorf.c
 //! https://git.musl-libc.org/cgit/musl/tree/src/math/floor.c
+//! https://git.musl-libc.org/cgit/musl/tree/src/math/floorl.c
 
 const std = @import("std");
 const builtin = @import("builtin");
@@ -31,32 +32,17 @@ pub fn __floorh(x: f16) callconv(.c) f16 {
     const e = @as(i16, @intCast((u >> 10) & 31)) - 15;
     var m: u16 = undefined;
 
-    // TODO: Shouldn't need this explicit check.
+    if (e >= 10) return x;
-    if (x == 0.0) {
-        return x;
-    }
-
-    if (e >= 10) {
-        return x;
-    }
 
     if (e >= 0) {
-        m = @as(u16, 1023) >> @intCast(e);
+        m = @as(u16, 0x03FF) >> @intCast(e);
-        if (u & m == 0) {
+        if (u & m == 0) return x;
-            return x;
-        }
         if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1.0p120);
-        if (u >> 15 != 0) {
+        if (u >> 15 != 0) u += m;
-            u += m;
-        }
         return @bitCast(u & ~m);
     } else {
         if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1.0p120);
-        if (u >> 15 == 0) {
+        return if (u >> 15 == 0) 0.0 else if (u << 1 != 0) -1.0 else x;
-            return 0.0;
-        } else {
-            return -1.0;
-        }
     }
 }
 
@@ -65,32 +51,17 @@ pub fn floorf(x: f32) callconv(.c) f32 {
     const e = @as(i32, @intCast((u >> 23) & 0xFF)) - 0x7F;
     var m: u32 = undefined;
 
-    // TODO: Shouldn't need this explicit check.
+    if (e >= 23) return x;
-    if (x == 0.0) {
-        return x;
-    }
-
-    if (e >= 23) {
-        return x;
-    }
 
     if (e >= 0) {
         m = @as(u32, 0x007FFFFF) >> @intCast(e);
-        if (u & m == 0) {
+        if (u & m == 0) return x;
-            return x;
-        }
         if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1.0p120);
-        if (u >> 31 != 0) {
+        if (u >> 31 != 0) u += m;
-            u += m;
-        }
         return @bitCast(u & ~m);
     } else {
         if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1.0p120);
-        if (u >> 31 == 0) {
+        return if (u >> 31 == 0) 0.0 else if (u << 1 != 0) -1.0 else x;
-            return 0.0;
-        } else {
-            return -1.0;
-        }
     }
 }
 
@@ -126,8 +97,34 @@ pub fn floor(x: f64) callconv(.c) f64 {
 }
 
 pub fn __floorx(x: f80) callconv(.c) f80 {
-    // TODO: more efficient implementation
+    const f80_toint = 1.0 / math.floatEps(f80);
-    return @floatCast(floorq(x));
+
+    const u: u80 = @bitCast(x);
+    const e = (u >> 64) & 0x7FFF;
+    var y: f80 = undefined;
+
+    if (e >= 0x3FFF + 64 or x == 0) {
+        return x;
+    }
+
+    if (u >> 79 != 0) {
+        y = x - f80_toint + f80_toint - x;
+    } else {
+        y = x + f80_toint - f80_toint - x;
+    }
+
+    if (e <= 0x3FFF - 1) {
+        if (common.want_float_exceptions) mem.doNotOptimizeAway(y);
+        if (u >> 79 != 0) {
+            return -1.0;
+        } else {
+            return 0.0;
+        }
+    } else if (y > 0) {
+        return x + y - 1;
+    } else {
+        return x + y;
+    }
 }
 
 pub fn floorq(x: f128) callconv(.c) f128 {
@@ -188,6 +185,12 @@ test "floor64" {
     try expect(floor(0.2) == 0.0);
 }
 
+test "floor80" {
+    try expect(__floorx(1.3) == 1.0);
+    try expect(__floorx(-1.3) == -2.0);
+    try expect(__floorx(0.2) == 0.0);
+}
+
 test "floor128" {
     try expect(floorq(1.3) == 1.0);
     try expect(floorq(-1.3) == -2.0);
@@ -218,6 +221,14 @@ test "floor64.special" {
     try expect(math.isNan(floor(math.nan(f64))));
 }
 
+test "floor80.special" {
+    try expect(__floorx(0.0) == 0.0);
+    try expect(__floorx(-0.0) == -0.0);
+    try expect(math.isPositiveInf(__floorx(math.inf(f80))));
+    try expect(math.isNegativeInf(__floorx(-math.inf(f80))));
+    try expect(math.isNan(__floorx(math.nan(f80))));
+}
+
 test "floor128.special" {
     try expect(floorq(0.0) == 0.0);
     try expect(floorq(-0.0) == -0.0);