rustymath_core/
hyperbolic.rs

1//! The hyperbolic functions and their inverses. The functions build on `exp` and `expm1`, the inverses on `ln`'s reduction. Each computes at |x| and restores the sign of x where the function is odd.
2use crate::{float::Native, ln::ln_parts};
3
4/// sinh |x| = h(1 + 1/(2h + 1)), with h = (e^|x| - 1)/2. The halving joins `expm1`'s scaling, so sinh stays finite where e^|x| would overflow. Where h is ∞, 1/(2h + 1) is 0, and the result is ∞.
5#[inline]
6pub(crate) fn sinh<R: Native, const FUSED: bool>(x: R) -> R { sinh_formula(&crate::recipe::NativeOps::<FUSED>::DEFAULT, x) }
7#[inline]
8pub(crate) fn cosh<R: Native, const FUSED: bool>(x: R) -> R { cosh_formula(&crate::recipe::NativeOps::<FUSED>::DEFAULT, x) }
9#[inline]
10pub(crate) fn tanh<R: Native, const FUSED: bool>(x: R) -> R { tanh_formula(&crate::recipe::NativeOps::<FUSED>::DEFAULT, x) }
11
12#[inline]
13pub(crate) fn sinh_formula<F: crate::Float, O: crate::recipe::Ops<F> + ?Sized>(o: &O, x: O::V) -> O::V {
14    let one = o.number(F::one());
15    let h = crate::exp::minus_one_formula(o, o.magnitude(x), o.negate(one));
16    let result = o.multiply(h, o.add(one, o.divide(one, o.add(o.add(h, h), one))));
17    o.select(o.less(x, o.number(F::zero())), o.negate(result), result)
18}
19#[inline]
20pub(crate) fn cosh_formula<F: crate::Float, O: crate::recipe::Ops<F> + ?Sized>(o: &O, x: O::V) -> O::V {
21    let e = crate::exp::scaled_formula(o, o.magnitude(x), o.number(-F::one()));
22    o.add(e, o.divide(o.number(F::narrow(0.25)), e))
23}
24#[inline]
25pub(crate) fn tanh_formula<F: crate::Float, O: crate::recipe::Ops<F> + ?Sized>(o: &O, x: O::V) -> O::V {
26    let top = o.number(F::narrow(if o.width() == crate::buffer::FloatWidth::F64 { 20.0 } else { 9.0 }));
27    let a = o.magnitude(x);
28    let a = o.select(o.less(top, a), top, a);
29    let e = crate::exp::minus_one_formula(o, o.add(a, a), o.number(F::zero()));
30    let result = o.divide(e, o.add(e, o.number(F::narrow(2.0))));
31    o.select(o.less(x, o.number(F::zero())), o.negate(result), result)
32}
33
34/// Beyond this magnitude, `asinh` and `acosh` take v = x in place of x + √(x² ± 1), which rounds to 2x and may overflow, and add 1 to e. x² hasn't yet overflowed below it at either width.
35#[inline]
36fn large<R: Native>() -> R { R::narrow(4503599627370496.0) }
37
38/// ln v for v = m × 2^e, with e raised by 1 when `far`. While e is 0, `near` gives the quotient s = (m-1)/(m+1) as (num, den) in a form that keeps its accuracy. An infinite or NaN v gives itself. Every caller's v is 1 or more, ∞, NaN or negative, so none needs the subnormal step of `split_positive`.
39#[inline]
40fn ln_of<R: Native, const FUSED: bool>(v: R, far: bool, near: (R, R)) -> R {
41    let (m, e) = v.split();
42    let e = if far { e + R::one() } else { e };
43    let (num, den) = if e == R::zero() { near } else { (m - R::one(), m + R::one()) };
44    let y = ln_parts::<R, FUSED>(num / den, e);
45    if v.is_nan() | (v == R::infinity()) { v } else { y }
46}
47
48/// asinh |x| = ln v with v = |x| + r and r = √(x² + 1). Near 0, s = (v-1)/(v+1) = (|x|(1+r) + x²)/((2+|x|)(1+r) + x²).
49#[inline]
50pub(crate) fn asinh<R: Native, const FUSED: bool>(x: R) -> R {
51    let a = x.abs();
52    let far = a > large();
53    let r = a.mla::<FUSED>(a, R::one()).sqrt();
54    let (q, a2) = (R::one() + r, a * a);
55    let near = (a.mla::<FUSED>(q, a2), (a + R::narrow(2.0)).mla::<FUSED>(q, a2));
56    ln_of::<R, FUSED>(if far { a } else { a + r }, far, near).copysign(x)
57}
58
59/// acosh x = ln v with v = x + √(x² - 1). Near 1, s = t/(2 + t) with t = v - 1 = d + √(d² + 2d), where d = x - 1 is exact. x below 1 gives NaN.
60#[inline]
61pub(crate) fn acosh<R: Native, const FUSED: bool>(x: R) -> R {
62    let far = x > large();
63    let d = x - R::one();
64    let t = d + d.mla::<FUSED>(d, d + d).sqrt();
65    let y = ln_of::<R, FUSED>(if far { x } else { R::one() + t }, far, (t, t + R::narrow(2.0)));
66    if x < R::one() { R::nan() } else { y }
67}
68
69/// atanh |x| = ln v / 2 with v = (1 + |x|)/(1 - |x|). Near 0, s = (v-1)/(v+1) is |x| itself. ±1 gives ±∞, and a magnitude beyond 1 gives NaN.
70#[inline]
71pub(crate) fn atanh<R: Native, const FUSED: bool>(x: R) -> R {
72    let a = x.abs();
73    let y = ln_of::<R, FUSED>((R::one() + a) / (R::one() - a), false, (a, R::one())) * R::narrow(0.5);
74    (if a <= R::one() { y } else { R::nan() }).copysign(x)
75}
76
77/// The real part of acosh x, or of its principal value at a complex z, as a function object: acosh |x| for |x| at least 1, and 0 between -1 and 1. The real part of `6○⁻¹`.
78#[derive(Clone, Copy, Debug)]
79pub struct RealAcosh;
80
81/// The real part of atanh x, or of its principal value at a complex z, as a function object: atanh x for x between -1 and 1, and atanh (1 ÷ x) beyond, which is ½ ln |(1 + x) ÷ (1 - x)|. The real part of `7○⁻¹`.
82#[derive(Clone, Copy, Debug)]
83pub struct RealAtanh;
rustymath_core/hyperbolic.rs

rustymath_core/
hyperbolic.rs
