From 470bce6c6d243009076c647b52d860c630c759ec Mon Sep 17 00:00:00 2001 From: bumbread Date: Thu, 9 Jun 2022 20:04:39 +1100 Subject: [PATCH] Some ieee754 funcs of math --- code/math.c | 42 ------ code/math/basic.h | 25 ---- code/math/{exponential.h => exp.h} | 0 code/math/floats.h | 45 ------ code/math/ieee754.h | 134 ++++++++++++++++++ code/math/math.c | 50 +++++++ code/math/{hyperbolic.h => near.h} | 0 code/math/{nearest-int.h => pow.h} | 0 code/math/{power.h => stat.h} | 0 .../math/{probability-statistics.h => trig.h} | 0 code/math/{trigonometry.h => trigh.h} | 0 code/stdlib.c | 36 +++++ inc/math.h | 30 +++- inc/stdio.h | 4 + inc/stdlib.h | 24 ++-- 15 files changed, 261 insertions(+), 129 deletions(-) delete mode 100644 code/math.c rename code/math/{exponential.h => exp.h} (100%) delete mode 100644 code/math/floats.h create mode 100644 code/math/ieee754.h create mode 100644 code/math/math.c rename code/math/{hyperbolic.h => near.h} (100%) rename code/math/{nearest-int.h => pow.h} (100%) rename code/math/{power.h => stat.h} (100%) rename code/math/{probability-statistics.h => trig.h} (100%) rename code/math/{trigonometry.h => trigh.h} (100%) diff --git a/code/math.c b/code/math.c deleted file mode 100644 index 05de3f2..0000000 --- a/code/math.c +++ /dev/null @@ -1,42 +0,0 @@ - -#include - -#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 diff --git a/code/math/basic.h b/code/math/basic.h index 73b7390..e69de29 100644 --- a/code/math/basic.h +++ b/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; -} - - diff --git a/code/math/exponential.h b/code/math/exp.h similarity index 100% rename from code/math/exponential.h rename to code/math/exp.h diff --git a/code/math/floats.h b/code/math/floats.h deleted file mode 100644 index 79e2e72..0000000 --- a/code/math/floats.h +++ /dev/null @@ -1,45 +0,0 @@ - -#include - -#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; - } -} - diff --git a/code/math/ieee754.h b/code/math/ieee754.h new file mode 100644 index 0000000..a013427 --- /dev/null +++ b/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; +} \ No newline at end of file diff --git a/code/math/math.c b/code/math/math.c new file mode 100644 index 0000000..fa1fbe2 --- /dev/null +++ b/code/math/math.c @@ -0,0 +1,50 @@ + +#include +#include + +// 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 diff --git a/code/math/hyperbolic.h b/code/math/near.h similarity index 100% rename from code/math/hyperbolic.h rename to code/math/near.h diff --git a/code/math/nearest-int.h b/code/math/pow.h similarity index 100% rename from code/math/nearest-int.h rename to code/math/pow.h diff --git a/code/math/power.h b/code/math/stat.h similarity index 100% rename from code/math/power.h rename to code/math/stat.h diff --git a/code/math/probability-statistics.h b/code/math/trig.h similarity index 100% rename from code/math/probability-statistics.h rename to code/math/trig.h diff --git a/code/math/trigonometry.h b/code/math/trigh.h similarity index 100% rename from code/math/trigonometry.h rename to code/math/trigh.h diff --git a/code/stdlib.c b/code/stdlib.c index 39ff5b9..c853f49 100644 --- a/code/stdlib.c +++ b/code/stdlib.c @@ -1,5 +1,41 @@ #include +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; diff --git a/inc/math.h b/inc/math.h index 0a954f4..92da1b2 100644 --- a/inc/math.h +++ b/inc/math.h @@ -17,21 +17,41 @@ typedef double double_t; // FP_ILOGB0 // FP_ILOGBNAN -#define MATH_ERRNO 0 -#define MATH_ERREXCEPT 1 +#define MATH_ERRNO 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); -float acosf(float x); +// Ordering +#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); diff --git a/inc/stdio.h b/inc/stdio.h index 1dffee6..b5e6289 100644 --- a/inc/stdio.h +++ b/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 diff --git a/inc/stdlib.h b/inc/stdlib.h index 672cfbd..0e587a9 100644 --- a/inc/stdlib.h +++ b/inc/stdlib.h @@ -15,20 +15,20 @@ #define EXIT_SUCCESS 0 #define EXIT_FAILURE 1 -// typedef struct div_t { -// int quot; -// int rem; -// } div_t; +typedef struct div_t { + int quot; + int rem; +} div_t; -// typedef struct ldiv_t { -// long quot; -// long rem; -// } ldiv_t; +typedef struct ldiv_t { + long quot; + long rem; +} ldiv_t; -// typedef struct lldiv_t { -// long long quot; -// long long rem; -// } lldiv_t; +typedef struct lldiv_t { + long long quot; + long long rem; +} lldiv_t; // #define EXIT_FAILURE 1 // #define EXIT_SUCCESS 0