From f533aa779c2bf6ee7b677935e4cadca1b2af5c16 Mon Sep 17 00:00:00 2001 From: bumbread Date: Sun, 12 Jun 2022 09:34:36 +1100 Subject: [PATCH] Float division functions --- bake.cmd | 2 +- code/math/_util.h | 10 ++ code/math/basic.c | 112 +++++++++++++++++++ code/math/division.c | 250 +++++++++++++++++++++++++++++++++++++++++++ code/math/ieee754.c | 6 -- code/math/round.c | 8 +- 6 files changed, 374 insertions(+), 14 deletions(-) create mode 100644 code/math/_util.h create mode 100644 code/math/basic.c create mode 100644 code/math/division.c diff --git a/bake.cmd b/bake.cmd index fd28549..a243b17 100644 --- a/bake.cmd +++ b/bake.cmd @@ -13,7 +13,7 @@ if "%~1" neq "_start_" ( ) shift /1 -set CIABATTA_OPTIONS=-Iinc -Wall -g -gcodeview -nodefaultlibs -D_CRT_SECURE_NO_WARNINGS +set CIABATTA_OPTIONS=-Iinc -Wall -g -gcodeview -nodefaultlibs -D_CRT_SECURE_NO_WARNINGS -mfma set PLATFORM=win if "%1"=="test" ( diff --git a/code/math/_util.h b/code/math/_util.h new file mode 100644 index 0000000..e640e7a --- /dev/null +++ b/code/math/_util.h @@ -0,0 +1,10 @@ + +#define asuint64(x) ((union {double f; uint64_t i;}){x}).i +#define asdouble(x) ((union {double f; uint64_t i;}){x}).f + +#include <_compiler.h> +#if defined(_compiler_clang) || defined(_compiler_gnu) + #define just_do_it(t) __attribute__((unused)) volatile t +#else + #define just_do_it(t) volatile t +#endif diff --git a/code/math/basic.c b/code/math/basic.c new file mode 100644 index 0000000..27e9d68 --- /dev/null +++ b/code/math/basic.c @@ -0,0 +1,112 @@ + +#include +#include + +#if !defined(__FMA__) + #error "Get a better CPU (the kind that supports FMA) or enable -mfma" +#endif + +#include + +double fabs(double x) { + union {double f; uint64_t i;} u = {x}; + u.i &= -1ULL/2; + return u.f; +} + +float fabsf(float x) { + union {float f; uint32_t i;} u = {x}; + u.i &= 0x7fffffff; + return u.f; +} + +long double fabsl(long double x) { + return fabs(x); +} + +double fdim(double x, double y) { + if (isnan(x)) return x; + if (isnan(y)) return y; + return x > y ? x - y : 0; +} + +float fdimf(float x, float y) { + if (isnan(x)) return x; + if (isnan(y)) return y; + return x > y ? x - y : 0; +} + +long double fdiml(long double x, long double y) { + return fdim(x, y); +} + +double fmax(double x, double y) { + if (isnan(x)) return y; + if (isnan(y)) return x; + /* handle signed zeros, see C99 Annex F.9.9.2 */ + if (signbit(x) != signbit(y)) + return signbit(x) ? y : x; + return x < y ? y : x; +} + +float fmaxf(float x, float y) { + if (isnan(x)) return y; + if (isnan(y)) return x; + /* handle signed zeroes, see C99 Annex F.9.9.2 */ + if (signbit(x) != signbit(y)) + return signbit(x) ? y : x; + return x < y ? y : x; +} + +long double fmaxl(long double x, long double y) { + return fmax(x, y); +} + +double fmin(double x, double y) { + if (isnan(x)) return y; + if (isnan(y)) return x; + /* handle signed zeros, see C99 Annex F.9.9.2 */ + if (signbit(x) != signbit(y)) + return signbit(x) ? x : y; + return x < y ? x : y; +} + +float fminf(float x, float y) { + if (isnan(x)) return y; + if (isnan(y)) return x; + /* handle signed zeros, see C99 Annex F.9.9.2 */ + if (signbit(x) != signbit(y)) + return signbit(x) ? x : y; + return x < y ? x : y; +} + +long double fminl(long double x, long double y) { + return fmin(x, y); +} + +double fma(double x, double y, double z) { + __m128d xd = _mm_set_sd(x); + __m128d yd = _mm_set_sd(y); + __m128d zd = _mm_set_sd(z); + __m128d rd = _mm_fmadd_sd(xd, yd, zd); + double res = _mm_cvtsd_f64(rd); + return res; +} + +float fmaf(float x, float y, float z) { + __m128 xd = _mm_set_ss(x); + __m128 yd = _mm_set_ss(y); + __m128 zd = _mm_set_ss(z); + __m128 rd = _mm_fmadd_sd(xd, yd, zd); + float res = _mm_cvtsd_f64(rd); + return res; +} + +long double fmal(long double x, long double y, long double z) { + return fma(x, y, z); +} + +float remainderf(float x, float y) { + int q; + return remquof(x, y, &q); +} diff --git a/code/math/division.c b/code/math/division.c new file mode 100644 index 0000000..5e173ef --- /dev/null +++ b/code/math/division.c @@ -0,0 +1,250 @@ + +#include +#include + +double remquo(double x, double y, int *quo) { + union {double f; uint64_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>52 & 0x7ff; + int ey = uy.i>>52 & 0x7ff; + int sx = ux.i>>63; + int sy = uy.i>>63; + uint32_t q; + uint64_t i; + uint64_t uxi = ux.i; + + *quo = 0; + if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff) + return (x*y)/(x*y); + if (ux.i<<1 == 0) + return x; + + /* normalize x and y */ + if (!ex) { + for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1ULL >> 12; + uxi |= 1ULL << 52; + } + if (!ey) { + for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1ULL >> 12; + uy.i |= 1ULL << 52; + } + + q = 0; + if (ex < ey) { + if (ex+1 == ey) + goto end; + return x; + } + + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 63 == 0) { + uxi = i; + q++; + } + uxi <<= 1; + q <<= 1; + } + i = uxi - uy.i; + if (i >> 63 == 0) { + uxi = i; + q++; + } + if (uxi == 0) + ex = -60; + else + for (; uxi>>52 == 0; uxi <<= 1, ex--); +end: + /* scale result and decide between |x| and |x|-|y| */ + if (ex > 0) { + uxi -= 1ULL << 52; + uxi |= (uint64_t)ex << 52; + } else { + uxi >>= -ex + 1; + } + ux.i = uxi; + x = ux.f; + if (sy) + y = -y; + if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) { + x -= y; + q++; + } + q &= 0x7fffffff; + *quo = sx^sy ? -(int)q : (int)q; + return sx ? -x : x; +} + +float remquof(float x, float y, int *quo) { + union {float f; uint32_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>23 & 0xff; + int ey = uy.i>>23 & 0xff; + int sx = ux.i>>31; + int sy = uy.i>>31; + uint32_t q; + uint32_t i; + uint32_t uxi = ux.i; + + *quo = 0; + if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) + return (x*y)/(x*y); + if (ux.i<<1 == 0) + return x; + + /* normalize x and y */ + if (!ex) { + for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1U >> 9; + uxi |= 1U << 23; + } + if (!ey) { + for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1U >> 9; + uy.i |= 1U << 23; + } + + q = 0; + if (ex < ey) { + if (ex+1 == ey) + goto end; + return x; + } + + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 31 == 0) { + uxi = i; + q++; + } + uxi <<= 1; + q <<= 1; + } + i = uxi - uy.i; + if (i >> 31 == 0) { + uxi = i; + q++; + } + if (uxi == 0) + ex = -30; + else + for (; uxi>>23 == 0; uxi <<= 1, ex--); +end: + /* scale result and decide between |x| and |x|-|y| */ + if (ex > 0) { + uxi -= 1U << 23; + uxi |= (uint32_t)ex << 23; + } else { + uxi >>= -ex + 1; + } + ux.i = uxi; + x = ux.f; + if (sy) + y = -y; + if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) { + x -= y; + q++; + } + q &= 0x7fffffff; + *quo = sx^sy ? -(int)q : (int)q; + return sx ? -x : x; +} + +long double remquol(long double x, long double y, int *quo) { + return remquo(x, y, quo); +} + +double remainder(double x, double y) { + int q; + return remquo(x, y, &q); +} + +float remainderf(float x, float y) { + int q; + return remquof(x, y, &q); +} + +long double remainderl(long double x, long double y) { + return remainder(x, y); +} + +double modf(double x, double *iptr) { + union {double f; uint64_t i;} u = {x}; + uint64_t mask; + int e = (int)(u.i>>52 & 0x7ff) - 0x3ff; + + /* no fractional part */ + if (e >= 52) { + *iptr = x; + if (e == 0x400 && u.i<<12 != 0) /* nan */ + return x; + u.i &= 1ULL<<63; + return u.f; + } + + /* no integral part*/ + if (e < 0) { + u.i &= 1ULL<<63; + *iptr = u.f; + return x; + } + + mask = -1ULL>>12>>e; + if ((u.i & mask) == 0) { + *iptr = x; + u.i &= 1ULL<<63; + return u.f; + } + u.i &= ~mask; + *iptr = u.f; + return x - u.f; +} + +float modff(float x, float *iptr) { + union {float f; uint32_t i;} u = {x}; + uint32_t mask; + int e = (int)(u.i>>23 & 0xff) - 0x7f; + + /* no fractional part */ + if (e >= 23) { + *iptr = x; + if (e == 0x80 && u.i<<9 != 0) { /* nan */ + return x; + } + u.i &= 0x80000000; + return u.f; + } + /* no integral part */ + if (e < 0) { + u.i &= 0x80000000; + *iptr = u.f; + return x; + } + + mask = 0x007fffff>>e; + if ((u.i & mask) == 0) { + *iptr = x; + u.i &= 0x80000000; + return u.f; + } + u.i &= ~mask; + *iptr = u.f; + return x - u.f; +} + +long double modfl(long double x, long double *iptr) { + double d; + long double r = modf(x, &d); + *iptr = d; + return r; +} diff --git a/code/math/ieee754.c b/code/math/ieee754.c index 6727503..bc93923 100644 --- a/code/math/ieee754.c +++ b/code/math/ieee754.c @@ -4,12 +4,6 @@ #include #include -#include <_compiler.h> -#if defined(_compiler_clang) || defined(_compiler_gnu) - #define just_do_it(t) __attribute__((unused)) volatile t -#endif - - int _fpclassify(double x) { union {double f; uint64_t i;} u = {x}; int e = u.i>>52 & 0x7ff; diff --git a/code/math/round.c b/code/math/round.c index 4d1fe6b..4b7f1f5 100644 --- a/code/math/round.c +++ b/code/math/round.c @@ -5,13 +5,7 @@ #include #include -#include <_compiler.h> -#if defined(_compiler_clang) || defined(_compiler_gnu) - #define just_do_it(t) __attribute__((unused)) volatile t -#endif - -#define asuint64(x) ((union {double f; uint64_t i;}){x}).i -#define asdouble(x) ((union {double f; uint64_t i;}){x}).f +#include "_util.h" double nearbyint(double x) { #pragma STDC FENV_ACCESS ON