From 440a1b7273e663ce0c7680aba62efe744d2ab556 Mon Sep 17 00:00:00 2001
From: C Jones <code@calebjones.net>
Date: Fri, 13 Jul 2018 22:10:30 -0400
Subject: [PATCH 1/3] Implement part of sin/cos with quadrant selection
 unimplemented

---
 library/compiler-builtins/libm/src/lib.rs     |   4 -
 .../compiler-builtins/libm/src/math/cos.rs    |  73 +++++
 .../compiler-builtins/libm/src/math/mod.rs    |  14 +-
 .../compiler-builtins/libm/src/math/sin.rs    |  80 +++++
 .../libm/src/math/trig_common.rs              | 289 ++++++++++++++++++
 5 files changed, 452 insertions(+), 8 deletions(-)
 create mode 100644 library/compiler-builtins/libm/src/math/cos.rs
 create mode 100644 library/compiler-builtins/libm/src/math/sin.rs
 create mode 100644 library/compiler-builtins/libm/src/math/trig_common.rs

diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs
index df7ee813fbaf..0b9efeeb338a 100644
--- a/library/compiler-builtins/libm/src/lib.rs
+++ b/library/compiler-builtins/libm/src/lib.rs
@@ -386,10 +386,8 @@ pub trait F64Ext: private::Sealed {
 
     fn hypot(self, other: Self) -> Self;
 
-    #[cfg(todo)]
     fn sin(self) -> Self;
 
-    #[cfg(todo)]
     fn cos(self) -> Self;
 
     #[cfg(todo)]
@@ -548,13 +546,11 @@ impl F64Ext for f64 {
         hypot(self, other)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn sin(self) -> Self {
         sin(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn cos(self) -> Self {
         cos(self)
diff --git a/library/compiler-builtins/libm/src/math/cos.rs b/library/compiler-builtins/libm/src/math/cos.rs
new file mode 100644
index 000000000000..91412c6489a7
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/cos.rs
@@ -0,0 +1,73 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cos.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 math::trig_common::{_cos, _sin, rem_pio2};
+
+// cos(x)
+// Return cosine function of x.
+//
+// kernel function:
+//      __sin           ... sine function on [-pi/4,pi/4]
+//      __cos           ... cosine function on [-pi/4,pi/4]
+//      __rem_pio2      ... argument reduction routine
+//
+// Method.
+//      Let S,C and T denote the sin, cos and tan respectively on
+//      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+//      in [-pi/4 , +pi/4], and let n = k mod 4.
+//      We have
+//
+//          n        sin(x)      cos(x)        tan(x)
+//     ----------------------------------------------------------
+//          0          S           C             T
+//          1          C          -S            -1/T
+//          2         -S          -C             T
+//          3         -C           S            -1/T
+//     ----------------------------------------------------------
+//
+// Special cases:
+//      Let trig be any of sin, cos, or tan.
+//      trig(+-INF)  is NaN, with signals;
+//      trig(NaN)    is that NaN;
+//
+// Accuracy:
+//      TRIG(x) returns trig(x) nearly rounded
+//
+pub fn cos(x: f64) -> f64 {
+    let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
+
+    /* |x| ~< pi/4 */
+    if ix <= 0x3fe921fb {
+        if ix < 0x3e46a09e {
+            /* if x < 2**-27 * sqrt(2) */
+            /* raise inexact if x != 0 */
+            if x as i32 == 0 {
+                return 1.0;
+            }
+        }
+        return _cos(x, 0.0);
+    }
+
+    /* cos(Inf or NaN) is NaN */
+    if ix >= 0x7ff00000 {
+        return x - x;
+    }
+
+    /* argument reduction needed */
+    let (n, y0, y1) = rem_pio2(x);
+    match n & 3 {
+        0 => _cos(y0, y1),
+        1 => -_sin(y0, y1, 1),
+        2 => -_cos(y0, y1),
+        _ => _sin(y0, y1, 1),
+    }
+}
diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs
index 42c596857d7e..a657da810e33 100644
--- a/library/compiler-builtins/libm/src/math/mod.rs
+++ b/library/compiler-builtins/libm/src/math/mod.rs
@@ -15,6 +15,7 @@ mod cbrt;
 mod cbrtf;
 mod ceil;
 mod ceilf;
+mod cos;
 mod cosf;
 mod exp;
 mod exp2;
@@ -46,6 +47,7 @@ mod round;
 mod roundf;
 mod scalbn;
 mod scalbnf;
+mod sin;
 mod sinf;
 mod sqrt;
 mod sqrtf;
@@ -63,6 +65,7 @@ pub use self::cbrt::cbrt;
 pub use self::cbrtf::cbrtf;
 pub use self::ceil::ceil;
 pub use self::ceilf::ceilf;
+pub use self::cos::cos;
 pub use self::cosf::cosf;
 pub use self::exp::exp;
 pub use self::exp2::exp2;
@@ -94,6 +97,7 @@ pub use self::round::round;
 pub use self::roundf::roundf;
 pub use self::scalbn::scalbn;
 pub use self::scalbnf::scalbnf;
+pub use self::sin::sin;
 pub use self::sinf::sinf;
 pub use self::sqrt::sqrt;
 pub use self::sqrtf::sqrtf;
@@ -106,11 +110,13 @@ mod k_sinf;
 mod k_tanf;
 mod rem_pio2_large;
 mod rem_pio2f;
+mod trig_common;
 
-use self::{
-    k_cosf::k_cosf, k_sinf::k_sinf, k_tanf::k_tanf, rem_pio2_large::rem_pio2_large,
-    rem_pio2f::rem_pio2f,
-};
+use self::k_cosf::k_cosf;
+use self::k_sinf::k_sinf;
+use self::k_tanf::k_tanf;
+use self::rem_pio2_large::rem_pio2_large;
+use self::rem_pio2f::rem_pio2f;
 
 #[inline]
 pub fn get_high_word(x: f64) -> u32 {
diff --git a/library/compiler-builtins/libm/src/math/sin.rs b/library/compiler-builtins/libm/src/math/sin.rs
new file mode 100644
index 000000000000..0717c1ae3edf
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/sin.rs
@@ -0,0 +1,80 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_sin.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;
+
+use math::trig_common::{_cos, _sin, rem_pio2};
+
+// sin(x)
+// Return sine function of x.
+//
+// kernel function:
+//      __sin            ... sine function on [-pi/4,pi/4]
+//      __cos            ... cose function on [-pi/4,pi/4]
+//      __rem_pio2       ... argument reduction routine
+//
+// Method.
+//      Let S,C and T denote the sin, cos and tan respectively on
+//      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+//      in [-pi/4 , +pi/4], and let n = k mod 4.
+//      We have
+//
+//          n        sin(x)      cos(x)        tan(x)
+//     ----------------------------------------------------------
+//          0          S           C             T
+//          1          C          -S            -1/T
+//          2         -S          -C             T
+//          3         -C           S            -1/T
+//     ----------------------------------------------------------
+//
+// Special cases:
+//      Let trig be any of sin, cos, or tan.
+//      trig(+-INF)  is NaN, with signals;
+//      trig(NaN)    is that NaN;
+//
+// Accuracy:
+//      TRIG(x) returns trig(x) nearly rounded
+pub fn sin(x: f64) -> f64 {
+    let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120f === 2 ^ 120
+
+    /* High word of x. */
+    let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
+
+    /* |x| ~< pi/4 */
+    if ix <= 0x3fe921fb {
+        if ix < 0x3e500000 {
+            /* |x| < 2**-26 */
+            /* raise inexact if x != 0 and underflow if subnormal*/
+            if ix < 0x00100000 {
+                force_eval!(x / x1p120);
+            } else {
+                force_eval!(x + x1p120);
+            }
+            return x;
+        }
+        return _sin(x, 0.0, 0);
+    }
+
+    /* sin(Inf or NaN) is NaN */
+    if ix >= 0x7ff00000 {
+        return x - x;
+    }
+
+    /* argument reduction needed */
+    let (n, y0, y1) = rem_pio2(x);
+    match n & 3 {
+        0 => _sin(y0, y1, 1),
+        1 => _cos(y0, y1),
+        2 => -_sin(y0, y1, 1),
+        _ => -_cos(y0, y1),
+    }
+}
diff --git a/library/compiler-builtins/libm/src/math/trig_common.rs b/library/compiler-builtins/libm/src/math/trig_common.rs
new file mode 100644
index 000000000000..59e75e2dd6b9
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/trig_common.rs
@@ -0,0 +1,289 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* origin: FreeBSD /usr/src/lib/msun/src/k_sin.c */
+
+const S1: f64 = -1.66666666666666324348e-01; /* 0xBFC55555, 0x55555549 */
+const S2: f64 = 8.33333333332248946124e-03; /* 0x3F811111, 0x1110F8A6 */
+const S3: f64 = -1.98412698298579493134e-04; /* 0xBF2A01A0, 0x19C161D5 */
+const S4: f64 = 2.75573137070700676789e-06; /* 0x3EC71DE3, 0x57B1FE7D */
+const S5: f64 = -2.50507602534068634195e-08; /* 0xBE5AE5E6, 0x8A2B9CEB */
+const S6: f64 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
+
+// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
+// Input x is assumed to be bounded by ~pi/4 in magnitude.
+// Input y is the tail of x.
+// Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
+//
+// Algorithm
+//      1. Since sin(-x) = -sin(x), we need only to consider positive x.
+//      2. Callers must return sin(-0) = -0 without calling here since our
+//         odd polynomial is not evaluated in a way that preserves -0.
+//         Callers may do the optimization sin(x) ~ x for tiny x.
+//      3. sin(x) is approximated by a polynomial of degree 13 on
+//         [0,pi/4]
+//                               3            13
+//              sin(x) ~ x + S1*x + ... + S6*x
+//         where
+//
+//      |sin(x)         2     4     6     8     10     12  |     -58
+//      |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
+//      |  x                                               |
+//
+//      4. sin(x+y) = sin(x) + sin'(x')*y
+//                  ~ sin(x) + (1-x*x/2)*y
+//         For better accuracy, let
+//                   3      2      2      2      2
+//              r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
+//         then                   3    2
+//              sin(x) = x + (S1*x + (x *(r-y/2)+y))
+pub fn _sin(x: f64, y: f64, iy: i32) -> f64 {
+    let z = x * x;
+    let w = z * z;
+    let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
+    let v = z * x;
+    if iy == 0 {
+        x + v * (S1 + z * r)
+    } else {
+        x - ((z * (0.5 * y - v * r) - y) - v * S1)
+    }
+}
+
+/* origin: FreeBSD /usr/src/lib/msun/src/k_cos.c */
+const C1: f64 = 4.16666666666666019037e-02; /* 0x3FA55555, 0x5555554C */
+const C2: f64 = -1.38888888888741095749e-03; /* 0xBF56C16C, 0x16C15177 */
+const C3: f64 = 2.48015872894767294178e-05; /* 0x3EFA01A0, 0x19CB1590 */
+const C4: f64 = -2.75573143513906633035e-07; /* 0xBE927E4F, 0x809C52AD */
+const C5: f64 = 2.08757232129817482790e-09; /* 0x3E21EE9E, 0xBDB4B1C4 */
+const C6: f64 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
+
+// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
+// Input x is assumed to be bounded by ~pi/4 in magnitude.
+// Input y is the tail of x.
+//
+// Algorithm
+//      1. Since cos(-x) = cos(x), we need only to consider positive x.
+//      2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
+//      3. cos(x) is approximated by a polynomial of degree 14 on
+//         [0,pi/4]
+//                                       4            14
+//              cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
+//         where the remez error is
+//
+//      |              2     4     6     8     10    12     14 |     -58
+//      |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
+//      |                                                      |
+//
+//                     4     6     8     10    12     14
+//      4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
+//             cos(x) ~ 1 - x*x/2 + r
+//         since cos(x+y) ~ cos(x) - sin(x)*y
+//                        ~ cos(x) - x*y,
+//         a correction term is necessary in cos(x) and hence
+//              cos(x+y) = 1 - (x*x/2 - (r - x*y))
+//         For better accuracy, rearrange to
+//              cos(x+y) ~ w + (tmp + (r-x*y))
+//         where w = 1 - x*x/2 and tmp is a tiny correction term
+//         (1 - x*x/2 == w + tmp exactly in infinite precision).
+//         The exactness of w + tmp in infinite precision depends on w
+//         and tmp having the same precision as x.  If they have extra
+//         precision due to compiler bugs, then the extra precision is
+//         only good provided it is retained in all terms of the final
+//         expression for cos().  Retention happens in all cases tested
+//         under FreeBSD, so don't pessimize things by forcibly clipping
+//         any extra precision in w.
+pub fn _cos(x: f64, y: f64) -> f64 {
+    let z = x * x;
+    let w = z * z;
+    let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
+    let hz = 0.5 * z;
+    let w = 1.0 - hz;
+    w + (((1.0 - w) - hz) + (z * r - x * y))
+}
+
+// origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c
+// Optimized by Bruce D. Evans. */
+
+// #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+// #define EPS DBL_EPSILON
+const EPS: f64 = 2.2204460492503131e-16;
+// #elif FLT_EVAL_METHOD==2
+// #define EPS LDBL_EPSILON
+// #endif
+
+// TODO: Support FLT_EVAL_METHOD?
+
+#[allow(unused, non_upper_case_globals)]
+const toint: f64 = 1.5 / EPS;
+/// 53 bits of 2/pi
+#[allow(unused, non_upper_case_globals)]
+const invpio2: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
+/// first 33 bits of pi/2
+#[allow(unused, non_upper_case_globals)]
+const pio2_1: f64 = 1.57079632673412561417e+00; /* 0x3FF921FB, 0x54400000 */
+/// pi/2 - pio2_1
+#[allow(unused, non_upper_case_globals)]
+const pio2_1t: f64 = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */
+/// second 33 bits of pi/2
+#[allow(unused, non_upper_case_globals)]
+const pio2_2: f64 = 6.07710050630396597660e-11; /* 0x3DD0B461, 0x1A600000 */
+/// pi/2 - (pio2_1+pio2_2)
+#[allow(unused, non_upper_case_globals)]
+const pio2_2t: f64 = 2.02226624879595063154e-21; /* 0x3BA3198A, 0x2E037073 */
+/// third 33 bits of pi/2
+#[allow(unused, non_upper_case_globals)]
+const pio2_3: f64 = 2.02226624871116645580e-21; /* 0x3BA3198A, 0x2E000000 */
+/// pi/2 - (pio2_1+pio2_2+pio2_3)
+#[allow(unused, non_upper_case_globals)]
+const pio2_3t: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
+
+/* __rem_pio2(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __rem_pio2_large() for large x
+ */
+
+/* caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
+/*
+{
+	union {double f; uint64_t i;} u = {x};
+	double_t z,w,t,r,fn;
+	double tx[3],ty[2];
+	uint32_t ix;
+	int sign, n, ex, ey, i;
+
+	sign = u.i>>63;
+	ix = u.i>>32 & 0x7fffffff;
+	if (ix <= 0x400f6a7a) {  /* |x| ~<= 5pi/4 */
+		if ((ix & 0xfffff) == 0x921fb)  /* |x| ~= pi/2 or 2pi/2 */
+			goto medium;  /* cancellation -- use medium case */
+		if (ix <= 0x4002d97c) {  /* |x| ~<= 3pi/4 */
+			if (!sign) {
+				z = x - pio2_1;  /* one round good to 85 bits */
+				y[0] = z - pio2_1t;
+				y[1] = (z-y[0]) - pio2_1t;
+				return 1;
+			} else {
+				z = x + pio2_1;
+				y[0] = z + pio2_1t;
+				y[1] = (z-y[0]) + pio2_1t;
+				return -1;
+			}
+		} else {
+			if (!sign) {
+				z = x - 2*pio2_1;
+				y[0] = z - 2*pio2_1t;
+				y[1] = (z-y[0]) - 2*pio2_1t;
+				return 2;
+			} else {
+				z = x + 2*pio2_1;
+				y[0] = z + 2*pio2_1t;
+				y[1] = (z-y[0]) + 2*pio2_1t;
+				return -2;
+			}
+		}
+	}
+	if (ix <= 0x401c463b) {  /* |x| ~<= 9pi/4 */
+		if (ix <= 0x4015fdbc) {  /* |x| ~<= 7pi/4 */
+			if (ix == 0x4012d97c)  /* |x| ~= 3pi/2 */
+				goto medium;
+			if (!sign) {
+				z = x - 3*pio2_1;
+				y[0] = z - 3*pio2_1t;
+				y[1] = (z-y[0]) - 3*pio2_1t;
+				return 3;
+			} else {
+				z = x + 3*pio2_1;
+				y[0] = z + 3*pio2_1t;
+				y[1] = (z-y[0]) + 3*pio2_1t;
+				return -3;
+			}
+		} else {
+			if (ix == 0x401921fb)  /* |x| ~= 4pi/2 */
+				goto medium;
+			if (!sign) {
+				z = x - 4*pio2_1;
+				y[0] = z - 4*pio2_1t;
+				y[1] = (z-y[0]) - 4*pio2_1t;
+				return 4;
+			} else {
+				z = x + 4*pio2_1;
+				y[0] = z + 4*pio2_1t;
+				y[1] = (z-y[0]) + 4*pio2_1t;
+				return -4;
+			}
+		}
+	}
+	if (ix < 0x413921fb) {  /* |x| ~< 2^20*(pi/2), medium size */
+medium:
+		/* rint(x/(pi/2)), Assume round-to-nearest. */
+		fn = (double_t)x*invpio2 + toint - toint;
+		n = (int32_t)fn;
+		r = x - fn*pio2_1;
+		w = fn*pio2_1t;  /* 1st round, good to 85 bits */
+		y[0] = r - w;
+		u.f = y[0];
+		ey = u.i>>52 & 0x7ff;
+		ex = ix>>20;
+		if (ex - ey > 16) { /* 2nd round, good to 118 bits */
+			t = r;
+			w = fn*pio2_2;
+			r = t - w;
+			w = fn*pio2_2t - ((t-r)-w);
+			y[0] = r - w;
+			u.f = y[0];
+			ey = u.i>>52 & 0x7ff;
+			if (ex - ey > 49) {  /* 3rd round, good to 151 bits, covers all cases */
+				t = r;
+				w = fn*pio2_3;
+				r = t - w;
+				w = fn*pio2_3t - ((t-r)-w);
+				y[0] = r - w;
+			}
+		}
+		y[1] = (r - y[0]) - w;
+		return n;
+	}
+	/*
+	 * all other (large) arguments
+	 */
+	if (ix >= 0x7ff00000) {  /* x is inf or NaN */
+		y[0] = y[1] = x - x;
+		return 0;
+	}
+	/* set z = scalbn(|x|,-ilogb(x)+23) */
+	u.f = x;
+	u.i &= (uint64_t)-1>>12;
+	u.i |= (uint64_t)(0x3ff + 23)<<52;
+	z = u.f;
+	for (i=0; i < 2; i++) {
+		tx[i] = (double)(int32_t)z;
+		z     = (z-tx[i])*0x1p24;
+	}
+	tx[i] = z;
+	/* skip zero terms, first term is non-zero */
+	while (tx[i] == 0.0)
+		i--;
+	n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
+	if (sign) {
+		y[0] = -ty[0];
+		y[1] = -ty[1];
+		return -n;
+	}
+	y[0] = ty[0];
+	y[1] = ty[1];
+	return n;
+}
+*/
+
+pub fn rem_pio2(_x: f64) -> (i32, f64, f64) {
+    unimplemented!()
+}

From b2751897d2fde5c7aa0ec30b7d45b714011cd0dd Mon Sep 17 00:00:00 2001
From: C Jones <code@calebjones.net>
Date: Fri, 13 Jul 2018 23:03:21 -0400
Subject: [PATCH 2/3] Convert rem_pio2 code, split up modules

---
 .../compiler-builtins/libm/src/math/cos.rs    |  39 ++-
 .../compiler-builtins/libm/src/math/exp2.rs   |   2 +-
 .../compiler-builtins/libm/src/math/k_cos.rs  |  62 ++++
 .../compiler-builtins/libm/src/math/k_sin.rs  |  57 ++++
 .../compiler-builtins/libm/src/math/mod.rs    |  10 +-
 .../libm/src/math/rem_pio2.rs                 | 187 ++++++++++++
 .../compiler-builtins/libm/src/math/sin.rs    |  41 ++-
 .../libm/src/math/trig_common.rs              | 289 ------------------
 .../libm/test-generator/src/main.rs           |   4 +-
 9 files changed, 356 insertions(+), 335 deletions(-)
 create mode 100644 library/compiler-builtins/libm/src/math/k_cos.rs
 create mode 100644 library/compiler-builtins/libm/src/math/k_sin.rs
 create mode 100644 library/compiler-builtins/libm/src/math/rem_pio2.rs
 delete mode 100644 library/compiler-builtins/libm/src/math/trig_common.rs

diff --git a/library/compiler-builtins/libm/src/math/cos.rs b/library/compiler-builtins/libm/src/math/cos.rs
index 91412c6489a7..e6e9b373674c 100644
--- a/library/compiler-builtins/libm/src/math/cos.rs
+++ b/library/compiler-builtins/libm/src/math/cos.rs
@@ -1,24 +1,23 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_cos.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.
- * ====================================================
- */
+// origin: FreeBSD /usr/src/lib/msun/src/s_cos.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 math::trig_common::{_cos, _sin, rem_pio2};
+use super::{k_cos, k_sin, rem_pio2};
 
 // cos(x)
 // Return cosine function of x.
 //
 // kernel function:
-//      __sin           ... sine function on [-pi/4,pi/4]
-//      __cos           ... cosine function on [-pi/4,pi/4]
-//      __rem_pio2      ... argument reduction routine
+//      k_sin           ... sine function on [-pi/4,pi/4]
+//      k_cos           ... cosine function on [-pi/4,pi/4]
+//      rem_pio2        ... argument reduction routine
 //
 // Method.
 //      Let S,C and T denote the sin, cos and tan respectively on
@@ -54,7 +53,7 @@ pub fn cos(x: f64) -> f64 {
                 return 1.0;
             }
         }
-        return _cos(x, 0.0);
+        return k_cos(x, 0.0);
     }
 
     /* cos(Inf or NaN) is NaN */
@@ -65,9 +64,9 @@ pub fn cos(x: f64) -> f64 {
     /* argument reduction needed */
     let (n, y0, y1) = rem_pio2(x);
     match n & 3 {
-        0 => _cos(y0, y1),
-        1 => -_sin(y0, y1, 1),
-        2 => -_cos(y0, y1),
-        _ => _sin(y0, y1, 1),
+        0 => k_cos(y0, y1),
+        1 => -k_sin(y0, y1, 1),
+        2 => -k_cos(y0, y1),
+        _ => k_sin(y0, y1, 1),
     }
 }
diff --git a/library/compiler-builtins/libm/src/math/exp2.rs b/library/compiler-builtins/libm/src/math/exp2.rs
index 0ab11989679b..3952e93007e1 100644
--- a/library/compiler-builtins/libm/src/math/exp2.rs
+++ b/library/compiler-builtins/libm/src/math/exp2.rs
@@ -24,7 +24,7 @@
 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 // SUCH DAMAGE.
 
-use super::scalbn::scalbn;
+use super::scalbn;
 
 const TBLSIZE: usize = 256;
 
diff --git a/library/compiler-builtins/libm/src/math/k_cos.rs b/library/compiler-builtins/libm/src/math/k_cos.rs
new file mode 100644
index 000000000000..693950d1d418
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/k_cos.rs
@@ -0,0 +1,62 @@
+// origin: FreeBSD /usr/src/lib/msun/src/k_cos.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.
+// ====================================================
+
+const C1: f64 = 4.16666666666666019037e-02; /* 0x3FA55555, 0x5555554C */
+const C2: f64 = -1.38888888888741095749e-03; /* 0xBF56C16C, 0x16C15177 */
+const C3: f64 = 2.48015872894767294178e-05; /* 0x3EFA01A0, 0x19CB1590 */
+const C4: f64 = -2.75573143513906633035e-07; /* 0xBE927E4F, 0x809C52AD */
+const C5: f64 = 2.08757232129817482790e-09; /* 0x3E21EE9E, 0xBDB4B1C4 */
+const C6: f64 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
+
+// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
+// Input x is assumed to be bounded by ~pi/4 in magnitude.
+// Input y is the tail of x.
+//
+// Algorithm
+//      1. Since cos(-x) = cos(x), we need only to consider positive x.
+//      2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
+//      3. cos(x) is approximated by a polynomial of degree 14 on
+//         [0,pi/4]
+//                                       4            14
+//              cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
+//         where the remez error is
+//
+//      |              2     4     6     8     10    12     14 |     -58
+//      |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
+//      |                                                      |
+//
+//                     4     6     8     10    12     14
+//      4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
+//             cos(x) ~ 1 - x*x/2 + r
+//         since cos(x+y) ~ cos(x) - sin(x)*y
+//                        ~ cos(x) - x*y,
+//         a correction term is necessary in cos(x) and hence
+//              cos(x+y) = 1 - (x*x/2 - (r - x*y))
+//         For better accuracy, rearrange to
+//              cos(x+y) ~ w + (tmp + (r-x*y))
+//         where w = 1 - x*x/2 and tmp is a tiny correction term
+//         (1 - x*x/2 == w + tmp exactly in infinite precision).
+//         The exactness of w + tmp in infinite precision depends on w
+//         and tmp having the same precision as x.  If they have extra
+//         precision due to compiler bugs, then the extra precision is
+//         only good provided it is retained in all terms of the final
+//         expression for cos().  Retention happens in all cases tested
+//         under FreeBSD, so don't pessimize things by forcibly clipping
+//         any extra precision in w.
+#[inline]
+pub fn k_cos(x: f64, y: f64) -> f64 {
+    let z = x * x;
+    let w = z * z;
+    let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
+    let hz = 0.5 * z;
+    let w = 1.0 - hz;
+    w + (((1.0 - w) - hz) + (z * r - x * y))
+}
diff --git a/library/compiler-builtins/libm/src/math/k_sin.rs b/library/compiler-builtins/libm/src/math/k_sin.rs
new file mode 100644
index 000000000000..3e07c3594475
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/k_sin.rs
@@ -0,0 +1,57 @@
+// origin: FreeBSD /usr/src/lib/msun/src/k_sin.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.
+// ====================================================
+
+const S1: f64 = -1.66666666666666324348e-01; /* 0xBFC55555, 0x55555549 */
+const S2: f64 = 8.33333333332248946124e-03; /* 0x3F811111, 0x1110F8A6 */
+const S3: f64 = -1.98412698298579493134e-04; /* 0xBF2A01A0, 0x19C161D5 */
+const S4: f64 = 2.75573137070700676789e-06; /* 0x3EC71DE3, 0x57B1FE7D */
+const S5: f64 = -2.50507602534068634195e-08; /* 0xBE5AE5E6, 0x8A2B9CEB */
+const S6: f64 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
+
+// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
+// Input x is assumed to be bounded by ~pi/4 in magnitude.
+// Input y is the tail of x.
+// Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
+//
+// Algorithm
+//      1. Since sin(-x) = -sin(x), we need only to consider positive x.
+//      2. Callers must return sin(-0) = -0 without calling here since our
+//         odd polynomial is not evaluated in a way that preserves -0.
+//         Callers may do the optimization sin(x) ~ x for tiny x.
+//      3. sin(x) is approximated by a polynomial of degree 13 on
+//         [0,pi/4]
+//                               3            13
+//              sin(x) ~ x + S1*x + ... + S6*x
+//         where
+//
+//      |sin(x)         2     4     6     8     10     12  |     -58
+//      |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
+//      |  x                                               |
+//
+//      4. sin(x+y) = sin(x) + sin'(x')*y
+//                  ~ sin(x) + (1-x*x/2)*y
+//         For better accuracy, let
+//                   3      2      2      2      2
+//              r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
+//         then                   3    2
+//              sin(x) = x + (S1*x + (x *(r-y/2)+y))
+#[inline]
+pub fn k_sin(x: f64, y: f64, iy: i32) -> f64 {
+    let z = x * x;
+    let w = z * z;
+    let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
+    let v = z * x;
+    if iy == 0 {
+        x + v * (S1 + z * r)
+    } else {
+        x - ((z * (0.5 * y - v * r) - y) - v * S1)
+    }
+}
diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs
index a657da810e33..7b3d9abee697 100644
--- a/library/compiler-builtins/libm/src/math/mod.rs
+++ b/library/compiler-builtins/libm/src/math/mod.rs
@@ -6,6 +6,7 @@ macro_rules! force_eval {
     };
 }
 
+// Public modules
 mod acos;
 mod acosf;
 mod asin;
@@ -105,16 +106,23 @@ pub use self::tanf::tanf;
 pub use self::trunc::trunc;
 pub use self::truncf::truncf;
 
+// Private modules
+mod k_cos;
 mod k_cosf;
+mod k_sin;
 mod k_sinf;
 mod k_tanf;
+mod rem_pio2;
 mod rem_pio2_large;
 mod rem_pio2f;
-mod trig_common;
 
+// Private re-imports
+use self::k_cos::k_cos;
 use self::k_cosf::k_cosf;
+use self::k_sin::k_sin;
 use self::k_sinf::k_sinf;
 use self::k_tanf::k_tanf;
+use self::rem_pio2::rem_pio2;
 use self::rem_pio2_large::rem_pio2_large;
 use self::rem_pio2f::rem_pio2f;
 
diff --git a/library/compiler-builtins/libm/src/math/rem_pio2.rs b/library/compiler-builtins/libm/src/math/rem_pio2.rs
new file mode 100644
index 000000000000..68db7056b63c
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/rem_pio2.rs
@@ -0,0 +1,187 @@
+// origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.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.
+// ====================================================
+//
+// Optimized by Bruce D. Evans. */
+
+use super::rem_pio2_large;
+
+// #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+// #define EPS DBL_EPSILON
+const EPS: f64 = 2.2204460492503131e-16;
+// #elif FLT_EVAL_METHOD==2
+// #define EPS LDBL_EPSILON
+// #endif
+
+// TODO: Support FLT_EVAL_METHOD?
+
+const TO_INT: f64 = 1.5 / EPS;
+/// 53 bits of 2/pi
+const INV_PIO2: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
+/// first 33 bits of pi/2
+const PIO2_1: f64 = 1.57079632673412561417e+00; /* 0x3FF921FB, 0x54400000 */
+/// pi/2 - PIO2_1
+const PIO2_1T: f64 = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */
+/// second 33 bits of pi/2
+const PIO2_2: f64 = 6.07710050630396597660e-11; /* 0x3DD0B461, 0x1A600000 */
+/// pi/2 - (PIO2_1+PIO2_2)
+const PIO2_2T: f64 = 2.02226624879595063154e-21; /* 0x3BA3198A, 0x2E037073 */
+/// third 33 bits of pi/2
+const PIO2_3: f64 = 2.02226624871116645580e-21; /* 0x3BA3198A, 0x2E000000 */
+/// pi/2 - (PIO2_1+PIO2_2+PIO2_3)
+const PIO2_3T: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
+
+// return the remainder of x rem pi/2 in y[0]+y[1]
+// use rem_pio2_large() for large x
+//
+// caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
+pub fn rem_pio2(x: f64) -> (i32, f64, f64) {
+    let x1p24 = f64::from_bits(0x7041);
+
+    let sign = (f64::to_bits(x) >> 63) as i32;
+    let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
+
+    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;
+        let n = f_n as i32;
+        let mut r = x - f_n * PIO2_1;
+        let mut w = f_n * PIO2_1T; /* 1st round, good to 85 bits */
+        let mut y0 = r - w;
+        let ui = f64::to_bits(y0);
+        let ey = (ui >> 52) as i32 & 0x7ff;
+        let ex = (ix >> 20) as i32;
+        if ex - ey > 16 {
+            /* 2nd round, good to 118 bits */
+            let t = r;
+            w = f_n * PIO2_2;
+            r = t - w;
+            w = f_n * PIO2_2T - ((t - r) - w);
+            y0 = r - w;
+            let ey = (f64::to_bits(y0) >> 52) as i32 & 0x7ff;
+            if ex - ey > 49 {
+                /* 3rd round, good to 151 bits, covers all cases */
+                let t = r;
+                w = f_n * PIO2_3;
+                r = t - w;
+                w = f_n * PIO2_3T - ((t - r) - w);
+                y0 = r - w;
+            }
+        }
+        let y1 = (r - y0) - w;
+        return (n, y0, y1);
+    }
+
+    if ix <= 0x400f6a7a {
+        /* |x| ~<= 5pi/4 */
+        if (ix & 0xfffff) == 0x921fb {
+            /* |x| ~= pi/2 or 2pi/2 */
+            medium(x, ix); /* cancellation -- use medium case */
+        }
+        if ix <= 0x4002d97c {
+            /* |x| ~<= 3pi/4 */
+            if sign == 0 {
+                let z = x - PIO2_1; /* one round good to 85 bits */
+                let y0 = z - PIO2_1T;
+                let y1 = (z - y0) - PIO2_1T;
+                return (1, y0, y1);
+            } else {
+                let z = x + PIO2_1;
+                let y0 = z + PIO2_1T;
+                let y1 = (z - y0) + PIO2_1T;
+                return (-1, y0, y1);
+            }
+        } else {
+            if sign == 0 {
+                let z = x - 2.0 * PIO2_1;
+                let y0 = z - 2.0 * PIO2_1T;
+                let y1 = (z - y0) - 2.0 * PIO2_1T;
+                return (2, y0, y1);
+            } else {
+                let z = x + 2.0 * PIO2_1;
+                let y0 = z + 2.0 * PIO2_1T;
+                let y1 = (z - y0) + 2.0 * PIO2_1T;
+                return (-2, y0, y1);
+            }
+        }
+    }
+    if ix <= 0x401c463b {
+        /* |x| ~<= 9pi/4 */
+        if ix <= 0x4015fdbc {
+            /* |x| ~<= 7pi/4 */
+            if ix == 0x4012d97c {
+                /* |x| ~= 3pi/2 */
+                return medium(x, ix);
+            }
+            if sign == 0 {
+                let z = x - 3.0 * PIO2_1;
+                let y0 = z - 3.0 * PIO2_1T;
+                let y1 = (z - y0) - 3.0 * PIO2_1T;
+                return (3, y0, y1);
+            } else {
+                let z = x + 3.0 * PIO2_1;
+                let y0 = z + 3.0 * PIO2_1T;
+                let y1 = (z - y0) + 3.0 * PIO2_1T;
+                return (-3, y0, y1);
+            }
+        } else {
+            if ix == 0x401921fb {
+                /* |x| ~= 4pi/2 */
+                return medium(x, ix);
+            }
+            if sign == 0 {
+                let z = x - 4.0 * PIO2_1;
+                let y0 = z - 4.0 * PIO2_1T;
+                let y1 = (z - y0) - 4.0 * PIO2_1T;
+                return (4, y0, y1);
+            } else {
+                let z = x + 4.0 * PIO2_1;
+                let y0 = z + 4.0 * PIO2_1T;
+                let y1 = (z - y0) + 4.0 * PIO2_1T;
+                return (-4, y0, y1);
+            }
+        }
+    }
+    if ix < 0x413921fb {
+        /* |x| ~< 2^20*(pi/2), medium size */
+        return medium(x, ix);
+    }
+    /*
+     * all other (large) arguments
+     */
+    if ix >= 0x7ff00000 {
+        /* x is inf or NaN */
+        let y0 = x - x;
+        let y1 = y0;
+        return (0, y0, y1);
+    }
+    /* set z = scalbn(|x|,-ilogb(x)+23) */
+    let mut ui = f64::to_bits(x);
+    ui &= (!1) >> 12;
+    ui |= (0x3ff + 23) << 52;
+    let mut z = f64::from_bits(ui);
+    let mut tx = [0.0; 3];
+    for i in 0..2 {
+        tx[i] = z as i32 as f64;
+        z = (z - tx[i]) * x1p24;
+    }
+    tx[2] = z;
+    /* skip zero terms, first term is non-zero */
+    let mut i = 2;
+    while tx[i] == 0.0 {
+        i -= 1;
+    }
+    let mut ty = [0.0; 3];
+    let n = rem_pio2_large(&tx[..=i], &mut ty, ((ix >> 20) - (0x3ff + 23)) as i32, 1);
+    if sign != 0 {
+        return (-n, -ty[0], -ty[1]);
+    }
+    return (n, ty[0], ty[1]);
+}
diff --git a/library/compiler-builtins/libm/src/math/sin.rs b/library/compiler-builtins/libm/src/math/sin.rs
index 0717c1ae3edf..13eb30248ab7 100644
--- a/library/compiler-builtins/libm/src/math/sin.rs
+++ b/library/compiler-builtins/libm/src/math/sin.rs
@@ -1,26 +1,23 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_sin.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.
- * ====================================================
- */
+// origin: FreeBSD /usr/src/lib/msun/src/s_sin.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;
-
-use math::trig_common::{_cos, _sin, rem_pio2};
+use super::{k_cos, k_sin, rem_pio2};
 
 // sin(x)
 // Return sine function of x.
 //
 // kernel function:
