From 888014c83724008194b4234d96dc3af3c86dc38a Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Wed, 11 Mar 2026 17:58:42 +0100
Subject: [PATCH 01/10] `libzigc`: `std.math` -> `math`

Remove the two uses of `std.math` since that namespace is
already assigned to the `math` constant and used as such.
---
 lib/c/math.zig | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 501df810d6..15340ef622 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -70,7 +70,7 @@ fn acos(x: f64) callconv(.c) f64 {
 }
 
 fn acosf(x: f32) callconv(.c) f32 {
-    return std.math.acos(x);
+    return math.acos(x);
 }
 
 fn atan(x: f64) callconv(.c) f64 {
@@ -177,7 +177,7 @@ fn pow10f(x: f32) callconv(.c) f32 {
 }
 
 fn rint(x: f64) callconv(.c) f64 {
-    const toint: f64 = 1.0 / @as(f64, std.math.floatEps(f64));
+    const toint: f64 = 1.0 / @as(f64, math.floatEps(f64));
     const a: u64 = @bitCast(x);
     const e = a >> 52 & 0x7ff;
     const s = a >> 63;

From 69cbe8ea3c54b181a42612973ae25aada7bd8b6e Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Wed, 11 Mar 2026 18:51:39 +0100
Subject: [PATCH 02/10] `libzigc`: Implement `asin`

The changes were tested by running:

```
$ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib

$ stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=asin -fqemu -fwasmtime --summary line
Build Summary: 921/921 steps succeeded
```
---
 lib/c/math.zig                     |   5 ++
 lib/libc/musl/src/math/asin.c      | 107 -----------------------------
 lib/libc/musl/src/math/i386/asin.s |  21 ------
 src/libs/musl.zig                  |   2 -
 src/libs/wasi_libc.zig             |   1 -
 5 files changed, 5 insertions(+), 131 deletions(-)
 delete mode 100644 lib/libc/musl/src/math/asin.c
 delete mode 100644 lib/libc/musl/src/math/i386/asin.s

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 15340ef622..0ee97ea674 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -42,6 +42,7 @@ comptime {
     if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
         symbol(&acos, "acos");
         symbol(&acosf, "acosf");
+        symbol(&asin, "asin");
         symbol(&atan, "atan");
         symbol(&atanf, "atanf");
         symbol(&atanl, "atanl");
@@ -73,6 +74,10 @@ fn acosf(x: f32) callconv(.c) f32 {
     return math.acos(x);
 }
 
+fn asin(x: f64) callconv(.c) f64 {
+    return math.asin(x);
+}
+
 fn atan(x: f64) callconv(.c) f64 {
     return math.atan(x);
 }
diff --git a/lib/libc/musl/src/math/asin.c b/lib/libc/musl/src/math/asin.c
deleted file mode 100644
index c926b18855..0000000000
--- a/lib/libc/musl/src/math/asin.c
+++ /dev/null
@@ -1,107 +0,0 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/e_asin.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* asin(x)
- * Method :
- *      Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
- *      we approximate asin(x) on [0,0.5] by
- *              asin(x) = x + x*x^2*R(x^2)
- *      where
- *              R(x^2) is a rational approximation of (asin(x)-x)/x^3
- *      and its remez error is bounded by
- *              |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
- *
- *      For x in [0.5,1]
- *              asin(x) = pi/2-2*asin(sqrt((1-x)/2))
- *      Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
- *      then for x>0.98
- *              asin(x) = pi/2 - 2*(s+s*z*R(z))
- *                      = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
- *      For x<=0.98, let pio4_hi = pio2_hi/2, then
- *              f = hi part of s;
- *              c = sqrt(z) - f = (z-f*f)/(s+f)         ...f+c=sqrt(z)
- *      and
- *              asin(x) = pi/2 - 2*(s+s*z*R(z))
- *                      = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
- *                      = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
- *
- * Special cases:
- *      if x is NaN, return x itself;
- *      if |x|>1, return NaN with invalid signal.
- *
- */
-
-#include "libm.h"
-
-static const double
-pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
-pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
-/* coefficients for R(x^2) */
-pS0 =  1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
-pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
-pS2 =  2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
-pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
-pS4 =  7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
-pS5 =  3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
-qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
-qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
-qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
-qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
-
-static double R(double z)
-{
-	double_t p, q;
-	p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
-	q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
-	return p/q;
-}
-
-double asin(double x)
-{
-	double z,r,s;
-	uint32_t hx,ix;
-
-	GET_HIGH_WORD(hx, x);
-	ix = hx & 0x7fffffff;
-	/* |x| >= 1 or nan */
-	if (ix >= 0x3ff00000) {
-		uint32_t lx;
-		GET_LOW_WORD(lx, x);
-		if ((ix-0x3ff00000 | lx) == 0)
-			/* asin(1) = +-pi/2 with inexact */
-			return x*pio2_hi + 0x1p-120f;
-		return 0/(x-x);
-	}
-	/* |x| < 0.5 */
-	if (ix < 0x3fe00000) {
-		/* if 0x1p-1022 <= |x| < 0x1p-26, avoid raising underflow */
-		if (ix < 0x3e500000 && ix >= 0x00100000)
-			return x;
-		return x + x*R(x*x);
-	}
-	/* 1 > |x| >= 0.5 */
-	z = (1 - fabs(x))*0.5;
-	s = sqrt(z);
-	r = R(z);
-	if (ix >= 0x3fef3333) {  /* if |x| > 0.975 */
-		x = pio2_hi-(2*(s+s*r)-pio2_lo);
-	} else {
-		double f,c;
-		/* f+c = sqrt(z) */
-		f = s;
-		SET_LOW_WORD(f,0);
-		c = (z-f*f)/(s+f);
-		x = 0.5*pio2_hi - (2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
-	}
-	if (hx >> 31)
-		return -x;
-	return x;
-}
diff --git a/lib/libc/musl/src/math/i386/asin.s b/lib/libc/musl/src/math/i386/asin.s
deleted file mode 100644
index 2bc8356f49..0000000000
--- a/lib/libc/musl/src/math/i386/asin.s
+++ /dev/null
@@ -1,21 +0,0 @@
-.global asin
-.type asin,@function
-asin:
-	fldl 4(%esp)
-	mov 8(%esp),%eax
-	add %eax,%eax
-	cmp $0x00200000,%eax
-	jb 1f
-	fld %st(0)
-	fld1
-	fsub %st(0),%st(1)
-	fadd %st(2)
-	fmulp
-	fsqrt
-	fpatan
-	fstpl 4(%esp)
-	fldl 4(%esp)
-	ret
-		# subnormal x, return x with underflow
-1:	fsts 4(%esp)
-	ret
diff --git a/src/libs/musl.zig b/src/libs/musl.zig
index d3d471459a..31391b0736 100644
--- a/src/libs/musl.zig
+++ b/src/libs/musl.zig
@@ -801,7 +801,6 @@ const src_files = [_][]const u8{
     "musl/src/math/acosl.c",
     "musl/src/math/arm/fma.c",
     "musl/src/math/arm/fmaf.c",
-    "musl/src/math/asin.c",
     "musl/src/math/asinf.c",
     "musl/src/math/asinh.c",
     "musl/src/math/asinhf.c",
@@ -849,7 +848,6 @@ const src_files = [_][]const u8{
     "musl/src/math/i386/acosl.s",
     "musl/src/math/i386/asinf.s",
     "musl/src/math/i386/asinl.s",
-    "musl/src/math/i386/asin.s",
     "musl/src/math/i386/atan2f.s",
     "musl/src/math/i386/atan2l.s",
     "musl/src/math/i386/atan2.s",
diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig
index 67c1b0e74c..5afb9e5382 100644
--- a/src/libs/wasi_libc.zig
+++ b/src/libs/wasi_libc.zig
@@ -668,7 +668,6 @@ const libc_top_half_src_files = [_][]const u8{
     "musl/src/math/acoshf.c",
     "musl/src/math/acoshl.c",
     "musl/src/math/acosl.c",
-    "musl/src/math/asin.c",
     "musl/src/math/asinf.c",
     "musl/src/math/asinh.c",
     "musl/src/math/asinhf.c",

From c77103add76a320d1ce877216e26dddd390caadc Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Sat, 14 Mar 2026 10:55:04 +0100
Subject: [PATCH 03/10] `libzigc`: implement `acoshf`

The functions were tested by running:

```
$ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib

$ stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=acoshf -fqemu -fwasmtime --summary line
Build Summary: 369/369 steps succeeded
```
---
 lib/c/math.zig                  |  5 +++++
 lib/libc/musl/src/math/acoshf.c | 26 --------------------------
 src/libs/musl.zig               |  1 -
 src/libs/wasi_libc.zig          |  1 -
 4 files changed, 5 insertions(+), 28 deletions(-)
 delete mode 100644 lib/libc/musl/src/math/acoshf.c

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 0ee97ea674..367da8d48f 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -42,6 +42,7 @@ comptime {
     if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
         symbol(&acos, "acos");
         symbol(&acosf, "acosf");
+        symbol(&acoshf, "acoshf");
         symbol(&asin, "asin");
         symbol(&atan, "atan");
         symbol(&atanf, "atanf");
@@ -74,6 +75,10 @@ fn acosf(x: f32) callconv(.c) f32 {
     return math.acos(x);
 }
 
+fn acoshf(x: f32) callconv(.c) f32 {
+    return math.acosh(x);
+}
+
 fn asin(x: f64) callconv(.c) f64 {
     return math.asin(x);
 }
diff --git a/lib/libc/musl/src/math/acoshf.c b/lib/libc/musl/src/math/acoshf.c
deleted file mode 100644
index b773d48e2b..0000000000
--- a/lib/libc/musl/src/math/acoshf.c
+++ /dev/null
@@ -1,26 +0,0 @@
-#include "libm.h"
-
-#if FLT_EVAL_METHOD==2
-#undef sqrtf
-#define sqrtf sqrtl
-#elif FLT_EVAL_METHOD==1
-#undef sqrtf
-#define sqrtf sqrt
-#endif
-
-/* acosh(x) = log(x + sqrt(x*x-1)) */
-float acoshf(float x)
-{
-	union {float f; uint32_t i;} u = {x};
-	uint32_t a = u.i & 0x7fffffff;
-
-	if (a < 0x3f800000+(1<<23))
-		/* |x| < 2, invalid if x < 1 */
-		/* up to 2ulp error in [1,1.125] */
-		return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1)));
-	if (u.i < 0x3f800000+(12<<23))
-		/* 2 <= x < 0x1p12 */
-		return logf(2*x - 1/(x+sqrtf(x*x-1)));
-	/* x >= 0x1p12 or x <= -2 or nan */
-	return logf(x) + 0.693147180559945309417232121458176568f;
-}
diff --git a/src/libs/musl.zig b/src/libs/musl.zig
index 31391b0736..19ec8fc36e 100644
--- a/src/libs/musl.zig
+++ b/src/libs/musl.zig
@@ -796,7 +796,6 @@ const src_files = [_][]const u8{
     "musl/src/math/aarch64/nearbyintf.c",
     "musl/src/math/aarch64/rintf.c",
     "musl/src/math/acosh.c",
-    "musl/src/math/acoshf.c",
     "musl/src/math/acoshl.c",
     "musl/src/math/acosl.c",
     "musl/src/math/arm/fma.c",
diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig
index 5afb9e5382..be691070e0 100644
--- a/src/libs/wasi_libc.zig
+++ b/src/libs/wasi_libc.zig
@@ -665,7 +665,6 @@ const libc_top_half_src_files = [_][]const u8{
     "musl/src/locale/wcscoll.c",
     "musl/src/locale/wcsxfrm.c",
     "musl/src/math/acosh.c",
-    "musl/src/math/acoshf.c",
     "musl/src/math/acoshl.c",
     "musl/src/math/acosl.c",
     "musl/src/math/asinf.c",

From eab22ab2b88313f86158230ddcdcb22f28483d85 Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Sat, 14 Mar 2026 20:47:49 +0100
Subject: [PATCH 04/10] `libzigc`: Implement `modf`

The behaviour regarding special cases differs between `libc` and Zig's
`stdlib` for `modf`, so the implementation couldn't be a straightforward
calling of `stdlib` function.

Other than the obvious documented differences, I also had problems with
the `INVALID` flag being raised while running `libc-test` suite on riscv
arch through qemu. The solution was to test if the argument is `NaN`,
and then return a quiet `NaN` if so.

Passing tests, that should include all the special cases to be wary of,
were also added.

Test results:
```
$ stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=modf -fqemu -fwasmtime --summary line
Build Summary: 921/921 steps succeeded
```
---
 lib/c/math.zig                | 63 +++++++++++++++++++++++++++++++++++
 lib/libc/musl/src/math/modf.c | 34 -------------------
 src/libs/musl.zig             |  1 -
 src/libs/wasi_libc.zig        |  1 -
 4 files changed, 63 insertions(+), 36 deletions(-)
 delete mode 100644 lib/libc/musl/src/math/modf.c

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 367da8d48f..4c0f7c5978 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -53,6 +53,7 @@ comptime {
         symbol(&exp10, "exp10");
         symbol(&exp10f, "exp10f");
         symbol(&hypot, "hypot");
+        symbol(&modf, "modf");
         symbol(&pow, "pow");
         symbol(&pow10, "pow10");
         symbol(&pow10f, "pow10f");
@@ -162,6 +163,68 @@ fn isnanl(x: c_longdouble) callconv(.c) c_int {
     return if (math.isNan(x)) 1 else 0;
 }
 
+fn modf(x: f64, iptr: *f64) callconv(.c) f64 {
+    if (math.isNegativeInf(x)) {
+        iptr.* = -math.inf(f64);
+        return -0.0;
+    }
+
+    if (math.isPositiveInf(x)) {
+        iptr.* = math.inf(f64);
+        return 0.0;
+    }
+
+    // Avoids raising the INVALID flag on qemu-riscv
+    if (math.isNan(x)) {
+        iptr.* = math.nan(f64);
+        return math.nan(f64);
+    }
+
+    const r = math.modf(x);
+    iptr.* = r.ipart;
+
+    // If the result would be a negative zero, we must be explicit about
+    // returning a negative zero.
+    return if (math.isNegativeZero(x) or (x < 0.0 and x == r.ipart)) -0.0 else r.fpart;
+}
+
+test "modf" {
+    var int: f64 = undefined;
+    const iptr = &int;
+    const eps_val = 1e-6;
+
+    const normal_frac = modf(1234.5678, iptr);
+    try std.testing.expectApproxEqRel(0.5678, normal_frac, eps_val);
+    try std.testing.expectApproxEqRel(1234.0, iptr.*, eps_val);
+
+    // When `x` is a NaN, NaN is returned and `*iptr` is set to NaN
+    const nan_frac = modf(math.nan(f64), iptr);
+    try std.testing.expect(math.isNan(nan_frac));
+    try std.testing.expect(math.isNan(iptr.*));
+
+    // When `x` is positive infinity, +0 is returned and `*iptr` is set to
+    // positive infinity
+    const pos_zero_frac = modf(math.inf(f64), iptr);
+    try std.testing.expect(math.isPositiveZero(pos_zero_frac));
+    try std.testing.expect(math.isPositiveInf(iptr.*));
+
+    // When `x` is negative infinity, -0 is returned and `*iptr` is set to
+    // negative infinity
+    const neg_zero_frac = modf(-math.inf(f64), iptr);
+    try std.testing.expect(math.isNegativeZero(neg_zero_frac));
+    try std.testing.expect(math.isNegativeInf(iptr.*));
+
+    // Return -0 when `x` is a negative integer
+    const nz_frac = modf(-1000.0, iptr);
+    try std.testing.expect(math.isNegativeZero(nz_frac));
+    try std.testing.expectEqual(-1000.0, iptr.*);
+
+    // Return +0 when `x` is a positive integer
+    const pz_frac = modf(1000.0, iptr);
+    try std.testing.expect(math.isPositiveZero(pz_frac));
+    try std.testing.expectEqual(1000.0, iptr.*);
+}
+
 fn nan(_: [*:0]const c_char) callconv(.c) f64 {
     return math.nan(f64);
 }
diff --git a/lib/libc/musl/src/math/modf.c b/lib/libc/musl/src/math/modf.c
deleted file mode 100644
index 1c8a1db90d..0000000000
--- a/lib/libc/musl/src/math/modf.c
+++ /dev/null
@@ -1,34 +0,0 @@
-#include "libm.h"
-
-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;
-}
diff --git a/src/libs/musl.zig b/src/libs/musl.zig
index 19ec8fc36e..3db43a1c76 100644
--- a/src/libs/musl.zig
+++ b/src/libs/musl.zig
@@ -934,7 +934,6 @@ const src_files = [_][]const u8{
     "musl/src/math/__math_uflowf.c",
     "musl/src/math/__math_xflow.c",
     "musl/src/math/__math_xflowf.c",
-    "musl/src/math/modf.c",
     "musl/src/math/modff.c",
     "musl/src/math/modfl.c",
     "musl/src/math/nearbyint.c",
diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig
index be691070e0..5845b4784c 100644
--- a/src/libs/wasi_libc.zig
+++ b/src/libs/wasi_libc.zig
@@ -755,7 +755,6 @@ const libc_top_half_src_files = [_][]const u8{
     "musl/src/math/__math_uflowf.c",
     "musl/src/math/__math_xflow.c",
     "musl/src/math/__math_xflowf.c",
-    "musl/src/math/modf.c",
     "musl/src/math/modff.c",
     "musl/src/math/modfl.c",
     "musl/src/math/nearbyintl.c",

From 1c280032aae151b184e9224e2e0d1414f7436a01 Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Sat, 14 Mar 2026 21:03:40 +0100
Subject: [PATCH 05/10] `libzigc`: Assign `std.testing` functions to consts

---
 lib/c/math.zig | 57 ++++++++++++++++++++++++++------------------------
 1 file changed, 30 insertions(+), 27 deletions(-)

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 4c0f7c5978..2593d5bf86 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -2,6 +2,9 @@ const builtin = @import("builtin");
 
 const std = @import("std");
 const math = std.math;
+const expect = std.testing.expect;
+const expectEqual = std.testing.expectEqual;
+const expectApproxEqRel = std.testing.expectApproxEqRel;
 
 const symbol = @import("../c.zig").symbol;
 
@@ -194,35 +197,35 @@ test "modf" {
     const eps_val = 1e-6;
 
     const normal_frac = modf(1234.5678, iptr);
-    try std.testing.expectApproxEqRel(0.5678, normal_frac, eps_val);
-    try std.testing.expectApproxEqRel(1234.0, iptr.*, eps_val);
+    try expectApproxEqRel(0.5678, normal_frac, eps_val);
+    try expectApproxEqRel(1234.0, iptr.*, eps_val);
 
     // When `x` is a NaN, NaN is returned and `*iptr` is set to NaN
     const nan_frac = modf(math.nan(f64), iptr);
-    try std.testing.expect(math.isNan(nan_frac));
-    try std.testing.expect(math.isNan(iptr.*));
+    try expect(math.isNan(nan_frac));
+    try expect(math.isNan(iptr.*));
 
     // When `x` is positive infinity, +0 is returned and `*iptr` is set to
     // positive infinity
     const pos_zero_frac = modf(math.inf(f64), iptr);
-    try std.testing.expect(math.isPositiveZero(pos_zero_frac));
-    try std.testing.expect(math.isPositiveInf(iptr.*));
+    try expect(math.isPositiveZero(pos_zero_frac));
+    try expect(math.isPositiveInf(iptr.*));
 
     // When `x` is negative infinity, -0 is returned and `*iptr` is set to
     // negative infinity
     const neg_zero_frac = modf(-math.inf(f64), iptr);
-    try std.testing.expect(math.isNegativeZero(neg_zero_frac));
-    try std.testing.expect(math.isNegativeInf(iptr.*));
+    try expect(math.isNegativeZero(neg_zero_frac));
+    try expect(math.isNegativeInf(iptr.*));
 
     // Return -0 when `x` is a negative integer
     const nz_frac = modf(-1000.0, iptr);
-    try std.testing.expect(math.isNegativeZero(nz_frac));
-    try std.testing.expectEqual(-1000.0, iptr.*);
+    try expect(math.isNegativeZero(nz_frac));
+    try expectEqual(-1000.0, iptr.*);
 
     // Return +0 when `x` is a positive integer
     const pz_frac = modf(1000.0, iptr);
-    try std.testing.expect(math.isPositiveZero(pz_frac));
-    try std.testing.expectEqual(1000.0, iptr.*);
+    try expect(math.isPositiveZero(pz_frac));
+    try expectEqual(1000.0, iptr.*);
 }
 
 fn nan(_: [*:0]const c_char) callconv(.c) f64 {
@@ -272,35 +275,35 @@ fn rint(x: f64) callconv(.c) f64 {
 
 test "rint" {
     // Positive numbers round correctly
-    try std.testing.expectEqual(@as(f64, 42.0), rint(42.2));
-    try std.testing.expectEqual(@as(f64, 42.0), rint(41.8));
+    try expectEqual(@as(f64, 42.0), rint(42.2));
+    try expectEqual(@as(f64, 42.0), rint(41.8));
 
     // Negative numbers round correctly
-    try std.testing.expectEqual(@as(f64, -6.0), rint(-5.9));
-    try std.testing.expectEqual(@as(f64, -6.0), rint(-6.1));
+    try expectEqual(@as(f64, -6.0), rint(-5.9));
+    try expectEqual(@as(f64, -6.0), rint(-6.1));
 
     // No rounding needed test
-    try std.testing.expectEqual(@as(f64, 5.0), rint(5.0));
-    try std.testing.expectEqual(@as(f64, -10.0), rint(-10.0));
-    try std.testing.expectEqual(@as(f64, 0.0), rint(0.0));
+    try expectEqual(@as(f64, 5.0), rint(5.0));
+    try expectEqual(@as(f64, -10.0), rint(-10.0));
+    try expectEqual(@as(f64, 0.0), rint(0.0));
 
     // Very large numbers return unchanged
     const large: f64 = 9007199254740992.0; // 2^53
-    try std.testing.expectEqual(large, rint(large));
-    try std.testing.expectEqual(-large, rint(-large));
+    try expectEqual(large, rint(large));
+    try expectEqual(-large, rint(-large));
 
     // Small positive numbers round to zero
     const pos_result = rint(0.3);
-    try std.testing.expectEqual(@as(f64, 0.0), pos_result);
-    try std.testing.expect(@as(u64, @bitCast(pos_result)) == 0);
+    try expectEqual(@as(f64, 0.0), pos_result);
+    try expect(@as(u64, @bitCast(pos_result)) == 0);
 
     // Small negative numbers round to negative zero
     const neg_result = rint(-0.3);
-    try std.testing.expectEqual(@as(f64, 0.0), neg_result);
+    try expectEqual(@as(f64, 0.0), neg_result);
     const bits: u64 = @bitCast(neg_result);
-    try std.testing.expect((bits >> 63) == 1);
+    try expect((bits >> 63) == 1);
 
     // Exact half rounds to nearest even (banker's rounding)
-    try std.testing.expectEqual(@as(f64, 2.0), rint(2.5));
-    try std.testing.expectEqual(@as(f64, 4.0), rint(3.5));
+    try expectEqual(@as(f64, 2.0), rint(2.5));
+    try expectEqual(@as(f64, 4.0), rint(3.5));
 }

From 014178725744eda46db04ccfe44187ef616ceaf7 Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Sun, 15 Mar 2026 20:24:28 +0100
Subject: [PATCH 06/10] `libzigc`: implement `modff`

`modf` function was generalized and renamed to `modfGeneric`, `modf` and
`modff` provide the appropriate type while calling that function. The
unit tests were also generalized so they can be reused for different
float types.

Both `modf` and `modff` were tested after making these changes:
```
$ stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=modf -fqemu -fwasmtime --summary line
Build Summary: 921/921 steps succeeded
```

```
stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=modff -fqemu -fwasmtime --summary line
Build Summary: 369/369 steps succeeded
```
---
 lib/c/math.zig                 | 64 +++++++++++++++++++++++-----------
 lib/libc/mingw/math/modff.c    | 42 ----------------------
 lib/libc/musl/src/math/modff.c | 34 ------------------
 src/libs/mingw.zig             |  1 -
 src/libs/musl.zig              |  1 -
 src/libs/wasi_libc.zig         |  1 -
 6 files changed, 44 insertions(+), 99 deletions(-)
 delete mode 100644 lib/libc/mingw/math/modff.c
 delete mode 100644 lib/libc/musl/src/math/modff.c

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 2593d5bf86..497154ea04 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -4,6 +4,7 @@ const std = @import("std");
 const math = std.math;
 const expect = std.testing.expect;
 const expectEqual = std.testing.expectEqual;
+const expectApproxEqAbs = std.testing.expectApproxEqAbs;
 const expectApproxEqRel = std.testing.expectApproxEqRel;
 
 const symbol = @import("../c.zig").symbol;
@@ -37,6 +38,7 @@ comptime {
         symbol(&coshf, "coshf");
         symbol(&hypotf, "hypotf");
         symbol(&hypotl, "hypotl");
+        symbol(&modff, "modff");
         symbol(&nan, "nan");
         symbol(&nanf, "nanf");
         symbol(&nanl, "nanl");
@@ -166,66 +168,88 @@ fn isnanl(x: c_longdouble) callconv(.c) c_int {
     return if (math.isNan(x)) 1 else 0;
 }
 
-fn modf(x: f64, iptr: *f64) callconv(.c) f64 {
+fn modfGeneric(comptime T: type, x: T, iptr: *T) T {
     if (math.isNegativeInf(x)) {
-        iptr.* = -math.inf(f64);
+        iptr.* = -math.inf(T);
         return -0.0;
     }
 
     if (math.isPositiveInf(x)) {
-        iptr.* = math.inf(f64);
+        iptr.* = math.inf(T);
         return 0.0;
     }
 
     // Avoids raising the INVALID flag on qemu-riscv
     if (math.isNan(x)) {
-        iptr.* = math.nan(f64);
-        return math.nan(f64);
+        iptr.* = math.nan(T);
+        return math.nan(T);
     }
 
     const r = math.modf(x);
     iptr.* = r.ipart;
 
-    // If the result would be a negative zero, we must be explicit about
+    // If the result is a negative zero, we must be explicit about
     // returning a negative zero.
     return if (math.isNegativeZero(x) or (x < 0.0 and x == r.ipart)) -0.0 else r.fpart;
 }
 
-test "modf" {
-    var int: f64 = undefined;
-    const iptr = &int;
-    const eps_val = 1e-6;
+fn modf(x: f64, iptr: *f64) callconv(.c) f64 {
+    return modfGeneric(f64, x, iptr);
+}
 
-    const normal_frac = modf(1234.5678, iptr);
-    try expectApproxEqRel(0.5678, normal_frac, eps_val);
-    try expectApproxEqRel(1234.0, iptr.*, eps_val);
+fn modff(x: f32, iptr: *f32) callconv(.c) f32 {
+    return modfGeneric(f32, x, iptr);
+}
+
+fn testModf(comptime T: type) !void {
+    // Choose the appropriate `modf` impl to test based on type
+    const f = switch (T) {
+        f64 => modf,
+        f32 => modff,
+        else => @compileError("modf not implemented for " ++ @typeName(T)),
+    };
+
+    var int: T = undefined;
+    const iptr = &int;
+    const eps_val: comptime_float = @max(1e-6, math.floatEps(T));
+
+    const normal_frac = f(@as(T, 1234.567), iptr);
+    // Account for precision error
+    const expected = 1234.567 - @as(T, 1234);
+    try expectApproxEqAbs(expected, normal_frac, eps_val);
+    try expectApproxEqRel(@as(T, 1234.0), iptr.*, eps_val);
 
     // When `x` is a NaN, NaN is returned and `*iptr` is set to NaN
-    const nan_frac = modf(math.nan(f64), iptr);
+    const nan_frac = f(math.nan(T), iptr);
     try expect(math.isNan(nan_frac));
     try expect(math.isNan(iptr.*));
 
     // When `x` is positive infinity, +0 is returned and `*iptr` is set to
     // positive infinity
-    const pos_zero_frac = modf(math.inf(f64), iptr);
+    const pos_zero_frac = f(math.inf(T), iptr);
     try expect(math.isPositiveZero(pos_zero_frac));
     try expect(math.isPositiveInf(iptr.*));
 
     // When `x` is negative infinity, -0 is returned and `*iptr` is set to
     // negative infinity
-    const neg_zero_frac = modf(-math.inf(f64), iptr);
+    const neg_zero_frac = f(-math.inf(T), iptr);
     try expect(math.isNegativeZero(neg_zero_frac));
     try expect(math.isNegativeInf(iptr.*));
 
     // Return -0 when `x` is a negative integer
-    const nz_frac = modf(-1000.0, iptr);
+    const nz_frac = f(@as(T, -1000.0), iptr);
     try expect(math.isNegativeZero(nz_frac));
-    try expectEqual(-1000.0, iptr.*);
+    try expectEqual(@as(T, -1000.0), iptr.*);
 
     // Return +0 when `x` is a positive integer
-    const pz_frac = modf(1000.0, iptr);
+    const pz_frac = f(@as(T, 1000.0), iptr);
     try expect(math.isPositiveZero(pz_frac));
-    try expectEqual(1000.0, iptr.*);
+    try expectEqual(@as(T, 1000.0), iptr.*);
+}
+
+test "modf" {
+    try testModf(f64);
+    try testModf(f32);
 }
 
 fn nan(_: [*:0]const c_char) callconv(.c) f64 {
diff --git a/lib/libc/mingw/math/modff.c b/lib/libc/mingw/math/modff.c
deleted file mode 100644
index d2f8514d4d..0000000000
--- a/lib/libc/mingw/math/modff.c
+++ /dev/null
@@ -1,42 +0,0 @@
-/**
- * This file has no copyright assigned and is placed in the Public Domain.
- * This file is part of the mingw-w64 runtime package.
- * No warranty is given; refer to the file DISCLAIMER.PD within this package.
- */
-#include <fenv.h>
-#include <math.h>
-#include <errno.h>
-
-float
-modff (float value, float* iptr)
-{
-  float int_part = 0.0F;
-  /* truncate */ 
-  /* truncate */
-#if (defined(_AMD64_) && !defined(_ARM64EC_)) || (defined(__x86_64__) && !defined(__arm64ec__))
-  asm volatile ("subq $8, %%rsp\n"
-    "fnstcw 4(%%rsp)\n"
-    "movzwl 4(%%rsp), %%eax\n"
-    "orb $12, %%ah\n"
-    "movw %%ax, (%%rsp)\n"
-    "fldcw (%%rsp)\n"
-    "frndint\n"
-    "fldcw 4(%%rsp)\n"
-    "addq $8, %%rsp\n" : "=t" (int_part) : "0" (value) : "eax"); /* round */
-#elif defined(_X86_) || defined(__i386__)
-  asm volatile ("push %%eax\n\tsubl $8, %%esp\n"
-    "fnstcw 4(%%esp)\n"
-    "movzwl 4(%%esp), %%eax\n"
-    "orb $12, %%ah\n"
-    "movw %%ax, (%%esp)\n"
-    "fldcw (%%esp)\n"
-    "frndint\n"
-    "fldcw 4(%%esp)\n"
-    "addl $8, %%esp\n\tpop %%eax\n" : "=t" (int_part) : "0" (value) : "eax"); /* round */
-#else
-  int_part = truncf(value);
-#endif
-  if (iptr)
-    *iptr = int_part;
-  return (isinf (value) ?  0.0F : value - int_part);
-}
diff --git a/lib/libc/musl/src/math/modff.c b/lib/libc/musl/src/math/modff.c
deleted file mode 100644
index 639514effa..0000000000
--- a/lib/libc/musl/src/math/modff.c
+++ /dev/null
@@ -1,34 +0,0 @@
-#include "libm.h"
-
-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;
-}
diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig
index 70e0d6e8c0..8168febe6d 100644
--- a/src/libs/mingw.zig
+++ b/src/libs/mingw.zig
@@ -977,7 +977,6 @@ const mingw32_x86_src = [_][]const u8{
 
 const mingw32_x86_32_src = [_][]const u8{
     // ucrtbase
-    "math" ++ path.sep_str ++ "modff.c",
     "math" ++ path.sep_str ++ "powf.c",
     "math" ++ path.sep_str ++ "sinhf.c",
     "math" ++ path.sep_str ++ "tanhf.c",
diff --git a/src/libs/musl.zig b/src/libs/musl.zig
index 3db43a1c76..ee7f63bf56 100644
--- a/src/libs/musl.zig
+++ b/src/libs/musl.zig
@@ -934,7 +934,6 @@ const src_files = [_][]const u8{
     "musl/src/math/__math_uflowf.c",
     "musl/src/math/__math_xflow.c",
     "musl/src/math/__math_xflowf.c",
-    "musl/src/math/modff.c",
     "musl/src/math/modfl.c",
     "musl/src/math/nearbyint.c",
     "musl/src/math/nearbyintf.c",
diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig
index 5845b4784c..3077aad385 100644
--- a/src/libs/wasi_libc.zig
+++ b/src/libs/wasi_libc.zig
@@ -755,7 +755,6 @@ const libc_top_half_src_files = [_][]const u8{
     "musl/src/math/__math_uflowf.c",
     "musl/src/math/__math_xflow.c",
     "musl/src/math/__math_xflowf.c",
-    "musl/src/math/modff.c",
     "musl/src/math/modfl.c",
     "musl/src/math/nearbyintl.c",
     "musl/src/math/nextafter.c",

From 4cad6f5372f6f17e249c24fba4d23466fa853fc0 Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Sun, 15 Mar 2026 21:22:01 +0100
Subject: [PATCH 07/10] `libzigc`: Implement `modfl`

The changes were tested by running:

```
$ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib

$ stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=modfl -fqemu -fwasmtime --summary line
Build Summary: 369/369 steps succeeded
```
---
 lib/c/math.zig                 | 11 +++++--
 lib/libc/mingw/math/modfl.c    | 41 --------------------------
 lib/libc/musl/src/math/modfl.c | 53 ----------------------------------
 src/libs/mingw.zig             |  1 -
 src/libs/musl.zig              |  1 -
 src/libs/wasi_libc.zig         |  1 -
 6 files changed, 9 insertions(+), 99 deletions(-)
 delete mode 100644 lib/libc/mingw/math/modfl.c
 delete mode 100644 lib/libc/musl/src/math/modfl.c

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 497154ea04..6f86c925df 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -39,6 +39,7 @@ comptime {
         symbol(&hypotf, "hypotf");
         symbol(&hypotl, "hypotl");
         symbol(&modff, "modff");
+        symbol(&modfl, "modfl");
         symbol(&nan, "nan");
         symbol(&nanf, "nanf");
         symbol(&nanl, "nanl");
@@ -201,11 +202,16 @@ fn modff(x: f32, iptr: *f32) callconv(.c) f32 {
     return modfGeneric(f32, x, iptr);
 }
 
+fn modfl(x: c_longdouble, iptr: *c_longdouble) callconv(.c) c_longdouble {
+    return modfGeneric(c_longdouble, x, iptr);
+}
+
 fn testModf(comptime T: type) !void {
     // Choose the appropriate `modf` impl to test based on type
     const f = switch (T) {
-        f64 => modf,
         f32 => modff,
+        f64 => modf,
+        c_longdouble => modfl,
         else => @compileError("modf not implemented for " ++ @typeName(T)),
     };
 
@@ -248,8 +254,9 @@ fn testModf(comptime T: type) !void {
 }
 
 test "modf" {
-    try testModf(f64);
     try testModf(f32);
+    try testModf(f64);
+    try testModf(c_longdouble);
 }
 
 fn nan(_: [*:0]const c_char) callconv(.c) f64 {
diff --git a/lib/libc/mingw/math/modfl.c b/lib/libc/mingw/math/modfl.c
deleted file mode 100644
index 60b9b81b92..0000000000
--- a/lib/libc/mingw/math/modfl.c
+++ /dev/null
@@ -1,41 +0,0 @@
-/**
- * This file has no copyright assigned and is placed in the Public Domain.
- * This file is part of the mingw-w64 runtime package.
- * No warranty is given; refer to the file DISCLAIMER.PD within this package.
- */
-#include <fenv.h>
-#include <math.h>
-#include <errno.h>
-
-long double
-modfl (long double value, long double* iptr)
-{
-  long double int_part = 0.0L;
-  /* truncate */
-#if (defined(_AMD64_) && !defined(_ARM64EC_)) || (defined(__x86_64__) && !defined(__arm64ec__))
-  asm volatile ("subq $8, %%rsp\n"
-    "fnstcw 4(%%rsp)\n"
-    "movzwl 4(%%rsp), %%eax\n"
-    "orb $12, %%ah\n"
-    "movw %%ax, (%%rsp)\n"
-    "fldcw (%%rsp)\n"
-    "frndint\n"
-    "fldcw 4(%%rsp)\n"
-    "addq $8, %%rsp\n" : "=t" (int_part) : "0" (value) : "eax"); /* round */
-#elif defined(_X86_) || defined(__i386__)
-  asm volatile ("push %%eax\n\tsubl $8, %%esp\n"
-    "fnstcw 4(%%esp)\n"
-    "movzwl 4(%%esp), %%eax\n"
-    "orb $12, %%ah\n"
-    "movw %%ax, (%%esp)\n"
-    "fldcw (%%esp)\n"
-    "frndint\n"
-    "fldcw 4(%%esp)\n"
-    "addl $8, %%esp\n\tpop %%eax\n" : "=t" (int_part) : "0" (value) : "eax"); /* round */
-#else
-  int_part = truncl(value);
-#endif
-  if (iptr)
-    *iptr = int_part;
-  return (isinf (value) ?  0.0L : value - int_part);
-}
diff --git a/lib/libc/musl/src/math/modfl.c b/lib/libc/musl/src/math/modfl.c
deleted file mode 100644
index a47b1924f7..0000000000
--- a/lib/libc/musl/src/math/modfl.c
+++ /dev/null
@@ -1,53 +0,0 @@
-#include "libm.h"
-
-#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
-long double modfl(long double x, long double *iptr)
-{
-	double d;
-	long double r;
-
-	r = modf(x, &d);
-	*iptr = d;
-	return r;
-}
-#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-
-static const long double toint = 1/LDBL_EPSILON;
-
-long double modfl(long double x, long double *iptr)
-{
-	union ldshape u = {x};
-	int e = (u.i.se & 0x7fff) - 0x3fff;
-	int s = u.i.se >> 15;
-	long double absx;
-	long double y;
-
-	/* no fractional part */
-	if (e >= LDBL_MANT_DIG-1) {
-		*iptr = x;
-		if (isnan(x))
-			return x;
-		return s ? -0.0 : 0.0;
-	}
-
-	/* no integral part*/
-	if (e < 0) {
-		*iptr = s ? -0.0 : 0.0;
-		return x;
-	}
-
-	/* raises spurious inexact */
-	absx = s ? -x : x;
-	y = absx + toint - toint - absx;
-	if (y == 0) {
-		*iptr = x;
-		return s ? -0.0 : 0.0;
-	}
-	if (y > 0)
-		y -= 1;
-	if (s)
-		y = -y;
-	*iptr = x + y;
-	return -y;
-}
-#endif
diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig
index 8168febe6d..312a9a1e26 100644
--- a/src/libs/mingw.zig
+++ b/src/libs/mingw.zig
@@ -619,7 +619,6 @@ const mingw32_generic_src = [_][]const u8{
     "math" ++ path.sep_str ++ "lgamma.c",
     "math" ++ path.sep_str ++ "lgammaf.c",
     "math" ++ path.sep_str ++ "lgammal.c",
-    "math" ++ path.sep_str ++ "modfl.c",
     "math" ++ path.sep_str ++ "powi.c",
     "math" ++ path.sep_str ++ "powif.c",
     "math" ++ path.sep_str ++ "powil.c",
diff --git a/src/libs/musl.zig b/src/libs/musl.zig
index ee7f63bf56..9ac12e26e5 100644
--- a/src/libs/musl.zig
+++ b/src/libs/musl.zig
@@ -934,7 +934,6 @@ const src_files = [_][]const u8{
     "musl/src/math/__math_uflowf.c",
     "musl/src/math/__math_xflow.c",
     "musl/src/math/__math_xflowf.c",
-    "musl/src/math/modfl.c",
     "musl/src/math/nearbyint.c",
     "musl/src/math/nearbyintf.c",
     "musl/src/math/nearbyintl.c",
diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig
index 3077aad385..35640264ef 100644
--- a/src/libs/wasi_libc.zig
+++ b/src/libs/wasi_libc.zig
@@ -755,7 +755,6 @@ const libc_top_half_src_files = [_][]const u8{
     "musl/src/math/__math_uflowf.c",
     "musl/src/math/__math_xflow.c",
     "musl/src/math/__math_xflowf.c",
-    "musl/src/math/modfl.c",
     "musl/src/math/nearbyintl.c",
     "musl/src/math/nextafter.c",
     "musl/src/math/nextafterf.c",

From fd680ba5fbe99dde71b6ec28ecae1cddf4fed71b Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Sun, 15 Mar 2026 23:13:24 +0100
Subject: [PATCH 08/10] `libzigc`: Implement `tanh`

The changes were tested by running following commands:

```
$ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib

$ stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=tanh -fqemu -fwasmtime --summary line
Build Summary: 1289/1289 steps succeeded
```
---
 lib/c/math.zig                |  5 ++++
 lib/libc/musl/src/math/tanh.c | 45 -----------------------------------
 src/libs/musl.zig             |  1 -
 src/libs/wasi_libc.zig        |  1 -
 4 files changed, 5 insertions(+), 47 deletions(-)
 delete mode 100644 lib/libc/musl/src/math/tanh.c

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 6f86c925df..3520c63457 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -63,6 +63,7 @@ comptime {
         symbol(&pow, "pow");
         symbol(&pow10, "pow10");
         symbol(&pow10f, "pow10f");
+        symbol(&tanh, "tanh");
     }
 
     if (builtin.target.isMuslLibC()) {
@@ -338,3 +339,7 @@ test "rint" {
     try expectEqual(@as(f64, 2.0), rint(2.5));
     try expectEqual(@as(f64, 4.0), rint(3.5));
 }
+
+fn tanh(x: f64) callconv(.c) f64 {
+    return math.tanh(x);
+}
diff --git a/lib/libc/musl/src/math/tanh.c b/lib/libc/musl/src/math/tanh.c
deleted file mode 100644
index 20d6dbcf41..0000000000
--- a/lib/libc/musl/src/math/tanh.c
+++ /dev/null
@@ -1,45 +0,0 @@
-#include "libm.h"
-
-/* tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x))
- *         = (exp(2*x) - 1)/(exp(2*x) - 1 + 2)
- *         = (1 - exp(-2*x))/(exp(-2*x) - 1 + 2)
- */
-double tanh(double x)
-{
-	union {double f; uint64_t i;} u = {.f = x};
-	uint32_t w;
-	int sign;
-	double_t t;
-
-	/* x = |x| */
-	sign = u.i >> 63;
-	u.i &= (uint64_t)-1/2;
-	x = u.f;
-	w = u.i >> 32;
-
-	if (w > 0x3fe193ea) {
-		/* |x| > log(3)/2 ~= 0.5493 or nan */
-		if (w > 0x40340000) {
-			/* |x| > 20 or nan */
-			/* note: this branch avoids raising overflow */
-			t = 1 - 0/x;
-		} else {
-			t = expm1(2*x);
-			t = 1 - 2/(t+2);
-		}
-	} else if (w > 0x3fd058ae) {
-		/* |x| > log(5/3)/2 ~= 0.2554 */
-		t = expm1(2*x);
-		t = t/(t+2);
-	} else if (w >= 0x00100000) {
-		/* |x| >= 0x1p-1022, up to 2ulp error in [0.1,0.2554] */
-		t = expm1(-2*x);
-		t = -t/(t+2);
-	} else {
-		/* |x| is subnormal */
-		/* note: the branch above would not raise underflow in [0x1p-1023,0x1p-1022) */
-		FORCE_EVAL((float)x);
-		t = x;
-	}
-	return sign ? -t : t;
-}
diff --git a/src/libs/musl.zig b/src/libs/musl.zig
index 9ac12e26e5..4e5408beca 100644
--- a/src/libs/musl.zig
+++ b/src/libs/musl.zig
@@ -1003,7 +1003,6 @@ const src_files = [_][]const u8{
     "musl/src/math/sinl.c",
     "musl/src/math/__tan.c",
     "musl/src/math/__tandf.c",
-    "musl/src/math/tanh.c",
     "musl/src/math/tanhf.c",
     "musl/src/math/tanhl.c",
     "musl/src/math/__tanl.c",
diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig
index 35640264ef..1d35b0bf19 100644
--- a/src/libs/wasi_libc.zig
+++ b/src/libs/wasi_libc.zig
@@ -793,7 +793,6 @@ const libc_top_half_src_files = [_][]const u8{
     "musl/src/math/sinl.c",
     "musl/src/math/__tan.c",
     "musl/src/math/__tandf.c",
-    "musl/src/math/tanh.c",
     "musl/src/math/tanhf.c",
     "musl/src/math/tanhl.c",
     "musl/src/math/__tanl.c",

From 691416c24d3c99574f1b168931f8fcd6227760b8 Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Sun, 15 Mar 2026 23:28:13 +0100
Subject: [PATCH 09/10] `libzigc`: Implement `tanhf`

The changes were tested by running:

```
$ ./build/stage3/bin/zig build -p stage4 -Denable-llvm -Dno-lib

$ stage4/bin/zig build test-libc -Dlibc-test-path=<LIBC-TEST-PATH> -Dtest-filter=tanhf -fqemu -fwasmtime --summary line
Build Summary: 553/553 steps succeeded
```
---
 lib/c/math.zig                 |  5 +++++
 lib/libc/mingw/math/tanhf.c    | 10 ---------
 lib/libc/musl/src/math/tanhf.c | 39 ----------------------------------
 src/libs/mingw.zig             |  1 -
 src/libs/musl.zig              |  1 -
 src/libs/wasi_libc.zig         |  1 -
 6 files changed, 5 insertions(+), 52 deletions(-)
 delete mode 100644 lib/libc/mingw/math/tanhf.c
 delete mode 100644 lib/libc/musl/src/math/tanhf.c

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 3520c63457..15022cb4a6 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -43,6 +43,7 @@ comptime {
         symbol(&nan, "nan");
         symbol(&nanf, "nanf");
         symbol(&nanl, "nanl");
+        symbol(&tanhf, "tanhf");
     }
 
     if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
@@ -343,3 +344,7 @@ test "rint" {
 fn tanh(x: f64) callconv(.c) f64 {
     return math.tanh(x);
 }
+
+fn tanhf(x: f32) callconv(.c) f32 {
+    return math.tanh(x);
+}
diff --git a/lib/libc/mingw/math/tanhf.c b/lib/libc/mingw/math/tanhf.c
deleted file mode 100644
index cb6f4c75a4..0000000000
--- a/lib/libc/mingw/math/tanhf.c
+++ /dev/null
@@ -1,10 +0,0 @@
-/**
- * This file has no copyright assigned and is placed in the Public Domain.
- * This file is part of the mingw-w64 runtime package.
- * No warranty is given; refer to the file DISCLAIMER.PD within this package.
- */
-#include <math.h>
-float tanhf (float x)
-{
-  return (float) tanh (x);
-}
diff --git a/lib/libc/musl/src/math/tanhf.c b/lib/libc/musl/src/math/tanhf.c
deleted file mode 100644
index 10636fbd7b..0000000000
--- a/lib/libc/musl/src/math/tanhf.c
+++ /dev/null
@@ -1,39 +0,0 @@
-#include "libm.h"
-
-float tanhf(float x)
-{
-	union {float f; uint32_t i;} u = {.f = x};
-	uint32_t w;
-	int sign;
-	float t;
-
-	/* x = |x| */
-	sign = u.i >> 31;
-	u.i &= 0x7fffffff;
-	x = u.f;
-	w = u.i;
-
-	if (w > 0x3f0c9f54) {
-		/* |x| > log(3)/2 ~= 0.5493 or nan */
-		if (w > 0x41200000) {
-			/* |x| > 10 */
-			t = 1 + 0/x;
-		} else {
-			t = expm1f(2*x);
-			t = 1 - 2/(t+2);
-		}
-	} else if (w > 0x3e82c578) {
-		/* |x| > log(5/3)/2 ~= 0.2554 */
-		t = expm1f(2*x);
-		t = t/(t+2);
-	} else if (w >= 0x00800000) {
-		/* |x| >= 0x1p-126 */
-		t = expm1f(-2*x);
-		t = -t/(t+2);
-	} else {
-		/* |x| is subnormal */
-		FORCE_EVAL(x*x);
-		t = x;
-	}
-	return sign ? -t : t;
-}
diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig
index 312a9a1e26..07e67ea468 100644
--- a/src/libs/mingw.zig
+++ b/src/libs/mingw.zig
@@ -978,7 +978,6 @@ const mingw32_x86_32_src = [_][]const u8{
     // ucrtbase
     "math" ++ path.sep_str ++ "powf.c",
     "math" ++ path.sep_str ++ "sinhf.c",
-    "math" ++ path.sep_str ++ "tanhf.c",
     "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "acosf.c",
     "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "asinf.c",
     "math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "atan2f.c",
diff --git a/src/libs/musl.zig b/src/libs/musl.zig
index 4e5408beca..b96ffbf18a 100644
--- a/src/libs/musl.zig
+++ b/src/libs/musl.zig
@@ -1003,7 +1003,6 @@ const src_files = [_][]const u8{
     "musl/src/math/sinl.c",
     "musl/src/math/__tan.c",
     "musl/src/math/__tandf.c",
-    "musl/src/math/tanhf.c",
     "musl/src/math/tanhl.c",
     "musl/src/math/__tanl.c",
     "musl/src/math/tanl.c",
diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig
index 1d35b0bf19..de7cabd542 100644
--- a/src/libs/wasi_libc.zig
+++ b/src/libs/wasi_libc.zig
@@ -793,7 +793,6 @@ const libc_top_half_src_files = [_][]const u8{
     "musl/src/math/sinl.c",
     "musl/src/math/__tan.c",
     "musl/src/math/__tandf.c",
-    "musl/src/math/tanhf.c",
     "musl/src/math/tanhl.c",
     "musl/src/math/__tanl.c",
     "musl/src/math/tanl.c",

From 153990d407e53d965bbb0c2c9ece333737c90363 Mon Sep 17 00:00:00 2001
From: mihael <hi@mihaelm.com>
Date: Tue, 17 Mar 2026 23:31:53 +0100
Subject: [PATCH 10/10] `libzigc`: Remove qemu-riscv comment

---
 lib/c/math.zig | 1 -
 1 file changed, 1 deletion(-)

diff --git a/lib/c/math.zig b/lib/c/math.zig
index 15022cb4a6..3f7b9ca237 100644
--- a/lib/c/math.zig
+++ b/lib/c/math.zig
@@ -182,7 +182,6 @@ fn modfGeneric(comptime T: type, x: T, iptr: *T) T {
         return 0.0;
     }
 
-    // Avoids raising the INVALID flag on qemu-riscv
     if (math.isNan(x)) {
         iptr.* = math.nan(T);
         return math.nan(T);