diff --git a/lib/compiler_rt/sqrt.zig b/lib/compiler_rt/sqrt.zig
index 2c321863682a547f405708cf97c52e0acf3f5a17..870b1ee560ae680fd731107577d0d21f3b948f38 100644
--- a/lib/compiler_rt/sqrt.zig
+++ b/lib/compiler_rt/sqrt.zig
@@ -54,7 +54,7 @@ pub fn __sqrth(x: f16) callconv(.c) f16 {
     //   m: 2.14 r: 0.16, s: 2.14, d: 2.14, u: 2.14, three: 2.14
     const three: u16 = 0xC000;
     const i: usize = @intCast((ix >> 4) & 0x7F);
-    const r = __rsqrt_tab[i];
+    const r = rsqrt_tab[i];
     // |r*sqrt(m) - 1| < 0x1p-8
     var s = mul16(m, r);
     // |s/sqrt(m) - 1| < 0x1p-8
@@ -92,74 +92,56 @@ pub fn __sqrth(x: f16) callconv(.c) f16 {
 
 pub fn sqrtf(x: f32) callconv(.c) f32 {
     var ix: u32 = @bitCast(x);
-    var top = ix >> 23;
 
-    // special case handling.
-    if (top -% 0x01 >= 0xFF - 0x01) {
+    if (ix < @as(u32, @bitCast(@as(f32, 0x1p-126))) or @as(u32, @bitCast(std.math.inf(f32))) <= ix) {
         @branchHint(.unlikely);
-        // x < 0x1p-126 or inf or nan.
-        if (ix & 0x7FFF_FFFF == 0) return x;
-        if (ix == 0x7F80_0000) return x;
-        if (ix > 0x7F80_0000) return math.nan(f32);
-        // x is subnormal, normalize it.
-        ix = @bitCast(x * 0x1p23);
-        top = (ix >> 23) -% 23;
+
+        if (ix & 0x7fffffff == 0)
+            return x;
+
+        if (ix == @as(u32, @bitCast(std.math.inf(f32))))
+            return x;
+
+        if (ix > @as(u32, @bitCast(std.math.inf(f32))))
+            return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f32);
+
+        ix = @as(u32, @bitCast(@as(i32, @bitCast(x * 0x1p23)) - (23 << 23)));
     }
 
-    // argument reduction:
-    // x = 4^e m; with integer e, and m in [1, 4)
-    // m: fixed point representation [2.30]
-    // 2^e is the exponent part of the result.
-    const even = (top & 1) != 0;
-    const m = if (even) (ix << 7) & 0x7FFF_FFFF else (ix << 8) | 0x8000_0000;
-    top = (top +% 0x7F) >> 1;
+    const m: u32 = if (ix & 0x00800000 != 0)
+        (ix << 7) & 0x7fffffff
+    else
+        (ix << 8) | 0x80000000;
+
+    const ey = ((ix >> 1) + (0x3f800000 >> 1)) & 0x7f800000;
+    // const ey = ((ix + 0x3f800000) & 0xff000000) >> 1;
+
+    const three = 0xc0000000;
+    const i = (ix >> 17) & 0x7f;
+    var r = @as(u32, rsqrt_tab[i]) << 16;
 
-    // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
-    // the fixed point representations are
-    //   m: 2.30 r: 0.32, s: 2.30, d: 2.30, u: 2.30, three: 2.30
-    const three: u32 = 0xC000_0000;
-    var i: usize = @intCast((ix >> 17) & 0x3F);
-    if (even) i += 64;
-    var r = @as(u32, @intCast(__rsqrt_tab[i])) << 16;
-    // |r*sqrt(m) - 1| < 0x1p-8
     var s = mul32(m, r);
-    // |s/sqrt(m) - 1| < 0x1p-8
     var d = mul32(s, r);
     var u = three - d;
     r = mul32(r, u) << 1;
-    // |r*sqrt(m) - 1| < 0x1.7bp-16
     s = mul32(s, u) << 1;
-    // |s/sqrt(m) - 1| < 0x1.7bp-16
     d = mul32(s, r);
     u = three - d;
-    s = mul32(s, u); // repr: 3.29
-    // -0x1.03p-28 < s/sqrt(m) - 1 < 0x1.fp-31
-    s = (s - 1) >> 6; // repr: 9.23
-    // s < sqrt(m) < s + 0x1.08p-23
+    s = mul32(s, u);
+    s = (s - 1) >> 6;
 
-    // compute nearest rounded result:
-    // the nearest result to 23 bits is either s or s+0x1p-23,
-    // we can decide by comparing (2^23 s + 0.5)^2 to 2^46 m.
     const d0 = (m << 16) -% s *% s;
     const d1 = s -% d0;
     const d2 = d1 +% s +% 1;
-    s += d1 >> 31;
-    s &= 0x007F_FFFF;
-    s |= top << 23;
-    const y: f32 = @bitCast(s);
+    const y: f32 = @bitCast(((s + (d1 >> 31)) & 0x007fffff) | ey);
 
-    // handle rounding modes and inexact exception:
-    // only (s+1)^2 == 2^16 m case is exact otherwise
-    // add a tiny value to cause the fenv effects.
-    if (d2 != 0) {
-        @branchHint(.likely);
-        var tiny: u32 = 0x0100_0000;
-        tiny |= (d1 ^ d2) & 0x8000_0000;
-        const t: f32 = @bitCast(tiny);
-        return y + t;
-    }
+    const tiny: u32 = if (d2 == 0) blk: {
+        @branchHint(.unlikely);
+        break :blk 0;
+    } else 0x01000000;
+    const t: f32 = @bitCast(tiny | ((d1 ^ d2) & 0x80000000));
 
-    return y;
+    return y + t;
 }
 
 pub fn sqrt(x: f64) callconv(.c) f64 {
@@ -172,7 +154,7 @@ pub fn sqrt(x: f64) callconv(.c) f64 {
         // x < 0x1p-1022 or inf or nan.
         if (ix & 0x7FFF_FFFF_FFFF_FFFF == 0) return x;
         if (ix == 0x7FF0_0000_0000_0000) return x;
-        if (ix > 0x7FF0_0000_0000_0000) return math.nan(f64);
+        if (ix > 0x7FF0_0000_0000_0000) return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f64);
         // x is subnormal, normalize it.
         ix = @bitCast(x * 0x1p52);
         top = (ix >> 52) -% 52;
@@ -248,7 +230,7 @@ pub fn sqrt(x: f64) callconv(.c) f64 {
     var d: struct { u32, u64 } = undefined;
     var u: struct { u32, u64 } = undefined;
     const i: usize = @intCast((ix >> 46) & 0x7F);
-    r[0] = @intCast(__rsqrt_tab[i]);
+    r[0] = @intCast(rsqrt_tab[i]);
     r[0] <<= 16;
     // |r sqrt(m) - 1| < 0x1.fdp-9
     s[0] = mul32(@intCast(m >> 32), r[0]);
@@ -309,7 +291,7 @@ pub fn __sqrtx(x: f80) callconv(.c) f80 {
         // x < 0x1p-16382 or inf or nan.
         if (ix & 0x7FFF_FFFF_FFFF_FFFF_FFFF == 0) return x;
         if (ix == 0x7FFF_8000_0000_0000_0000) return x;
-        if (ix > 0x7FFF_8000_0000_0000_0000) return math.nan(f80);
+        if (ix > 0x7FFF_8000_0000_0000_0000) return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f80);
         // x is subnormal, normalize it.
         ix = @bitCast(x * 0x1p63);
         top = (ix >> 64) -% 63;
@@ -341,7 +323,7 @@ pub fn __sqrtx(x: f80) callconv(.c) f80 {
     var u: struct { u32, u64, u80 } = undefined;
     var i: usize = @intCast((ix >> 57) & 0x3F);
     if (even) i += 64;
-    r[0] = @intCast(__rsqrt_tab[i]);
+    r[0] = @intCast(rsqrt_tab[i]);
     r[0] <<= 16;
     // |r sqrt(m) - 1| < 0x1p-8
     s[0] = mul32(@intCast(m >> 48), r[0]);
@@ -437,7 +419,7 @@ pub fn sqrtq(x: f128) callconv(.c) f128 {
     var d: struct { u32, u64, u128 } = undefined;
     var u: struct { u32, u64, u128 } = undefined;
     const i: usize = @intCast((ix >> 106) & 0x7F);
-    r[0] = @intCast(__rsqrt_tab[i]);
+    r[0] = @intCast(rsqrt_tab[i]);
     r[0] <<= 16;
     // |r sqrt(m) - 1| < 0x1p-8
     s[0] = mul32(@intCast(m >> 96), r[0]);
@@ -507,7 +489,7 @@ pub fn sqrtl(x: c_longdouble) callconv(.c) c_longdouble {
     }
 }
 
-const __rsqrt_tab: [128]u16 = .{
+const rsqrt_tab: [128]u16 = .{
     0xB451, 0xB2F0, 0xB196, 0xB044, 0xAEF9, 0xADB6, 0xAC79, 0xAB43,
     0xAA14, 0xA8EB, 0xA7C8, 0xA6AA, 0xA592, 0xA480, 0xA373, 0xA26B,
     0xA168, 0xA06A, 0x9F70, 0x9E7B, 0x9D8A, 0x9C9D, 0x9BB5, 0x9AD1,
@@ -527,15 +509,15 @@ const __rsqrt_tab: [128]u16 = .{
 };
 
 inline fn mul16(a: u16, b: u16) u16 {
-    return @intCast(@as(u32, @intCast(a)) * @as(u32, @intCast(b)) >> 16);
+    return @intCast(@as(u32, a) * b >> 16);
 }
 
 inline fn mul32(a: u32, b: u32) u32 {
-    return @intCast(@as(u64, @intCast(a)) * @as(u64, @intCast(b)) >> 32);
+    return @intCast(@as(u64, a) * b >> 32);
 }
 
 inline fn mul64(a: u64, b: u64) u64 {
-    return @intCast(@as(u128, @intCast(a)) * @as(u128, @intCast(b)) >> 64);
+    return @intCast(@as(u128, a) * b >> 64);
 }
 
 inline fn mul80(a: u80, b: u80) u80 {
diff --git a/lib/libc/mingw/math/sqrt.def.h b/lib/libc/mingw/math/sqrt.def.h
deleted file mode 100644
index 43bd5a37d1636a4ef8819599f731046c2eb28cb4..0000000000000000000000000000000000000000
--- a/lib/libc/mingw/math/sqrt.def.h
+++ /dev/null
@@ -1,92 +0,0 @@
-/*
- This Software is provided under the Zope Public License (ZPL) Version 2.1.
-
- Copyright (c) 2009, 2010 by the mingw-w64 project
-
- See the AUTHORS file for the list of contributors to the mingw-w64 project.
-
- This license has been certified as open source. It has also been designated
- as GPL compatible by the Free Software Foundation (FSF).
-
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions are met:
-
-   1. Redistributions in source code must retain the accompanying copyright
-      notice, this list of conditions, and the following disclaimer.
-   2. Redistributions in binary form must reproduce the accompanying
-      copyright notice, this list of conditions, and the following disclaimer
-      in the documentation and/or other materials provided with the
-      distribution.
-   3. Names of the copyright holders must not be used to endorse or promote
-      products derived from this software without prior written permission
-      from the copyright holders.
-   4. The right to distribute this software or to use it for any purpose does
-      not give you the right to use Servicemarks (sm) or Trademarks (tm) of
-      the copyright holders.  Use of them is covered by separate agreement
-      with the copyright holders.
-   5. If any files are modified, you must cause the modified files to carry
-      prominent notices stating that you changed the files and the date of
-      any change.
-
- Disclaimer
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
- OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
- EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
- INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
- LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, 
- OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
- EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-*/
-
-#include "../complex/complex_internal.h"
-#include <errno.h>
-
-__FLT_TYPE
-__FLT_ABI (sqrt) (__FLT_TYPE x)
-{
-  __FLT_TYPE res = __FLT_CST (0.0);
-  int x_class = fpclassify (x);
-  if (x_class == FP_NAN || signbit (x))
-    {
-      if (x_class == FP_ZERO)
-	return __FLT_CST (-0.0);
-
-      if (x_class == FP_NAN)
-        {
-          __FLT_RPT_DOMAIN ("sqrt", x, 0.0, x);
-          return x;
-        }
-
-      res = -__FLT_NAN;
-      __FLT_RPT_DOMAIN ("sqrt", x, 0.0, res);
-      return res;
-    }
-  else if (x_class == FP_ZERO)
-    return __FLT_CST (0.0);
-  else if (x_class == FP_INFINITE)
-    return __FLT_HUGE_VAL;
-  else if (x == __FLT_CST (1.0))
-   return __FLT_CST (1.0);
-#if defined(__arm__) || defined(_ARM_)
-#if _NEW_COMPLEX_FLOAT
-  asm volatile ("fsqrts %[dst], %[src];\n" : [dst] "=t" (res) : [src] "t" (x));
-#else
-  asm volatile ("fsqrtd %[dst], %[src];\n" : [dst] "=w" (res) : [src] "w" (x));
-#endif
-#elif defined(__aarch64__) || defined(_ARM64_) || defined(__arm64ec__) || defined(_ARM64EC_)
-#if _NEW_COMPLEX_FLOAT
-  asm volatile ("fsqrt %s[dst], %s[src]\n" : [dst] "=w" (res) : [src] "w" (x));
-#else
-  asm volatile ("fsqrt %d[dst], %d[src]\n" : [dst] "=w" (res) : [src] "w" (x));
-#endif
-#elif defined(_X86_) || defined(__i386__) || defined(_AMD64_) || defined(__x86_64__)
-  asm volatile ("fsqrt" : "=t" (res) : "0" (x));
-#else
-#error Not supported on your platform yet
-#endif
-  return res;
-}
diff --git a/lib/libc/mingw/math/sqrtf.c b/lib/libc/mingw/math/sqrtf.c
deleted file mode 100644
index 316fda408d7151c831e32442df4a22ef83d5aaec..0000000000000000000000000000000000000000
--- a/lib/libc/mingw/math/sqrtf.c
+++ /dev/null
@@ -1,46 +0,0 @@
-/*
- This Software is provided under the Zope Public License (ZPL) Version 2.1.
-
- Copyright (c) 2009, 2010 by the mingw-w64 project
-
- See the AUTHORS file for the list of contributors to the mingw-w64 project.
-
- This license has been certified as open source. It has also been designated
- as GPL compatible by the Free Software Foundation (FSF).
-
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions are met:
-
-   1. Redistributions in source code must retain the accompanying copyright
-      notice, this list of conditions, and the following disclaimer.
-   2. Redistributions in binary form must reproduce the accompanying
-      copyright notice, this list of conditions, and the following disclaimer
-      in the documentation and/or other materials provided with the
-      distribution.
-   3. Names of the copyright holders must not be used to endorse or promote
-      products derived from this software without prior written permission
-      from the copyright holders.
-   4. The right to distribute this software or to use it for any purpose does
-      not give you the right to use Servicemarks (sm) or Trademarks (tm) of
-      the copyright holders.  Use of them is covered by separate agreement
-      with the copyright holders.
-   5. If any files are modified, you must cause the modified files to carry
-      prominent notices stating that you changed the files and the date of
-      any change.
-
- Disclaimer
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
- OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
- EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
- INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
- LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, 
- OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
- EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-*/
-
-#define _NEW_COMPLEX_FLOAT 1
-#include "sqrt.def.h"
diff --git a/lib/libc/mingw/math/sqrtl.c b/lib/libc/mingw/math/sqrtl.c
deleted file mode 100644
index ffd8185911fdf6b50e264f43294bf816d184d8c4..0000000000000000000000000000000000000000
--- a/lib/libc/mingw/math/sqrtl.c
+++ /dev/null
@@ -1,46 +0,0 @@
-/*
- This Software is provided under the Zope Public License (ZPL) Version 2.1.
-
- Copyright (c) 2009, 2010 by the mingw-w64 project
-
- See the AUTHORS file for the list of contributors to the mingw-w64 project.
-
- This license has been certified as open source. It has also been designated
- as GPL compatible by the Free Software Foundation (FSF).
-
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions are met:
-
-   1. Redistributions in source code must retain the accompanying copyright
-      notice, this list of conditions, and the following disclaimer.
-   2. Redistributions in binary form must reproduce the accompanying
-      copyright notice, this list of conditions, and the following disclaimer
-      in the documentation and/or other materials provided with the
-      distribution.
-   3. Names of the copyright holders must not be used to endorse or promote
-      products derived from this software without prior written permission
-      from the copyright holders.
-   4. The right to distribute this software or to use it for any purpose does
-      not give you the right to use Servicemarks (sm) or Trademarks (tm) of
-      the copyright holders.  Use of them is covered by separate agreement
-      with the copyright holders.
-   5. If any files are modified, you must cause the modified files to carry
-      prominent notices stating that you changed the files and the date of
-      any change.
-
- Disclaimer
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
- OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
- EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
- INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
- LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, 
- OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
- EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-*/
-
-#define _NEW_COMPLEX_LDOUBLE 1
-#include "sqrt.def.h"
diff --git a/lib/libc/musl/src/math/aarch64/sqrt.c b/lib/libc/musl/src/math/aarch64/sqrt.c
deleted file mode 100644
index fe93c3e6ad554edd593325694bbc8669f640452c..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/aarch64/sqrt.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <math.h>
-
-double sqrt(double x)
-{
-	__asm__ ("fsqrt %d0, %d1" : "=w"(x) : "w"(x));
-	return x;
-}
diff --git a/lib/libc/musl/src/math/aarch64/sqrtf.c b/lib/libc/musl/src/math/aarch64/sqrtf.c
deleted file mode 100644
index 275c7f399c106f76616fb91b1e92e11a50e73387..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/aarch64/sqrtf.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <math.h>
-
-float sqrtf(float x)
-{
-	__asm__ ("fsqrt %s0, %s1" : "=w"(x) : "w"(x));
-	return x;
-}
diff --git a/lib/libc/musl/src/math/arm/sqrt.c b/lib/libc/musl/src/math/arm/sqrt.c
deleted file mode 100644
index 567e2e91015374af42b5ce05f9cde2b92497a81b..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/arm/sqrt.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if (__ARM_PCS_VFP || (__VFP_FP__ && !__SOFTFP__)) && (__ARM_FP&8)
-
-double sqrt(double x)
-{
-	__asm__ ("vsqrt.f64 %P0, %P1" : "=w"(x) : "w"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrt.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/arm/sqrtf.c b/lib/libc/musl/src/math/arm/sqrtf.c
deleted file mode 100644
index 32693293b49b777667923a5008bbb2145ddd930d..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/arm/sqrtf.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if (__ARM_PCS_VFP || (__VFP_FP__ && !__SOFTFP__)) && !BROKEN_VFP_ASM
-
-float sqrtf(float x)
-{
-	__asm__ ("vsqrt.f32 %0, %1" : "=t"(x) : "t"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrtf.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/i386/sqrt.c b/lib/libc/musl/src/math/i386/sqrt.c
deleted file mode 100644
index 934fbccab82c993b51b9ee0d81b19826f73f56ff..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/i386/sqrt.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include "libm.h"
-
-double sqrt(double x)
-{
-	union ldshape ux;
-	unsigned fpsr;
-	__asm__ ("fsqrt; fnstsw %%ax": "=t"(ux.f), "=a"(fpsr) : "0"(x));
-	if ((ux.i.m & 0x7ff) != 0x400)
-		return (double)ux.f;
-	/* Rounding to double would have encountered an exact halfway case.
-	   Adjust mantissa downwards if fsqrt rounded up, else upwards.
-	   (result of fsqrt could not have been exact) */
-	ux.i.m ^= (fpsr & 0x200) + 0x300;
-	return (double)ux.f;
-}
diff --git a/lib/libc/musl/src/math/i386/sqrtf.c b/lib/libc/musl/src/math/i386/sqrtf.c
deleted file mode 100644
index 41c65c2bdbbddaf1fb74ea8f467a547cc27a50ea..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/i386/sqrtf.c
+++ /dev/null
@@ -1,12 +0,0 @@
-#include <math.h>
-
-float sqrtf(float x)
-{
-	long double t;
-	/* The long double result has sufficient precision so that
-	 * second rounding to float still keeps the returned value
-	 * correctly rounded, see Pierre Roux, "Innocuous Double
-	 * Rounding of Basic Arithmetic Operations". */
-	__asm__ ("fsqrt" : "=t"(t) : "0"(x));
-	return (float)t;
-}
diff --git a/lib/libc/musl/src/math/i386/sqrtl.c b/lib/libc/musl/src/math/i386/sqrtl.c
deleted file mode 100644
index 864cfcc4f66eb5c0b8c87f38299fc727eb6e7eac..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/i386/sqrtl.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <math.h>
-
-long double sqrtl(long double x)
-{
-	__asm__ ("fsqrt" : "+t"(x));
-	return x;
-}
diff --git a/lib/libc/musl/src/math/m68k/sqrtl.c b/lib/libc/musl/src/math/m68k/sqrtl.c
deleted file mode 100644
index b1c303c7e2696bd21ece2f50c85857b9dcb4c978..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/m68k/sqrtl.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if __HAVE_68881__
-
-long double sqrtl(long double x)
-{
-	__asm__ ("fsqrt.x %1,%0" : "=f"(x) : "fm"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrtl.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/mips/sqrt.c b/lib/libc/musl/src/math/mips/sqrt.c
deleted file mode 100644
index 595c9dbc36232db1469359f36654e37dba2b534d..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/mips/sqrt.c
+++ /dev/null
@@ -1,16 +0,0 @@
-#if !defined(__mips_soft_float) && __mips >= 3
-
-#include <math.h>
-
-double sqrt(double x)
-{
-	double r;
-	__asm__("sqrt.d %0,%1" : "=f"(r) : "f"(x));
-	return r;
-}
-
-#else
-
-#include "../sqrt.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/mips/sqrtf.c b/lib/libc/musl/src/math/mips/sqrtf.c
deleted file mode 100644
index 84090d2d31f48b39c7c0a3c4ec39014764c7162c..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/mips/sqrtf.c
+++ /dev/null
@@ -1,16 +0,0 @@
-#if !defined(__mips_soft_float) && __mips >= 2
-
-#include <math.h>
-
-float sqrtf(float x)
-{
-	float r;
-	__asm__("sqrt.s %0,%1" : "=f"(r) : "f"(x));
-	return r;
-}
-
-#else
-
-#include "../sqrtf.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/powerpc/sqrt.c b/lib/libc/musl/src/math/powerpc/sqrt.c
deleted file mode 100644
index 8718dbd0ca7b4086fc1ebb2e7cd3b2e6eb3a20e2..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/powerpc/sqrt.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if !defined _SOFT_FLOAT && defined _ARCH_PPCSQ
-
-double sqrt(double x)
-{
-	__asm__ ("fsqrt %0, %1\n" : "=d" (x) : "d" (x));
-	return x;
-}
-
-#else
-
-#include "../sqrt.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/powerpc/sqrtf.c b/lib/libc/musl/src/math/powerpc/sqrtf.c
deleted file mode 100644
index 3431b672d065ffdf0472a8df8f981da16b117981..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/powerpc/sqrtf.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if !defined _SOFT_FLOAT && defined _ARCH_PPCSQ
-
-float sqrtf(float x)
-{
-	__asm__ ("fsqrts %0, %1\n" : "=f" (x) : "f" (x));
-	return x;
-}
-
-#else
-
-#include "../sqrtf.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/powerpc64/sqrt.c b/lib/libc/musl/src/math/powerpc64/sqrt.c
deleted file mode 100644
index 13bb98d91cf204caefe08a36455527c8edb04df5..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/powerpc64/sqrt.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <math.h>
-
-double sqrt(double x)
-{
-	__asm__ ("fsqrt %0, %1" : "=d"(x) : "d"(x));
-	return x;
-}
diff --git a/lib/libc/musl/src/math/powerpc64/sqrtf.c b/lib/libc/musl/src/math/powerpc64/sqrtf.c
deleted file mode 100644
index b6ecb106b687b790bb2228fbd6c1fe79ca820aa2..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/powerpc64/sqrtf.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <math.h>
-
-float sqrtf(float x)
-{
-	__asm__ ("fsqrts %0, %1" : "=f"(x) : "f"(x));
-	return x;
-}
diff --git a/lib/libc/musl/src/math/riscv32/sqrt.c b/lib/libc/musl/src/math/riscv32/sqrt.c
deleted file mode 100644
index 867a504c278c8088e46a18a5a496898e50ce90a6..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/riscv32/sqrt.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if __riscv_flen >= 64
-
-double sqrt(double x)
-{
-	__asm__ ("fsqrt.d %0, %1" : "=f"(x) : "f"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrt.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/riscv32/sqrtf.c b/lib/libc/musl/src/math/riscv32/sqrtf.c
deleted file mode 100644
index 610c2cf800c3f1a2f9e2646f6ffef660abc30346..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/riscv32/sqrtf.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if __riscv_flen >= 32
-
-float sqrtf(float x)
-{
-	__asm__ ("fsqrt.s %0, %1" : "=f"(x) : "f"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrtf.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/riscv64/sqrt.c b/lib/libc/musl/src/math/riscv64/sqrt.c
deleted file mode 100644
index 867a504c278c8088e46a18a5a496898e50ce90a6..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/riscv64/sqrt.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if __riscv_flen >= 64
-
-double sqrt(double x)
-{
-	__asm__ ("fsqrt.d %0, %1" : "=f"(x) : "f"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrt.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/riscv64/sqrtf.c b/lib/libc/musl/src/math/riscv64/sqrtf.c
deleted file mode 100644
index 610c2cf800c3f1a2f9e2646f6ffef660abc30346..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/riscv64/sqrtf.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if __riscv_flen >= 32
-
-float sqrtf(float x)
-{
-	__asm__ ("fsqrt.s %0, %1" : "=f"(x) : "f"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrtf.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/s390x/sqrt.c b/lib/libc/musl/src/math/s390x/sqrt.c
deleted file mode 100644
index a80dc4a7d9eb1f06c7fe5a889b56ede26101a298..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/s390x/sqrt.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if defined(__HTM__) || __ARCH__ >= 9
-
-double sqrt(double x)
-{
-	__asm__ ("sqdbr %0, %1" : "=f"(x) : "f"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrt.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/s390x/sqrtf.c b/lib/libc/musl/src/math/s390x/sqrtf.c
deleted file mode 100644
index 52f57eb99810c253768353ceec92f633b39c0854..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/s390x/sqrtf.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if defined(__HTM__) || __ARCH__ >= 9
-
-float sqrtf(float x)
-{
-	__asm__ ("sqebr %0, %1" : "=f"(x) : "f"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrtf.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/s390x/sqrtl.c b/lib/libc/musl/src/math/s390x/sqrtl.c
deleted file mode 100644
index 5702d54da0cf0ed3124eb1057852056630ecaa8d..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/s390x/sqrtl.c
+++ /dev/null
@@ -1,15 +0,0 @@
-#include <math.h>
-
-#if defined(__HTM__) || __ARCH__ >= 9
-
-long double sqrtl(long double x)
-{
-	__asm__ ("sqxbr %0, %1" : "=f"(x) : "f"(x));
-	return x;
-}
-
-#else
-
-#include "../sqrtl.c"
-
-#endif
diff --git a/lib/libc/musl/src/math/sqrt.c b/lib/libc/musl/src/math/sqrt.c
deleted file mode 100644
index 5ba26559621357018857a49e40b5745aaca4cc51..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/sqrt.c
+++ /dev/null
@@ -1,158 +0,0 @@
-#include <stdint.h>
-#include <math.h>
-#include "libm.h"
-#include "sqrt_data.h"
-
-#define FENV_SUPPORT 1
-
-/* returns a*b*2^-32 - e, with error 0 <= e < 1.  */
-static inline uint32_t mul32(uint32_t a, uint32_t b)
-{
-	return (uint64_t)a*b >> 32;
-}
-
-/* returns a*b*2^-64 - e, with error 0 <= e < 3.  */
-static inline uint64_t mul64(uint64_t a, uint64_t b)
-{
-	uint64_t ahi = a>>32;
-	uint64_t alo = a&0xffffffff;
-	uint64_t bhi = b>>32;
-	uint64_t blo = b&0xffffffff;
-	return ahi*bhi + (ahi*blo >> 32) + (alo*bhi >> 32);
-}
-
-double sqrt(double x)
-{
-	uint64_t ix, top, m;
-
-	/* special case handling.  */
-	ix = asuint64(x);
-	top = ix >> 52;
-	if (predict_false(top - 0x001 >= 0x7ff - 0x001)) {
-		/* x < 0x1p-1022 or inf or nan.  */
-		if (ix * 2 == 0)
-			return x;
-		if (ix == 0x7ff0000000000000)
-			return x;
-		if (ix > 0x7ff0000000000000)
-			return __math_invalid(x);
-		/* x is subnormal, normalize it.  */
-		ix = asuint64(x * 0x1p52);
-		top = ix >> 52;
-		top -= 52;
-	}
-
-	/* argument reduction:
-	   x = 4^e m; with integer e, and m in [1, 4)
-	   m: fixed point representation [2.62]
-	   2^e is the exponent part of the result.  */
-	int even = top & 1;
-	m = (ix << 11) | 0x8000000000000000;
-	if (even) m >>= 1;
-	top = (top + 0x3ff) >> 1;
-
-	/* approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
-
-	   initial estimate:
-	   7bit table lookup (1bit exponent and 6bit significand).
-
-	   iterative approximation:
-	   using 2 goldschmidt iterations with 32bit int arithmetics
-	   and a final iteration with 64bit int arithmetics.
-
-	   details:
-
-	   the relative error (e = r0 sqrt(m)-1) of a linear estimate
-	   (r0 = a m + b) is |e| < 0.085955 ~ 0x1.6p-4 at best,
-	   a table lookup is faster and needs one less iteration
-	   6 bit lookup table (128b) gives |e| < 0x1.f9p-8
-	   7 bit lookup table (256b) gives |e| < 0x1.fdp-9
-	   for single and double prec 6bit is enough but for quad
-	   prec 7bit is needed (or modified iterations). to avoid
-	   one more iteration >=13bit table would be needed (16k).
-
-	   a newton-raphson iteration for r is
-	     w = r*r
-	     u = 3 - m*w
-	     r = r*u/2
-	   can use a goldschmidt iteration for s at the end or
-	     s = m*r
-
-	   first goldschmidt iteration is
-	     s = m*r
-	     u = 3 - s*r
-	     r = r*u/2
-	     s = s*u/2
-	   next goldschmidt iteration is
-	     u = 3 - s*r
-	     r = r*u/2
-	     s = s*u/2
-	   and at the end r is not computed only s.
-
-	   they use the same amount of operations and converge at the
-	   same quadratic rate, i.e. if
-	     r1 sqrt(m) - 1 = e, then
-	     r2 sqrt(m) - 1 = -3/2 e^2 - 1/2 e^3
-	   the advantage of goldschmidt is that the mul for s and r
-	   are independent (computed in parallel), however it is not
-	   "self synchronizing": it only uses the input m in the
-	   first iteration so rounding errors accumulate. at the end
-	   or when switching to larger precision arithmetics rounding
-	   errors dominate so the first iteration should be used.
-
-	   the fixed point representations are
-	     m: 2.30 r: 0.32, s: 2.30, d: 2.30, u: 2.30, three: 2.30
-	   and after switching to 64 bit
-	     m: 2.62 r: 0.64, s: 2.62, d: 2.62, u: 2.62, three: 2.62  */
-
-	static const uint64_t three = 0xc0000000;
-	uint64_t r, s, d, u, i;
-
-	i = (ix >> 46) % 128;
-	r = (uint32_t)__rsqrt_tab[i] << 16;
-	/* |r sqrt(m) - 1| < 0x1.fdp-9 */
-	s = mul32(m>>32, r);
-	/* |s/sqrt(m) - 1| < 0x1.fdp-9 */
-	d = mul32(s, r);
-	u = three - d;
-	r = mul32(r, u) << 1;
-	/* |r sqrt(m) - 1| < 0x1.7bp-16 */
-	s = mul32(s, u) << 1;
-	/* |s/sqrt(m) - 1| < 0x1.7bp-16 */
-	d = mul32(s, r);
-	u = three - d;
-	r = mul32(r, u) << 1;
-	/* |r sqrt(m) - 1| < 0x1.3704p-29 (measured worst-case) */
-	r = r << 32;
-	s = mul64(m, r);
-	d = mul64(s, r);
-	u = (three<<32) - d;
-	s = mul64(s, u);  /* repr: 3.61 */
-	/* -0x1p-57 < s - sqrt(m) < 0x1.8001p-61 */
-	s = (s - 2) >> 9; /* repr: 12.52 */
-	/* -0x1.09p-52 < s - sqrt(m) < -0x1.fffcp-63 */
-
-	/* s < sqrt(m) < s + 0x1.09p-52,
-	   compute nearest rounded result:
-	   the nearest result to 52 bits is either s or s+0x1p-52,
-	   we can decide by comparing (2^52 s + 0.5)^2 to 2^104 m.  */
-	uint64_t d0, d1, d2;
-	double y, t;
-	d0 = (m << 42) - s*s;
-	d1 = s - d0;
-	d2 = d1 + s + 1;
-	s += d1 >> 63;
-	s &= 0x000fffffffffffff;
-	s |= top << 52;
-	y = asdouble(s);
-	if (FENV_SUPPORT) {
-		/* handle rounding modes and inexact exception:
-		   only (s+1)^2 == 2^42 m case is exact otherwise
-		   add a tiny value to cause the fenv effects.  */
-		uint64_t tiny = predict_false(d2==0) ? 0 : 0x0010000000000000;
-		tiny |= (d1^d2) & 0x8000000000000000;
-		t = asdouble(tiny);
-		y = eval_as_double(y + t);
-	}
-	return y;
-}
diff --git a/lib/libc/musl/src/math/sqrt_data.c b/lib/libc/musl/src/math/sqrt_data.c
deleted file mode 100644
index 61bc22f4309586e220da450aaff948a561071d54..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/sqrt_data.c
+++ /dev/null
@@ -1,19 +0,0 @@
-#include "sqrt_data.h"
-const uint16_t __rsqrt_tab[128] = {
-0xb451,0xb2f0,0xb196,0xb044,0xaef9,0xadb6,0xac79,0xab43,
-0xaa14,0xa8eb,0xa7c8,0xa6aa,0xa592,0xa480,0xa373,0xa26b,
-0xa168,0xa06a,0x9f70,0x9e7b,0x9d8a,0x9c9d,0x9bb5,0x9ad1,
-0x99f0,0x9913,0x983a,0x9765,0x9693,0x95c4,0x94f8,0x9430,
-0x936b,0x92a9,0x91ea,0x912e,0x9075,0x8fbe,0x8f0a,0x8e59,
-0x8daa,0x8cfe,0x8c54,0x8bac,0x8b07,0x8a64,0x89c4,0x8925,
-0x8889,0x87ee,0x8756,0x86c0,0x862b,0x8599,0x8508,0x8479,
-0x83ec,0x8361,0x82d8,0x8250,0x81c9,0x8145,0x80c2,0x8040,
-0xff02,0xfd0e,0xfb25,0xf947,0xf773,0xf5aa,0xf3ea,0xf234,
-0xf087,0xeee3,0xed47,0xebb3,0xea27,0xe8a3,0xe727,0xe5b2,
-0xe443,0xe2dc,0xe17a,0xe020,0xdecb,0xdd7d,0xdc34,0xdaf1,
-0xd9b3,0xd87b,0xd748,0xd61a,0xd4f1,0xd3cd,0xd2ad,0xd192,
-0xd07b,0xcf69,0xce5b,0xcd51,0xcc4a,0xcb48,0xca4a,0xc94f,
-0xc858,0xc764,0xc674,0xc587,0xc49d,0xc3b7,0xc2d4,0xc1f4,
-0xc116,0xc03c,0xbf65,0xbe90,0xbdbe,0xbcef,0xbc23,0xbb59,
-0xba91,0xb9cc,0xb90a,0xb84a,0xb78c,0xb6d0,0xb617,0xb560,
-};
diff --git a/lib/libc/musl/src/math/sqrt_data.h b/lib/libc/musl/src/math/sqrt_data.h
deleted file mode 100644
index 260c7f9c292beb48fd7a0985a58a7740d38c0d5e..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/sqrt_data.h
+++ /dev/null
@@ -1,13 +0,0 @@
-#ifndef _SQRT_DATA_H
-#define _SQRT_DATA_H
-
-#include <features.h>
-#include <stdint.h>
-
-/* if x in [1,2): i = (int)(64*x);
-   if x in [2,4): i = (int)(32*x-64);
-   __rsqrt_tab[i]*2^-16 is estimating 1/sqrt(x) with small relative error:
-   |__rsqrt_tab[i]*0x1p-16*sqrt(x) - 1| < -0x1.fdp-9 < 2^-8 */
-extern hidden const uint16_t __rsqrt_tab[128];
-
-#endif
diff --git a/lib/libc/musl/src/math/sqrtf.c b/lib/libc/musl/src/math/sqrtf.c
deleted file mode 100644
index 740d81cbab421707a090d2475523807fd27b1337..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/sqrtf.c
+++ /dev/null
@@ -1,83 +0,0 @@
-#include <stdint.h>
-#include <math.h>
-#include "libm.h"
-#include "sqrt_data.h"
-
-#define FENV_SUPPORT 1
-
-static inline uint32_t mul32(uint32_t a, uint32_t b)
-{
-	return (uint64_t)a*b >> 32;
-}
-
-/* see sqrt.c for more detailed comments.  */
-
-float sqrtf(float x)
-{
-	uint32_t ix, m, m1, m0, even, ey;
-
-	ix = asuint(x);
-	if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) {
-		/* x < 0x1p-126 or inf or nan.  */
-		if (ix * 2 == 0)
-			return x;
-		if (ix == 0x7f800000)
-			return x;
-		if (ix > 0x7f800000)
-			return __math_invalidf(x);
-		/* x is subnormal, normalize it.  */
-		ix = asuint(x * 0x1p23f);
-		ix -= 23 << 23;
-	}
-
-	/* x = 4^e m; with int e and m in [1, 4).  */
-	even = ix & 0x00800000;
-	m1 = (ix << 8) | 0x80000000;
-	m0 = (ix << 7) & 0x7fffffff;
-	m = even ? m0 : m1;
-
-	/* 2^e is the exponent part of the return value.  */
-	ey = ix >> 1;
-	ey += 0x3f800000 >> 1;
-	ey &= 0x7f800000;
-
-	/* compute r ~ 1/sqrt(m), s ~ sqrt(m) with 2 goldschmidt iterations.  */
-	static const uint32_t three = 0xc0000000;
-	uint32_t r, s, d, u, i;
-	i = (ix >> 17) % 128;
-	r = (uint32_t)__rsqrt_tab[i] << 16;
-	/* |r*sqrt(m) - 1| < 0x1p-8 */
-	s = mul32(m, r);
-	/* |s/sqrt(m) - 1| < 0x1p-8 */
-	d = mul32(s, r);
-	u = three - d;
-	r = mul32(r, u) << 1;
-	/* |r*sqrt(m) - 1| < 0x1.7bp-16 */
-	s = mul32(s, u) << 1;
-	/* |s/sqrt(m) - 1| < 0x1.7bp-16 */
-	d = mul32(s, r);
-	u = three - d;
-	s = mul32(s, u);
-	/* -0x1.03p-28 < s/sqrt(m) - 1 < 0x1.fp-31 */
-	s = (s - 1)>>6;
-	/* s < sqrt(m) < s + 0x1.08p-23 */
-
-	/* compute nearest rounded result.  */
-	uint32_t d0, d1, d2;
-	float y, t;
-	d0 = (m << 16) - s*s;
-	d1 = s - d0;
-	d2 = d1 + s + 1;
-	s += d1 >> 31;
-	s &= 0x007fffff;
-	s |= ey;
-	y = asfloat(s);
-	if (FENV_SUPPORT) {
-		/* handle rounding and inexact exception. */
-		uint32_t tiny = predict_false(d2==0) ? 0 : 0x01000000;
-		tiny |= (d1^d2) & 0x80000000;
-		t = asfloat(tiny);
-		y = eval_as_float(y + t);
-	}
-	return y;
-}
diff --git a/lib/libc/musl/src/math/sqrtl.c b/lib/libc/musl/src/math/sqrtl.c
deleted file mode 100644
index a231b3f2016b391c778c9cb5da04913bd73493c9..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/sqrtl.c
+++ /dev/null
@@ -1,259 +0,0 @@
-#include <stdint.h>
-#include <math.h>
-#include <float.h>
-#include "libm.h"
-
-#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
-long double sqrtl(long double x)
-{
-	return sqrt(x);
-}
-#elif (LDBL_MANT_DIG == 113 || LDBL_MANT_DIG == 64) && LDBL_MAX_EXP == 16384
-#include "sqrt_data.h"
-
-#define FENV_SUPPORT 1
-
-typedef struct {
-	uint64_t hi;
-	uint64_t lo;
-} u128;
-
-/* top: 16 bit sign+exponent, x: significand.  */
-static inline long double mkldbl(uint64_t top, u128 x)
-{
-	union ldshape u;
-#if LDBL_MANT_DIG == 113
-	u.i2.hi = x.hi;
-	u.i2.lo = x.lo;
-	u.i2.hi &= 0x0000ffffffffffff;
-	u.i2.hi |= top << 48;
-#elif LDBL_MANT_DIG == 64
-	u.i.se = top;
-	u.i.m = x.lo;
-	/* force the top bit on non-zero (and non-subnormal) results.  */
-	if (top & 0x7fff)
-		u.i.m |= 0x8000000000000000;
-#endif
-	return u.f;
-}
-
-/* return: top 16 bit is sign+exp and following bits are the significand.  */
-static inline u128 asu128(long double x)
-{
-	union ldshape u = {.f=x};
-	u128 r;
-#if LDBL_MANT_DIG == 113
-	r.hi = u.i2.hi;
-	r.lo = u.i2.lo;
-#elif LDBL_MANT_DIG == 64
-	r.lo = u.i.m<<49;
-	/* ignore the top bit: pseudo numbers are not handled. */
-	r.hi = u.i.m>>15;
-	r.hi &= 0x0000ffffffffffff;
-	r.hi |= (uint64_t)u.i.se << 48;
-#endif
-	return r;
-}
-
-/* returns a*b*2^-32 - e, with error 0 <= e < 1.  */
-static inline uint32_t mul32(uint32_t a, uint32_t b)
-{
-	return (uint64_t)a*b >> 32;
-}
-
-/* returns a*b*2^-64 - e, with error 0 <= e < 3.  */
-static inline uint64_t mul64(uint64_t a, uint64_t b)
-{
-	uint64_t ahi = a>>32;
-	uint64_t alo = a&0xffffffff;
-	uint64_t bhi = b>>32;
-	uint64_t blo = b&0xffffffff;
-	return ahi*bhi + (ahi*blo >> 32) + (alo*bhi >> 32);
-}
-
-static inline u128 add64(u128 a, uint64_t b)
-{
-	u128 r;
-	r.lo = a.lo + b;
-	r.hi = a.hi;
-	if (r.lo < a.lo)
-		r.hi++;
-	return r;
-}
-
-static inline u128 add128(u128 a, u128 b)
-{
-	u128 r;
-	r.lo = a.lo + b.lo;
-	r.hi = a.hi + b.hi;
-	if (r.lo < a.lo)
-		r.hi++;
-	return r;
-}
-
-static inline u128 sub64(u128 a, uint64_t b)
-{
-	u128 r;
-	r.lo = a.lo - b;
-	r.hi = a.hi;
-	if (a.lo < b)
-		r.hi--;
-	return r;
-}
-
-static inline u128 sub128(u128 a, u128 b)
-{
-	u128 r;
-	r.lo = a.lo - b.lo;
-	r.hi = a.hi - b.hi;
-	if (a.lo < b.lo)
-		r.hi--;
-	return r;
-}
-
-/* a<<n, 0 <= n <= 127 */
-static inline u128 lsh(u128 a, int n)
-{
-	if (n == 0)
-		return a;
-	if (n >= 64) {
-		a.hi = a.lo<<(n-64);
-		a.lo = 0;
-	} else {
-		a.hi = (a.hi<<n) | (a.lo>>(64-n));
-		a.lo = a.lo<<n;
-	}
-	return a;
-}
-
-/* a>>n, 0 <= n <= 127 */
-static inline u128 rsh(u128 a, int n)
-{
-	if (n == 0)
-		return a;
-	if (n >= 64) {
-		a.lo = a.hi>>(n-64);
-		a.hi = 0;
-	} else {
-		a.lo = (a.lo>>n) | (a.hi<<(64-n));
-		a.hi = a.hi>>n;
-	}
-	return a;
-}
-
-/* returns a*b exactly.  */
-static inline u128 mul64_128(uint64_t a, uint64_t b)
-{
-	u128 r;
-	uint64_t ahi = a>>32;
-	uint64_t alo = a&0xffffffff;
-	uint64_t bhi = b>>32;
-	uint64_t blo = b&0xffffffff;
-	uint64_t lo1 = ((ahi*blo)&0xffffffff) + ((alo*bhi)&0xffffffff) + (alo*blo>>32);
-	uint64_t lo2 = (alo*blo)&0xffffffff;
-	r.hi = ahi*bhi + (ahi*blo>>32) + (alo*bhi>>32) + (lo1>>32);
-	r.lo = (lo1<<32) + lo2;
-	return r;
-}
-
-/* returns a*b*2^-128 - e, with error 0 <= e < 7.  */
-static inline u128 mul128(u128 a, u128 b)
-{
-	u128 hi = mul64_128(a.hi, b.hi);
-	uint64_t m1 = mul64(a.hi, b.lo);
-	uint64_t m2 = mul64(a.lo, b.hi);
-	return add64(add64(hi, m1), m2);
-}
-
-/* returns a*b % 2^128.  */
-static inline u128 mul128_tail(u128 a, u128 b)
-{
-	u128 lo = mul64_128(a.lo, b.lo);
-	lo.hi += a.hi*b.lo + a.lo*b.hi;
-	return lo;
-}
-
-
-/* see sqrt.c for detailed comments.  */
-
-long double sqrtl(long double x)
-{
-	u128 ix, ml;
-	uint64_t top;
-
-	ix = asu128(x);
-	top = ix.hi >> 48;
-	if (predict_false(top - 0x0001 >= 0x7fff - 0x0001)) {
-		/* x < 0x1p-16382 or inf or nan.  */
-		if (2*ix.hi == 0 && ix.lo == 0)
-			return x;
-		if (ix.hi == 0x7fff000000000000 && ix.lo == 0)
-			return x;
-		if (top >= 0x7fff)
-			return __math_invalidl(x);
-		/* x is subnormal, normalize it.  */
-		ix = asu128(x * 0x1p112);
-		top = ix.hi >> 48;
-		top -= 112;
-	}
-
-	/* x = 4^e m; with int e and m in [1, 4) */
-	int even = top & 1;
-	ml = lsh(ix, 15);
-	ml.hi |= 0x8000000000000000;
-	if (even) ml = rsh(ml, 1);
-	top = (top + 0x3fff) >> 1;
-
-	/* r ~ 1/sqrt(m) */
-	const uint64_t three = 0xc0000000;
-	uint64_t r, s, d, u, i;
-	i = (ix.hi >> 42) % 128;
-	r = (uint32_t)__rsqrt_tab[i] << 16;
-	/* |r sqrt(m) - 1| < 0x1p-8 */
-	s = mul32(ml.hi>>32, r);
-	d = mul32(s, r);
-	u = three - d;
-	r = mul32(u, r) << 1;
-	/* |r sqrt(m) - 1| < 0x1.7bp-16, switch to 64bit */
-	r = r<<32;
-	s = mul64(ml.hi, r);
-	d = mul64(s, r);
-	u = (three<<32) - d;
-	r = mul64(u, r) << 1;
-	/* |r sqrt(m) - 1| < 0x1.a5p-31 */
-	s = mul64(u, s) << 1;
-	d = mul64(s, r);
-	u = (three<<32) - d;
-	r = mul64(u, r) << 1;
-	/* |r sqrt(m) - 1| < 0x1.c001p-59, switch to 128bit */
-
-	const u128 threel = {.hi=three<<32, .lo=0};
-	u128 rl, sl, dl, ul;
-	rl.hi = r;
-	rl.lo = 0;
-	sl = mul128(ml, rl);
-	dl = mul128(sl, rl);
-	ul = sub128(threel, dl);
-	sl = mul128(ul, sl); /* repr: 3.125 */
-	/* -0x1p-116 < s - sqrt(m) < 0x3.8001p-125 */
-	sl = rsh(sub64(sl, 4), 125-(LDBL_MANT_DIG-1));
-	/* s < sqrt(m) < s + 1 ULP + tiny */
-
-	long double y;
-	u128 d2, d1, d0;
-	d0 = sub128(lsh(ml, 2*(LDBL_MANT_DIG-1)-126), mul128_tail(sl,sl));
-	d1 = sub128(sl, d0);
-	d2 = add128(add64(sl, 1), d1);
-	sl = add64(sl, d1.hi >> 63);
-	y = mkldbl(top, sl);
-	if (FENV_SUPPORT) {
-		/* handle rounding modes and inexact exception.  */
-		top = predict_false((d2.hi|d2.lo)==0) ? 0 : 1;
-		top |= ((d1.hi^d2.hi)&0x8000000000000000) >> 48;
-		y += mkldbl(top, (u128){0});
-	}
-	return y;
-}
-#else
-#error unsupported long double format
-#endif
diff --git a/lib/libc/musl/src/math/x32/sqrt.s b/lib/libc/musl/src/math/x32/sqrt.s
deleted file mode 100644
index d3c609f9f8d540beac7cc5a436aab807ea511b1d..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/x32/sqrt.s
+++ /dev/null
@@ -1,4 +0,0 @@
-.global sqrt
-.type sqrt,@function
-sqrt:	sqrtsd %xmm0, %xmm0
-	ret
diff --git a/lib/libc/musl/src/math/x32/sqrtf.s b/lib/libc/musl/src/math/x32/sqrtf.s
deleted file mode 100644
index eec48c609412ff88ad576c56df8d4fe62c5ea448..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/x32/sqrtf.s
+++ /dev/null
@@ -1,4 +0,0 @@
-.global sqrtf
-.type sqrtf,@function
-sqrtf:  sqrtss %xmm0, %xmm0
-	ret
diff --git a/lib/libc/musl/src/math/x32/sqrtl.s b/lib/libc/musl/src/math/x32/sqrtl.s
deleted file mode 100644
index 8d70856ec035d79a2f5e38ba10d9bf8abcfd3983..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/x32/sqrtl.s
+++ /dev/null
@@ -1,5 +0,0 @@
-.global sqrtl
-.type sqrtl,@function
-sqrtl:	fldt 8(%esp)
-	fsqrt
-	ret
diff --git a/lib/libc/musl/src/math/x86_64/sqrt.c b/lib/libc/musl/src/math/x86_64/sqrt.c
deleted file mode 100644
index 657e09e3b4d44b1cc12bf3300d7e5c5dc06c11ec..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/x86_64/sqrt.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <math.h>
-
-double sqrt(double x)
-{
-	__asm__ ("sqrtsd %1, %0" : "=x"(x) : "x"(x));
-	return x;
-}
diff --git a/lib/libc/musl/src/math/x86_64/sqrtf.c b/lib/libc/musl/src/math/x86_64/sqrtf.c
deleted file mode 100644
index 720baec601289ed560a73b7a071bfb5b61e6a49a..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/x86_64/sqrtf.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <math.h>
-
-float sqrtf(float x)
-{
-	__asm__ ("sqrtss %1, %0" : "=x"(x) : "x"(x));
-	return x;
-}
diff --git a/lib/libc/musl/src/math/x86_64/sqrtl.c b/lib/libc/musl/src/math/x86_64/sqrtl.c
deleted file mode 100644
index 864cfcc4f66eb5c0b8c87f38299fc727eb6e7eac..0000000000000000000000000000000000000000
--- a/lib/libc/musl/src/math/x86_64/sqrtl.c
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <math.h>
-
-long double sqrtl(long double x)
-{
-	__asm__ ("fsqrt" : "+t"(x));
-	return x;
-}
diff --git a/src/libs/mingw.zig b/src/libs/mingw.zig
index ded4db0c3e8f2690479f406794de59493ec92eb5..eb9eb25a34f573a243ba0cf9698ab3970e8b2c45 100644
--- a/src/libs/mingw.zig
+++ b/src/libs/mingw.zig
@@ -630,7 +630,6 @@ const mingw32_generic_src = [_][]const u8{
     "math" ++ path.sep_str ++ "signbitl.c",
     "math" ++ path.sep_str ++ "signgam.c",
     "math" ++ path.sep_str ++ "sinhl.c",
-    "math" ++ path.sep_str ++ "sqrtl.c",
     "math" ++ path.sep_str ++ "tanhl.c",
     "misc" ++ path.sep_str ++ "alarm.c",
     "misc" ++ path.sep_str ++ "btowc.c",
@@ -995,7 +994,6 @@ const mingw32_x86_32_src = [_][]const u8{
     "math" ++ path.sep_str ++ "modff.c",
     "math" ++ path.sep_str ++ "powf.c",
     "math" ++ path.sep_str ++ "sinhf.c",
-    "math" ++ path.sep_str ++ "sqrtf.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",
diff --git a/src/libs/musl.zig b/src/libs/musl.zig
index 3393cecbb041159e92a73b7aeb173d3619a1f886..302d22b062aae5f960885ffe8f3107114db25ff4 100644
--- a/src/libs/musl.zig
+++ b/src/libs/musl.zig
@@ -826,8 +826,6 @@ const src_files = [_][]const u8{
     "musl/src/math/aarch64/rintf.c",
     "musl/src/math/aarch64/round.c",
     "musl/src/math/aarch64/roundf.c",
-    "musl/src/math/aarch64/sqrt.c",
-    "musl/src/math/aarch64/sqrtf.c",
     "musl/src/math/acos.c",
     "musl/src/math/acosf.c",
     "musl/src/math/acosh.c",
@@ -836,8 +834,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/arm/sqrt.c",
-    "musl/src/math/arm/sqrtf.c",
     "musl/src/math/asin.c",
     "musl/src/math/asinf.c",
     "musl/src/math/asinh.c",
@@ -944,9 +940,6 @@ const src_files = [_][]const u8{
     "musl/src/math/i386/scalbnf.s",
     "musl/src/math/i386/scalbnl.s",
     "musl/src/math/i386/scalbn.s",
-    "musl/src/math/i386/sqrt.c",
-    "musl/src/math/i386/sqrtf.c",
-    "musl/src/math/i386/sqrtl.c",
     "musl/src/math/ilogb.c",
     "musl/src/math/ilogbf.c",
     "musl/src/math/ilogbl.c",
@@ -986,7 +979,6 @@ const src_files = [_][]const u8{
     "musl/src/math/lround.c",
     "musl/src/math/lroundf.c",
     "musl/src/math/lroundl.c",
-    "musl/src/math/m68k/sqrtl.c",
     "musl/src/math/__math_divzero.c",
     "musl/src/math/__math_divzerof.c",
     "musl/src/math/__math_invalid.c",
@@ -998,8 +990,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/mips/sqrt.c",
-    "musl/src/math/mips/sqrtf.c",
     "musl/src/math/modf.c",
     "musl/src/math/modff.c",
     "musl/src/math/modfl.c",
@@ -1023,12 +1013,8 @@ const src_files = [_][]const u8{
     "musl/src/math/powerpc64/lroundf.c",
     "musl/src/math/powerpc64/round.c",
     "musl/src/math/powerpc64/roundf.c",
-    "musl/src/math/powerpc64/sqrt.c",
-    "musl/src/math/powerpc64/sqrtf.c",
     "musl/src/math/powerpc/fma.c",
     "musl/src/math/powerpc/fmaf.c",
-    "musl/src/math/powerpc/sqrt.c",
-    "musl/src/math/powerpc/sqrtf.c",
     "musl/src/math/powf.c",
     "musl/src/math/powf_data.c",
     "musl/src/math/powl.c",
@@ -1047,12 +1033,8 @@ const src_files = [_][]const u8{
     "musl/src/math/rintl.c",
     "musl/src/math/riscv32/fma.c",
     "musl/src/math/riscv32/fmaf.c",
-    "musl/src/math/riscv32/sqrt.c",
-    "musl/src/math/riscv32/sqrtf.c",
     "musl/src/math/riscv64/fma.c",
     "musl/src/math/riscv64/fmaf.c",
-    "musl/src/math/riscv64/sqrt.c",
-    "musl/src/math/riscv64/sqrtf.c",
     "musl/src/math/round.c",
     "musl/src/math/roundf.c",
     "musl/src/math/roundl.c",
@@ -1067,9 +1049,6 @@ const src_files = [_][]const u8{
     "musl/src/math/s390x/round.c",
     "musl/src/math/s390x/roundf.c",
     "musl/src/math/s390x/roundl.c",
-    "musl/src/math/s390x/sqrt.c",
-    "musl/src/math/s390x/sqrtf.c",
-    "musl/src/math/s390x/sqrtl.c",
     "musl/src/math/scalb.c",
     "musl/src/math/scalbf.c",
     "musl/src/math/scalbln.c",
@@ -1092,10 +1071,6 @@ const src_files = [_][]const u8{
     "musl/src/math/sinhl.c",
     "musl/src/math/__sinl.c",
     "musl/src/math/sinl.c",
-    "musl/src/math/sqrt.c",
-    "musl/src/math/sqrt_data.c",
-    "musl/src/math/sqrtf.c",
-    "musl/src/math/sqrtl.c",
     "musl/src/math/__tan.c",
     "musl/src/math/__tandf.c",
     "musl/src/math/tanh.c",
@@ -1128,9 +1103,6 @@ const src_files = [_][]const u8{
     "musl/src/math/x32/lrint.s",
     "musl/src/math/x32/remainderl.s",
     "musl/src/math/x32/rintl.s",
-    "musl/src/math/x32/sqrtf.s",
-    "musl/src/math/x32/sqrtl.s",
-    "musl/src/math/x32/sqrt.s",
     "musl/src/math/x86_64/acosl.s",
     "musl/src/math/x86_64/asinl.s",
     "musl/src/math/x86_64/atan2l.s",
@@ -1154,9 +1126,6 @@ const src_files = [_][]const u8{
     "musl/src/math/x86_64/remainderl.c",
     "musl/src/math/x86_64/remquol.c",
     "musl/src/math/x86_64/rintl.c",
-    "musl/src/math/x86_64/sqrt.c",
-    "musl/src/math/x86_64/sqrtf.c",
-    "musl/src/math/x86_64/sqrtl.c",
     "musl/src/misc/a64l.c",
     "musl/src/misc/basename.c",
     "musl/src/misc/dirname.c",
diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig
index 1d002e335313b65c4c40c119ecd2a1abb46d49d9..161b60d55271e0112e1aa28d05c8f33642a53937 100644
--- a/src/libs/wasi_libc.zig
+++ b/src/libs/wasi_libc.zig
@@ -834,8 +834,6 @@ const libc_top_half_src_files = [_][]const u8{
     "musl/src/math/sinhl.c",
     "musl/src/math/__sinl.c",
     "musl/src/math/sinl.c",
-    "musl/src/math/sqrt_data.c",
-    "musl/src/math/sqrtl.c",
     "musl/src/math/__tan.c",
     "musl/src/math/__tandf.c",
     "musl/src/math/tanh.c",