dashu_float/math/
hyper.rs

1//! Hyperbolic functions, built from the cancellation-free `exp_m1` / `ln_1p` primitives:
2//!
3//! - `sinh(x) = (exp_m1(x) - exp_m1(-x)) / 2`
4//! - `cosh(x) = (exp_m1(x) + exp_m1(-x)) / 2 + 1`
5//! - `tanh(x) = exp_m1(2x) / (exp_m1(2x) + 2)`
6//! - `asinh(x) = sign(x) · ln_1p(|x| + x²/(sqrt(x²+1)+1))`
7//! - `acosh(x) = ln_1p((x-1) + sqrt((x-1)(x+1)))`  (x ≥ 1)
8//! - `atanh(x) = ln_1p(2x/(1-x)) / 2`  (|x| < 1)
9//!
10//! The `exp_m1` / `ln_1p` forms avoid the catastrophic cancellation that the naive
11//! `(exp(x)-exp(-x))/2` and `ln(1+…)` formulas suffer for small arguments. Special
12//! values follow IEEE 754: infinities are values (not errors) for the forward functions
13//! and `asinh`; `acosh(x<1)` and `atanh(|x|>1)` are domain errors.
14
15use crate::{
16    error::{assert_limited_precision, FpError},
17    fbig::FBig,
18    math::{
19        cache::{reborrow_cache, ConstCache},
20        FpResult,
21    },
22    repr::{Context, Repr, Word},
23    round::Round,
24};
25use dashu_base::{Abs, AbsOrd, Approximation::Exact, Sign};
26
27impl<R: Round> Context<R> {
28    /// Hyperbolic sine.
29    pub fn sinh<const B: Word>(
30        &self,
31        x: &Repr<B>,
32        mut cache: Option<&mut ConstCache>,
33    ) -> FpResult<FBig<R, B>> {
34        if x.is_infinite() {
35            return Ok(Exact(FBig::new(Repr::infinity_with_sign(x.sign()), *self)));
36        }
37        assert_limited_precision(self.precision);
38        if x.significand.is_zero() {
39            // sinh(±0) = ±0
40            return Ok(Exact(FBig::new(signed_zero_repr(x), *self)));
41        }
42
43        // sinh(x) = (exp_m1(x) - exp_m1(-x)) / 2  (cancellation-free)
44        let work = Context::<R>::new(self.precision + 50);
45        let x_f = FBig::<R, B>::new(work.repr_round_ref(x).value(), work);
46        let neg_x = -x_f.clone();
47        let ep = work.exp_m1(&x_f.repr, reborrow_cache(&mut cache));
48        let em = work.exp_m1(&neg_x.repr, reborrow_cache(&mut cache));
49        match (ep, em) {
50            (Ok(ep), Ok(em)) => {
51                Ok(((ep.value() - em.value()) / 2i32).with_precision(self.precision))
52            }
53            // |x| large enough that exp_m1 overflowed: sinh(x) → ±inf (sign of x).
54            _ => Err(FpError::Overflow(x.sign())),
55        }
56    }
57
58    /// Hyperbolic cosine.
59    pub fn cosh<const B: Word>(
60        &self,
61        x: &Repr<B>,
62        mut cache: Option<&mut ConstCache>,
63    ) -> FpResult<FBig<R, B>> {
64        if x.is_infinite() {
65            // cosh(±inf) = +inf
66            return Ok(Exact(FBig::new(Repr::infinity(), *self)));
67        }
68        assert_limited_precision(self.precision);
69        if x.significand.is_zero() {
70            // cosh(±0) = 1
71            return Ok(Exact(FBig::new(Repr::one(), *self)));
72        }
73
74        // cosh(x) = (exp_m1(x) + exp_m1(-x)) / 2 + 1  (no cancellation: same-sign sum)
75        let work = Context::<R>::new(self.precision + 50);
76        let x_f = FBig::<R, B>::new(work.repr_round_ref(x).value(), work);
77        let neg_x = -x_f.clone();
78        let ep = work.exp_m1(&x_f.repr, reborrow_cache(&mut cache));
79        let em = work.exp_m1(&neg_x.repr, reborrow_cache(&mut cache));
80        match (ep, em) {
81            (Ok(ep), Ok(em)) => Ok(((ep.value() + em.value()) / 2i32 + FBig::<R, B>::ONE)
82                .with_precision(self.precision)),
83            // cosh(x) ≥ 0 always, so overflow → +inf regardless of x's sign.
84            _ => Err(FpError::Overflow(Sign::Positive)),
85        }
86    }
87
88    /// Simultaneously compute `sinh(x)` and `cosh(x)` (context layer). Returns
89    /// `(sinh_result, cosh_result)` where each is a [`FpResult`].
90    ///
91    /// This is more efficient than calling [`sinh`](Context::sinh) and [`cosh`](Context::cosh)
92    /// separately, since the two share the `exp_m1(±x)` sub-computations.
93    pub fn sinh_cosh<const B: Word>(
94        &self,
95        x: &Repr<B>,
96        mut cache: Option<&mut ConstCache>,
97    ) -> (FpResult<FBig<R, B>>, FpResult<FBig<R, B>>) {
98        if x.is_infinite() {
99            return (
100                Ok(Exact(FBig::new(Repr::infinity_with_sign(x.sign()), *self))),
101                Ok(Exact(FBig::new(Repr::infinity(), *self))),
102            );
103        }
104        assert_limited_precision(self.precision);
105        if x.significand.is_zero() {
106            return (
107                Ok(Exact(FBig::new(signed_zero_repr(x), *self))),
108                Ok(Exact(FBig::new(Repr::one(), *self))),
109            );
110        }
111
112        // sinh = (exp_m1(x) - exp_m1(-x)) / 2;  cosh = (exp_m1(x) + exp_m1(-x)) / 2 + 1
113        let work = Context::<R>::new(self.precision + 50);
114        let x_f = FBig::<R, B>::new(work.repr_round_ref(x).value(), work);
115        let neg_x = -x_f.clone();
116        let ep = work.exp_m1(&x_f.repr, reborrow_cache(&mut cache));
117        let em = work.exp_m1(&neg_x.repr, reborrow_cache(&mut cache));
118        match (ep, em) {
119            (Ok(ep), Ok(em)) => {
120                let ep = ep.value();
121                let em = em.value();
122                let sinh_val = ((ep.clone() - em.clone()) / 2i32).with_precision(self.precision);
123                let cosh_val =
124                    ((ep + em) / 2i32 + FBig::<R, B>::ONE).with_precision(self.precision);
125                (Ok(sinh_val), Ok(cosh_val))
126            }
127            // |x| large enough that exp_m1 overflowed:
128            //   sinh(x) → ±inf (sign of x), cosh(x) → +inf
129            _ => (Err(FpError::Overflow(x.sign())), Err(FpError::Overflow(Sign::Positive))),
130        }
131    }
132
133    /// Hyperbolic tangent.
134    pub fn tanh<const B: Word>(
135        &self,
136        x: &Repr<B>,
137        mut cache: Option<&mut ConstCache>,
138    ) -> FpResult<FBig<R, B>> {
139        if x.is_infinite() {
140            // tanh(±inf) = ±1
141            let one = FBig::new(Repr::one(), *self);
142            return Ok(Exact(if x.sign() == Sign::Negative {
143                -one
144            } else {
145                one
146            }));
147        }
148        assert_limited_precision(self.precision);
149        if x.significand.is_zero() {
150            // tanh(±0) = ±0
151            return Ok(Exact(FBig::new(signed_zero_repr(x), *self)));
152        }
153
154        // tanh(x) = exp_m1(2x) / (exp_m1(2x) + 2). For large negative x, exp_m1(2x) → -1,
155        // giving -1/(-1+2) = -1; for large positive x, exp_m1(2x) overflows, but
156        // tanh(+huge) = 1, so short-circuit before the division would yield +inf/+inf = NaN.
157        let work = Context::<R>::new(self.precision + 50);
158        let x_f = FBig::<R, B>::new(work.repr_round_ref(x).value(), work);
159        let two_x = x_f.clone() * 2i32;
160        match work.exp_m1(&two_x.repr, reborrow_cache(&mut cache)) {
161            Err(FpError::Overflow(_)) => Ok(FBig::ONE.with_precision(self.precision)),
162            Ok(e) => {
163                let e = e.value();
164                Ok((e.clone() / (e.clone() + 2i32)).with_precision(self.precision))
165            }
166            Err(other) => Err(other),
167        }
168    }
169
170    /// Inverse hyperbolic sine.
171    pub fn asinh<const B: Word>(
172        &self,
173        x: &Repr<B>,
174        mut cache: Option<&mut ConstCache>,
175    ) -> FpResult<FBig<R, B>> {
176        if x.is_infinite() {
177            return Ok(Exact(FBig::new(Repr::infinity_with_sign(x.sign()), *self)));
178        }
179        assert_limited_precision(self.precision);
180        if x.significand.is_zero() {
181            // asinh(±0) = ±0
182            return Ok(Exact(FBig::new(signed_zero_repr(x), *self)));
183        }
184
185        // asinh(x) = sign(x) · ln_1p(|x| + x²/(sqrt(x²+1)+1)).
186        // The x²/(sqrt+1) form avoids the `sqrt(x²+1) - 1` cancellation near 0.
187        let work = Context::<R>::new(self.precision + 50);
188        let x_f = FBig::<R, B>::new(work.repr_round_ref(x).value(), work);
189        let sign = x_f.sign();
190        let abs_x = x_f.abs();
191        let arg = match work.sqr(&abs_x.repr) {
192            Ok(x_sq) => {
193                let x_sq = x_sq.value();
194                let sqrt_plus_one = work.sqrt(&(x_sq.clone() + FBig::<R, B>::ONE).repr)?.value()
195                    + FBig::<R, B>::ONE;
196                abs_x.clone() + x_sq / sqrt_plus_one
197            }
198            // |x| so large that x² overflows: asinh(x) ≈ sign·ln(2|x|) (the √(1+1/x²)
199            // correction is far below representable precision here).
200            Err(FpError::Overflow(_)) => {
201                let ln_val = work
202                    .ln(&(abs_x.clone() * 2i32).repr, reborrow_cache(&mut cache))?
203                    .value();
204                return Ok(apply_sign(ln_val, sign).with_precision(self.precision));
205            }
206            Err(other) => return Err(other),
207        };
208        let res = work.ln_1p(&arg.repr, reborrow_cache(&mut cache))?.value();
209        Ok(apply_sign(res, sign).with_precision(self.precision))
210    }
211
212    /// Inverse hyperbolic cosine. Domain: `x ≥ 1`.
213    pub fn acosh<const B: Word>(
214        &self,
215        x: &Repr<B>,
216        mut cache: Option<&mut ConstCache>,
217    ) -> FpResult<FBig<R, B>> {
218        if x.is_infinite() {
219            if x.sign() == Sign::Negative {
220                return Err(FpError::OutOfDomain);
221            }
222            return Ok(Exact(FBig::new(Repr::infinity(), *self)));
223        }
224        assert_limited_precision(self.precision);
225        // domain x ≥ 1 (acosh(1) = 0 is handled below; x < 1 is an error)
226        if x.sign() == Sign::Negative
227            || FBig::<R, B>::new(x.clone(), *self)
228                .abs_cmp(&FBig::ONE)
229                .is_lt()
230        {
231            return Err(FpError::OutOfDomain);
232        }
233        if x.is_one() {
234            return Ok(Exact(FBig::new(Repr::zero(), *self)));
235        }
236
237        // acosh(x) = ln_1p((x-1) + sqrt((x-1)(x+1))). The (x-1)(x+1) form avoids the
238        // `x²-1` cancellation near x = 1.
239        let work = Context::<R>::new(self.precision + 50);
240        let x_f = FBig::<R, B>::new(work.repr_round_ref(x).value(), work);
241        let xm1 = &x_f - FBig::<R, B>::ONE;
242        let xp1 = &x_f + FBig::<R, B>::ONE;
243        let arg = match work.mul(&xm1.repr, &xp1.repr) {
244            Ok(prod) => xm1.clone() + work.sqrt(&prod.value().repr)?.value(),
245            // (x-1)(x+1) overflowed: acosh(x) ≈ ln(2x).
246            Err(FpError::Overflow(_)) => {
247                let ln_val = work
248                    .ln(&(x_f.clone() * 2i32).repr, reborrow_cache(&mut cache))?
249                    .value();
250                return Ok(ln_val.with_precision(self.precision));
251            }
252            Err(other) => return Err(other),
253        };
254        let res = work.ln_1p(&arg.repr, reborrow_cache(&mut cache))?.value();
255        Ok(res.with_precision(self.precision))
256    }
257
258    /// Inverse hyperbolic tangent. Domain: `-1 < x < 1` (`x = ±1` → ±∞, `|x| > 1` is an error).
259    pub fn atanh<const B: Word>(
260        &self,
261        x: &Repr<B>,
262        mut cache: Option<&mut ConstCache>,
263    ) -> FpResult<FBig<R, B>> {
264        if x.is_infinite() {
265            return Err(FpError::OutOfDomain);
266        }
267        assert_limited_precision(self.precision);
268        if x.significand.is_zero() {
269            // atanh(±0) = ±0
270            return Ok(Exact(FBig::new(signed_zero_repr(x), *self)));
271        }
272        // domain |x| < 1: |x| = 1 → ±∞ (value), |x| > 1 → error
273        match FBig::<R, B>::new(x.clone(), *self).abs_cmp(&FBig::ONE) {
274            core::cmp::Ordering::Greater => return Err(FpError::OutOfDomain),
275            core::cmp::Ordering::Equal => {
276                return Ok(Exact(FBig::new(Repr::infinity_with_sign(x.sign()), *self)));
277            }
278            _ => {}
279        }
280
281        // atanh(x) = ln_1p(2x/(1-x)) / 2.
282        let work = Context::<R>::new(self.precision + 50);
283        let x_f = FBig::<R, B>::new(work.repr_round_ref(x).value(), work);
284        let ratio = (x_f.clone() * 2i32) / (FBig::<R, B>::ONE - &x_f);
285        let res = work.ln_1p(&ratio.repr, reborrow_cache(&mut cache))?.value();
286        Ok((res / 2i32).with_precision(self.precision))
287    }
288}
289
290impl<R: Round, const B: Word> FBig<R, B> {
291    /// Calculate the hyperbolic sine of the floating point number.
292    ///
293    /// # Examples
294    ///
295    /// ```
296    /// # use core::str::FromStr;
297    /// # use dashu_base::ParseError;
298    /// # use dashu_float::DBig;
299    /// let a = DBig::from_str("0.5000000")?;
300    /// assert_eq!(a.sinh(), DBig::from_str("0.52109531")?);
301    /// # Ok::<(), ParseError>(())
302    /// ```
303    #[inline]
304    pub fn sinh(&self) -> Self {
305        self.context.unwrap_fp(self.context.sinh(&self.repr, None))
306    }
307
308    /// Calculate the hyperbolic cosine of the floating point number.
309    ///
310    /// # Examples
311    ///
312    /// ```
313    /// # use core::str::FromStr;
314    /// # use dashu_base::ParseError;
315    /// # use dashu_float::DBig;
316    /// let a = DBig::from_str("0.5000000")?;
317    /// assert_eq!(a.cosh(), DBig::from_str("1.127626")?);
318    /// # Ok::<(), ParseError>(())
319    /// ```
320    #[inline]
321    pub fn cosh(&self) -> Self {
322        self.context.unwrap_fp(self.context.cosh(&self.repr, None))
323    }
324
325    /// Simultaneously calculate the hyperbolic sine and cosine of the number.
326    ///
327    /// This is more efficient than calling [`sinh`](FBig::sinh) and [`cosh`](FBig::cosh)
328    /// separately, since the two share the `exp_m1(±x)` sub-computations.
329    ///
330    /// # Examples
331    ///
332    /// ```
333    /// # use core::str::FromStr;
334    /// # use dashu_base::ParseError;
335    /// # use dashu_float::DBig;
336    /// let a = DBig::from_str("0.5000000")?;
337    /// let (s, c) = a.sinh_cosh();
338    /// assert_eq!(s, DBig::from_str("0.52109531")?);
339    /// assert_eq!(c, DBig::from_str("1.127626")?);
340    /// # Ok::<(), ParseError>(())
341    /// ```
342    #[inline]
343    pub fn sinh_cosh(&self) -> (Self, Self) {
344        let (s, c) = self.context.sinh_cosh(&self.repr, None);
345        (self.context.unwrap_fp(s), self.context.unwrap_fp(c))
346    }
347
348    /// Calculate the hyperbolic tangent of the floating point number.
349    ///
350    /// # Examples
351    ///
352    /// ```
353    /// # use core::str::FromStr;
354    /// # use dashu_base::ParseError;
355    /// # use dashu_float::DBig;
356    /// let a = DBig::from_str("0.5000000")?;
357    /// assert_eq!(a.tanh(), DBig::from_str("0.46211716")?);
358    /// # Ok::<(), ParseError>(())
359    /// ```
360    #[inline]
361    pub fn tanh(&self) -> Self {
362        self.context.unwrap_fp(self.context.tanh(&self.repr, None))
363    }
364
365    /// Calculate the inverse hyperbolic sine of the floating point number.
366    ///
367    /// # Examples
368    ///
369    /// ```
370    /// # use core::str::FromStr;
371    /// # use dashu_base::ParseError;
372    /// # use dashu_float::DBig;
373    /// let a = DBig::from_str("0.5000000")?;
374    /// assert_eq!(a.asinh(), DBig::from_str("0.48121183")?);
375    /// # Ok::<(), ParseError>(())
376    /// ```
377    #[inline]
378    pub fn asinh(&self) -> Self {
379        self.context.unwrap_fp(self.context.asinh(&self.repr, None))
380    }
381
382    /// Calculate the inverse hyperbolic cosine of the floating point number.
383    ///
384    /// # Panics
385    ///
386    /// Panics if the number is less than 1 (out of domain).
387    ///
388    /// # Examples
389    ///
390    /// ```
391    /// # use core::str::FromStr;
392    /// # use dashu_base::ParseError;
393    /// # use dashu_float::DBig;
394    /// let a = DBig::from_str("2.000000")?;
395    /// assert_eq!(a.acosh(), DBig::from_str("1.316958")?);
396    /// # Ok::<(), ParseError>(())
397    /// ```
398    #[inline]
399    pub fn acosh(&self) -> Self {
400        self.context.unwrap_fp(self.context.acosh(&self.repr, None))
401    }
402
403    /// Calculate the inverse hyperbolic tangent of the floating point number.
404    ///
405    /// # Panics
406    ///
407    /// Panics if the absolute value is greater than or equal to 1 (out of domain;
408    /// `|x| = 1` is infinite and `|x| > 1` is not real).
409    ///
410    /// # Examples
411    ///
412    /// ```
413    /// # use core::str::FromStr;
414    /// # use dashu_base::ParseError;
415    /// # use dashu_float::DBig;
416    /// let a = DBig::from_str("0.5000000")?;
417    /// assert_eq!(a.atanh(), DBig::from_str("0.54930614")?);
418    /// # Ok::<(), ParseError>(())
419    /// ```
420    #[inline]
421    pub fn atanh(&self) -> Self {
422        self.context.unwrap_fp(self.context.atanh(&self.repr, None))
423    }
424}
425
426/// `±0` `Repr` carrying the sign of `x` (used by the odd hyperbolics at zero input).
427fn signed_zero_repr<const B: Word>(x: &Repr<B>) -> Repr<B> {
428    if x.is_neg_zero() {
429        Repr::neg_zero()
430    } else {
431        Repr::zero()
432    }
433}
434
435/// Negate `v` when `sign` is `Negative` (used to apply `sign(x)` in `asinh`).
436fn apply_sign<R: Round, const B: Word>(v: FBig<R, B>, sign: Sign) -> FBig<R, B> {
437    if sign == Sign::Negative {
438        -v
439    } else {
440        v
441    }
442}
dashu_float/math/hyper.rs

dashu_float/math/
hyper.rs
