From bd5969fc876a10b18613302db7087ef3c40f18e1 Mon Sep 17 00:00:00 2001 From: Ian Moffett Date: Thu, 7 Mar 2024 17:28:00 -0500 Subject: lib: Add mlibc Signed-off-by: Ian Moffett --- lib/mlibc/options/ansi/musl-generic-math/hypot.c | 67 ++++++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 lib/mlibc/options/ansi/musl-generic-math/hypot.c (limited to 'lib/mlibc/options/ansi/musl-generic-math/hypot.c') diff --git a/lib/mlibc/options/ansi/musl-generic-math/hypot.c b/lib/mlibc/options/ansi/musl-generic-math/hypot.c new file mode 100644 index 0000000..6071bf1 --- /dev/null +++ b/lib/mlibc/options/ansi/musl-generic-math/hypot.c @@ -0,0 +1,67 @@ +#include +#include +#include + +#if FLT_EVAL_METHOD > 1U && LDBL_MANT_DIG == 64 +#define SPLIT (0x1p32 + 1) +#else +#define SPLIT (0x1p27 + 1) +#endif + +static void sq(double_t *hi, double_t *lo, double x) +{ + double_t xh, xl, xc; + + xc = (double_t)x*SPLIT; + xh = x - xc + xc; + xl = x - xh; + *hi = (double_t)x*x; + *lo = xh*xh - *hi + 2*xh*xl + xl*xl; +} + +double hypot(double x, double y) +{ + union {double f; uint64_t i;} ux = {x}, uy = {y}, ut; + int ex, ey; + double_t hx, lx, hy, ly, z; + + /* arrange |x| >= |y| */ + ux.i &= -1ULL>>1; + uy.i &= -1ULL>>1; + if (ux.i < uy.i) { + ut = ux; + ux = uy; + uy = ut; + } + + /* special cases */ + ex = ux.i>>52; + ey = uy.i>>52; + x = ux.f; + y = uy.f; + /* note: hypot(inf,nan) == inf */ + if (ey == 0x7ff) + return y; + if (ex == 0x7ff || uy.i == 0) + return x; + /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */ + /* 64 difference is enough for ld80 double_t */ + if (ex - ey > 64) + return x + y; + + /* precise sqrt argument in nearest rounding mode without overflow */ + /* xh*xh must not overflow and xl*xl must not underflow in sq */ + z = 1; + if (ex > 0x3ff+510) { + z = 0x1p700; + x *= 0x1p-700; + y *= 0x1p-700; + } else if (ey < 0x3ff-450) { + z = 0x1p-700; + x *= 0x1p700; + y *= 0x1p700; + } + sq(&hx, &lx, x); + sq(&hy, &ly, y); + return z*sqrt(ly+lx+hy+hx); +} -- cgit v1.2.3