Implement tan

Also includes implementing the private k_tan function.

Closes rust-lang#36

Author: porglezomp

src/lib.rs

@@ -390,7 +390,6 @@ pub trait F64Ext: private::Sealed {
 
     fn cos(self) -> Self;
 
-    #[cfg(todo)]
     fn tan(self) -> Self;
 
     #[cfg(todo)]
@@ -556,7 +555,6 @@ impl F64Ext for f64 {
         cos(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn tan(self) -> Self {
         tan(self)

src/math/k_tan.rs

@@ -0,0 +1,105 @@
+// 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.
+// ====================================================
+
+// kernel tan 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 odd indicates whether tan (if odd = 0) or -1/tan (if odd = 1) is returned.
+//
+// Algorithm
+//      1. Since tan(-x) = -tan(x), we need only to consider positive x.
+//      2. Callers must return tan(-0) = -0 without calling here since our
+//         odd polynomial is not evaluated in a way that preserves -0.
+//         Callers may do the optimization tan(x) ~ x for tiny x.
+//      3. tan(x) is approximated by a odd polynomial of degree 27 on
+//         [0,0.67434]
+//                               3             27
+//              tan(x) ~ x + T1*x + ... + T13*x
+//         where
+//
+//              |tan(x)         2     4            26   |     -59.2
+//              |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
+//              |  x                                    |
+//
+//         Note: tan(x+y) = tan(x) + tan'(x)*y
+//                        ~ tan(x) + (1+x*x)*y
+//         Therefore, for better accuracy in computing tan(x+y), let
+//                   3      2      2       2       2
+//              r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
+//         then
+//                                  3    2
+//              tan(x+y) = x + (T1*x + (x *(r+y)+y))
+//
+//      4. For x in [0.67434,pi/4],  let y = pi/4 - x, then
+//              tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
+//                     = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
+static T: [f64; 13] = [
+    3.33333333333334091986e-01,  /* 3FD55555, 55555563 */
+    1.33333333333201242699e-01,  /* 3FC11111, 1110FE7A */
+    5.39682539762260521377e-02,  /* 3FABA1BA, 1BB341FE */
+    2.18694882948595424599e-02,  /* 3F9664F4, 8406D637 */
+    8.86323982359930005737e-03,  /* 3F8226E3, E96E8493 */
+    3.59207910759131235356e-03,  /* 3F6D6D22, C9560328 */
+    1.45620945432529025516e-03,  /* 3F57DBC8, FEE08315 */
+    5.88041240820264096874e-04,  /* 3F4344D8, F2F26501 */
+    2.46463134818469906812e-04,  /* 3F3026F7, 1A8D1068 */
+    7.81794442939557092300e-05,  /* 3F147E88, A03792A6 */
+    7.14072491382608190305e-05,  /* 3F12B80F, 32F0A7E9 */
+    -1.85586374855275456654e-05, /* BEF375CB, DB605373 */
+    2.59073051863633712884e-05,  /* 3EFB2A70, 74BF7AD4 */
+];
+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 {
+    let hx = (f64::to_bits(x) >> 32) as u32;
+    let big = (hx & 0x7fffffff) >= 0x3FE59428; /* |x| >= 0.6744 */
+    if big {
+        let sign = hx >> 31;
+        if sign != 0 {
+            x = -x;
+            y = -y;
+        }
+        x = (PIO4 - x) + (PIO4_LO - y);
+        y = 0.0;
+    }
+    let z = x * x;
+    let w = z * z;
+    /*
+     * Break x^5*(T[1]+x^2*T[2]+...) into
+     * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
+     * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
+     */
+    let r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11]))));
+    let v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12])))));
+    let s = z * x;
+    let r = y + z * (s * (r + v) + y) + s * T[0];
+    let w = x + r;
+    if big {
+        let sign = hx >> 31;
+        let s = 1.0 - 2.0 * odd as f64;
+        let v = s - 2.0 * (x + (r - w * w / (w + s)));
+        return if sign != 0 { -v } else { v };
+    }
+    if odd == 0 {
+        return w;
+    }
+    /* -1.0/(x+r) has up to 2ulp error, so compute it accurately */
+    let w0 = zero_low_word(w);
+    let v = r - (w0 - x); /* w0+v = r+x */
+    let a = -1.0 / w;
+    let a0 = zero_low_word(a);
+    a0 + a * (1.0 + a0 * w0 + a0 * v)
+}
+
+#[inline]
+fn zero_low_word(x: f64) -> f64 {
+    f64::from_bits(f64::to_bits(x) & 0xFFFF_FFFF_0000_0000)
+}

src/math/mod.rs

@@ -52,6 +52,7 @@ mod sin;
 mod sinf;
 mod sqrt;
 mod sqrtf;
+mod tan;
 mod tanf;
 mod trunc;
 mod truncf;
@@ -102,6 +103,7 @@ pub use self::sin::sin;
 pub use self::sinf::sinf;
 pub use self::sqrt::sqrt;
 pub use self::sqrtf::sqrtf;
+pub use self::tan::tan;
 pub use self::tanf::tanf;
 pub use self::trunc::trunc;
 pub use self::truncf::truncf;
@@ -111,6 +113,7 @@ mod k_cos;
 mod k_cosf;
 mod k_sin;
 mod k_sinf;
+mod k_tan;
 mod k_tanf;
 mod rem_pio2;
 mod rem_pio2_large;
@@ -121,6 +124,7 @@ 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_tan::k_tan;
 use self::k_tanf::k_tanf;
 use self::rem_pio2::rem_pio2;
 use self::rem_pio2_large::rem_pio2_large;

src/math/tan.rs

@@ -0,0 +1,69 @@
+// origin: FreeBSD /usr/src/lib/msun/src/s_tan.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_tan, rem_pio2};
+
+// tan(x)
+// Return tangent function of x.
+//
+// kernel function:
+//      k_tan           ... tangent 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 tan(x: f64) -> f64 {
+    let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
+
+    let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
+    /* |x| ~< pi/4 */
+    if ix <= 0x3fe921fb {
+        if ix < 0x3e400000 {
+            /* |x| < 2**-27 */
+            /* raise inexact if x!=0 and underflow if subnormal */
+            force_eval!(if ix < 0x00100000 {
+                x / x1p120 as f64
+            } else {
+                x + x1p120 as f64
+            });
+            return x;
+        }
+        return k_tan(x, 0.0, 0);
+    }
+
+    /* tan(Inf or NaN) is NaN */
+    if ix >= 0x7ff00000 {
+        return x - x;
+    }
+
+    /* argument reduction */
+    let (n, y0, y1) = rem_pio2(x);
+    k_tan(y0, y1, n & 1)
+}

test-generator/src/main.rs

@@ -716,7 +716,7 @@ f64_f64! {
     sin,
     // sinh,
     sqrt,
-    // tan,
+    tan,
     // tanh,
     trunc,
     fabs,