Merge rust-lang#80

80: Add sin and cos r=japaric a=porglezomp

Closes rust-lang#11 
Closes rust-lang#33

Co-authored-by: C Jones <[email protected]>

Author: bors[bot]

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)

src/math/cos.rs

@@ -0,0 +1,72 @@
+// 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 super::{k_cos, k_sin, rem_pio2};
+
+// cos(x)
+// Return cosine function of x.
+//
+// kernel function:
+//      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
+//      [-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 k_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 => k_cos(y0, y1),
+        1 => -k_sin(y0, y1, 1),
+        2 => -k_cos(y0, y1),
+        _ => k_sin(y0, y1, 1),
+    }
+}

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;
 

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))
+}

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)
+    }
+}

src/math/mod.rs

@@ -6,6 +6,7 @@ macro_rules! force_eval {
     };
 }
 
+// Public modules
 mod acos;
 mod acosf;
 mod asin;
@@ -15,6 +16,7 @@ mod cbrt;
 mod cbrtf;
 mod ceil;
 mod ceilf;
+mod cos;
 mod cosf;
 mod exp;
 mod exp2;
@@ -46,6 +48,7 @@ mod round;
 mod roundf;
 mod scalbn;
 mod scalbnf;
+mod sin;
 mod sinf;
 mod sqrt;
 mod sqrtf;
@@ -63,6 +66,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,23 +98,33 @@ 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;
 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;
 
-use self::{
-    k_cosf::k_cosf, k_sinf::k_sinf, k_tanf::k_tanf, rem_pio2_large::rem_pio2_large,
-    rem_pio2f::rem_pio2f,
-};
+// 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;
 
 #[inline]
 pub fn get_high_word(x: f64) -> u32 {