From 8c1bb593ce7aee9b64eb88453f173c59686f5f2b Mon Sep 17 00:00:00 2001
From: "Justin C. Miller" <justin@devjustinian.com>
Date: Thu, 12 Jan 2023 21:51:36 -0800
Subject: [PATCH] [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 |  6 +--
 src/libraries/libc/math/ceil.cpp         |  4 ++
 src/libraries/libc/math/ceilf.cpp        |  4 ++
 src/libraries/libc/math/ceill.cpp        |  4 ++
 src/libraries/libc/math/float_ops.h      | 59 ++++++++++++++++++++++++
 src/libraries/libc/math/frexp.cpp        |  4 ++
 src/libraries/libc/math/frexpf.cpp       |  4 ++
 src/libraries/libc/math/frexpl.cpp       |  4 ++
 8 files changed, 84 insertions(+), 5 deletions(-)
 create mode 100644 src/libraries/libc/math/ceil.cpp
 create mode 100644 src/libraries/libc/math/ceilf.cpp
 create mode 100644 src/libraries/libc/math/ceill.cpp
 create mode 100644 src/libraries/libc/math/float_ops.h
 create mode 100644 src/libraries/libc/math/frexp.cpp
 create mode 100644 src/libraries/libc/math/frexpf.cpp
 create mode 100644 src/libraries/libc/math/frexpl.cpp

diff --git a/src/libraries/libc/math/builtins.cpp.cog b/src/libraries/libc/math/builtins.cpp.cog
index 79fef92..674ae6b 100644
--- a/src/libraries/libc/math/builtins.cpp.cog
+++ b/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); }
diff --git a/src/libraries/libc/math/ceil.cpp b/src/libraries/libc/math/ceil.cpp
new file mode 100644
index 0000000..88db15a
--- /dev/null
+++ b/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); }
diff --git a/src/libraries/libc/math/ceilf.cpp b/src/libraries/libc/math/ceilf.cpp
new file mode 100644
index 0000000..dad1776
--- /dev/null
+++ b/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); }
diff --git a/src/libraries/libc/math/ceill.cpp b/src/libraries/libc/math/ceill.cpp
new file mode 100644
index 0000000..8d518fd
--- /dev/null
+++ b/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; }
diff --git a/src/libraries/libc/math/float_ops.h b/src/libraries/libc/math/float_ops.h
new file mode 100644
index 0000000..dc39d1a
--- /dev/null
+++ b/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);
+}
diff --git a/src/libraries/libc/math/frexp.cpp b/src/libraries/libc/math/frexp.cpp
new file mode 100644
index 0000000..a5406f3
--- /dev/null
+++ b/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); }
diff --git a/src/libraries/libc/math/frexpf.cpp b/src/libraries/libc/math/frexpf.cpp
new file mode 100644
index 0000000..51a8858
--- /dev/null
+++ b/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); }
diff --git a/src/libraries/libc/math/frexpl.cpp b/src/libraries/libc/math/frexpl.cpp
new file mode 100644
index 0000000..7f1f6d8
--- /dev/null
+++ b/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; }