branchless fpclassify

bake-test.cmd

@@ -1 +0,0 @@
-clang test\test_%test%.c utf8.obj -Iinc -g -lciabatta.lib

inc/math.h

@@ -23,11 +23,11 @@ typedef double double_t;
 #define math_errhandling (MATH_ERRNO | MATH_ERREXCEPT)
 
 // Classification
-#define FP_INFINITE  0
-#define FP_NAN       1
+#define FP_ZERO      0
+#define FP_SUBNORMAL 1
 #define FP_NORMAL    2
-#define FP_SUBNORMAL 3
-#define FP_ZERO      4
+#define FP_INFINITE  4
+#define FP_NAN       5
 int     _fpclassify(double);
 int     _fpclassifyf(float);
 int     _fpclassifyl(long double);

inc/stdio.h

@@ -27,7 +27,7 @@
 typedef struct FILE FILE;
 
 typedef struct {
-    int64_t   offset;
+    unsigned long long offset;
     mbstate_t mbstate;
 } fpos_t;
 

src/ciabatta.c

@@ -46,7 +46,7 @@
 #include "math/basic.c"
 #include "math/division.c"
 #include "math/gen_math.c"
-#include "math/ieee754.c"
+#include "math/bits.c"
 #include "math/round.c"
 #include "conv/digits.c"
 #include "conv/strpfx.c"

src/math/bits.c

@@ -0,0 +1,106 @@
+
+// This thing doesn't compile branchlessly on MSVC lmao
+// This is a branchless version of fpclassify
+// The way it works is by assuming the following structure of float class
+// bit 2 | bit 1 | bit 0
+//   e       n      m
+// e is 1 for infinities and nan, 0 otherwise
+// n is 1 for normal numbers, 0 otherwise
+// m is 1 for nan and subnormal numbers, 0 otherwise
+// This leaves the following values for classes:
+//  0 - FP_ZERO
+//  1 - FP_SUBNORMAL
+//  2 - FP_NORMAL
+//  4 - FP_INFINITE
+//  5 - FP_NAN
+int _fpclassify(f64 x) {
+    // First, extract bits
+    u64 bits = F64_BITS(x);
+    i64 exp  = F64_BEXP(bits);
+    i64 mant = F64_MANT(bits);
+    // Get the number that's only zero when exp = 0x7ff
+    i64 ee = F64_BEXP_MAX - exp;
+    // Invert it to get number that's 1 iff x is infinity or nan
+    i64 e = !(ee);
+    // The value (ee*exp) is zero in two cases:
+    //   exp = 0x7ff, (infinities, nan)
+    //   exp = 0 (zero, subnormal)
+    // We negate this to make it so that this value is 1 only if number is not
+    // normal
+    i64 nn = !(ee*exp);
+    // Negate the previous thing. Now n is 1 iff number is normal
+    i64 n = !nn;
+    // m is 1 if mantissa is nonzero and the number is not normal
+    i64 m = !!mant & nn;
+    // Construct the float class
+    return (e<<2) | (n<<1) | m;
+}
+
+// Same recipe as above, different constants
+int _fpclassifyf(float x) {
+    u64 bits = F32_BITS(x);
+    i64 exp  = F32_BEXP(bits);
+    i64 mant = F32_MANT(bits);
+    i64 ee = F32_BEXP_MAX - exp;
+    i64 e = !(ee);
+    i64 nn = !(ee*exp);
+    i64 n = !nn;
+    i64 m = !!mant & nn;
+    return (e<<2) | (n<<1) | m;
+}
+
+int _fpclassifyl(fl64 x) {
+    return _fpclassify((f64)x);
+}
+
+
+int _signbit(f64 x) {
+    union {
+        f64 d;
+        uint64_t i;
+    } y = { x };
+    return y.i>>63;
+}
+
+int _signbitf(float x) {
+    union {
+        float f;
+        uint32_t i;
+    } y = { x };
+    return y.i>>31;
+}
+
+int _signbitl(fl64 x) {
+    return _signbit(x);
+}
+
+f64 copysign(f64 x, f64 y) {
+    union {f64 f; uint64_t i;} ux={x}, uy={y};
+    ux.i &= ~(1ULL<<63);
+    ux.i |= uy.i & (1ULL<<63);
+    return ux.f;
+}
+
+float copysignf(float x, float y) {
+    union {float f; uint32_t i;} ux={x}, uy={y};
+    ux.i &= 0x7fffffff;
+    ux.i |= uy.i & 0x80000000;
+    return ux.f;
+}
+
+fl64 copysignl(fl64 x, fl64 y) {
+    return copysign(x, y);
+}
+
+f64 nan(const char *s) {
+    return NAN;
+}
+
+float nanf(const char *s) {
+    return NAN;
+}
+
+fl64 nanl(const char *s) {
+    return NAN;
+}
+