-//      __sin            ... sine function on [-pi/4,pi/4]
-//      __cos            ... cose function on [-pi/4,pi/4]
-//      __rem_pio2       ... argument reduction routine
+//      k_sin            ... sine function on [-pi/4,pi/4]
+//      k_cos            ... cose function on [-pi/4,pi/4]
+//      rem_pio2         ... argument reduction routine
 //
 // Method.
 //      Let S,C and T denote the sin, cos and tan respectively on
@@ -61,7 +58,7 @@ pub fn sin(x: f64) -> f64 {
             }
             return x;
         }
-        return _sin(x, 0.0, 0);
+        return k_sin(x, 0.0, 0);
     }
 
     /* sin(Inf or NaN) is NaN */
@@ -72,9 +69,9 @@ pub fn sin(x: f64) -> f64 {
     /* argument reduction needed */
     let (n, y0, y1) = rem_pio2(x);
     match n & 3 {
-        0 => _sin(y0, y1, 1),
-        1 => _cos(y0, y1),
-        2 => -_sin(y0, y1, 1),
-        _ => -_cos(y0, y1),
+        0 => k_sin(y0, y1, 1),
+        1 => k_cos(y0, y1),
+        2 => -k_sin(y0, y1, 1),
+        _ => -k_cos(y0, y1),
     }
 }
