Some ieee754 funcs of math

code/math.c

@@ -1,42 +0,0 @@
-
-#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,25 +0,0 @@
-
-ftype fabs(ftype x) {
-    if(x >= 0) {
-        return x;
-    }
-    return -x;
-}
-
-double fmod(double x, double y) {
-    int n = 0;
-    while(y < x*n) {
-        ++n;
-    }
-    return y - n*x;
-}
-
-float fmodf(float x, float y) {
-    int n = 0;
-    while(y < x*n) {
-        ++n;
-    }
-    return y - n*x;
-}
-
-

code/math/floats.h

@@ -1,45 +0,0 @@
-
-#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;
-    }
-}
-

code/math/ieee754.h

@@ -0,0 +1,134 @@
+
+// Instantiates 'template' macros for floating-point values
+// When including expecting the following parameters to be defined:
+//  ftype: Floating-point type to consider
+//  itype: Signed integer type corresponding to ftype's bitwidth
+//  f_ebits: Number of bits in the exponent
+//  f_mbits: Number of bits in the mantissa
+//  suffix(name): appends corresponding suffix to the given name,
+//                e.g. f for floats
+
+#define f_nbits (1+f_ebits+f_mbits)
+#define f_emax ((1ULL << (f_ebits-1)) - 1)
+#define f_emin (1 - f_emax)
+
+// Extracting fields from the float
+#define f_eoffs (f_mbits)
+#define f_soffs (f_mbits+f_ebits)
+#define f_emask ((1ULL << f_ebits) - 1)
+#define f_mmask ((1ULL << f_mbits) - 1)
+#define f_smask 1ULL
+#define f_eval(b) (((b) >> f_eoffs) & f_emask)
+#define f_sval(b) (((b) >> f_soffs) & f_smask)
+#define f_mval(b) (((b) >> 0) & f_mmask)
+#define f_abs(b)  ((b) & ~(f_smask << f_soffs))
+
+#define b_cons(s,e,m) ((s << f_soffs) | (e << f_eoffs) | (m))
+
+// Converting float to integer bits
+static inline itype suffix(f_bits)(ftype f) {
+    union _u {
+        ftype f;
+        itype b;
+    } u;
+    u.f = f;
+    return u.b;
+}
+
+static inline itype suffix(f_frombits)(itype b) {
+        union _u {
+        ftype f;
+        itype b;
+    } u;
+    u.b = b;
+    return u.f;
+}
+
+// Floating-point classification
+int suffix(_fpclassify)(ftype f) {
+    itype bits = suffix(f_bits)(f);
+    itype exp = f_eval(bits);
+    itype man = f_mval(bits);
+    if(exp == f_emax) {
+        if(man == 0) return FP_INFINITE;
+        else return FP_NAN;
+    }
+    else if(exp == 0) {
+        if(man == 0) return FP_ZERO;
+        else return FP_SUBNORMAL;
+    }
+    else return FP_NORMAL;
+}
+
+int suffix(_signbit)(ftype f) {
+    itype bits = suffix(f_bits)(f);
+    itype sign = f_sval(bits);
+    return sign;
+}
+
+ftype suffix(copysign)(ftype x, ftype y) {
+    itype xbits = suffix(f_bits)(x);
+    itype ybits = suffix(f_bits)(y);
+    itype exp = f_eval(xbits);
+    itype man = f_mval(xbits);
+    itype sgn = f_sval(ybits);
+    itype rbits = b_cons(sgn, exp, man);
+    return suffix(f_frombits)(rbits);
+}
+
+// Floating-point non-signaling comparison
+
+static Ordering suffix(_ordering)(ftype x, ftype y) {
+    itype xclass = suffix(_fpclassify)(x);
+    itype yclass = suffix(_fpclassify)(y);
+    if(xclass == FP_NAN || yclass == FP_NAN) {
+        return UN;
+    }
+    itype xbits = suffix(f_bits)(x);
+    itype ybits = suffix(f_bits)(y);
+    itype xsgn = f_sval(xbits);
+    itype ysgn = f_sval(ybits);
+    itype xabs = f_abs(xbits);
+    itype yabs = f_abs(ybits);
+    if(xsgn == ysgn) {
+        if(xabs > yabs) return GR;
+        if(xabs < yabs) return LS;
+        return EQ;
+    }
+    else {
+        if(xabs == 0 && yabs == 0) return EQ;
+        if(xsgn) return LS;
+        if(ysgn) return GR;
+    }
+    return UN; // I may be stupid
+}
+
+int suffix(_isgrt)(ftype x, ftype y) {
+    int ord = suffix(_ordering)(x, y);
+    return ord == GR;
+}
+
+int suffix(_isgeq)(ftype x, ftype y) {
+    int ord = suffix(_ordering)(x, y);
+    return ord == GR || ord == EQ;
+}
+
+int suffix(_isles)(ftype x, ftype y) {
+    int ord = suffix(_ordering)(x, y);
+    return ord == LS;
+}
+
+int suffix(_isleq)(ftype x, ftype y) {
+    int ord = suffix(_ordering)(x, y);
+    return ord == LS || ord == EQ;
+}
+
+int suffix(_isleg)(ftype x, ftype y) {
+    int ord = suffix(_ordering)(x, y);
+    return ord == LS || ord == GR;
+}
+
+int suffix(_isuno)(ftype x, ftype y) {
+    int ord = suffix(_ordering)(x, y);
+    return ord == UN;
+}

code/math/math.c

