dashu_cmplx/math/
trig.rs

1//! Complex trigonometric functions via the real–imaginary decomposition, reusing `dashu-float`'s
2//! real `sin`/`cos` and cancellation-free `sinh`/`cosh`.
3//!
4//! `sin(x+iy) = sin x·cosh y + i·cos x·sinh y`, `cos(x+iy) = cos x·cosh y − i·sin x·sinh y`. This
5//! form avoids the `exp(±iz)` identity's exponential blow-up for large `|Im z|`.
6
7use crate::cbig::CBig;
8use crate::repr::{combine_parts, reborrow_cache, CfpResult, Context};
9use dashu_float::round::Round;
10use dashu_float::{ConstCache, FBig, FpError, Repr};
11use dashu_int::{IBig, Word};
12
13/// Guard digits (base-B) for the forward trig. Composes real `sin_cos` + `sinh_cosh` + two
14/// products; the cancellation near the trig zeros is absorbed by the re-round.
15const TRIG_GUARD: usize = 16;
16
17impl<R: Round> Context<R> {
18    /// Simultaneously compute `sin z` and `cos z` (context layer). Returns `(sin, cos)` each as a
19    /// [`CfpResult`]. An infinite input maps to [`FpError::Indeterminate`] (the C99 NaN cases).
20    pub fn sin_cos<const B: Word>(
21        &self,
22        z: &CBig<R, B>,
23        mut cache: Option<&mut ConstCache>,
24    ) -> (CfpResult<R, B>, CfpResult<R, B>) {
25        if z.is_infinite() {
26            return (Err(FpError::Indeterminate), Err(FpError::Indeterminate));
27        }
28        if z.is_zero() {
29            let zero = Ok(crate::repr::exact(
30                FBig::from_repr(Repr::zero(), self.float()),
31                FBig::from_repr(Repr::zero(), self.float()),
32            ));
33            let one = Ok(crate::repr::exact(
34                FBig::from_repr(Repr::one(), self.float()),
35                FBig::from_repr(Repr::zero(), self.float()),
36            ));
37            return (zero, one);
38        }
39
40        let gctx = self.guard(TRIG_GUARD);
41        let p = self.precision();
42        let (sinx, cosx) = gctx.sin_cos(z.re(), reborrow_cache(&mut cache));
43        let sinx = match sinx {
44            Ok(v) => v.value(),
45            Err(e) => return (Err(e), Err(FpError::Indeterminate)),
46        };
47        let cosx = match cosx {
48            Ok(v) => v.value(),
49            Err(e) => return (Err(FpError::Indeterminate), Err(e)),
50        };
51        let (sinhy_res, coshy_res) = gctx.sinh_cosh(z.im(), reborrow_cache(&mut cache));
52        let sinhy = match sinhy_res {
53            Ok(v) => v.value(),
54            Err(e) => return (Err(e), Err(FpError::Indeterminate)),
55        };
56        let coshy = match coshy_res {
57            Ok(v) => v.value(),
58            Err(e) => return (Err(FpError::Indeterminate), Err(e)),
59        };
60
61        // sin z = (sinx·coshy) + i·(cosx·sinhy); cos z = (cosx·coshy) − i·(sinx·sinhy).
62        // `sin_cos` returns a tuple, so the products are matched explicitly (no `?`).
63        let prod = |a: &FBig<R, B>, b: &FBig<R, B>| -> Result<_, FpError> {
64            Ok(gctx.mul(a.repr(), b.repr())?.value().with_precision(p))
65        };
66        let sin_re = match prod(&sinx, &coshy) {
67            Ok(v) => v,
68            Err(e) => return (Err(e), Err(FpError::Indeterminate)),
69        };
70        let sin_im = match prod(&cosx, &sinhy) {
71            Ok(v) => v,
72            Err(e) => return (Err(e), Err(FpError::Indeterminate)),
73        };
74        let cos_re = match prod(&cosx, &coshy) {
75            Ok(v) => v,
76            Err(e) => return (Err(FpError::Indeterminate), Err(e)),
77        };
78        let neg_sinx = -sinx;
79        let cos_im = match prod(&neg_sinx, &sinhy) {
80            Ok(v) => v,
81            Err(e) => return (Err(FpError::Indeterminate), Err(e)),
82        };
83        (Ok(combine_parts(sin_re, sin_im)), Ok(combine_parts(cos_re, cos_im)))
84    }
85
86    /// Complex sine (context layer).
87    #[inline]
88    pub fn sin<const B: Word>(
89        &self,
90        z: &CBig<R, B>,
91        cache: Option<&mut ConstCache>,
92    ) -> CfpResult<R, B> {
93        self.sin_cos(z, cache).0
94    }
95
96    /// Complex cosine (context layer).
97    #[inline]
98    pub fn cos<const B: Word>(
99        &self,
100        z: &CBig<R, B>,
101        cache: Option<&mut ConstCache>,
102    ) -> CfpResult<R, B> {
103        self.sin_cos(z, cache).1
104    }
105
106    /// Complex tangent `sin z / cos z` (context layer).
107    pub fn tan<const B: Word>(
108        &self,
109        z: &CBig<R, B>,
110        cache: Option<&mut ConstCache>,
111    ) -> CfpResult<R, B> {
112        let (sin_z, cos_z) = self.sin_cos(z, cache);
113        let sin_z = sin_z?;
114        let cos_z = cos_z?;
115        self.div(&sin_z.value(), &cos_z.value())
116    }
117
118    /// Inverse sine `asin z = -i·log(iz + sqrt(1-z²))` (context layer, Kahan form). The argument of
119    /// the inner `log` always has positive real part, so the branch cut comes entirely from the
120    /// `sqrt`; an infinite input maps to [`FpError::Indeterminate`].
121    pub fn asin<const B: Word>(
122        &self,
123        z: &CBig<R, B>,
124        mut cache: Option<&mut ConstCache>,
125    ) -> CfpResult<R, B> {
126        if z.is_infinite() {
127            return Err(FpError::Indeterminate);
128        }
129        let gctx = Context::new(self.precision() + ITRIG_GUARD);
130        let p = self.precision();
131        let one = CBig::ONE;
132        let z2 = gctx.sqr(z)?.value();
133        let one_m_z2 = gctx.sub(&one, &z2)?.value();
134        let sqrt_term = gctx.sqrt(&one_m_z2)?.value();
135        let iz = z.mul_i(false); // exact rotation
136        let w = gctx.add(&iz, &sqrt_term)?.value();
137        let log_w = gctx.log(&w, reborrow_cache(&mut cache))?.value();
138        let asin_z = log_w.mul_i(true); // -i·log(w)
139        let (re, im) = asin_z.into_parts();
140        Ok(combine_parts(re.with_precision(p), im.with_precision(p)))
141    }
142
143    /// Inverse cosine `acos z = -i·log(z + i·sqrt(1-z²))` (context layer, Kahan form).
144    pub fn acos<const B: Word>(
145        &self,
146        z: &CBig<R, B>,
147        mut cache: Option<&mut ConstCache>,
148    ) -> CfpResult<R, B> {
149        if z.is_infinite() {
150            return Err(FpError::Indeterminate);
151        }
152        let gctx = Context::new(self.precision() + ITRIG_GUARD);
153        let p = self.precision();
154        let one = CBig::ONE;
155        let z2 = gctx.sqr(z)?.value();
156        let one_m_z2 = gctx.sub(&one, &z2)?.value();
157        let sqrt_term = gctx.sqrt(&one_m_z2)?.value();
158        let i_sqrt = sqrt_term.mul_i(false); // i·sqrt(1-z²)
159        let w = gctx.add(z, &i_sqrt)?.value();
160        let log_w = gctx.log(&w, reborrow_cache(&mut cache))?.value();
161        let acos_z = log_w.mul_i(true); // -i·log(w)
162        let (re, im) = acos_z.into_parts();
163        Ok(combine_parts(re.with_precision(p), im.with_precision(p)))
164    }
165
166    /// Inverse tangent `atan z = (i/2)·(log(1-iz) - log(1+iz))` (context layer).
167    pub fn atan<const B: Word>(
168        &self,
169        z: &CBig<R, B>,
170        mut cache: Option<&mut ConstCache>,
171    ) -> CfpResult<R, B> {
172        if z.is_infinite() {
173            // atan(±∞) = ±π/2; defer the exact constant to the formula via the limit, but the
174            // 1±iz terms become infinite and the log diverges — report Indeterminate for now.
175            return Err(FpError::Indeterminate);
176        }
177        let gctx = Context::new(self.precision() + ITRIG_GUARD);
178        let p = self.precision();
179        let one = CBig::ONE;
180        let iz = z.mul_i(false);
181        let a = gctx.sub(&one, &iz)?.value(); // 1 - iz
182        let b = gctx.add(&one, &iz)?.value(); // 1 + iz
183        let log_a = gctx.log(&a, reborrow_cache(&mut cache))?.value();
184        let log_b = gctx.log(&b, reborrow_cache(&mut cache))?.value();
185        let diff = gctx.sub(&log_a, &log_b)?.value();
186        let i_half_diff = diff.mul_i(false); // i·diff, then /2 below
187        let two: CBig<R, B> = IBig::from(2).into();
188        let atan_z = gctx.div(&i_half_diff, &two)?.value();
189        let (re, im) = atan_z.into_parts();
190        Ok(combine_parts(re.with_precision(p), im.with_precision(p)))
191    }
192}
193
194/// Guard digits (base-B) for the inverse trig (squares, a sqrt, logs, and a divide).
195const ITRIG_GUARD: usize = 18;
196
197impl<R: Round, const B: Word> CBig<R, B> {
198    /// Complex sine (convenience layer). Panics on an indeterminate special value.
199    #[inline]
200    pub fn sin(&self) -> Self {
201        self.context().unwrap_cfp(self.context().sin(self, None))
202    }
203
204    /// Complex cosine (convenience layer). Panics on an indeterminate special value.
205    #[inline]
206    pub fn cos(&self) -> Self {
207        self.context().unwrap_cfp(self.context().cos(self, None))
208    }
209
210    /// Simultaneously compute `(sin z, cos z)` (convenience layer).
211    #[inline]
212    pub fn sin_cos(&self) -> (Self, Self) {
213        let (s, c) = self.context().sin_cos(self, None);
214        (self.context().unwrap_cfp(s), self.context().unwrap_cfp(c))
215    }
216
217    /// Complex tangent (convenience layer).
218    #[inline]
219    pub fn tan(&self) -> Self {
220        self.context().unwrap_cfp(self.context().tan(self, None))
221    }
222
223    /// Inverse sine (convenience layer).
224    #[inline]
225    pub fn asin(&self) -> Self {
226        self.context().unwrap_cfp(self.context().asin(self, None))
227    }
228
229    /// Inverse cosine (convenience layer).
230    #[inline]
231    pub fn acos(&self) -> Self {
232        self.context().unwrap_cfp(self.context().acos(self, None))
233    }
234
235    /// Inverse tangent (convenience layer).
236    #[inline]
237    pub fn atan(&self) -> Self {
238        self.context().unwrap_cfp(self.context().atan(self, None))
239    }
240}
241
242#[cfg(test)]
243mod tests {
244    use super::*;
245    use dashu_float::round::mode;
246
247    type C = CBig<mode::HalfAway, 10>;
248    type F = FBig<mode::HalfAway, 10>;
249
250    fn c(re: i32, im: i32) -> C {
251        let mk = |v: i32| -> F { F::from(v).with_precision(53).value() };
252        CBig::from_parts(mk(re), mk(im))
253    }
254
255    #[test]
256    fn sin_zero_is_zero() {
257        assert!(C::ZERO.sin() == C::ZERO);
258    }
259
260    #[test]
261    fn cos_zero_is_one() {
262        assert!(C::ZERO.cos() == C::ONE);
263    }
264
265    #[test]
266    fn pythagorean_identity() {
267        // sin²z + cos²z = 1
268        let z = c(1, 1);
269        let s = z.sin();
270        let co = z.cos();
271        let sum = &s.sqr() + &co.sqr();
272        // purely real ≈ 1, imaginary ≈ 0
273        let (re, im) = sum.into_parts();
274        use dashu_base::{Abs, AbsOrd};
275        assert!((re.clone() - F::ONE)
276            .abs()
277            .abs_cmp(&F::from_parts(1.into(), -12))
278            .is_le());
279        assert!(im.abs_cmp(&F::from_parts(1.into(), -12)).is_le());
280    }
281
282    #[test]
283    fn sin_i_is_i_sinh_one() {
284        // sin(i) = i·sinh(1) = i·1.1752… ; purely imaginary
285        let s = C::I.sin();
286        assert!(s.re().significand().is_zero());
287        assert!(!s.im().significand().is_zero());
288    }
289
290    #[test]
291    fn asin_zero_is_zero() {
292        assert!(C::ZERO.asin() == C::ZERO);
293    }
294
295    #[test]
296    fn asin_one_is_half_pi() {
297        use dashu_base::{Abs, AbsOrd};
298        // asin(1) = π/2
299        let (re, im) = C::ONE.asin().into_parts();
300        let half_pi = F::from_parts(15707963267948966i64.into(), -16)
301            .with_precision(60)
302            .value();
303        assert!((re.clone() - half_pi)
304            .abs()
305            .abs_cmp(&F::from_parts(1.into(), -12))
306            .is_le());
307        assert!(im.abs_cmp(&F::from_parts(1.into(), -12)).is_le());
308    }
309
310    #[test]
311    fn acos_zero_is_half_pi() {
312        use dashu_base::{Abs, AbsOrd};
313        let (re, _im) = C::ZERO.acos().into_parts();
314        let half_pi = F::from_parts(15707963267948966i64.into(), -16)
315            .with_precision(60)
316            .value();
317        assert!((re - half_pi)
318            .abs()
319            .abs_cmp(&F::from_parts(1.into(), -12))
320            .is_le());
321    }
322
323    #[test]
324    fn atan_one_is_quarter_pi() {
325        use dashu_base::{Abs, AbsOrd};
326        // atan(1) = π/4
327        let (re, _im) = C::ONE.atan().into_parts();
328        let quarter_pi = F::from_parts(7853981633974483i64.into(), -16)
329            .with_precision(60)
330            .value();
331        assert!((re - quarter_pi)
332            .abs()
333            .abs_cmp(&F::from_parts(1.into(), -12))
334            .is_le());
335    }
336
337    #[test]
338    fn sin_asin_roundtrip() {
339        // asin(sin z) ≈ z for a small z (within the principal range)
340        let z = c(1, 1);
341        let r = z.sin().asin();
342        assert!(r == z);
343    }
344}
dashu_cmplx/math/trig.rs

dashu_cmplx/math/
trig.rs