src/math/ieee754.c

@@ -1,77 +0,0 @@
-
-#include <math.h>
-#include <fenv.h>
-#include <stdint.h>
-#include <float.h>
-
-int _fpclassify(double x) {
-    union {double f; uint64_t i;} u = {x};
-    int e = u.i>>52 & 0x7ff;
-    if (!e) return u.i<<1 ? FP_SUBNORMAL : FP_ZERO;
-    if (e==0x7ff) return u.i<<12 ? FP_NAN : FP_INFINITE;
-    return FP_NORMAL;
-}
-
-int _fpclassifyf(float x) {
-    union {float f; uint32_t i;} u = {x};
-    int e = u.i>>23 & 0xff;
-    if (!e) return u.i<<1 ? FP_SUBNORMAL : FP_ZERO;
-    if (e==0xff) return u.i<<9 ? FP_NAN : FP_INFINITE;
-    return FP_NORMAL;
-}
-
-int _fpclassifyl(long double x) {
-    return _fpclassify(x);
-}
-
-
-int _signbit(double x) {
-    union {
-        double d;
-        uint64_t i;
-    } y = { x };
-    return y.i>>63;
-}
-
-int _signbitf(float x) {
-    union {
-        float f;
-        uint32_t i;
-    } y = { x };
-    return y.i>>31;
-}
-
-int _signbitl(long double x) {
-    return _signbit(x);
-}
-
-double copysign(double x, double y) {
-    union {double f; uint64_t i;} ux={x}, uy={y};
-    ux.i &= ~(1ULL<<63);
-    ux.i |= uy.i & (1ULL<<63);
-    return ux.f;
-}
-
-float copysignf(float x, float y) {
-    union {float f; uint32_t i;} ux={x}, uy={y};
-    ux.i &= 0x7fffffff;
-    ux.i |= uy.i & 0x80000000;
-    return ux.f;
-}
-
-long double copysignl(long double x, long double y) {
-    return copysign(x, y);
-}
-
-double nan(const char *s) {
-    return NAN;
-}
-
-float nanf(const char *s) {
-    return NAN;
-}
-
-long double nanl(const char *s) {
-    return NAN;
-}
-

src/util.c

@@ -44,8 +44,30 @@ typedef wchar_t wchar;
 #define STR_(a) #a
 #define STR(a) STR_(a)
 
-#define DOUBLE_BITS(x) ((union {double f; u64 i;}){x}).i
-#define DOUBLE_CONS(x) ((union {double f; u64 i;}){x}).f
-#define FLOAT_BITS(x)  ((union {float f;  u32 i;}){x}).i
-#define FLOAT_CONS(x)  ((union {float f;  u32 i;}){x}).f
+#define F64_BITS(x) ((union {f64 f; u64 i;}){x}).i
+#define F64_CONS(x) ((union {f64 f; u64 i;}){x}).f
+#define F32_BITS(x) ((union {f32 f; u32 i;}){x}).i
+#define F32_CONS(x) ((union {f32 f; u32 i;}){x}).f
 
+#define F64_MANT_MASK UINT64_C(0xfffffffffffff)
+#define F64_MANT_MAX  UINT64_C(0xfffffffffffff)
+#define F64_MANT_BITS 52
+#define F64_BEXP_BITS 11
+#define F64_BEXP_MASK 0x7ff
+#define F64_BEXP_MAX  0x7ff
+
+#define F32_MANT_MASK 0x7fffff
+#define F32_MANT_MAX  0x7fffff
+#define F32_MANT_BITS 23
+#define F32_BEXP_BITS 8
+#define F32_BEXP_MASK 0xff
+#define F32_BEXP_MAX  0xff
+
+
+#define F64_SIGN(bits) (bits >> (F64_MANT_BITS + F64_BEXP_BITS))
+#define F64_BEXP(bits) ((bits >> F64_MANT_BITS) & F64_BEXP_MASK)
+#define F64_MANT(bits) ((bits) & F64_MANT_MASK)
+
+#define F32_SIGN(bits) (bits >> (F32_MANT_BITS + F32_BEXP_BITS))
+#define F32_BEXP(bits) ((bits >> F32_MANT_BITS) & F32_BEXP_MASK)
+#define F32_MANT(bits) ((bits) & F32_MANT_MASK)

test.cmd

@@ -0,0 +1,2 @@
+
+clang test\%1 utf8.obj -Iinc -g -lciabatta.lib -nostdlib

test/math.c

@@ -5,6 +5,7 @@
 #include <fenv.h>
 #include <inttypes.h>
 #include <errno.h>
+#include <limits.h>
 
 const char *show_classification(double x) {
     switch(fpclassify(x)) {