diff --git a/library/compiler-builtins/libm/src/math/trig_common.rs b/library/compiler-builtins/libm/src/math/trig_common.rs
deleted file mode 100644
index 59e75e2dd6b9..000000000000
--- a/library/compiler-builtins/libm/src/math/trig_common.rs
+++ /dev/null
@@ -1,289 +0,0 @@
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-
-/* origin: FreeBSD /usr/src/lib/msun/src/k_sin.c */
-
-const S1: f64 = -1.66666666666666324348e-01; /* 0xBFC55555, 0x55555549 */
-const S2: f64 = 8.33333333332248946124e-03; /* 0x3F811111, 0x1110F8A6 */
-const S3: f64 = -1.98412698298579493134e-04; /* 0xBF2A01A0, 0x19C161D5 */
-const S4: f64 = 2.75573137070700676789e-06; /* 0x3EC71DE3, 0x57B1FE7D */
-const S5: f64 = -2.50507602534068634195e-08; /* 0xBE5AE5E6, 0x8A2B9CEB */
-const S6: f64 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
-
-// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
-// Input x is assumed to be bounded by ~pi/4 in magnitude.
-// Input y is the tail of x.
-// Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
-//
-// Algorithm
-//      1. Since sin(-x) = -sin(x), we need only to consider positive x.
-//      2. Callers must return sin(-0) = -0 without calling here since our
-//         odd polynomial is not evaluated in a way that preserves -0.
-//         Callers may do the optimization sin(x) ~ x for tiny x.
-//      3. sin(x) is approximated by a polynomial of degree 13 on
-//         [0,pi/4]
-//                               3            13
-//              sin(x) ~ x + S1*x + ... + S6*x
-//         where
-//
-//      |sin(x)         2     4     6     8     10     12  |     -58
-//      |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
-//      |  x                                               |
-//
-//      4. sin(x+y) = sin(x) + sin'(x')*y
-//                  ~ sin(x) + (1-x*x/2)*y
-//         For better accuracy, let
-//                   3      2      2      2      2
-//              r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
-//         then                   3    2
-//              sin(x) = x + (S1*x + (x *(r-y/2)+y))
-pub fn _sin(x: f64, y: f64, iy: i32) -> f64 {
-    let z = x * x;
-    let w = z * z;
-    let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
-    let v = z * x;
-    if iy == 0 {
-        x + v * (S1 + z * r)
-    } else {
-        x - ((z * (0.5 * y - v * r) - y) - v * S1)
-    }
-}
-
-/* origin: FreeBSD /usr/src/lib/msun/src/k_cos.c */
-const C1: f64 = 4.16666666666666019037e-02; /* 0x3FA55555, 0x5555554C */
-const C2: f64 = -1.38888888888741095749e-03; /* 0xBF56C16C, 0x16C15177 */
-const C3: f64 = 2.48015872894767294178e-05; /* 0x3EFA01A0, 0x19CB1590 */
-const C4: f64 = -2.75573143513906633035e-07; /* 0xBE927E4F, 0x809C52AD */
-const C5: f64 = 2.08757232129817482790e-09; /* 0x3E21EE9E, 0xBDB4B1C4 */
-const C6: f64 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
-
-// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
-// Input x is assumed to be bounded by ~pi/4 in magnitude.
-// Input y is the tail of x.
-//
-// Algorithm
-//      1. Since cos(-x) = cos(x), we need only to consider positive x.
-//      2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
-//      3. cos(x) is approximated by a polynomial of degree 14 on
-//         [0,pi/4]
-//                                       4            14
-//              cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
-//         where the remez error is
-//
-//      |              2     4     6     8     10    12     14 |     -58
-//      |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
-//      |                                                      |
-//
-//                     4     6     8     10    12     14
-//      4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
-//             cos(x) ~ 1 - x*x/2 + r
-//         since cos(x+y) ~ cos(x) - sin(x)*y
-//                        ~ cos(x) - x*y,
-//         a correction term is necessary in cos(x) and hence
-//              cos(x+y) = 1 - (x*x/2 - (r - x*y))
-//         For better accuracy, rearrange to
-//              cos(x+y) ~ w + (tmp + (r-x*y))
-//         where w = 1 - x*x/2 and tmp is a tiny correction term
-//         (1 - x*x/2 == w + tmp exactly in infinite precision).
-//         The exactness of w + tmp in infinite precision depends on w
-//         and tmp having the same precision as x.  If they have extra
-//         precision due to compiler bugs, then the extra precision is
-//         only good provided it is retained in all terms of the final
-//         expression for cos().  Retention happens in all cases tested
-//         under FreeBSD, so don't pessimize things by forcibly clipping
-//         any extra precision in w.
-pub fn _cos(x: f64, y: f64) -> f64 {
-    let z = x * x;
-    let w = z * z;
-    let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
-    let hz = 0.5 * z;
-    let w = 1.0 - hz;
-    w + (((1.0 - w) - hz) + (z * r - x * y))
-}
-
-// origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c
-// Optimized by Bruce D. Evans. */
-
-// #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
-// #define EPS DBL_EPSILON
-const EPS: f64 = 2.2204460492503131e-16;
-// #elif FLT_EVAL_METHOD==2
-// #define EPS LDBL_EPSILON
-// #endif
-
-// TODO: Support FLT_EVAL_METHOD?
-
-#[allow(unused, non_upper_case_globals)]
-const toint: f64 = 1.5 / EPS;
-/// 53 bits of 2/pi
-#[allow(unused, non_upper_case_globals)]
-const invpio2: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
-/// first 33 bits of pi/2
-#[allow(unused, non_upper_case_globals)]
-const pio2_1: f64 = 1.57079632673412561417e+00; /* 0x3FF921FB, 0x54400000 */
-/// pi/2 - pio2_1
-#[allow(unused, non_upper_case_globals)]
-const pio2_1t: f64 = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */
-/// second 33 bits of pi/2
-#[allow(unused, non_upper_case_globals)]
-const pio2_2: f64 = 6.07710050630396597660e-11; /* 0x3DD0B461, 0x1A600000 */
-/// pi/2 - (pio2_1+pio2_2)
-#[allow(unused, non_upper_case_globals)]
-const pio2_2t: f64 = 2.02226624879595063154e-21; /* 0x3BA3198A, 0x2E037073 */
-/// third 33 bits of pi/2
-#[allow(unused, non_upper_case_globals)]
-const pio2_3: f64 = 2.02226624871116645580e-21; /* 0x3BA3198A, 0x2E000000 */
-/// pi/2 - (pio2_1+pio2_2+pio2_3)
-#[allow(unused, non_upper_case_globals)]
-const pio2_3t: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
-
-/* __rem_pio2(x,y)
- *
- * return the remainder of x rem pi/2 in y[0]+y[1]
- * use __rem_pio2_large() for large x
- */
-
-/* caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
-/*
-{
-	union {double f; uint64_t i;} u = {x};
-	double_t z,w,t,r,fn;
-	double tx[3],ty[2];
-	uint32_t ix;
-	int sign, n, ex, ey, i;
-
-	sign = u.i>>63;
-	ix = u.i>>32 & 0x7fffffff;
-	if (ix <= 0x400f6a7a) {  /* |x| ~<= 5pi/4 */
-		if ((ix & 0xfffff) == 0x921fb)  /* |x| ~= pi/2 or 2pi/2 */
-			goto medium;  /* cancellation -- use medium case */
-		if (ix <= 0x4002d97c) {  /* |x| ~<= 3pi/4 */
-			if (!sign) {
-				z = x - pio2_1;  /* one round good to 85 bits */
-				y[0] = z - pio2_1t;
-				y[1] = (z-y[0]) - pio2_1t;
-				return 1;
-			} else {
-				z = x + pio2_1;
-				y[0] = z + pio2_1t;
-				y[1] = (z-y[0]) + pio2_1t;
-				return -1;
-			}
-		} else {
-			if (!sign) {
-				z = x - 2*pio2_1;
-				y[0] = z - 2*pio2_1t;
-				y[1] = (z-y[0]) - 2*pio2_1t;
-				return 2;
-			} else {
-				z = x + 2*pio2_1;
-				y[0] = z + 2*pio2_1t;
-				y[1] = (z-y[0]) + 2*pio2_1t;
-				return -2;
-			}
-		}
-	}
-	if (ix <= 0x401c463b) {  /* |x| ~<= 9pi/4 */
-		if (ix <= 0x4015fdbc) {  /* |x| ~<= 7pi/4 */
-			if (ix == 0x4012d97c)  /* |x| ~= 3pi/2 */
-				goto medium;
-			if (!sign) {
-				z = x - 3*pio2_1;
-				y[0] = z - 3*pio2_1t;
-				y[1] = (z-y[0]) - 3*pio2_1t;
-				return 3;
-			} else {
-				z = x + 3*pio2_1;
-				y[0] = z + 3*pio2_1t;
-				y[1] = (z-y[0]) + 3*pio2_1t;
-				return -3;
-			}
-		} else {
-			if (ix == 0x401921fb)  /* |x| ~= 4pi/2 */
-				goto medium;
-			if (!sign) {
-				z = x - 4*pio2_1;
-				y[0] = z - 4*pio2_1t;
-				y[1] = (z-y[0]) - 4*pio2_1t;
-				return 4;
-			} else {
-				z = x + 4*pio2_1;
-				y[0] = z + 4*pio2_1t;
-				y[1] = (z-y[0]) + 4*pio2_1t;
-				return -4;
-			}
-		}
-	}
-	if (ix < 0x413921fb) {  /* |x| ~< 2^20*(pi/2), medium size */
-medium:
-		/* rint(x/(pi/2)), Assume round-to-nearest. */
-		fn = (double_t)x*invpio2 + toint - toint;
-		n = (int32_t)fn;
-		r = x - fn*pio2_1;
-		w = fn*pio2_1t;  /* 1st round, good to 85 bits */
-		y[0] = r - w;
-		u.f = y[0];
-		ey = u.i>>52 & 0x7ff;
-		ex = ix>>20;
-		if (ex - ey > 16) { /* 2nd round, good to 118 bits */
-			t = r;
-			w = fn*pio2_2;
-			r = t - w;
-			w = fn*pio2_2t - ((t-r)-w);
-			y[0] = r - w;
-			u.f = y[0];
-			ey = u.i>>52 & 0x7ff;
-			if (ex - ey > 49) {  /* 3rd round, good to 151 bits, covers all cases */
-				t = r;
-				w = fn*pio2_3;
-				r = t - w;
-				w = fn*pio2_3t - ((t-r)-w);
-				y[0] = r - w;
-			}
-		}
-		y[1] = (r - y[0]) - w;
-		return n;
-	}
-	/*
-	 * all other (large) arguments
-	 */
-	if (ix >= 0x7ff00000) {  /* x is inf or NaN */
-		y[0] = y[1] = x - x;
-		return 0;
-	}
-	/* set z = scalbn(|x|,-ilogb(x)+23) */
-	u.f = x;
-	u.i &= (uint64_t)-1>>12;
-	u.i |= (uint64_t)(0x3ff + 23)<<52;
-	z = u.f;
-	for (i=0; i < 2; i++) {
-		tx[i] = (double)(int32_t)z;
-		z     = (z-tx[i])*0x1p24;
-	}
-	tx[i] = z;
-	/* skip zero terms, first term is non-zero */
-	while (tx[i] == 0.0)
-		i--;
-	n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
-	if (sign) {
-		y[0] = -ty[0];
-		y[1] = -ty[1];
-		return -n;
-	}
-	y[0] = ty[0];
-	y[1] = ty[1];
-	return n;
-}
-*/
-
-pub fn rem_pio2(_x: f64) -> (i32, f64, f64) {
-    unimplemented!()
-}
diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs
index b639cf11b19a..c8e35cea7b00 100644
--- a/library/compiler-builtins/libm/test-generator/src/main.rs
+++ b/library/compiler-builtins/libm/test-generator/src/main.rs
@@ -702,7 +702,7 @@ f64_f64! {
     // atan,
     cbrt,
     ceil,
-    // cos,
+    cos,
     // cosh,
     exp,
     exp2,
@@ -713,7 +713,7 @@ f64_f64! {
     log1p,
     log2,
     round,
-    // sin,
+    sin,
     // sinh,
     sqrt,
     // tan,

From 75428a13cf93a710deabaa2ab475984e1f26e18c Mon Sep 17 00:00:00 2001
From: C Jones <code@calebjones.net>
Date: Sat, 14 Jul 2018 18:40:52 -0400
Subject: [PATCH 3/3] Fix x1p24 constant

---
 library/compiler-builtins/libm/src/math/rem_pio2.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/library/compiler-builtins/libm/src/math/rem_pio2.rs b/library/compiler-builtins/libm/src/math/rem_pio2.rs
index 68db7056b63c..47bab6e65f19 100644
--- a/library/compiler-builtins/libm/src/math/rem_pio2.rs
+++ b/library/compiler-builtins/libm/src/math/rem_pio2.rs
@@ -43,7 +43,7 @@ const PIO2_3T: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 //
 // caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
 pub fn rem_pio2(x: f64) -> (i32, f64, f64) {
-    let x1p24 = f64::from_bits(0x7041);
+    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;