@@ -0,0 +1,50 @@
+
+#include <math.h>
+#include <stdint.h>
+
+// Used for float comparisons in ieee754.h
+enum Ordering {
+    LS,
+    EQ,
+    GR,
+    UN,
+} typedef Ordering;
+
+#define ftype   float
+#define itype   int32_t
+#define f_ebits 8
+#define f_mbits 23
+#define suffix(n) n ## f
+#include "ieee754.h"
+#undef suffix
+#undef f_mbits
+#undef f_ebits
+#undef itype
+#undef ftype
+
+#define ftype   double
+#define itype   int64_t
+#define f_ebits 11
+#define f_mbits 52
+#define suffix(n) n
+#include "ieee754.h"
+#undef suffix
+#undef f_mbits
+#undef f_ebits
+#undef itype
+#undef ftype
+
+_Static_assert(sizeof(long double) == sizeof(double),
+    "Are these 'long doubles' in the same room with us right now?");
+
+#define ftype   long double
+#define itype   int64_t
+#define f_ebits 11
+#define f_mbits 52
+#define suffix(n) n ## l
+//#include "ieee754.h"
+#undef suffix
+#undef f_mbits
+#undef f_ebits
+#undef itype
+#undef ftype

code/stdlib.c

@@ -1,5 +1,41 @@
 #include <stdlib.h>
 
+int abs(int n) {
+    if(n < 0) return -n;
+    return n;
+}
+
+long absl(long n) {
+    if(n < 0) return -n;
+    return n;
+}
+
+long long absll(long long n) {
+    if(n < 0) return -n;
+    return n;
+}
+
+div_t div(int x, int y) {
+    div_t res;
+    res.quot = x / y;
+    res.rem = x % y;
+    return res;
+}
+
+ldiv_t ldiv(long x, long y) {
+    ldiv_t res;
+    res.quot = x / y;
+    res.rem = x % y;
+    return res;
+}
+
+lldiv_t lldiv(long long x, long long y) {
+    lldiv_t res;
+    res.quot = x / y;
+    res.rem = x % y;
+    return res;
+}
+
 const void *bsearch(const void *key, const void *base, size_t nmemb, size_t size, int (*compar)(const void *, const void *)) {
     size_t left = 0;
     size_t right = nmemb;

inc/math.h

@@ -17,21 +17,41 @@ typedef double double_t;
 // FP_ILOGB0
 // FP_ILOGBNAN
 
-#define MATH_ERRNO     0
+#define MATH_ERRNO     1
-#define MATH_ERREXCEPT 1
+#define MATH_ERREXCEPT 2
 
 #define math_errhandling MATH_ERRNO
 
+// Classification
 #define FP_INFINITE  0
 #define FP_NAN       1
 #define FP_NORMAL    2
 #define FP_SUBNORMAL 3
 #define FP_ZERO      4
+int     _fpclassify(double f);
+int     _fpclassifyf(float f);
+#define fpclassify(x) (sizeof(x)==4?_fpclassifyf(x):_fpclassify(x))
+#define isfinite(x)   (fpclassify(x) != FP_INFINITE && fpclassify(x) != FP_NAN)
+#define isinf(x)      (fpclassify(x) == FP_INFINITE)
+#define isnan(x)      (fpclassify(x) == FP_NAN)
+#define isnormal(x)   (fpclassify(x) == FP_NORMAL)
 
-#define fpclassify(x) (_is_float(x)  ? _fpclassifyf(x) : _fpclassify(x))
+// Sign bit shit
+int     _signbit(double f);
+int     _signbitf(float f);
+#define signbit(x)    (sizeof(x)==4?_signbitf(x):signbit(x))
+float   copysignf(float x, float y);
 
-double acos(double x);
+// Ordering
-float acosf(float x);
+#define isgreater(x)      (sizeof(x)==4?_isgrtf(x):_isgrt(x))
+#define isgreaterequal(x) (sizeof(x)==4?_isgeqf(x):_isgeq(x))
+#define isless(x)         (sizeof(x)==4?_islesf(x):_isles(x))
+#define islessequal(x)    (sizeof(x)==4?_isleqf(x):_isleq(x))
+#define islessgreater(x)  (sizeof(x)==4?_islegf(x):_isleg(x))
+#define isunordered(x)    (sizeof(x)==4?_isunof(x):_isuno(x))
+
+double  acos(double x);
+float   acosf(float x);
 
 double asin(double x);
 float asinf(float x);

inc/stdio.h

@@ -6,6 +6,10 @@
 typedef struct FILE FILE;
 typedef int64_t fpos_t;
 
+typedef struct div_t div_t;
+typedef struct ldiv_t ldiv_t;
+typedef struct lldiv_t lldiv_t;
+
 #if !defined(NULL)
     #define NULL ((void *)0)
 #endif

inc/stdlib.h

@@ -15,20 +15,20 @@
 #define EXIT_SUCCESS 0
 #define EXIT_FAILURE 1
 
-// typedef struct div_t {
+typedef struct div_t {
-//     int quot;
+    int quot;
-//     int rem;
+    int rem;
-// } div_t;
+} div_t;
 
-// typedef struct ldiv_t {
+typedef struct ldiv_t {
-//     long quot;
+    long quot;
-//     long rem;
+    long rem;
-// } ldiv_t;
+} ldiv_t;
 
-// typedef struct lldiv_t {
+typedef struct lldiv_t {
-//     long long quot;
+    long long quot;
-//     long long rem;
+    long long rem;
-// } lldiv_t;
+} lldiv_t;
 
 // #define EXIT_FAILURE 1
 // #define EXIT_SUCCESS 0