From ef5fd0562f0b8f9736d8f63776246b6c7a7bec60 Mon Sep 17 00:00:00 2001 From: Kat Richey Date: Sat, 28 May 2022 11:15:26 -0500 Subject: [PATCH] Some rough floating point macros! --- include/math.h | 47 ++++++++++++++++++++++- makefile | 4 +- math/__fpclassify.c | 91 +++++++++++++++++++++++++++++++++++++++++++++ math/__signbit.c | 34 +++++++++++++++++ 4 files changed, 173 insertions(+), 3 deletions(-) create mode 100644 math/__fpclassify.c create mode 100644 math/__signbit.c diff --git a/include/math.h b/include/math.h index 3cec274..a952a34 100755 --- a/include/math.h +++ b/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 + #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); diff --git a/makefile b/makefile index 8c1be70..7b64263 100755 --- a/makefile +++ b/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 \ diff --git a/math/__fpclassify.c b/math/__fpclassify.c new file mode 100644 index 0000000..6a9ecec --- /dev/null +++ b/math/__fpclassify.c @@ -0,0 +1,91 @@ +#include + +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; +} \ No newline at end of file diff --git a/math/__signbit.c b/math/__signbit.c new file mode 100644 index 0000000..d990c0b --- /dev/null +++ b/math/__signbit.c @@ -0,0 +1,34 @@ +#include + +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; +} \ No newline at end of file