From 885f5974cdf65b59415837ae97f5a14ef1350670 Mon Sep 17 00:00:00 2001 From: Kacper Date: Tue, 9 Dec 2025 19:20:15 +0100 Subject: feat: add gzip and new headers --- lib/libm/fmodl.c | 107 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 107 insertions(+) create mode 100644 lib/libm/fmodl.c (limited to 'lib/libm/fmodl.c') diff --git a/lib/libm/fmodl.c b/lib/libm/fmodl.c new file mode 100644 index 00000000..72111dc1 --- /dev/null +++ b/lib/libm/fmodl.c @@ -0,0 +1,107 @@ +#include "libm.h" + +#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 +long double fmodl(long double x, long double y) +{ + return fmod(x, y); +} +#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +long double fmodl(long double x, long double y) +{ + union ldshape ux = { x }, uy = { y }; + int ex = ux.i.se & 0x7fff; + int ey = uy.i.se & 0x7fff; + int sx = ux.i.se & 0x8000; + + if (y == 0 || isnan(y) || ex == 0x7fff) + return (x * y) / (x * y); + ux.i.se = ex; + uy.i.se = ey; + if (ux.f <= uy.f) { + if (ux.f == uy.f) + return 0 * x; + return x; + } + + /* normalize x and y */ + if (!ex) { + ux.f *= 0x1p120f; + ex = ux.i.se - 120; + } + if (!ey) { + uy.f *= 0x1p120f; + ey = uy.i.se - 120; + } + + /* x mod y */ +#if LDBL_MANT_DIG == 64 + uint64_t i, mx, my; + mx = ux.i.m; + my = uy.i.m; + for (; ex > ey; ex--) { + i = mx - my; + if (mx >= my) { + if (i == 0) + return 0 * x; + mx = 2 * i; + } else if (2 * mx < mx) { + mx = 2 * mx - my; + } else { + mx = 2 * mx; + } + } + i = mx - my; + if (mx >= my) { + if (i == 0) + return 0 * x; + mx = i; + } + for (; mx >> 63 == 0; mx *= 2, ex--) + ; + ux.i.m = mx; +#elif LDBL_MANT_DIG == 113 + uint64_t hi, lo, xhi, xlo, yhi, ylo; + xhi = (ux.i2.hi & -1ULL >> 16) | 1ULL << 48; + yhi = (uy.i2.hi & -1ULL >> 16) | 1ULL << 48; + xlo = ux.i2.lo; + ylo = uy.i2.lo; + for (; ex > ey; ex--) { + hi = xhi - yhi; + lo = xlo - ylo; + if (xlo < ylo) + hi -= 1; + if (hi >> 63 == 0) { + if ((hi | lo) == 0) + return 0 * x; + xhi = 2 * hi + (lo >> 63); + xlo = 2 * lo; + } else { + xhi = 2 * xhi + (xlo >> 63); + xlo = 2 * xlo; + } + } + hi = xhi - yhi; + lo = xlo - ylo; + if (xlo < ylo) + hi -= 1; + if (hi >> 63 == 0) { + if ((hi | lo) == 0) + return 0 * x; + xhi = hi; + xlo = lo; + } + for (; xhi >> 48 == 0; xhi = 2 * xhi + (xlo >> 63), xlo = 2 * xlo, ex--) + ; + ux.i2.hi = xhi; + ux.i2.lo = xlo; +#endif + + /* scale result */ + if (ex <= 0) { + ux.i.se = (ex + 120) | sx; + ux.f *= 0x1p-120f; + } else + ux.i.se = ex | sx; + return ux.f; +} +#endif -- cgit v1.2.3