aljsndlkjasn

code/math.c

@@ -0,0 +1,42 @@
+
+#include <stdint.h>
+
+#define ftype        float
+#define itype        uint32_t
+#define f_bits       32
+#define f_sig_bits   24
+#define f_exp_max    127
+#define suffix(name) #name ## "f"
+#include "math/floats.h"
+#include "math/basic.h"
+#include "math/exponential.h"
+#include "math/hyperbolic.h"
+#include "math/nearest-int.h"
+#include "math/power.h"
+#include "math/probability-statistics.h"
+#undef ftype
+#undef itype
+#undef f_bits
+#undef f_sig_bits
+#undef f_exp_max
+#undef suffix
+
+#define ftype        double
+#define itype        uint64_t
+#define f_bits       64
+#define f_sig_bits   53
+#define f_exp_max    1023
+#define suffix(name) #name
+#include "math/basic.h"
+#include "math/exponential.h"
+#include "math/floats.h"
+#include "math/hyperbolic.h"
+#include "math/nearest-int.h"
+#include "math/power.h"
+#include "math/probability-statistics.h"
+#undef ftype
+#undef itype
+#undef f_bits
+#undef f_sig_bits
+#undef f_exp_max
+#undef suffix

code/math/basic.h

@@ -1,14 +1,5 @@
 
-#include <math.h>
-
-double fabs(double x) {
-    if(x >= 0) {
-        return x;
-    }
-    return -x;
-}
-
-float fabsf(float x) {
+ftype fabs(ftype x) {
     if(x >= 0) {
         return x;
     }

code/math/floats.h

@@ -0,0 +1,45 @@
+
+#include <math.h>
+
+#define f_man_bits (f_sig_bits - 1)
+#define f_exp_bits (f_bits - f_man_bits - 1)
+
+#define b_sign_val(b) ((b) >> (f_exp_bits + f_man_bits))
+#define b_bexp_val(b) (((b) >> f_man_bits) & f_exp_mask)
+#define b_mant_val(b) ((b) & f_man_mask)
+#define b_mant_mask   ((1 << f_man_bits) - 1)
+#define b_bexp_mask   ((1 << f_exp_bits) - 1)
+#define b_mant_max    ((1 << f_man_bits) - 1)
+#define b_bexp_max    ((1 << f_exp_bits) - 1)
+
+#define b_exp_bias    ((1 << (f_exp_bits-1)) - 1)
+#define b_exp_val(b)  (b_bexp_val(b) - b_exp_bias)
+
+static inline itype suffix(getbits)(ftype f) {
+    union _f {
+        ftype f;
+        itype i;
+    } u;
+    u.f = f;
+    int bits = u.i;
+    return bits;
+}
+
+int suffix(_fpclassify)(ftype f) {
+    itype bits = suffix(getbits)(f);
+    itype sign = b_sign_val(bits);
+    itype bexp = b_bexp_val(bits);
+    itype mant = b_mant_val(bits);
+    if(bexp == b_bexp_max) {
+        if(mant == 0) return FP_INFINITE;
+        else          return FP_NAN;
+    }
+    else if(bexp == 0) {
+        if(mant == 0) return FP_ZERO;
+        else          return FP_SUBNORMAL;
+    }
+    else {
+        return FP_NORMAL;
+    }
+}
+

inc/math.h

@@ -27,10 +27,7 @@ typedef double double_t;
 #define FP_NORMAL    2
 #define FP_SUBNORMAL 3
 #define FP_ZERO      4
-#define fpclassify(x) ( \
-    _is_float(x)  ? _fpclassify_f(x) :
-    _is_double(x) ? _fpclassify_d(x) :
-    FP_NAN)
+#define fpclassify(x) (_is_float(x)  ? _fpclassifyf(x) : _fpclassify(x))
 
 double acos(double x);
 float acosf(float x);