summaryrefslogtreecommitdiff
path: root/lib/libm/fmodl.c
diff options
context:
space:
mode:
authorKacper <kacper@mail.openlinux.dev>2025-12-09 19:20:15 +0100
committerKacper <kacper@mail.openlinux.dev>2025-12-09 19:20:15 +0100
commit885f5974cdf65b59415837ae97f5a14ef1350670 (patch)
tree66ac13de29c7f4932c5fcae11773df574e4e256a /lib/libm/fmodl.c
parent8f9e448b2ef6db7cd905540c21f3c5b190e7a1e7 (diff)
feat: add gzip and new headers
Diffstat (limited to 'lib/libm/fmodl.c')
-rw-r--r--lib/libm/fmodl.c107
1 files changed, 107 insertions, 0 deletions
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