Some rough floating point macros!

include/math.h

@@ -3,7 +3,7 @@
  * 
  * This header is a part of the FENIX C Library and is free software.
  * You can redistribute and/or modify it subject to the terms of the
- * Clumsy Wolf Public License v4. For more details, see the file COPYING.
+ * Clumsy Wolf Public License v5. For more details, see the file COPYING.
  *
  * The FENIX C Library is distributed WITH NO WARRANTY WHATSOEVER. See
  * The CWPL for more details.
@@ -12,6 +12,8 @@
 #ifndef _MATH_H
 #define _MATH_H
 
+#include <limits.h>
+
 #if FLT_EVAL_METHOD == 1
 typedef double float_t;
 typedef double double_t;
@@ -23,6 +25,44 @@ typedef float float_t;
 typedef double double_t;
 #endif
 
+/* Both of these are basically lifted from musl. */
+#define INFINITY (1e5000)
+#define NAN (0.0f / 0.0f)
+
+#define FP_INFINITE 1
+#define FP_NAN 2
+#define FP_NORMAL 3
+#define FP_SUBNORMAL 4
+#define FP_ZERO 0
+
+int __fpclassifyd(double x);
+int __fpclassifyf(float x);
+int __fpclassifyl(long double x);
+
+/* We're just gonna use the approach in the C standard here */
+#define fpclassify(x) ((sizeof(x) == sizeof(float)) ? __fpclassifyf(x) : \
+                      (sizeof(x) == sizeof(double)) ? __fpclassifyd(x) : \
+                      __fpclassifyl(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)
+
+int __signbitd(double x);
+int __signbitf(float x);
+int __signbitl(long double x);
+
+/*
+  In theory, I could just use (x < 0 || x == -0.0f), but would that work for
+  infinity and NaN? I feel like no. So fuck it. We'll just write some small
+  functions to pull the sign bit straight from the number.
+  -Kat
+*/
+#define signbit(x) ((sizeof(x) == sizeof(float)) ? __signbitf(x) : \
+                    (sizeof(x) == sizeof(double)) ? __signbitd(x) : \
+                    __signbitl(x))
+
 #define M_E 2.718281828
 #define M_PI 3.14159265
 #define PI 3.14159265
@@ -40,7 +80,10 @@ typedef double double_t;
 
 #define MATH_ERRNO 1
 #define MATH_ERREXCEPT 2
-#define math_errhandling (int) (MATH_ERRNO | MATH_ERREXCEPT)
+#define math_errhandling MATH_ERRNO
+
+#define FP_ILOGB0 INT_MIN
+#define FP_ILOGBNAN INT_MAX
 
 double acos(double);
 float acosf(float);

makefile

@@ -56,7 +56,7 @@ string/memset.o \
 string/strcmp.o \
 string/strcspn.o \
 string/strlen.o \
-string/strtok.o \
+string/strtok.o 
 
 HOSTEDOBJS=\
 $(ARCH_HOSTEDOBJS) \
@@ -66,6 +66,8 @@ math/cosl.o \
 math/fabs.o \
 math/fabsf.o \
 math/fabsl.o \
+math/__fpclassify.o \
+math/__signbit.o \
 stdlib/abs.o \
 stdlib/div.o \
 stdlib/labs.o \

math/__fpclassify.c

@@ -0,0 +1,91 @@
+#include <math.h>
+
+int __fpclassifyd(double x) {
+  union fp {
+    double d;
+    long i;
+  } f;
+
+  f.d = x;
+
+  if((f.i & 0x7fffffffffffffff) == 0) {
+    return FP_ZERO;
+  }
+
+  if((f.i & 0x7ff0000000000000) == 0 && (f.i & 0xfffffffffffff != 0)) {
+    return FP_SUBNORMAL;
+  }
+
+  if((f.i & 0x7ff0000000000000) == 0x7ff0000000000000) {
+    if((f.i & 0xfffffffffffff) == 0) {
+      return FP_INFINITE;
+    }
+    else {
+      return FP_NAN;
+    }
+  }
+
+  return FP_NORMAL;
+}
+
+int __fpclassifyf(float x) {
+  union fp {
+    float f;
+    int i;
+  } f;
+
+  f.f = x;
+
+  if((f.i & 0x7fffffff) == 0) {
+    return FP_ZERO;
+  }
+
+  if((f.i & 0x7f800000) == 0 && (f.i & 0x7fffff != 0)) {
+    return FP_SUBNORMAL;
+  }
+
+  if((f.i & 0x7f800000) == 0x7f800000) {
+    if((f.i & 0x7fffff) == 0) {
+      return FP_INFINITE;
+    }
+    else {
+      return FP_NAN;
+    }
+  }
+
+  return FP_NORMAL;
+}
+
+/*
+  This is based off IEEE 754 binary128 (A.K.A quad-precision) FP numbers. 
+  Apparently, C tends to actually use 80-bit extended precision, which would
+  mean this doesn't actually work. We'll worry about that later, though.
+  -Kat
+*/
+int __fpclassifyl(long double x) {
+  union fp {
+    long double d;
+    long long i;
+  } f;
+
+  f.d = x;
+
+  if((f.i & 0x7fffffffffffffffffffffffffffffff) == 0) {
+    return FP_ZERO;
+  }
+
+  if((f.i & 0x7fff0000000000000000000000000000) == 0 && (f.i & 0xffffffffffffffffffffffffffff != 0)) {
+    return FP_SUBNORMAL;
+  }
+
+  if((f.i & 0x7fff0000000000000000000000000000) == 0x7fff0000000000000000000000000000) {
+    if((f.i & 0xffffffffffffffffffffffffffff) == 0) {
+      return FP_INFINITE;
+    }
+    else {
+      return FP_NAN;
+    }
+  }
+
+  return FP_NORMAL;
+}

math/__signbit.c

@@ -0,0 +1,34 @@
+#include <math.h>
+
+int __signbitd(double x) {
+  union fp {
+    double d;
+    long i;
+  } f;
+
+  f.d = x;
+
+  return f.i & 0x8000000000000000;
+}
+
+int __signbitf(float x) {
+  union fp {
+    float f;
+    int i;
+  } f;
+
+  f.f = x;
+
+  return f.i & 0x80000000;
+}
+
+int __signbitl(long double x) {
+  union fp {
+    long double d;
+    long long i;
+  } f;
+
+  f.d = x;
+
+  return f.i & 0x80000000000000000000000000000000;
+}