[libc] Add ciel, frexpr implementations

The clang __builtin_* functions cannot be relied upon, as they may just
emit a call to the stdlib version. So this commit adds an implementation
for ceil and frexpr, as well as their float versions.

src/libraries/libc/math/builtins.cpp.cog

@@ -4,7 +4,7 @@
 /*[[[cog code generation
 
 builtins_single = [
-    "acos", "asin", "atan", "ceil", "cos", "cosh", "exp", "fabs",
+    "acos", "asin", "atan", "cos", "cosh", "exp", "fabs",
     "floor", "log10", "log", "sin", "sinh", "sqrt", "tan", "tanh",
     ]
 
@@ -29,10 +29,6 @@ for fname in builtins_double:
 ]]]*/
 /// [[[end]]]
 
-double frexp(double value, int *exp)                    { return __builtin_frexp(value, exp); }
-float frexpf(float value, int *exp)                     { return __builtin_frexpf(value, exp); }
-long double frexpl(long double value, int *exp)         { return __builtin_frexpl(value, exp); }
-
 double ldexp(double x, int exp)                         { return __builtin_ldexp(x, exp); }
 float ldexpf(float x, int exp)                          { return __builtin_ldexpf(x, exp); }
 long double ldexpl(long double x, int exp)              { return __builtin_ldexpl(x, exp); }

src/libraries/libc/math/ceil.cpp

@@ -0,0 +1,4 @@
+#include "float_ops.h"
+
+extern "C"
+double ceil(double f) { return __ceil<double_precision>(f); }

src/libraries/libc/math/ceilf.cpp

@@ -0,0 +1,4 @@
+#include "float_ops.h"
+
+extern "C"
+float ceilf(float f) { return __ceil<single_precision>(f); }

src/libraries/libc/math/ceill.cpp

@@ -0,0 +1,4 @@
+#include <assert.h>
+
+extern "C"
+long double ceill(long double f) { assert(false && "NYI"); return 0; }

src/libraries/libc/math/float_ops.h

@@ -0,0 +1,59 @@
+#pragma once
+#include <stddef.h>
+#include <stdint.h>
+
+struct single_precision
+{
+    using type = float;
+    static constexpr size_t sig_bits = 23;
+    static constexpr size_t exp_bits =  8;
+};
+
+struct double_precision
+{
+    using type = double;
+    static constexpr size_t sig_bits = 52;
+    static constexpr size_t exp_bits = 11;
+};
+
+template<typename Traits>
+inline typename Traits::type __ceil(typename Traits::type f) {
+    static constexpr int64_t exp_mask = (1ll<<Traits::exp_bits) - 1;
+    static constexpr int64_t exp_mid = (1ll<<(Traits::exp_bits-1)) - 1;
+
+    uint64_t bits = *reinterpret_cast<uint64_t*>(&f);
+    int64_t exponent = static_cast<int64_t>((bits >> Traits::sig_bits) & exp_mask) - exp_mid;
+    if (exponent < 0)
+        return f > 0 ? 1 : 0;
+
+    int64_t fraction_bits = Traits::sig_bits - exponent;
+    if (fraction_bits <= 0)
+        return f;
+
+    uint64_t nonfraction_mask = -1ull << fraction_bits;
+    uint64_t nonfraction = bits & nonfraction_mask;
+    typename Traits::type rounded = *reinterpret_cast<typename Traits::type*>(&nonfraction);
+
+    if (rounded > 0 && bits != nonfraction)
+        rounded += 1;
+
+    return rounded;
+}
+
+template<typename Traits>
+inline typename Traits::type __frexp(typename Traits::type f, int *exp) {
+    static constexpr int64_t exp_mask = ((1ll<<Traits::exp_bits) - 1) << Traits::sig_bits;
+    static constexpr int64_t exp_mid = (1ll<<(Traits::exp_bits-1)) - 1;
+
+    uint64_t bits = *reinterpret_cast<uint64_t*>(&f);
+
+    if (bits == 0) {
+        *exp = 0;
+        return 0;
+    }
+
+    *exp = ((bits & exp_mask) >> Traits::sig_bits) - exp_mid;
+
+    uint64_t result = (bits & ~exp_mask) | (exp_mid << Traits::sig_bits);
+    return *reinterpret_cast<typename Traits::type*>(&result);
+}

src/libraries/libc/math/frexp.cpp

@@ -0,0 +1,4 @@
+#include "float_ops.h"
+
+extern "C"
+double frexp(double f, int *exp) { return __frexp<double_precision>(f, exp); }

src/libraries/libc/math/frexpf.cpp

@@ -0,0 +1,4 @@
+#include "float_ops.h"
+
+extern "C"
+float frexpf(float f, int *exp) { return __frexp<single_precision>(f, exp); }

src/libraries/libc/math/frexpl.cpp

@@ -0,0 +1,4 @@
+#include <assert.h>
+
+extern "C"
+long double frexpl(long double f, int *exp) { assert(false && "NYI"); *exp = 0; return 0; }