Merge branch 'master' into master

library/compiler-builtins/libm/src/lib.rs

@@ -384,7 +384,6 @@ pub trait F64Ext: private::Sealed {
 
     fn cos(self) -> Self;
 
-    #[cfg(todo)]
     fn tan(self) -> Self;
 
     #[cfg(todo)]
@@ -408,7 +407,6 @@ pub trait F64Ext: private::Sealed {
 
     fn ln_1p(self) -> Self;
 
-    #[cfg(todo)]
     fn sinh(self) -> Self;
 
     #[cfg(todo)]
@@ -550,7 +548,6 @@ impl F64Ext for f64 {
         cos(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn tan(self) -> Self {
         tan(self)
@@ -589,7 +586,6 @@ impl F64Ext for f64 {
         log1p(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn sinh(self) -> Self {
         sinh(self)

library/compiler-builtins/libm/src/math/expo2.rs

@@ -0,0 +1,13 @@
+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 {
+    /* 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);
+
+    /* note that k is odd and scale*scale overflows */
+    let scale = combine_words(((0x3ff + K / 2) as u32) << 20, 0);
+    /* exp(x - k ln2) * 2**(k-1) */
+    return exp(x - kln2) * scale * scale;
+}

library/compiler-builtins/libm/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)
+}

library/compiler-builtins/libm/src/math/mod.rs

@@ -52,8 +52,10 @@ mod scalbn;
 mod scalbnf;
 mod sin;
 mod sinf;
+mod sinh;
 mod sqrt;
 mod sqrtf;
+mod tan;
 mod tanf;
 mod tanhf;
 mod trunc;
@@ -105,29 +107,36 @@ pub use self::scalbn::scalbn;
 pub use self::scalbnf::scalbnf;
 pub use self::sin::sin;
 pub use self::sinf::sinf;
+pub use self::sinh::sinh;
 pub use self::sqrt::sqrt;
 pub use self::sqrtf::sqrtf;
+pub use self::tan::tan;
 pub use self::tanf::tanf;
 pub use self::tanhf::tanhf;
 pub use self::trunc::trunc;
 pub use self::truncf::truncf;
 
 // Private modules
+mod expo2;
 mod k_cos;
 mod k_cosf;
 mod k_expo2f;
 mod k_sin;
 mod k_sinf;
+mod k_tan;
 mod k_tanf;
 mod rem_pio2;
 mod rem_pio2_large;
 mod rem_pio2f;
 
+// Private re-imports
+use self::expo2::expo2;
 use self::k_cos::k_cos;
 use self::k_cosf::k_cosf;
 use self::k_expo2f::k_expo2f;
 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;
@@ -158,3 +167,8 @@ pub fn with_set_low_word(f: f64, lo: u32) -> f64 {
     tmp |= lo as u64;
     f64::from_bits(tmp)
 }
+
+#[inline]
+fn combine_words(hi: u32, lo: u32) -> f64 {
+    f64::from_bits((hi as u64) << 32 | lo as u64)
+}

library/compiler-builtins/libm/src/math/sinh.rs

@@ -0,0 +1,48 @@
+use super::{expm1, expo2};
+
+// sinh(x) = (exp(x) - 1/exp(x))/2
+//         = (exp(x)-1 + (exp(x)-1)/exp(x))/2
+//         = x + x^3/6 + o(x^5)
+//
+pub fn sinh(x: f64) -> f64 {
+    // union {double f; uint64_t i;} u = {.f = x};
+    // uint32_t w;
+    // double t, h, absx;
+
+    let mut uf: f64 = x;
+    let mut ui: u64 = f64::to_bits(uf);
+    let w: u32;
+    let t: f64;
+    let mut h: f64;
+    let absx: f64;
+
+    h = 0.5;
+    if ui >> 63 != 0 {
+        h = -h;
+    }
+    /* |x| */
+    ui &= !1 / 2;
+    uf = f64::from_bits(ui);
+    absx = uf;
+    w = (ui >> 32) as u32;
+
+    /* |x| < log(DBL_MAX) */
+    if w < 0x40862e42 {
+        t = expm1(absx);
+        if w < 0x3ff00000 {
+            if w < 0x3ff00000 - (26 << 20) {
+                /* note: inexact and underflow are raised by expm1 */
+                /* note: this branch avoids spurious underflow */
+                return x;
+            }
+            return h * (2.0 * t - t * t / (t + 1.0));
+        }
+        /* note: |x|>log(0x1p26)+eps could be just h*exp(x) */
+        return h * (t + t / (t + 1.0));
+    }
+
+    /* |x| > log(DBL_MAX) or nan */
+    /* note: the result is stored to handle overflow */
+    t = 2.0 * h * expo2(absx);
+    return t;
+}

library/compiler-builtins/libm/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)
+}

library/compiler-builtins/libm/test-generator/src/main.rs

@@ -714,9 +714,9 @@ f64_f64! {
     log2,
     round,
     sin,
-    // sinh,
+    sinh,
     sqrt,
-    // tan,
+    tan,
     // tanh,
     trunc,
     fabs,