diff --git a/library/compiler-builtins/libm/CHANGELOG.md b/library/compiler-builtins/libm/CHANGELOG.md
index 3a496527bc04..107813b08136 100644
--- a/library/compiler-builtins/libm/CHANGELOG.md
+++ b/library/compiler-builtins/libm/CHANGELOG.md
@@ -5,6 +5,19 @@ This project adheres to [Semantic Versioning](http://semver.org/).
 
 ## [Unreleased]
 
+### Added
+
+- atan2f
+- cos
+- coshf
+- fmaf
+- sin
+- sinh
+- sinhf
+- tan
+- tanh
+- tanhf
+
 ## [v0.1.1] - 2018-07-14
 
 ### Added
diff --git a/library/compiler-builtins/libm/ci/script.sh b/library/compiler-builtins/libm/ci/script.sh
index f2a294b488f5..cf37ac1ca64b 100644
--- a/library/compiler-builtins/libm/ci/script.sh
+++ b/library/compiler-builtins/libm/ci/script.sh
@@ -15,11 +15,6 @@ main() {
     # generate tests
     cargo run --package test-generator --target x86_64-unknown-linux-musl
 
-    if cargo fmt --version >/dev/null 2>&1; then
-        # nicer syntax error messages (if any)
-        cargo fmt
-    fi
-
     # run tests
     cross test --target $TARGET --release
 
diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs
index e7abce540af4..1837b9c1f9ed 100644
--- a/library/compiler-builtins/libm/src/lib.rs
+++ b/library/compiler-builtins/libm/src/lib.rs
@@ -14,18 +14,19 @@
 
 mod math;
 
-#[cfg(todo)]
 use core::{f32, f64};
 
 pub use math::*;
 
 /// Approximate equality with 1 ULP of tolerance
 #[doc(hidden)]
+#[inline]
 pub fn _eqf(a: u32, b: u32) -> bool {
     (a as i32).wrapping_sub(b as i32).abs() <= 1
 }
 
 #[doc(hidden)]
+#[inline]
 pub fn _eq(a: u64, b: u64) -> bool {
     (a as i64).wrapping_sub(b as i64).abs() <= 1
 }
@@ -33,7 +34,7 @@ pub fn _eq(a: u64, b: u64) -> bool {
 /// Math support for `f32`
 ///
 /// This trait is sealed and cannot be implemented outside of `libm`.
-pub trait F32Ext: private::Sealed {
+pub trait F32Ext: private::Sealed + Sized {
     fn floor(self) -> Self;
 
     fn ceil(self) -> Self;
@@ -44,21 +45,17 @@ pub trait F32Ext: private::Sealed {
 
     fn fdim(self, rhs: Self) -> Self;
 
-    #[cfg(todo)]
     fn fract(self) -> Self;
 
     fn abs(self) -> Self;
 
-    #[cfg(todo)]
-    fn signum(self) -> Self;
+    // NOTE depends on unstable intrinsics::copysignf32
+    // fn signum(self) -> Self;
 
-    #[cfg(todo)]
     fn mul_add(self, a: Self, b: Self) -> Self;
 
-    #[cfg(todo)]
     fn div_euc(self, rhs: Self) -> Self;
 
-    #[cfg(todo)]
     fn mod_euc(self, rhs: Self) -> Self;
 
     // NOTE depends on unstable intrinsics::powif32
@@ -98,9 +95,11 @@ pub trait F32Ext: private::Sealed {
 
     fn atan2(self, other: Self) -> Self;
 
-    #[cfg(todo)]
     #[inline]
-    fn sin_cos(self) -> (Self, Self) {
+    fn sin_cos(self) -> (Self, Self)
+    where
+        Self: Copy,
+    {
         (self.sin(), self.cos())
     }
 
@@ -114,13 +113,10 @@ pub trait F32Ext: private::Sealed {
 
     fn tanh(self) -> Self;
 
-    #[cfg(todo)]
     fn asinh(self) -> Self;
 
-    #[cfg(todo)]
     fn acosh(self) -> Self;
 
-    #[cfg(todo)]
     fn atanh(self) -> Self;
 }
 
@@ -150,7 +146,6 @@ impl F32Ext for f32 {
         fdimf(self, rhs)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn fract(self) -> Self {
         self - self.trunc()
@@ -161,13 +156,11 @@ impl F32Ext for f32 {
         fabsf(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn mul_add(self, a: Self, b: Self) -> Self {
         fmaf(self, a, b)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn div_euc(self, rhs: Self) -> Self {
         let q = (self / rhs).trunc();
@@ -177,7 +170,6 @@ impl F32Ext for f32 {
         q
     }
 
-    #[cfg(todo)]
     #[inline]
     fn mod_euc(self, rhs: f32) -> f32 {
         let r = self % rhs;
@@ -298,7 +290,6 @@ impl F32Ext for f32 {
         tanhf(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn asinh(self) -> Self {
         if self == f32::NEG_INFINITY {
@@ -308,7 +299,6 @@ impl F32Ext for f32 {
         }
     }
 
-    #[cfg(todo)]
     #[inline]
     fn acosh(self) -> Self {
         match self {
@@ -317,7 +307,6 @@ impl F32Ext for f32 {
         }
     }
 
-    #[cfg(todo)]
     #[inline]
     fn atanh(self) -> Self {
         0.5 * ((2.0 * self) / (1.0 - self)).ln_1p()
@@ -327,7 +316,7 @@ impl F32Ext for f32 {
 /// Math support for `f64`
 ///
 /// This trait is sealed and cannot be implemented outside of `libm`.
-pub trait F64Ext: private::Sealed {
+pub trait F64Ext: private::Sealed + Sized {
     fn floor(self) -> Self;
 
     fn ceil(self) -> Self;
@@ -338,20 +327,17 @@ pub trait F64Ext: private::Sealed {
 
     fn fdim(self, rhs: Self) -> Self;
 
-    #[cfg(todo)]
     fn fract(self) -> Self;
 
     fn abs(self) -> Self;
 
-    #[cfg(todo)]
-    fn signum(self) -> Self;
+    // NOTE depends on unstable intrinsics::copysignf64
+    // fn signum(self) -> Self;
 
     fn mul_add(self, a: Self, b: Self) -> Self;
 
-    #[cfg(todo)]
     fn div_euc(self, rhs: Self) -> Self;
 
-    #[cfg(todo)]
     fn mod_euc(self, rhs: Self) -> Self;
 
     // NOTE depends on unstable intrinsics::powif64
@@ -384,19 +370,19 @@ pub trait F64Ext: private::Sealed {
 
     fn tan(self) -> Self;
 
-    #[cfg(todo)]
     fn asin(self) -> Self;
 
     fn acos(self) -> Self;
 
-    #[cfg(todo)]
     fn atan(self) -> Self;
 
     fn atan2(self, other: Self) -> Self;
 
-    #[cfg(todo)]
     #[inline]
-    fn sin_cos(self) -> (Self, Self) {
+    fn sin_cos(self) -> (Self, Self)
+    where
+        Self: Copy,
+    {
         (self.sin(), self.cos())
     }
 
@@ -406,19 +392,14 @@ pub trait F64Ext: private::Sealed {
 
     fn sinh(self) -> Self;
 
-    #[cfg(todo)]
     fn cosh(self) -> Self;
 
-    #[cfg(todo)]
     fn tanh(self) -> Self;
 
-    #[cfg(todo)]
     fn asinh(self) -> Self;
 
-    #[cfg(todo)]
     fn acosh(self) -> Self;
 
-    #[cfg(todo)]
     fn atanh(self) -> Self;
 }
 
@@ -447,7 +428,7 @@ impl F64Ext for f64 {
     fn fdim(self, rhs: Self) -> Self {
         fdim(self, rhs)
     }
-    #[cfg(todo)]
+
     #[inline]
     fn fract(self) -> Self {
         self - self.trunc()
@@ -463,7 +444,6 @@ impl F64Ext for f64 {
         fma(self, a, b)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn div_euc(self, rhs: Self) -> Self {
         let q = (self / rhs).trunc();
@@ -473,9 +453,8 @@ impl F64Ext for f64 {
         q
     }
 
-    #[cfg(todo)]
     #[inline]
-    fn mod_euc(self, rhs: f32) -> f32 {
+    fn mod_euc(self, rhs: f64) -> f64 {
         let r = self % rhs;
         if r < 0.0 {
             r + rhs.abs()
@@ -550,7 +529,6 @@ impl F64Ext for f64 {
         tan(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn asin(self) -> Self {
         asin(self)
@@ -561,7 +539,6 @@ impl F64Ext for f64 {
         acos(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn atan(self) -> Self {
         atan(self)
@@ -587,19 +564,16 @@ impl F64Ext for f64 {
         sinh(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn cosh(self) -> Self {
         cosh(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn tanh(self) -> Self {
         tanh(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn asinh(self) -> Self {
         if self == f64::NEG_INFINITY {
@@ -609,7 +583,6 @@ impl F64Ext for f64 {
         }
     }
 
-    #[cfg(todo)]
     #[inline]
     fn acosh(self) -> Self {
         match self {
@@ -618,7 +591,6 @@ impl F64Ext for f64 {
         }
     }
 
-    #[cfg(todo)]
     #[inline]
     fn atanh(self) -> Self {
         0.5 * ((2.0 * self) / (1.0 - self)).ln_1p()
diff --git a/library/compiler-builtins/libm/src/math/acosf.rs b/library/compiler-builtins/libm/src/math/acosf.rs
index bbe29c17c131..469601caba4d 100644
--- a/library/compiler-builtins/libm/src/math/acosf.rs
+++ b/library/compiler-builtins/libm/src/math/acosf.rs
@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::sqrtf::sqrtf;
 
 const PIO2_HI: f32 = 1.5707962513e+00; /* 0x3fc90fda */
@@ -7,6 +22,7 @@ const P_S1: f32 = -4.2743422091e-02;
 const P_S2: f32 = -8.6563630030e-03;
 const Q_S1: f32 = -7.0662963390e-01;
 
+#[inline]
 fn r(z: f32) -> f32 {
     let p = z * (P_S0 + z * (P_S1 + z * P_S2));
     let q = 1. + z * Q_S1;
diff --git a/library/compiler-builtins/libm/src/math/asin.rs b/library/compiler-builtins/libm/src/math/asin.rs
index 720169bdcf0b..a0bb4918c5a4 100644
--- a/library/compiler-builtins/libm/src/math/asin.rs
+++ b/library/compiler-builtins/libm/src/math/asin.rs
@@ -55,12 +55,14 @@ const Q_S2: f64 = 2.02094576023350569471e+00; /* 0x40002AE5, 0x9C598AC8 */
 const Q_S3: f64 = -6.88283971605453293030e-01; /* 0xBFE6066C, 0x1B8D0159 */
 const Q_S4: f64 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 
+#[inline]
 fn comp_r(z: f64) -> f64 {
     let p = z * (P_S0 + z * (P_S1 + z * (P_S2 + z * (P_S3 + z * (P_S4 + z * P_S5)))));
     let q = 1.0 + z * (Q_S1 + z * (Q_S2 + z * (Q_S3 + z * Q_S4)));
     return p / q;
 }
 
+#[inline]
 pub fn asin(mut x: f64) -> f64 {
     let z: f64;
     let r: f64;
diff --git a/library/compiler-builtins/libm/src/math/asinf.rs b/library/compiler-builtins/libm/src/math/asinf.rs
index 597be4cb7f98..79c85d81d006 100644
--- a/library/compiler-builtins/libm/src/math/asinf.rs
+++ b/library/compiler-builtins/libm/src/math/asinf.rs
@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::fabsf::fabsf;
 use super::sqrt::sqrt;
 
@@ -9,6 +24,7 @@ const P_S1: f32 = -4.2743422091e-02;
 const P_S2: f32 = -8.6563630030e-03;
 const Q_S1: f32 = -7.0662963390e-01;
 
+#[inline]
 fn r(z: f32) -> f32 {
     let p = z * (P_S0 + z * (P_S1 + z * P_S2));
     let q = 1. + z * Q_S1;
diff --git a/library/compiler-builtins/libm/src/math/atan.rs b/library/compiler-builtins/libm/src/math/atan.rs
new file mode 100644
index 000000000000..d057af6d6633
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/atan.rs
@@ -0,0 +1,170 @@
+/* atan(x)
+ * Method
+ *   1. Reduce x to positive by atan(x) = -atan(-x).
+ *   2. According to the integer k=4t+0.25 chopped, t=x, the argument
+ *      is further reduced to one of the following intervals and the
+ *      arctangent of t is evaluated by the corresponding formula:
+ *
+ *      [0,7/16]      atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
+ *      [7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
+ *      [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
+ *      [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
+ *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+use super::fabs;
+use core::f64;
+
+const ATANHI: [f64; 4] = [
+    4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
+    7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
+    9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
+    1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
+];
+
+const ATANLO: [f64; 4] = [
+    2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
+    3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
+    1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
+    6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
+];
+
+const AT: [f64; 11] = [
+    3.33333333333329318027e-01,  /* 0x3FD55555, 0x5555550D */
+    -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
+    1.42857142725034663711e-01,  /* 0x3FC24924, 0x920083FF */
+    -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
+    9.09088713343650656196e-02,  /* 0x3FB745CD, 0xC54C206E */
+    -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
+    6.66107313738753120669e-02,  /* 0x3FB10D66, 0xA0D03D51 */
+    -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
+    4.97687799461593236017e-02,  /* 0x3FA97B4B, 0x24760DEB */
+    -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
+    1.62858201153657823623e-02,  /* 0x3F90AD3A, 0xE322DA11 */
+];
+
+#[inline]
+pub fn atan(x: f64) -> f64 {
+    let mut x = x;
+    let mut ix = (x.to_bits() >> 32) as u32;
+    let sign = ix >> 31;
+    ix &= 0x7fff_ffff;
+    if ix >= 0x4410_0000 {
+        if x.is_nan() {
+            return x;
+        }
+
+        let z = ATANHI[3] + f64::from_bits(0x0380_0000); // 0x1p-120f
+        return if sign != 0 { -z } else { z };
+    }
+
+    let id = if ix < 0x3fdc_0000 {
+        /* |x| < 0.4375 */
+        if ix < 0x3e40_0000 {
+            /* |x| < 2^-27 */
+            if ix < 0x0010_0000 {
+                /* raise underflow for subnormal x */
+                force_eval!(x as f32);
+            }
+
+            return x;
+        }
+
+        -1
+    } else {
+        x = fabs(x);
+        if ix < 0x3ff30000 {
+            /* |x| < 1.1875 */
+            if ix < 0x3fe60000 {
+                /* 7/16 <= |x| < 11/16 */
+                x = (2. * x - 1.) / (2. + x);
+                0
+            } else {
+                /* 11/16 <= |x| < 19/16 */
+                x = (x - 1.) / (x + 1.);
+                1
+            }
+        } else {
+            if ix < 0x40038000 {
+                /* |x| < 2.4375 */
+                x = (x - 1.5) / (1. + 1.5 * x);
+                2
+            } else {
+                /* 2.4375 <= |x| < 2^66 */
+                x = -1. / x;
+                3
+            }
+        }
+    };
+
+    let z = x * x;
+    let w = z * z;
+    /* break sum from i=0 to 10 AT[i]z**(i+1) into odd and even poly */
+    let s1 = z * (AT[0] + w * (AT[2] + w * (AT[4] + w * (AT[6] + w * (AT[8] + w * AT[10])))));
+    let s2 = w * (AT[1] + w * (AT[3] + w * (AT[5] + w * (AT[7] + w * AT[9]))));
+
+    if id < 0 {
+        return x - x * (s1 + s2);
+    }
+
+    let z = ATANHI[id as usize] - (x * (s1 + s2) - ATANLO[id as usize] - x);
+
+    if sign != 0 {
+        -z
+    } else {
+        z
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::atan;
+    use core::f64;
+
+    #[test]
+    fn sanity_check() {
+        for (input, answer) in [
+            (3.0_f64.sqrt() / 3.0, f64::consts::FRAC_PI_6),
+            (1.0, f64::consts::FRAC_PI_4),
+            (3.0_f64.sqrt(), f64::consts::FRAC_PI_3),
+            (-3.0_f64.sqrt() / 3.0, -f64::consts::FRAC_PI_6),
+            (-1.0, -f64::consts::FRAC_PI_4),
+            (-3.0_f64.sqrt(), -f64::consts::FRAC_PI_3),
+        ].iter()
+        {
+            assert!(
+                (atan(*input) - answer) / answer < 1e-5,
+                "\natan({:.4}/16) = {:.4}, actual: {}",
+                input * 16.0,
+                answer,
+                atan(*input)
+            );
+        }
+    }
+
+    #[test]
+    fn zero() {
+        assert_eq!(atan(0.0), 0.0);
+    }
+
+    #[test]
+    fn infinity() {
+        assert_eq!(atan(f64::INFINITY), f64::consts::FRAC_PI_2);
+    }
+
+    #[test]
+    fn minus_infinity() {
+        assert_eq!(atan(f64::NEG_INFINITY), -f64::consts::FRAC_PI_2);
+    }
+
+    #[test]
+    fn nan() {
+        assert!(atan(f64::NAN).is_nan());
+    }
+}
diff --git a/library/compiler-builtins/libm/src/math/atanf.rs b/library/compiler-builtins/libm/src/math/atanf.rs
index 01c41f4ce178..b05152e2bc3e 100644
--- a/library/compiler-builtins/libm/src/math/atanf.rs
+++ b/library/compiler-builtins/libm/src/math/atanf.rs
@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::fabsf;
 
 const ATAN_HI: [f32; 4] = [
diff --git a/library/compiler-builtins/libm/src/math/cos.rs b/library/compiler-builtins/libm/src/math/cos.rs
index e6e9b373674c..df16b5c36a84 100644
--- a/library/compiler-builtins/libm/src/math/cos.rs
+++ b/library/compiler-builtins/libm/src/math/cos.rs
@@ -41,6 +41,7 @@ use super::{k_cos, k_sin, rem_pio2};
 // Accuracy:
 //      TRIG(x) returns trig(x) nearly rounded
 //
+#[inline]
 pub fn cos(x: f64) -> f64 {
     let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
 
diff --git a/library/compiler-builtins/libm/src/math/cosf.rs b/library/compiler-builtins/libm/src/math/cosf.rs
index 79df97e354a3..23faacdc26cf 100644
--- a/library/compiler-builtins/libm/src/math/cosf.rs
+++ b/library/compiler-builtins/libm/src/math/cosf.rs
@@ -1,3 +1,19 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cosf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::{k_cosf, k_sinf, rem_pio2f};
 
 use core::f64::consts::FRAC_PI_2;
diff --git a/library/compiler-builtins/libm/src/math/cosh.rs b/library/compiler-builtins/libm/src/math/cosh.rs
new file mode 100644
index 000000000000..f3f7fbfbeb16
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/cosh.rs
@@ -0,0 +1,33 @@
+use super::exp;
+use super::expm1;
+use super::k_expo2;
+
+#[inline]
+pub fn cosh(mut x: f64) -> f64 {
+    /* |x| */
+    let mut ix = x.to_bits();
+    ix &= 0x7fffffffffffffff;
+    x = f64::from_bits(ix);
+    let w = ix >> 32;
+
+    /* |x| < log(2) */
+    if w < 0x3fe62e42 {
+        if w < 0x3ff00000 - (26 << 20) {
+            let x1p120 = f64::from_bits(0x4770000000000000);
+            force_eval!(x + x1p120);
+            return 1.;
+        }
+        let t = expm1(x); // exponential minus 1
+        return 1. + t * t / (2. * (1. + t));
+    }
+
+    /* |x| < log(DBL_MAX) */
+    if w < 0x40862e42 {
+        let t = exp(x);
+        /* note: if x>log(0x1p26) then the 1/t is not needed */
+        return 0.5 * (t + 1. / t);
+    }
+
+    /* |x| > log(DBL_MAX) or nan */
+    k_expo2(x)
+}
diff --git a/library/compiler-builtins/libm/src/math/expm1.rs b/library/compiler-builtins/libm/src/math/expm1.rs
index 48082cd745cf..9da064ee7736 100644
--- a/library/compiler-builtins/libm/src/math/expm1.rs
+++ b/library/compiler-builtins/libm/src/math/expm1.rs
@@ -1,3 +1,15 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use core::f64;
 
 const O_THRESHOLD: f64 = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
@@ -11,7 +23,7 @@ const Q3: f64 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */
 const Q4: f64 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */
 const Q5: f64 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
 
-#[allow(warnings)]
+#[inline]
 pub fn expm1(mut x: f64) -> f64 {
     let hi: f64;
     let lo: f64;
diff --git a/library/compiler-builtins/libm/src/math/expm1f.rs b/library/compiler-builtins/libm/src/math/expm1f.rs
index 011e09b693d9..8f581733acd6 100644
--- a/library/compiler-builtins/libm/src/math/expm1f.rs
+++ b/library/compiler-builtins/libm/src/math/expm1f.rs
@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 const O_THRESHOLD: f32 = 8.8721679688e+01; /* 0x42b17180 */
 const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */
 const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */
diff --git a/library/compiler-builtins/libm/src/math/expo2.rs b/library/compiler-builtins/libm/src/math/expo2.rs
index 8a121b0a287a..39f9815c4047 100644
--- a/library/compiler-builtins/libm/src/math/expo2.rs
+++ b/library/compiler-builtins/libm/src/math/expo2.rs
@@ -1,7 +1,8 @@
 use super::{combine_words, exp};
 
 /* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */
-pub(crate) fn expo2(x: f64) -> f64 {
+#[inline]
+pub fn expo2(x: f64) -> f64 {
     /* k is such that k*ln2 has minimal relative error and x - kln2 > log(DBL_MIN) */
     const K: i32 = 2043;
     let kln2 = f64::from_bits(0x40962066151add8b);
diff --git a/library/compiler-builtins/libm/src/math/fenv.rs b/library/compiler-builtins/libm/src/math/fenv.rs
new file mode 100644
index 000000000000..63bb20368206
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/fenv.rs
@@ -0,0 +1,33 @@
+// src: musl/src/fenv/fenv.c
+/* Dummy functions for archs lacking fenv implementation */
+
+pub const FE_UNDERFLOW: i32 = 0;
+pub const FE_INEXACT: i32 = 0;
+
+pub const FE_TONEAREST: i32 = 0;
+pub const FE_TOWARDZERO: i32 = 0;
+
+#[inline]
+pub fn feclearexcept(_mask: i32) -> i32 {
+    0
+}
+
+#[inline]
+pub fn feraiseexcept(_mask: i32) -> i32 {
+    0
+}
+
+#[inline]
+pub fn fetestexcept(_mask: i32) -> i32 {
+    0
+}
+
+#[inline]
+pub fn fegetround() -> i32 {
+    FE_TONEAREST
+}
+
+#[inline]
+pub fn fesetround(_r: i32) -> i32 {
+    0
+}
diff --git a/library/compiler-builtins/libm/src/math/fmaf.rs b/library/compiler-builtins/libm/src/math/fmaf.rs
new file mode 100644
index 000000000000..25b04fc23966
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/fmaf.rs
@@ -0,0 +1,100 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_fmaf.c */
+/*-
+ * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS 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 AUTHOR OR CONTRIBUTORS 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.
+ */
+
+use core::f32;
+use core::ptr::read_volatile;
+
+use super::fenv::{
+    feclearexcept, fegetround, feraiseexcept, fesetround, fetestexcept, FE_INEXACT, FE_TONEAREST,
+    FE_TOWARDZERO, FE_UNDERFLOW,
+};
+
+/*
+ * Fused multiply-add: Compute x * y + z with a single rounding error.
+ *
+ * A double has more than twice as much precision than a float, so
+ * direct double-precision arithmetic suffices, except where double
+ * rounding occurs.
+ */
+#[inline]
+pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 {
+    let xy: f64;
+    let mut result: f64;
+    let mut ui: u64;
+    let e: i32;
+
+    xy = x as f64 * y as f64;
+    result = xy + z as f64;
+    ui = result.to_bits();
+    e = (ui >> 52) as i32 & 0x7ff;
+    /* Common case: The double precision result is fine. */
+    if (
+        /* not a halfway case */
+        ui & 0x1fffffff) != 0x10000000 ||
+        /* NaN */
+        e == 0x7ff ||
+        /* exact */
+        (result - xy == z as f64 && result - z as f64 == xy) ||
+        /* not round-to-nearest */
+        fegetround() != FE_TONEAREST
+    {
+        /*
+            underflow may not be raised correctly, example:
+            fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
+        */
+        if e < 0x3ff - 126 && e >= 0x3ff - 149 && fetestexcept(FE_INEXACT) != 0 {
+            feclearexcept(FE_INEXACT);
+            // prevent `xy + vz` from being CSE'd with `xy + z` above
+            let vz: f32 = unsafe { read_volatile(&z) };
+            result = xy + vz as f64;
+            if fetestexcept(FE_INEXACT) != 0 {
+                feraiseexcept(FE_UNDERFLOW);
+            } else {
+                feraiseexcept(FE_INEXACT);
+            }
+        }
+        z = result as f32;
+        return z;
+    }
+
+    /*
+     * If result is inexact, and exactly halfway between two float values,
+     * we need to adjust the low-order bit in the direction of the error.
+     */
+    fesetround(FE_TOWARDZERO);
+    // prevent `vxy + z` from being CSE'd with `xy + z` above
+    let vxy: f64 = unsafe { read_volatile(&xy) };
+    let mut adjusted_result: f64 = vxy + z as f64;
+    fesetround(FE_TONEAREST);
+    if result == adjusted_result {
+        ui = adjusted_result.to_bits();
+        ui += 1;
+        adjusted_result = f64::from_bits(ui);
+    }
+    z = adjusted_result as f32;
+    z
+}
diff --git a/library/compiler-builtins/libm/src/math/k_cosf.rs b/library/compiler-builtins/libm/src/math/k_cosf.rs
index 83d13b2e998a..4aa10c0f0894 100644
--- a/library/compiler-builtins/libm/src/math/k_cosf.rs
+++ b/library/compiler-builtins/libm/src/math/k_cosf.rs
@@ -1,3 +1,19 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/k_cosf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Debugged and optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 /* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
 const C0: f64 = -0.499999997251031003120; /* -0x1ffffffd0c5e81.0p-54 */
 const C1: f64 = 0.0416666233237390631894; /*  0x155553e1053a42.0p-57 */
@@ -5,7 +21,7 @@ const C2: f64 = -0.00138867637746099294692; /* -0x16c087e80f1e27.0p-62 */
 const C3: f64 = 0.0000243904487962774090654; /*  0x199342e0ee5069.0p-68 */
 
 #[inline]
-pub(crate) fn k_cosf(x: f64) -> f32 {
+pub fn k_cosf(x: f64) -> f32 {
     let z = x * x;
     let w = z * z;
     let r = C2 + z * C3;
diff --git a/library/compiler-builtins/libm/src/math/k_expo2.rs b/library/compiler-builtins/libm/src/math/k_expo2.rs
new file mode 100644
index 000000000000..e295c7a53462
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/k_expo2.rs
@@ -0,0 +1,14 @@
+use super::exp;
+
+/* k is such that k*ln2 has minimal relative error and x - kln2 > log(FLT_MIN) */
+const K: i32 = 2043;
+
+/* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
+#[inline]
+pub(crate) fn k_expo2(x: f64) -> f64 {
+    let k_ln2 = f64::from_bits(0x40962066151add8b);
+    /* note that k is odd and scale*scale overflows */
+    let scale = f64::from_bits(((((0x3ff + K / 2) as u32) << 20) as u64) << 32);
+    /* exp(x - k ln2) * 2**(k-1) */
+    exp(x - k_ln2) * scale * scale
+}
diff --git a/library/compiler-builtins/libm/src/math/k_expo2f.rs b/library/compiler-builtins/libm/src/math/k_expo2f.rs
index e2eaa4e6bb02..ec2a2c5e2b8a 100644
--- a/library/compiler-builtins/libm/src/math/k_expo2f.rs
+++ b/library/compiler-builtins/libm/src/math/k_expo2f.rs
@@ -5,7 +5,7 @@ const K: i32 = 235;
 
 /* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
 #[inline]
-pub(crate) fn k_expo2f(x: f32) -> f32 {
+pub fn k_expo2f(x: f32) -> f32 {
     let k_ln2 = f32::from_bits(0x4322e3bc);
     /* note that k is odd and scale*scale overflows */
     let scale = f32::from_bits(((0x7f + K / 2) as u32) << 23);
diff --git a/library/compiler-builtins/libm/src/math/k_sinf.rs b/library/compiler-builtins/libm/src/math/k_sinf.rs
index bb2183afc7ca..1c5f5f98a48f 100644
--- a/library/compiler-builtins/libm/src/math/k_sinf.rs
+++ b/library/compiler-builtins/libm/src/math/k_sinf.rs
@@ -1,3 +1,19 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/k_sinf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 /* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
 const S1: f64 = -0.166666666416265235595; /* -0x15555554cbac77.0p-55 */
 const S2: f64 = 0.0083333293858894631756; /*  0x111110896efbb2.0p-59 */
@@ -5,7 +21,7 @@ const S3: f64 = -0.000198393348360966317347; /* -0x1a00f9e2cae774.0p-65 */
 const S4: f64 = 0.0000027183114939898219064; /*  0x16cd878c3b46a7.0p-71 */
 
 #[inline]
-pub(crate) fn k_sinf(x: f64) -> f32 {
+pub fn k_sinf(x: f64) -> f32 {
     let z = x * x;
     let w = z * z;
     let r = S3 + z * S4;
diff --git a/library/compiler-builtins/libm/src/math/k_tan.rs b/library/compiler-builtins/libm/src/math/k_tan.rs
index b0dd317a2e72..e9ba21499d01 100644
--- a/library/compiler-builtins/libm/src/math/k_tan.rs
+++ b/library/compiler-builtins/libm/src/math/k_tan.rs
@@ -58,7 +58,8 @@ static T: [f64; 13] = [
 const PIO4: f64 = 7.85398163397448278999e-01; /* 3FE921FB, 54442D18 */
 const PIO4_LO: f64 = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
 
-pub(crate) fn k_tan(mut x: f64, mut y: f64, odd: i32) -> f64 {
+#[inline]
+pub fn k_tan(mut x: f64, mut y: f64, odd: i32) -> f64 {
     let hx = (f64::to_bits(x) >> 32) as u32;
     let big = (hx & 0x7fffffff) >= 0x3FE59428; /* |x| >= 0.6744 */
     if big {
diff --git a/library/compiler-builtins/libm/src/math/k_tanf.rs b/library/compiler-builtins/libm/src/math/k_tanf.rs
index db2e0caa7a5e..b9ccf2570f35 100644
--- a/library/compiler-builtins/libm/src/math/k_tanf.rs
+++ b/library/compiler-builtins/libm/src/math/k_tanf.rs
@@ -1,3 +1,14 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/k_tan.c */
+/*
+ * ====================================================
+ * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
 /* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */
 const T: [f64; 6] = [
     0.333331395030791399758,   /* 0x15554d3418c99f.0p-54 */
@@ -9,7 +20,7 @@ const T: [f64; 6] = [
 ];
 
 #[inline]
-pub(crate) fn k_tanf(x: f64, odd: bool) -> f32 {
+pub fn k_tanf(x: f64, odd: bool) -> f32 {
     let z = x * x;
     /*
      * Split up the polynomial into small independent terms to give
diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs
index 07032dcfbca5..3f6fa3211493 100644
--- a/library/compiler-builtins/libm/src/math/mod.rs
+++ b/library/compiler-builtins/libm/src/math/mod.rs
@@ -10,7 +10,6 @@ macro_rules! force_eval {
 mod acos;
 mod acosf;
 mod asin;
-mod asinf;
 mod atan2;
 mod atan2f;
 mod atanf;
@@ -20,6 +19,7 @@ mod ceil;
 mod ceilf;
 mod cos;
 mod cosf;
+mod cosh;
 mod coshf;
 mod exp;
 mod exp2;
@@ -34,6 +34,7 @@ mod fdimf;
 mod floor;
 mod floorf;
 mod fma;
+mod fmaf;
 mod fmod;
 mod fmodf;
 mod hypot;
@@ -59,6 +60,7 @@ mod sqrt;
 mod sqrtf;
 mod tan;
 mod tanf;
+mod tanh;
 mod tanhf;
 mod trunc;
 mod truncf;
@@ -69,6 +71,7 @@ pub use self::acosf::acosf;
 pub use self::asin::asin;
 pub use self::asinf::asinf;
 pub use self::atan2::atan2;
+pub use self::atan::atan;
 pub use self::atan2f::atan2f;
 pub use self::atanf::atanf;
 pub use self::cbrt::cbrt;
@@ -77,6 +80,7 @@ pub use self::ceil::ceil;
 pub use self::ceilf::ceilf;
 pub use self::cos::cos;
 pub use self::cosf::cosf;
+pub use self::cosh::cosh;
 pub use self::coshf::coshf;
 pub use self::exp::exp;
 pub use self::exp2::exp2;
@@ -91,6 +95,7 @@ pub use self::fdimf::fdimf;
 pub use self::floor::floor;
 pub use self::floorf::floorf;
 pub use self::fma::fma;
+pub use self::fmaf::fmaf;
 pub use self::fmod::fmod;
 pub use self::fmodf::fmodf;
 pub use self::hypot::hypot;
@@ -116,14 +121,17 @@ pub use self::sqrt::sqrt;
 pub use self::sqrtf::sqrtf;
 pub use self::tan::tan;
 pub use self::tanf::tanf;
+pub use self::tanh::tanh;
 pub use self::tanhf::tanhf;
 pub use self::trunc::trunc;
 pub use self::truncf::truncf;
 
 // Private modules
 mod expo2;
+mod fenv;
 mod k_cos;
 mod k_cosf;
+mod k_expo2;
 mod k_expo2f;
 mod k_sin;
 mod k_sinf;
@@ -137,6 +145,7 @@ mod rem_pio2f;
 use self::expo2::expo2;
 use self::k_cos::k_cos;
 use self::k_cosf::k_cosf;
+use self::k_expo2::k_expo2;
 use self::k_expo2f::k_expo2f;
 use self::k_sin::k_sin;
 use self::k_sinf::k_sinf;
@@ -147,17 +156,18 @@ use self::rem_pio2_large::rem_pio2_large;
 use self::rem_pio2f::rem_pio2f;
 
 #[inline]
-pub fn get_high_word(x: f64) -> u32 {
+fn get_high_word(x: f64) -> u32 {
     (x.to_bits() >> 32) as u32
 }
 
 #[inline]
-pub fn get_low_word(x: f64) -> u32 {
+fn get_low_word(x: f64) -> u32 {
     x.to_bits() as u32
 }
 
+#[allow(dead_code)]
 #[inline]
-pub fn with_set_high_word(f: f64, hi: u32) -> f64 {
+fn with_set_high_word(f: f64, hi: u32) -> f64 {
     let mut tmp = f.to_bits();
     tmp &= 0x00000000_ffffffff;
     tmp |= (hi as u64) << 32;
@@ -165,7 +175,7 @@ pub fn with_set_high_word(f: f64, hi: u32) -> f64 {
 }
 
 #[inline]
-pub fn with_set_low_word(f: f64, lo: u32) -> f64 {
+fn with_set_low_word(f: f64, lo: u32) -> f64 {
     let mut tmp = f.to_bits();
     tmp &= 0xffffffff_00000000;
     tmp |= lo as u64;
diff --git a/library/compiler-builtins/libm/src/math/rem_pio2.rs b/library/compiler-builtins/libm/src/math/rem_pio2.rs
index 47bab6e65f19..6e655e7d4192 100644
--- a/library/compiler-builtins/libm/src/math/rem_pio2.rs
+++ b/library/compiler-builtins/libm/src/math/rem_pio2.rs
@@ -42,12 +42,14 @@ const PIO2_3T: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 // use rem_pio2_large() for large x
 //
 // caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
+#[inline]
 pub fn rem_pio2(x: f64) -> (i32, f64, f64) {
     let x1p24 = f64::from_bits(0x4170000000000000);
 
     let sign = (f64::to_bits(x) >> 63) as i32;
     let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
 
+    #[inline]
     fn medium(x: f64, ix: u32) -> (i32, f64, f64) {
         /* rint(x/(pi/2)), Assume round-to-nearest. */
         let f_n = x as f64 * INV_PIO2 + TO_INT - TO_INT;
diff --git a/library/compiler-builtins/libm/src/math/rem_pio2_large.rs b/library/compiler-builtins/libm/src/math/rem_pio2_large.rs
index 52b47279cd91..745b700a5eb2 100644
--- a/library/compiler-builtins/libm/src/math/rem_pio2_large.rs
+++ b/library/compiler-builtins/libm/src/math/rem_pio2_large.rs
@@ -1,3 +1,15 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/k_rem_pio2.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.
+ * ====================================================
+ */
+
 use super::floor;
 use super::scalbn;
 
@@ -210,7 +222,7 @@ const PIO2: [f64; 8] = [
 /// more accurately, = 0 mod 8 ). Thus the number of operations are
 /// independent of the exponent of the input.
 #[inline]
-pub(crate) fn rem_pio2_large(x: &[f64], y: &mut [f64], e0: i32, prec: usize) -> i32 {
+pub fn rem_pio2_large(x: &[f64], y: &mut [f64], e0: i32, prec: usize) -> i32 {
     let x1p24 = f64::from_bits(0x4170000000000000); // 0x1p24 === 2 ^ 24
     let x1p_24 = f64::from_bits(0x3e70000000000000); // 0x1p_24 === 2 ^ (-24)
 
diff --git a/library/compiler-builtins/libm/src/math/rem_pio2f.rs b/library/compiler-builtins/libm/src/math/rem_pio2f.rs
index 73ec3775d3dc..5e7a7d43904d 100644
--- a/library/compiler-builtins/libm/src/math/rem_pio2f.rs
+++ b/library/compiler-builtins/libm/src/math/rem_pio2f.rs
@@ -1,3 +1,19 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2f.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Debugged and optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::rem_pio2_large;
 
 use core::f64;
@@ -16,7 +32,7 @@ const PIO2_1T: f64 = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
 /// use double precision for everything except passing x
 /// use __rem_pio2_large() for large x
 #[inline]
-pub(crate) fn rem_pio2f(x: f32) -> (i32, f64) {
+pub fn rem_pio2f(x: f32) -> (i32, f64) {
     let x64 = x as f64;
 
     let mut tx: [f64; 1] = [0.];
diff --git a/library/compiler-builtins/libm/src/math/sin.rs b/library/compiler-builtins/libm/src/math/sin.rs
index 13eb30248ab7..e749094e6681 100644
--- a/library/compiler-builtins/libm/src/math/sin.rs
+++ b/library/compiler-builtins/libm/src/math/sin.rs
@@ -40,6 +40,7 @@ use super::{k_cos, k_sin, rem_pio2};
 //
 // Accuracy:
 //      TRIG(x) returns trig(x) nearly rounded
+#[inline]
 pub fn sin(x: f64) -> f64 {
     let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120f === 2 ^ 120
 
diff --git a/library/compiler-builtins/libm/src/math/sinf.rs b/library/compiler-builtins/libm/src/math/sinf.rs
index 09f62cddcae3..c9b02bcdc513 100644
--- a/library/compiler-builtins/libm/src/math/sinf.rs
+++ b/library/compiler-builtins/libm/src/math/sinf.rs
@@ -1,3 +1,19 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_sinf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::{k_cosf, k_sinf, rem_pio2f};
 
 use core::f64::consts::FRAC_PI_2;
diff --git a/library/compiler-builtins/libm/src/math/sinh.rs b/library/compiler-builtins/libm/src/math/sinh.rs
index 0579871719a8..684e8e30984a 100644
--- a/library/compiler-builtins/libm/src/math/sinh.rs
+++ b/library/compiler-builtins/libm/src/math/sinh.rs
@@ -4,6 +4,7 @@ use super::{expm1, expo2};
 //         = (exp(x)-1 + (exp(x)-1)/exp(x))/2
 //         = x + x^3/6 + o(x^5)
 //
+#[inline]
 pub fn sinh(x: f64) -> f64 {
     // union {double f; uint64_t i;} u = {.f = x};
     // uint32_t w;
diff --git a/library/compiler-builtins/libm/src/math/tan.rs b/library/compiler-builtins/libm/src/math/tan.rs
index 92bbb6221229..5a5f178a51cc 100644
--- a/library/compiler-builtins/libm/src/math/tan.rs
+++ b/library/compiler-builtins/libm/src/math/tan.rs
@@ -39,6 +39,7 @@ use super::{k_tan, rem_pio2};
 //
 // Accuracy:
 //      TRIG(x) returns trig(x) nearly rounded
+#[inline]
 pub fn tan(x: f64) -> f64 {
     let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
 
diff --git a/library/compiler-builtins/libm/src/math/tanf.rs b/library/compiler-builtins/libm/src/math/tanf.rs
index 6bfbe06c1b4b..15a462d4e19f 100644
--- a/library/compiler-builtins/libm/src/math/tanf.rs
+++ b/library/compiler-builtins/libm/src/math/tanf.rs
@@ -1,3 +1,19 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_tanf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::{k_tanf, rem_pio2f};
 
 use core::f64::consts::FRAC_PI_2;
diff --git a/library/compiler-builtins/libm/src/math/tanh.rs b/library/compiler-builtins/libm/src/math/tanh.rs
new file mode 100644
index 000000000000..1c3dd0be438e
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/tanh.rs
@@ -0,0 +1,53 @@
+use super::expm1;
+
+/* 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)
+ */
+#[inline]
+pub fn tanh(mut x: f64) -> f64 {
+    let mut uf: f64 = x;
+    let mut ui: u64 = f64::to_bits(uf);
+
+    let w: u32;
+    let sign: bool;
+    let mut t: f64;
+
+    /* x = |x| */
+    sign = ui >> 63 != 0;
+    ui &= !1 / 2;
+    uf = f64::from_bits(ui);
+    x = uf;
+    w = (ui >> 32) as u32;
+
+    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 - 0.0 / x;
+        } else {
+            t = expm1(2.0 * x);
+            t = 1.0 - 2.0 / (t + 2.0);
+        }
+    } else if w > 0x3fd058ae {
+        /* |x| > log(5/3)/2 ~= 0.2554 */
+        t = expm1(2.0 * x);
+        t = t / (t + 2.0);
+    } else if w >= 0x00100000 {
+        /* |x| >= 0x1p-1022, up to 2ulp error in [0.1,0.2554] */
+        t = expm1(-2.0 * x);
+        t = -t / (t + 2.0);
+    } else {
+        /* |x| is subnormal */
+        /* note: the branch above would not raise underflow in [0x1p-1023,0x1p-1022) */
+        force_eval!(x as f32);
+        t = x;
+    }
+
+    if sign {
+        -t
+    } else {
+        t
+    }
+}
diff --git a/library/compiler-builtins/libm/test-generator/Cargo.toml b/library/compiler-builtins/libm/test-generator/Cargo.toml
index f45d173b4516..b810d9daf630 100644
--- a/library/compiler-builtins/libm/test-generator/Cargo.toml
+++ b/library/compiler-builtins/libm/test-generator/Cargo.toml
@@ -6,3 +6,4 @@ publish = false
 
 [dependencies]
 rand = "0.5.3"
+itertools = "0.7.8"
diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs
index e65f0f69bb04..09a9cb8f5887 100644
--- a/library/compiler-builtins/libm/test-generator/src/main.rs
+++ b/library/compiler-builtins/libm/test-generator/src/main.rs
@@ -4,13 +4,15 @@
 // NOTE usually the only thing you need to do to test a new math function is to add it to one of the
 // macro invocations found in the bottom of this file.
 
+#[macro_use]
+extern crate itertools;
 extern crate rand;
 
 use std::error::Error;
 use std::fmt::Write as _0;
 use std::fs::{self, File};
 use std::io::Write as _1;
-use std::{i16, u16, u32, u64, u8};
+use std::{f32, f64, i16, u16, u32, u64, u8};
 
 use rand::{Rng, SeedableRng, XorShiftRng};
 
@@ -34,6 +36,30 @@ fn f64(rng: &mut XorShiftRng) -> f64 {
     f64::from_bits(sign + exponent + mantissa)
 }
 
+const EDGE_CASES32: &[f32] = &[
+    -0.,
+    0.,
+    f32::EPSILON,
+    f32::INFINITY,
+    f32::MAX,
+    f32::MIN,
+    f32::MIN_POSITIVE,
+    f32::NAN,
+    f32::NEG_INFINITY,
+];
+
+const EDGE_CASES64: &[f64] = &[
+    -0.,
+    0.,
+    f64::EPSILON,
+    f64::INFINITY,
+    f64::MAX,
+    f64::MIN,
+    f64::MIN_POSITIVE,
+    f64::NAN,
+    f64::NEG_INFINITY,
+];
+
 // fn(f32) -> f32
 macro_rules! f32_f32 {
     ($($intr:ident,)*) => {
@@ -45,8 +71,9 @@ macro_rules! f32_f32 {
 
             $(
                 let mut cases = String::new();
-                for _ in 0..NTESTS {
-                    let inp = f32(rng);
+
+                // random inputs
+                for inp in EDGE_CASES32.iter().cloned().chain((0..NTESTS).map(|_| f32(rng))) {
                     let out = unsafe { $intr(inp) };
 
                     let inp = inp.to_bits();
@@ -112,11 +139,17 @@ macro_rules! f32f32_f32 {
                 $(fn $intr(_: f32, _: f32) -> f32;)*
             }
 
+            let mut rng2 = rng.clone();
+            let mut rng3 = rng.clone();
             $(
                 let mut cases = String::new();
-                for _ in 0..NTESTS {
-                    let i1 = f32(rng);
-                    let i2 = f32(rng);
+                for (i1, i2) in iproduct!(
+                    EDGE_CASES32.iter().cloned(),
+                    EDGE_CASES32.iter().cloned()
+                ).chain(EDGE_CASES32.iter().map(|i1| (*i1, f32(rng))))
+                    .chain(EDGE_CASES32.iter().map(|i2| (f32(&mut rng2), *i2)))
+                    .chain((0..NTESTS).map(|_| (f32(&mut rng3), f32(&mut rng3))))
+                {
                     let out = unsafe { $intr(i1, i2) };
 
                     let i1 = i1.to_bits();
@@ -186,12 +219,16 @@ macro_rules! f32f32f32_f32 {
                 $(fn $intr(_: f32, _: f32, _: f32) -> f32;)*
             }
 
+            let mut rng2 = rng.clone();
             $(
                 let mut cases = String::new();
-                for _ in 0..NTESTS {
-                    let i1 = f32(rng);
-                    let i2 = f32(rng);
-                    let i3 = f32(rng);
+                for (i1, i2, i3) in iproduct!(
+                    EDGE_CASES32.iter().cloned(),
+                    EDGE_CASES32.iter().cloned(),
+                    EDGE_CASES32.iter().cloned()
+                ).chain(EDGE_CASES32.iter().map(|i1| (*i1, f32(rng), f32(rng))))
+                    .chain((0..NTESTS).map(|_| (f32(&mut rng2), f32(&mut rng2), f32(&mut rng2))))
+                {
                     let out = unsafe { $intr(i1, i2, i3) };
 
                     let i1 = i1.to_bits();
@@ -266,10 +303,10 @@ macro_rules! f32i32_f32 {
                 $(fn $intr(_: f32, _: i32) -> f32;)*
             }
 
+            let mut rng2 = rng.clone();
             $(
                 let mut cases = String::new();
-                for _ in 0..NTESTS {
-                    let i1 = f32(rng);
+                for i1 in EDGE_CASES32.iter().cloned().chain((0..NTESTS).map(|_| f32(&mut rng2))) {
                     let i2 = rng.gen_range(i16::MIN, i16::MAX);
                     let out = unsafe { $intr(i1, i2 as i32) };
 
@@ -342,8 +379,7 @@ macro_rules! f64_f64 {
 
             $(
                 let mut cases = String::new();
-                for _ in 0..NTESTS {
-                    let inp = f64(rng);
+                for inp in EDGE_CASES64.iter().cloned().chain((0..NTESTS).map(|_| f64(rng))) {
                     let out = unsafe { $intr(inp) };
 
                     let inp = inp.to_bits();
@@ -412,11 +448,17 @@ macro_rules! f64f64_f64 {
                 $(fn $intr(_: f64, _: f64) -> f64;)*
             }
 
+            let mut rng2 = rng.clone();
+            let mut rng3 = rng.clone();
             $(
                 let mut cases = String::new();
-                for _ in 0..NTESTS {
-                    let i1 = f64(rng);
-                    let i2 = f64(rng);
+                for (i1, i2) in iproduct!(
+                    EDGE_CASES64.iter().cloned(),
+                    EDGE_CASES64.iter().cloned()
+                ).chain(EDGE_CASES64.iter().map(|i1| (*i1, f64(rng))))
+                    .chain(EDGE_CASES64.iter().map(|i2| (f64(&mut rng2), *i2)))
+                    .chain((0..NTESTS).map(|_| (f64(&mut rng3), f64(&mut rng3))))
+                {
                     let out = unsafe { $intr(i1, i2) };
 
                     let i1 = i1.to_bits();
@@ -485,12 +527,16 @@ macro_rules! f64f64f64_f64 {
                 $(fn $intr(_: f64, _: f64, _: f64) -> f64;)*
             }
 
+            let mut rng2 = rng.clone();
             $(
                 let mut cases = String::new();
-                for _ in 0..NTESTS {
-                    let i1 = f64(rng);
-                    let i2 = f64(rng);
-                    let i3 = f64(rng);
+                for (i1, i2, i3) in iproduct!(
+                    EDGE_CASES64.iter().cloned(),
+                    EDGE_CASES64.iter().cloned(),
+                    EDGE_CASES64.iter().cloned()
+                ).chain(EDGE_CASES64.iter().map(|i1| (*i1, f64(rng), f64(rng))))
+                    .chain((0..NTESTS).map(|_| (f64(&mut rng2), f64(&mut rng2), f64(&mut rng2))))
+                {
                     let out = unsafe { $intr(i1, i2, i3) };
 
                     let i1 = i1.to_bits();
@@ -565,10 +611,10 @@ macro_rules! f64i32_f64 {
                 $(fn $intr(_: f64, _: i32) -> f64;)*
             }
 
+            let mut rng2 = rng.clone();
             $(
                 let mut cases = String::new();
-                for _ in 0..NTESTS {
-                    let i1 = f64(rng);
+                for i1 in EDGE_CASES64.iter().cloned().chain((0..NTESTS).map(|_| f64(&mut rng2))) {
                     let i2 = rng.gen_range(i16::MIN, i16::MAX);
                     let out = unsafe { $intr(i1, i2 as i32) };
 
@@ -687,7 +733,7 @@ f32f32_f32! {
 
 // With signature `fn(f32, f32, f32) -> f32`
 f32f32f32_f32! {
-    // fmaf,
+    fmaf,
 }
 
 // With signature `fn(f32, i32) -> f32`
@@ -699,11 +745,11 @@ f32i32_f32! {
 f64_f64! {
     acos,
     asin,
-    // atan,
+    atan,
     cbrt,
     ceil,
     cos,
-    // cosh,
+    cosh,
     exp,
     exp2,
     expm1,
@@ -717,7 +763,7 @@ f64_f64! {
     sinh,
     sqrt,
     tan,
-    // tanh,
+    tanh,
     trunc,
     fabs,
 }