Skip to main content

oxinum_complex/native/
trig.rs

1//! Trigonometric and hyperbolic functions for native [`BigComplex`].
2//!
3//! The real [`BigFloat`] foundation provides `sin`, `cos`, `tan`, but no
4//! hyperbolic functions, so the private helpers [`bf_sinh`], [`bf_cosh`], and
5//! [`bf_tanh`] derive them from the exponential:
6//!
7//! ```text
8//! sinh x = (eˣ − e⁻ˣ) / 2
9//! cosh x = (eˣ + e⁻ˣ) / 2
10//! tanh x = sinh x / cosh x
11//! ```
12//!
13//! Halving is done by multiplying with the exact constant `0.5` so the
14//! panicking real `/` operator is never touched; the genuine division in
15//! `tanh` goes through [`BigFloat::div_ref_with_mode`].
16//!
17//! With `z = a + b·i` the complex identities are:
18//!
19//! ```text
20//! sin z  = sin a · cosh b + i · cos a · sinh b
21//! cos z  = cos a · cosh b − i · sin a · sinh b
22//! sinh z = sinh a · cos b + i · cosh a · sin b
23//! cosh z = cosh a · cos b + i · sinh a · sin b
24//! tan z  = sin z / cos z          (via checked_div)
25//! tanh z = sinh z / cosh z        (via checked_div)
26//! ```
27//!
28//! `tan`/`tanh` propagate [`OxiNumError::DivByZero`](oxinum_core::OxiNumError::DivByZero)
29//! at their poles (where the complex denominator collapses to zero).
30
31use oxinum_core::OxiNumResult;
32use oxinum_float::native::{BigFloat, RoundingMode};
33
34use super::BigComplex;
35
36/// Working-precision headroom added on top of the requested `prec`.
37const GUARD: u32 = 10;
38
39// ---------------------------------------------------------------------------
40// Private real hyperbolic helpers, derived from BigFloat::exp.
41// ---------------------------------------------------------------------------
42
43/// `sinh(x) = (eˣ − e⁻ˣ) / 2` at `prec` bits.
44fn bf_sinh(x: &BigFloat, prec: u32, mode: RoundingMode) -> OxiNumResult<BigFloat> {
45    let ex = x.exp(prec, mode)?;
46    let e_neg = (-x).exp(prec, mode)?;
47    let half = BigFloat::from_f64(0.5, prec)?;
48    let diff = &ex - &e_neg;
49    Ok((&diff * &half).with_precision(prec, mode))
50}
51
52/// `cosh(x) = (eˣ + e⁻ˣ) / 2` at `prec` bits.
53fn bf_cosh(x: &BigFloat, prec: u32, mode: RoundingMode) -> OxiNumResult<BigFloat> {
54    let ex = x.exp(prec, mode)?;
55    let e_neg = (-x).exp(prec, mode)?;
56    let half = BigFloat::from_f64(0.5, prec)?;
57    let sum = &ex + &e_neg;
58    Ok((&sum * &half).with_precision(prec, mode))
59}
60
61/// `tanh(x) = sinh(x) / cosh(x)` at `prec` bits.
62///
63/// `cosh(x) ≥ 1` for all real `x`, so the division never hits zero. Used by the
64/// real-axis fast path of the complex [`BigComplex::tanh`] (where `b == 0`,
65/// `tanh z` collapses to the real scalar `tanh a`), keeping that case off the
66/// general complex-division route.
67fn bf_tanh(x: &BigFloat, prec: u32, mode: RoundingMode) -> OxiNumResult<BigFloat> {
68    let sinh = bf_sinh(x, prec, mode)?;
69    let cosh = bf_cosh(x, prec, mode)?;
70    Ok(sinh
71        .div_ref_with_mode(&cosh, mode)?
72        .with_precision(prec, mode))
73}
74
75impl BigComplex {
76    /// The complex sine
77    /// `sin z = sin a · cosh b + i · cos a · sinh b` at `prec` bits.
78    ///
79    /// # Errors
80    ///
81    /// Propagates errors from the real `sin`/`cos`/`exp` routines.
82    ///
83    /// # Examples
84    ///
85    /// ```
86    /// use oxinum_complex::native::{BigComplex, RoundingMode};
87    /// let s = BigComplex::zero(80).sin(80, RoundingMode::HalfEven).expect("sin");
88    /// assert!(s.re().to_f64().abs() < 1e-12);
89    /// assert!(s.im().to_f64().abs() < 1e-12);
90    /// ```
91    pub fn sin(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
92        let guard = prec.saturating_add(GUARD);
93        let a = self.re.clone().with_precision(guard, mode);
94        let b = self.im.clone().with_precision(guard, mode);
95
96        let sin_a = a.sin(guard, mode)?;
97        let cos_a = a.cos(guard, mode)?;
98        let cosh_b = bf_cosh(&b, guard, mode)?;
99        let sinh_b = bf_sinh(&b, guard, mode)?;
100
101        let re = (&sin_a * &cosh_b).with_precision(prec, mode);
102        let im = (&cos_a * &sinh_b).with_precision(prec, mode);
103        Ok(BigComplex { re, im })
104    }
105
106    /// The complex cosine
107    /// `cos z = cos a · cosh b − i · sin a · sinh b` at `prec` bits.
108    ///
109    /// # Errors
110    ///
111    /// Propagates errors from the real `sin`/`cos`/`exp` routines.
112    ///
113    /// # Examples
114    ///
115    /// ```
116    /// use oxinum_complex::native::{BigComplex, RoundingMode};
117    /// let c = BigComplex::zero(80).cos(80, RoundingMode::HalfEven).expect("cos");
118    /// assert!((c.re().to_f64() - 1.0).abs() < 1e-12);
119    /// assert!(c.im().to_f64().abs() < 1e-12);
120    /// ```
121    pub fn cos(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
122        let guard = prec.saturating_add(GUARD);
123        let a = self.re.clone().with_precision(guard, mode);
124        let b = self.im.clone().with_precision(guard, mode);
125
126        let cos_a = a.cos(guard, mode)?;
127        let sin_a = a.sin(guard, mode)?;
128        let cosh_b = bf_cosh(&b, guard, mode)?;
129        let sinh_b = bf_sinh(&b, guard, mode)?;
130
131        let re = (&cos_a * &cosh_b).with_precision(prec, mode);
132        // im = −(sin a · sinh b)
133        let prod = &sin_a * &sinh_b;
134        let im = (-&prod).with_precision(prec, mode);
135        Ok(BigComplex { re, im })
136    }
137
138    /// The complex hyperbolic sine
139    /// `sinh z = sinh a · cos b + i · cosh a · sin b` at `prec` bits.
140    ///
141    /// # Errors
142    ///
143    /// Propagates errors from the real `sin`/`cos`/`exp` routines.
144    pub fn sinh(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
145        let guard = prec.saturating_add(GUARD);
146        let a = self.re.clone().with_precision(guard, mode);
147        let b = self.im.clone().with_precision(guard, mode);
148
149        let sinh_a = bf_sinh(&a, guard, mode)?;
150        let cosh_a = bf_cosh(&a, guard, mode)?;
151        let cos_b = b.cos(guard, mode)?;
152        let sin_b = b.sin(guard, mode)?;
153
154        let re = (&sinh_a * &cos_b).with_precision(prec, mode);
155        let im = (&cosh_a * &sin_b).with_precision(prec, mode);
156        Ok(BigComplex { re, im })
157    }
158
159    /// The complex hyperbolic cosine
160    /// `cosh z = cosh a · cos b + i · sinh a · sin b` at `prec` bits.
161    ///
162    /// # Errors
163    ///
164    /// Propagates errors from the real `sin`/`cos`/`exp` routines.
165    ///
166    /// # Examples
167    ///
168    /// ```
169    /// use oxinum_complex::native::{BigComplex, RoundingMode};
170    /// let c = BigComplex::zero(80).cosh(80, RoundingMode::HalfEven).expect("cosh");
171    /// assert!((c.re().to_f64() - 1.0).abs() < 1e-12);
172    /// assert!(c.im().to_f64().abs() < 1e-12);
173    /// ```
174    pub fn cosh(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
175        let guard = prec.saturating_add(GUARD);
176        let a = self.re.clone().with_precision(guard, mode);
177        let b = self.im.clone().with_precision(guard, mode);
178
179        let cosh_a = bf_cosh(&a, guard, mode)?;
180        let sinh_a = bf_sinh(&a, guard, mode)?;
181        let cos_b = b.cos(guard, mode)?;
182        let sin_b = b.sin(guard, mode)?;
183
184        let re = (&cosh_a * &cos_b).with_precision(prec, mode);
185        let im = (&sinh_a * &sin_b).with_precision(prec, mode);
186        Ok(BigComplex { re, im })
187    }
188
189    /// The complex tangent `tan z = sin z / cos z` at `prec` bits.
190    ///
191    /// # Errors
192    ///
193    /// - [`OxiNumError::DivByZero`](oxinum_core::OxiNumError::DivByZero) at the
194    ///   poles where `cos z = 0`.
195    /// - Propagates errors from the real `sin`/`cos`/`exp` routines.
196    pub fn tan(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
197        let guard = prec.saturating_add(GUARD);
198        let sin_z = self.sin(guard, mode)?;
199        let cos_z = self.cos(guard, mode)?;
200        sin_z.checked_div(&cos_z, prec, mode)
201    }
202
203    /// The complex hyperbolic tangent `tanh z = sinh z / cosh z` at `prec` bits.
204    ///
205    /// # Errors
206    ///
207    /// - [`OxiNumError::DivByZero`](oxinum_core::OxiNumError::DivByZero) at the
208    ///   poles where `cosh z = 0`.
209    /// - Propagates errors from the real `sin`/`cos`/`exp` routines.
210    pub fn tanh(&self, prec: u32, mode: RoundingMode) -> OxiNumResult<BigComplex> {
211        let guard = prec.saturating_add(GUARD);
212
213        // Real-axis fast path: b == 0 ⇒ tanh z = tanh(a) (purely real).
214        // Routes through the scalar `bf_tanh` instead of the general complex
215        // division, avoiding any spurious imaginary component from rounding.
216        if self.im.is_zero() {
217            let a = self.re.clone().with_precision(guard, mode);
218            let t = bf_tanh(&a, guard, mode)?.with_precision(prec, mode);
219            return Ok(BigComplex {
220                re: t,
221                im: BigFloat::zero(prec),
222            });
223        }
224
225        let sinh_z = self.sinh(guard, mode)?;
226        let cosh_z = self.cosh(guard, mode)?;
227        sinh_z.checked_div(&cosh_z, prec, mode)
228    }
229}
230
231// ---------------------------------------------------------------------------
232// Tests
233// ---------------------------------------------------------------------------
234
235#[cfg(test)]
236mod tests {
237    use super::*;
238
239    const PREC: u32 = 80;
240    const MODE: RoundingMode = RoundingMode::HalfEven;
241
242    fn c(re: f64, im: f64) -> BigComplex {
243        BigComplex::from_f64(re, im, PREC).expect("finite parts")
244    }
245
246    #[test]
247    fn sin_zero_is_zero() {
248        let s = BigComplex::zero(PREC).sin(PREC, MODE).expect("sin");
249        assert!(s.re().to_f64().abs() < 1e-12, "re = {}", s.re().to_f64());
250        assert!(s.im().to_f64().abs() < 1e-12, "im = {}", s.im().to_f64());
251    }
252
253    #[test]
254    fn cos_zero_is_one() {
255        let cz = BigComplex::zero(PREC).cos(PREC, MODE).expect("cos");
256        assert!(
257            (cz.re().to_f64() - 1.0).abs() < 1e-12,
258            "re = {}",
259            cz.re().to_f64()
260        );
261        assert!(cz.im().to_f64().abs() < 1e-12, "im = {}", cz.im().to_f64());
262    }
263
264    #[test]
265    fn cosh_zero_is_one() {
266        let cz = BigComplex::zero(PREC).cosh(PREC, MODE).expect("cosh");
267        assert!(
268            (cz.re().to_f64() - 1.0).abs() < 1e-12,
269            "re = {}",
270            cz.re().to_f64()
271        );
272        assert!(cz.im().to_f64().abs() < 1e-12, "im = {}", cz.im().to_f64());
273    }
274
275    #[test]
276    fn sinh_zero_is_zero() {
277        let s = BigComplex::zero(PREC).sinh(PREC, MODE).expect("sinh");
278        assert!(s.re().to_f64().abs() < 1e-12);
279        assert!(s.im().to_f64().abs() < 1e-12);
280    }
281
282    #[test]
283    fn pythagorean_identity_general() {
284        // sin²z + cos²z ≈ 1 at z = 0.5 + 0.3i.
285        let z = c(0.5, 0.3);
286        let s = z.sin(PREC, MODE).expect("sin");
287        let co = z.cos(PREC, MODE).expect("cos");
288        let s2 = s.mul_core(&s);
289        let c2 = co.mul_core(&co);
290        let sum = &s2 + &c2;
291        assert!(
292            (sum.re().to_f64() - 1.0).abs() < 1e-9,
293            "re(sum) = {}",
294            sum.re().to_f64()
295        );
296        assert!(
297            sum.im().to_f64().abs() < 1e-9,
298            "im(sum) = {}",
299            sum.im().to_f64()
300        );
301    }
302
303    #[test]
304    fn tan_zero_is_zero() {
305        let t = BigComplex::zero(PREC).tan(PREC, MODE).expect("tan");
306        assert!(t.re().to_f64().abs() < 1e-12);
307        assert!(t.im().to_f64().abs() < 1e-12);
308    }
309
310    #[test]
311    fn tanh_zero_is_zero() {
312        let t = BigComplex::zero(PREC).tanh(PREC, MODE).expect("tanh");
313        assert!(t.re().to_f64().abs() < 1e-12);
314        assert!(t.im().to_f64().abs() < 1e-12);
315    }
316
317    #[test]
318    fn tan_matches_real_axis() {
319        // tan(0.4 + 0i) ≈ tan(0.4) real.
320        let t = c(0.4, 0.0).tan(PREC, MODE).expect("tan");
321        let expected = 0.4_f64.tan();
322        assert!(
323            (t.re().to_f64() - expected).abs() < 1e-9,
324            "re = {}",
325            t.re().to_f64()
326        );
327        assert!(t.im().to_f64().abs() < 1e-9, "im = {}", t.im().to_f64());
328    }
329}