From dde8131c1e8ecb9c55d1b5f6888ce70d55119cd3 Mon Sep 17 00:00:00 2001
From: C Jones <code@calebjones.net>
Date: Sat, 14 Jul 2018 19:23:01 -0400
Subject: [PATCH] Implement sinh

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

Closes rust-lang/libm#35
---
 library/compiler-builtins/libm/src/lib.rs     |  2 -
 .../compiler-builtins/libm/src/math/expo2.rs  | 13 +++++
 .../compiler-builtins/libm/src/math/mod.rs    |  9 ++++
 .../compiler-builtins/libm/src/math/sinh.rs   | 48 +++++++++++++++++++
 .../libm/test-generator/src/main.rs           |  2 +-
 5 files changed, 71 insertions(+), 3 deletions(-)
 create mode 100644 library/compiler-builtins/libm/src/math/expo2.rs
 create mode 100644 library/compiler-builtins/libm/src/math/sinh.rs

diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs
index 0b9efeeb338a..ee864a36cf2e 100644
--- a/library/compiler-builtins/libm/src/lib.rs
+++ b/library/compiler-builtins/libm/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)
diff --git a/library/compiler-builtins/libm/src/math/expo2.rs b/library/compiler-builtins/libm/src/math/expo2.rs
new file mode 100644
index 000000000000..8a121b0a287a
--- /dev/null
+++ b/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;
+}
diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs
index 7b3d9abee697..19d8ad58d866 100644
--- a/library/compiler-builtins/libm/src/math/mod.rs
+++ b/library/compiler-builtins/libm/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)
+}
diff --git a/library/compiler-builtins/libm/src/math/sinh.rs b/library/compiler-builtins/libm/src/math/sinh.rs
new file mode 100644
index 000000000000..0579871719a8
--- /dev/null
+++ b/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;
+}
diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs
index c8e35cea7b00..758db9f527fb 100644
--- a/library/compiler-builtins/libm/test-generator/src/main.rs
+++ b/library/compiler-builtins/libm/test-generator/src/main.rs
@@ -714,7 +714,7 @@ f64_f64! {
     log2,
     round,
     sin,
-    // sinh,
+    sinh,
     sqrt,
     // tan,
     // tanh,