// Method : 
// mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
// 1. Replace x by |x| (cosh(x) = cosh(-x)). 
// 2. 
// 	                                        [ exp(x) - 1 ]^2 
//     0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
// 		       			           2*exp(x)
// 
// 	                                  exp(x) +  1/exp(x)
//     ln2/2    <= x <= 22     :  cosh(x) := -------------------
// 		       			          2
//     22       <= x <= lnovft :  cosh(x) := exp(x)/2 
//     lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
//     ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
// 
// Special cases:
// cosh(x) is |x| if x is +INF, -INF, or NaN.
// only cosh(0)=1 is exact for finite x.

///|
/// Compute hyperbolic cosine function of a double-precision floating point number.
///
/// # Examples
///
/// ```moonbit nocheck
/// assert_eq(cosh(0), 1)
/// assert_eq(cosh(1), 1.5430806348152437)
/// assert_eq(cosh(2), 3.7621956910836314)
/// assert_eq(cosh(3), 10.067661995777765)
/// assert_eq(cosh(-1), 1.5430806348152437)
/// ```
///
/// # Special Cases
///
/// 1. cosh(x) is |x| if x is +INF, -INF, or NaN.
/// 2. Only cosh(0)=1 is exact for finite x.
///
/// # Accuracy
///
/// 0 ulp.
pub fn cosh(x : Double) -> Double {
  if isnan(x) {
    return x
  }
  if isinf(x) {
    return @double.infinity
  }
  let ix = __hi(x).reinterpret_as_int() & 0x7fffffff
  if ix < 0x3fd62e43 {
    let t = expm1(fabs(x))
    let w = 1.0 + t
    if ix < 0x3c800000 {
      return w
    }
    return 1.0 + t * t / (w + w)
  }
  if ix < 0x40360000 {
    let t = exp(fabs(x))
    return 0.5 * t + 0.5 / t
  }
  if ix < 0x40862E42 {
    return 0.5 * exp(fabs(x))
  }
  let lx = __low(x).reinterpret_as_int()
  if ix < 0x408633ce || (ix == 0x408633ce && lx <= 0x8fb9f87d) {
    let w = exp(0.5 * fabs(x))
    let t = 0.5 * w
    return t * w
  }
  @double.infinity
}

///|
test "cosh" {
  fn assert_cosh_ulp(input, expect) raise {
    assert_ulp(expect, cosh(input), COSH_MAX_ULP)
  }

  assert_cosh_ulp(-0.8, 1.3374349463048447)
  assert_cosh_ulp(-0.7, 1.255169005630943)
  assert_cosh_ulp(-0.6, 1.1854652182422676)
  assert_cosh_ulp(-0.5, 1.1276259652063807)
  assert_cosh_ulp(-0.4, 1.081072371838455)
  assert_cosh_ulp(-0.3, 1.0453385141288605)
  assert_cosh_ulp(-0.2, 1.020066755619076)
  assert_cosh_ulp(-0.1, 1.0050041680558035)
  assert_cosh_ulp(-0, 1)
  assert_cosh_ulp(-3.141592653589793, 11.591953275521519)
  assert_cosh_ulp(-1.5707963267948966, 2.5091784786580567)
  assert_cosh_ulp(-0.7853981633974483, 1.3246090892520057)
  assert_cosh_ulp(0, 1)
  assert_cosh_ulp(0.1, 1.0050041680558035)
  assert_cosh_ulp(0.2, 1.020066755619076)
  assert_cosh_ulp(0.3, 1.0453385141288605)
  assert_cosh_ulp(0.4, 1.081072371838455)
  assert_cosh_ulp(0.5, 1.1276259652063807)
  assert_cosh_ulp(0.6, 1.1854652182422676)
  assert_cosh_ulp(0.7, 1.255169005630943)
  assert_cosh_ulp(0.8, 1.3374349463048447)
  assert_cosh_ulp(0.9, 1.4330863854487745)
  assert_cosh_ulp(1, 1.5430806348152437)
  assert_cosh_ulp(3.141592653589793, 11.591953275521519)
  assert_cosh_ulp(1.5707963267948966, 2.5091784786580567)
  assert_cosh_ulp(0.7853981633974483, 1.3246090892520057)
  assert_cosh_ulp(-1, 1.5430806348152437)
  assert_cosh_ulp(-2, 3.7621956910836314)
  assert_cosh_ulp(-3, 10.067661995777765)
  assert_cosh_ulp(-4, 27.308232836016487)
  assert_cosh_ulp(-5, 74.20994852478785)
  assert_cosh_ulp(-6, 201.7156361224559)
  assert_cosh_ulp(-7, 548.317035155212)
  assert_cosh_ulp(-8, 1490.479161252178)
  assert_cosh_ulp(-9, 4051.5420254925943)
  assert_cosh_ulp(1, 1.5430806348152437)
  assert_cosh_ulp(2, 3.7621956910836314)
  assert_cosh_ulp(3, 10.067661995777765)
  assert_cosh_ulp(4, 27.308232836016487)
  assert_cosh_ulp(5, 74.20994852478785)
  assert_cosh_ulp(6, 201.7156361224559)
  assert_cosh_ulp(7, 548.317035155212)
  assert_cosh_ulp(8, 1490.479161252178)
  assert_cosh_ulp(9, 4051.5420254925943)
  assert_cosh_ulp(10, 11013.232920103324)
  assert_cosh_ulp(100, 13440585709080678000000000000000000000000000)
  assert_cosh_ulp(1000, @double.infinity)
  assert_cosh_ulp(10000, @double.infinity)
  assert_cosh_ulp(2.5, 6.132289479663686)
  assert_cosh_ulp(3.4, 14.998736658678668)
  assert_cosh_ulp(5.3, 100.17090078434929)
  assert_cosh_ulp(6.2, 246.3755352619463)
  assert_cosh_ulp(7.1, 605.9839497987498)
  assert_cosh_ulp(8.9, 3665.986837772461)
  assert_cosh_ulp(9.8, 9016.872491640062)
  assert_cosh_ulp(10.7, 22177.92757642139)
  assert_cosh_ulp(101.6, 66571656819937820000000000000000000000000000)
  assert_cosh_ulp(1.542, 2.443940777316919)
  assert_cosh_ulp(2.846, 8.638422496385937)
  assert_cosh_ulp(7.881, 1323.2590662471498)
  assert_cosh_ulp(3.772, 21.74495889437803)
  assert_cosh_ulp(-1.542, 2.443940777316919)
  assert_cosh_ulp(-2.846, 8.638422496385937)
  assert_cosh_ulp(-7.881, 1323.2590662471498)
  assert_cosh_ulp(-3.772, 21.74495889437803)
  assert_cosh_ulp(-1, 1.5430806348152437)
  assert_cosh_ulp(0, 1)
  assert_cosh_ulp(-0, 1)
  assert_cosh_ulp(@double.not_a_number, @double.not_a_number)
  assert_cosh_ulp(@double.infinity, @double.infinity)
  assert_cosh_ulp(@double.neg_infinity, @double.infinity)
}