diff --git a/src/math/mod.rs b/src/math/mod.rs
index 48b400a92..3c7d05111 100644
--- a/src/math/mod.rs
+++ b/src/math/mod.rs
@@ -344,3 +344,93 @@ fn with_set_low_word(f: f64, lo: u32) -> f64 {
 fn combine_words(hi: u32, lo: u32) -> f64 {
     f64::from_bits((hi as u64) << 32 | lo as u64)
 }
+
+mod fdlibm {
+    #![allow(dead_code)] // FIXME: remove it after other dependent codes merged
+    #![allow(non_snake_case)]
+
+    /* @(#)fdlibm.h 5.1 93/09/24 */
+    /*
+     * ====================================================
+     * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+     *
+     * Developed at SunPro, a Sun Microsystems, Inc. business.
+     * Permission to use, copy, modify, and distribute this
+     * software is freely granted, provided that this notice
+     * is preserved.
+     * ====================================================
+     */
+
+    /* Most routines need to check whether a float is finite, infinite, or not a
+    number, and many need to know whether the result of an operation will
+    overflow.  These conditions depend on whether the largest exponent is
+    used for NaNs & infinities, or whether it's used for finite numbers. */
+
+    /// True if a positive float with bitmask X is finite.
+    #[inline]
+    pub(crate) fn FLT_UWORD_IS_FINITE(x: u32) -> bool {
+        x < INFINITE
+    }
+
+    /// True if a positive float with bitmask X is not a number.
+    #[inline]
+    pub(crate) fn FLT_UWORD_IS_NAN(x: u32) -> bool {
+        x > INFINITE
+    }
+
+    /// True if a positive float with bitmask X is +infinity.
+    #[inline]
+    pub(crate) fn FLT_UWORD_IS_INFINITE(x: u32) -> bool {
+        x == INFINITE
+    }
+
+    const INFINITE: u32 = 0x7f80_0000;
+
+    /// The bitmask of FLT_MAX.
+    pub(crate) const FLT_UWORD_MAX: u32 = 0x7f7f_ffff;
+    /// The bitmask of FLT_MAX/2.
+    pub(crate) const FLT_UWORD_HALF_MAX: u32 = FLT_UWORD_MAX - (1 << 23);
+    /// The bitmask of the largest finite exponent (129 if the largest
+    /// exponent is used for finite numbers, 128 otherwise).
+    pub(crate) const FLT_UWORD_EXP_MAX: u32 = 0x4300_0000;
+    /// The bitmask of log(FLT_MAX), rounded down.  This value is the largest
+    /// input that can be passed to exp() without producing overflow.
+    pub(crate) const FLT_UWORD_LOG_MAX: u32 = 0x42b1_7217;
+    /// The bitmask of log(2*FLT_MAX), rounded down.  This value is the
+    /// largest input than can be passed to cosh() without producing
+    /// overflow.
+    pub(crate) const FLT_UWORD_LOG_2MAX: u32 = 0x42b2_d4fc;
+    /// The largest biased exponent that can be used for finite numbers
+    /// (255 if the largest exponent is used for finite numbers, 254
+    /// otherwise)
+    pub(crate) const FLT_LARGEST_EXP: u32 = FLT_UWORD_MAX >> 23;
+    pub(crate) const HUGE: f32 = 3.402_823_466_385_288_60e+38;
+
+    /* Many routines check for zero and subnormal numbers.  Such things depend
+    on whether the target supports denormals or not */
+
+    /// True if a positive float with bitmask X is +0.  Without denormals,
+    /// any float with a zero exponent is a +0 representation.  With
+    /// denormals, the only +0 representation is a 0 bitmask.
+    #[inline]
+    pub(crate) fn FLT_UWORD_IS_ZERO(x: u32) -> bool {
+        x == 0
+    }
+
+    /// True if a non-zero positive float with bitmask X is subnormal.
+    /// (Routines should check for zeros first.)
+    #[inline]
+    pub(crate) fn FLT_UWORD_IS_SUBNORMAL(x: u32) -> bool {
+        x < 0x0080_0000
+    }
+
+    /// The bitmask of the smallest float above +0.  Call this number REAL_FLT_MIN...
+    pub(crate) const FLT_UWORD_MIN: u32 = 0x0000_0001;
+    /// The bitmask of the float representation of REAL_FLT_MIN's exponent.
+    pub(crate) const FLT_UWORD_EXP_MIN: u32 = 0x4316_0000;
+    /// The bitmask of |log(REAL_FLT_MIN)|, rounding down.
+    pub(crate) const FLT_UWORD_LOG_MIN: u32 = 0x42cf_f1b5;
+    /// REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
+    /// -22 if they are).
+    pub(crate) const FLT_SMALLEST_EXP: i32 = -22;
+}
diff --git a/src/math/sqrtf.rs b/src/math/sqrtf.rs
index b9365c617..9a4d50205 100644
--- a/src/math/sqrtf.rs
+++ b/src/math/sqrtf.rs
@@ -1,7 +1,7 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */
-/*
+/* ef_sqrtf.c -- float version of e_sqrt.c.
  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
  */
+
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -13,65 +13,69 @@
  * ====================================================
  */
 
+const ONE: f32 = 1.0;
 const TINY: f32 = 1.0e-30;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn sqrtf(x: f32) -> f32 {
+pub extern "C" fn sqrtf(x: f32) -> f32 {
+    use super::fdlibm::{FLT_UWORD_IS_FINITE, FLT_UWORD_IS_SUBNORMAL, FLT_UWORD_IS_ZERO};
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f32.sqrt` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
         #[cfg(target_arch = "wasm32")] {
             return if x < 0.0 {
-                ::core::f32::NAN
+                core::f32::NAN
             } else {
-                unsafe { ::core::intrinsics::sqrtf32(x) }
+                unsafe { core::intrinsics::sqrtf32(x) }
             }
         }
     }
+
     let mut z: f32;
-    let sign: i32 = 0x80000000u32 as i32;
+
+    let mut r: u32;
+    let hx: u32;
+
     let mut ix: i32;
     let mut s: i32;
     let mut q: i32;
     let mut m: i32;
     let mut t: i32;
     let mut i: i32;
-    let mut r: u32;
 
     ix = x.to_bits() as i32;
+    hx = ix as u32 & 0x7fff_ffff;
 
     /* take care of Inf and NaN */
-    if (ix as u32 & 0x7f800000) == 0x7f800000 {
+    if !FLT_UWORD_IS_FINITE(hx) {
         return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
     }
 
-    /* take care of zero */
-    if ix <= 0 {
-        if (ix & !sign) == 0 {
-            return x; /* sqrt(+-0) = +-0 */
-        }
-        if ix < 0 {
-            return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
-        }
+    /* take care of zero and -ves */
+    if FLT_UWORD_IS_ZERO(hx) {
+        return x; /* sqrt(+-0) = +-0 */
+    }
+    if ix < 0 {
+        return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
     }
 
     /* normalize x */
     m = ix >> 23;
-    if m == 0 {
+    if FLT_UWORD_IS_SUBNORMAL(hx) {
         /* subnormal x */
         i = 0;
-        while ix & 0x00800000 == 0 {
+        while ix & 0x0080_0000 == 0 {
             ix <<= 1;
-            i = i + 1;
+            i += 1;
         }
         m -= i - 1;
     }
     m -= 127; /* unbias exponent */
-    ix = (ix & 0x007fffff) | 0x00800000;
+    ix = (ix & 0x007f_ffff) | 0x0080_0000;
+    /* odd m, double x to make it even */
     if m & 1 == 1 {
-        /* odd m, double x to make it even */
         ix += ix;
     }
     m >>= 1; /* m = [m/2] */
@@ -79,8 +83,8 @@ pub fn sqrtf(x: f32) -> f32 {
     /* generate sqrt(x) bit by bit */
     ix += ix;
     q = 0;
-    s = 0;
-    r = 0x01000000; /* r = moving bit from right to left */
+    s = 0; /* q = sqrt(x) */
+    r = 0x0100_0000; /* r = moving bit from right to left */
 
     while r != 0 {
         t = s + r as i32;
@@ -95,18 +99,17 @@ pub fn sqrtf(x: f32) -> f32 {
 
     /* use floating add to find out rounding direction */
     if ix != 0 {
-        z = 1.0 - TINY; /* raise inexact flag */
-        if z >= 1.0 {
-            z = 1.0 + TINY;
-            if z > 1.0 {
+        z = ONE - TINY; /* trigger inexact flag */
+        if z >= ONE {
+            z = ONE + TINY;
+            if z > ONE {
                 q += 2;
             } else {
                 q += q & 1;
             }
         }
     }
-
-    ix = (q >> 1) + 0x3f000000;
+    ix = (q >> 1) + 0x3f00_0000;
     ix += m << 23;
     f32::from_bits(ix as u32)
 }