Merge #109

109: Implement sinh r=japaric a=porglezomp

This also adds `expo2` for the `__expo2` function, and `combine_words` to replace
the `INSERT_WORDS` macro.

Closes #35

Co-authored-by: C Jones <code@calebjones.net>

src/lib.rs

@@ -414,7 +414,6 @@ pub trait F64Ext: private::Sealed {
 
     fn ln_1p(self) -> Self;
 
-    #[cfg(todo)]
     fn sinh(self) -> Self;
 
     #[cfg(todo)]
@@ -595,7 +594,6 @@ impl F64Ext for f64 {
         log1p(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn sinh(self) -> Self {
         sinh(self)

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

src/math/mod.rs

@@ -50,6 +50,7 @@ mod scalbn;
 mod scalbnf;
 mod sin;
 mod sinf;
+mod sinh;
 mod sqrt;
 mod sqrtf;
 mod tanf;
@@ -100,6 +101,7 @@ 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::tanf::tanf;
@@ -107,6 +109,7 @@ pub use self::trunc::trunc;
 pub use self::truncf::truncf;
 
 // Private modules
+mod expo2;
 mod k_cos;
 mod k_cosf;
 mod k_sin;
@@ -117,6 +120,7 @@ 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_sin::k_sin;
@@ -151,3 +155,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)
+}

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

test-generator/src/main.rs

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