diff options
author | Ian Moffett <ian@osmora.org> | 2024-03-07 17:28:00 -0500 |
---|---|---|
committer | Ian Moffett <ian@osmora.org> | 2024-03-07 17:28:32 -0500 |
commit | bd5969fc876a10b18613302db7087ef3c40f18e1 (patch) | |
tree | 7c2b8619afe902abf99570df2873fbdf40a4d1a1 /lib/mlibc/options/ansi/musl-generic-math/remquo.c | |
parent | a95b38b1b92b172e6cc4e8e56a88a30cc65907b0 (diff) |
lib: Add mlibc
Signed-off-by: Ian Moffett <ian@osmora.org>
Diffstat (limited to 'lib/mlibc/options/ansi/musl-generic-math/remquo.c')
-rw-r--r-- | lib/mlibc/options/ansi/musl-generic-math/remquo.c | 82 |
1 files changed, 82 insertions, 0 deletions
diff --git a/lib/mlibc/options/ansi/musl-generic-math/remquo.c b/lib/mlibc/options/ansi/musl-generic-math/remquo.c new file mode 100644 index 0000000..59d5ad5 --- /dev/null +++ b/lib/mlibc/options/ansi/musl-generic-math/remquo.c @@ -0,0 +1,82 @@ +#include <math.h> +#include <stdint.h> + +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; +} |