dashu_float/math/
trig.rs

1use crate::{
2    error::assert_limited_precision,
3    fbig::FBig,
4    math::FpResult,
5    repr::{Context, Repr, Word},
6    round::Round,
7};
8use core::cmp::Ordering;
9use core::convert::TryFrom;
10use dashu_base::{AbsOrd, RemEuclid, Sign};
11use dashu_int::IBig;
12
13#[derive(Debug, Clone, Copy, PartialEq, Eq)]
14enum Quadrant {
15    First,
16    Second,
17    Third,
18    Fourth,
19}
20
21impl<R: Round> Context<R> {
22    /// Calculate the internal work context for trigonometric functions based on input magnitude.
23    ///
24    /// This ensures we have enough guard digits to prevent catastrophic cancellation
25    /// during range reduction for large inputs.
26    fn compute_work_context_trig<const B: Word>(self, x: &Repr<B>) -> Self {
27        // x_mag estimates m = floor(log_BASE(|x|))
28        let x_mag = (x.exponent.saturating_add(x.digits_ub() as isize)).max(0) as usize;
29
30        // We need precision + log10(x) digits to maintain 'precision' digits after reduction.
31        // We add a base of 50 guard digits, plus 10% of x_mag for very large arguments
32        // to account for cumulative errors in division and multiplication during reduction.
33        let extra_guards = 50 + x_mag / 10;
34        let work_precision = self
35            .precision
36            .saturating_add(x_mag)
37            .saturating_add(extra_guards);
38        Self::new(work_precision)
39    }
40
41    /// Reduces the argument to the first quadrant for trigonometric evaluation.
42    /// Returns the internal work context, the reduced argument `r`, and the quadrant `k % 4`.
43    fn reduce_to_quadrant<const B: Word>(self, x: &Repr<B>) -> (Self, FBig<R, B>, Quadrant) {
44        let work_context = self.compute_work_context_trig(x);
45        let x_f = FBig::<R, B>::new(work_context.repr_round(x.clone()).value(), work_context);
46
47        let pi = work_context.pi::<B>().value();
48        let half_pi = &pi / 2;
49        let x_scaled: FBig<R, B> = &x_f / &half_pi;
50        let k_f = x_scaled.round();
51        let r = x_f - &k_f * half_pi;
52        let Ok(k) = IBig::try_from(k_f) else {
53            unreachable!(
54                "round() always returns an integer and trig functions ensure input is finite"
55            );
56        };
57
58        let k_mod_4_big = k.rem_euclid(IBig::from(4));
59        let Ok(k_mod_4_int) = i8::try_from(k_mod_4_big) else {
60            unreachable!("k % 4 is always in [0, 3]");
61        };
62        let quadrant = match k_mod_4_int {
63            0 => Quadrant::First,
64            1 => Quadrant::Second,
65            2 => Quadrant::Third,
66            3 => Quadrant::Fourth,
67            _ => unreachable!(),
68        };
69
70        (work_context, r, quadrant)
71    }
72
73    /// Calculate the sine of the floating point representation.
74    #[must_use]
75    pub fn sin<const B: Word>(&self, x: &Repr<B>) -> FpResult<B> {
76        if x.is_infinite() {
77            return FpResult::NaN;
78        }
79        assert_limited_precision(self.precision);
80
81        if x.is_zero() {
82            let res = FBig::<R, B>::ZERO.with_precision(self.precision);
83            return FpResult::Normal(res.map(|v| v.repr));
84        }
85
86        let (work_context, r, quadrant) = self.reduce_to_quadrant(x);
87
88        // 3. Evaluate the reduced series based on the quadrant
89        let res = match quadrant {
90            Quadrant::First => work_context.sin_internal(&r),
91            Quadrant::Second => work_context.cos_internal(&r),
92            Quadrant::Third => -work_context.sin_internal(&r),
93            Quadrant::Fourth => -work_context.cos_internal(&r),
94        };
95        FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr))
96    }
97
98    /// Internal Taylor series for sine: S(x) = x - x^3/3! + x^5/5! - ...
99    fn sin_internal<const B: Word>(self, x: &FBig<R, B>) -> FBig<R, B> {
100        if x.repr.significand.is_zero() {
101            return FBig::ZERO;
102        }
103        let x2 = x.sqr();
104        let mut sum = x.clone();
105        let mut term = x.clone();
106        let mut k = 1usize;
107        let threshold = sum.sub_ulp();
108        loop {
109            term *= &x2;
110            term /= (2 * k) * (2 * k + 1);
111            if term.abs_cmp(&threshold).is_le() {
112                break;
113            }
114            if k % 2 == 1 {
115                sum -= &term;
116            } else {
117                sum += &term;
118            }
119            k += 1;
120        }
121        sum
122    }
123
124    /// Calculate the cosine of the floating point representation.
125    #[must_use]
126    pub fn cos<const B: Word>(&self, x: &Repr<B>) -> FpResult<B> {
127        if x.is_infinite() {
128            return FpResult::NaN;
129        }
130        assert_limited_precision(self.precision);
131
132        if x.is_zero() {
133            let res = FBig::<R, B>::ONE.with_precision(self.precision);
134            return FpResult::Normal(res.map(|v| v.repr));
135        }
136
137        let (work_context, r, quadrant) = self.reduce_to_quadrant(x);
138
139        // 3. Evaluate the reduced series based on the quadrant
140        let res = match quadrant {
141            Quadrant::First => work_context.cos_internal(&r),
142            Quadrant::Second => -work_context.sin_internal(&r),
143            Quadrant::Third => -work_context.cos_internal(&r),
144            Quadrant::Fourth => work_context.sin_internal(&r),
145        };
146        FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr))
147    }
148
149    /// Internal Taylor series for cosine: C(x) = 1 - x^2/2! + x^4/4! - ...
150    fn cos_internal<const B: Word>(self, x: &FBig<R, B>) -> FBig<R, B> {
151        if x.repr.significand.is_zero() {
152            return FBig::ONE.with_precision(self.precision).value();
153        }
154        let x2 = x.sqr();
155        let mut sum = FBig::<R, B>::ONE.with_precision(self.precision).value();
156        let mut term = sum.clone();
157        let mut k = 1usize;
158        let threshold = sum.sub_ulp();
159        loop {
160            term *= &x2;
161            term /= (2 * k) * (2 * k - 1);
162            if term.abs_cmp(&threshold).is_le() {
163                break;
164            }
165            if k % 2 == 1 {
166                sum -= &term;
167            } else {
168                sum += &term;
169            }
170            k += 1;
171        }
172        sum
173    }
174
175    /// Calculate both the sine and cosine of the floating point representation.
176    ///
177    /// This is more efficient than calling `sin` and `cos` separately.
178    #[must_use]
179    pub fn sin_cos<const B: Word>(&self, x: &Repr<B>) -> (FpResult<B>, FpResult<B>) {
180        if x.is_infinite() {
181            return (FpResult::NaN, FpResult::NaN);
182        }
183        assert_limited_precision(self.precision);
184
185        if x.is_zero() {
186            let s = FBig::<R, B>::ZERO.with_precision(self.precision);
187            let c = FBig::<R, B>::ONE.with_precision(self.precision);
188            return (FpResult::Normal(s.map(|v| v.repr)), FpResult::Normal(c.map(|v| v.repr)));
189        }
190
191        let (work_context, r, quadrant) = self.reduce_to_quadrant(x);
192
193        let (sin_r, cos_r) = work_context.sin_cos_internal(&r);
194
195        let (s, c) = match quadrant {
196            Quadrant::First => (sin_r, cos_r),
197            Quadrant::Second => (cos_r, -sin_r),
198            Quadrant::Third => (-sin_r, -cos_r),
199            Quadrant::Fourth => (-cos_r, sin_r),
200        };
201
202        (
203            FpResult::Normal(s.with_precision(self.precision).map(|v| v.repr)),
204            FpResult::Normal(c.with_precision(self.precision).map(|v| v.repr)),
205        )
206    }
207
208    /// Simultaneously evaluate Taylor series for sine and cosine.
209    pub(crate) fn sin_cos_internal<const B: Word>(
210        self,
211        x: &FBig<R, B>,
212    ) -> (FBig<R, B>, FBig<R, B>) {
213        if x.repr.significand.is_zero() {
214            return (FBig::ZERO, FBig::ONE.with_precision(self.precision).value());
215        }
216        let x2 = x.sqr();
217        let mut sin_sum = x.clone();
218        let mut cos_sum = FBig::<R, B>::ONE.with_precision(self.precision).value();
219        let mut sin_term = x.clone();
220        let mut cos_term = cos_sum.clone();
221        let mut k = 1usize;
222        let sin_threshold = sin_sum.sub_ulp();
223        let cos_threshold = cos_sum.sub_ulp();
224        loop {
225            cos_term *= &x2;
226            cos_term /= (2 * k) * (2 * k - 1);
227            sin_term *= &x2;
228            sin_term /= (2 * k) * (2 * k + 1);
229
230            if sin_term.abs_cmp(&sin_threshold).is_le() && cos_term.abs_cmp(&cos_threshold).is_le()
231            {
232                break;
233            }
234
235            if k % 2 == 1 {
236                cos_sum -= &cos_term;
237                sin_sum -= &sin_term;
238            } else {
239                cos_sum += &cos_term;
240                sin_sum += &sin_term;
241            }
242            k += 1;
243        }
244        (sin_sum, cos_sum)
245    }
246
247    /// Calculate the tangent of the floating point representation.
248    ///
249    /// # Note
250    /// Near odd multiples of π/2, the result returns `Infinite`.
251    #[must_use]
252    pub fn tan<const B: Word>(&self, x: &Repr<B>) -> FpResult<B> {
253        if x.is_infinite() {
254            return FpResult::NaN;
255        }
256        assert_limited_precision(self.precision);
257
258        if x.is_zero() {
259            let res = FBig::<R, B>::ZERO.with_precision(self.precision);
260            return FpResult::Normal(res.map(|v| v.repr));
261        }
262
263        let (work_context, r, quadrant) = self.reduce_to_quadrant(x);
264        let (sin_r, cos_r) = work_context.sin_cos_internal(&r);
265
266        let (s_f, c_f) = match quadrant {
267            Quadrant::First => (sin_r, cos_r),
268            Quadrant::Second => (cos_r, -sin_r),
269            Quadrant::Third => (-sin_r, -cos_r),
270            Quadrant::Fourth => (-cos_r, sin_r),
271        };
272
273        if c_f.repr.is_zero() {
274            return FpResult::Infinite;
275        }
276        FpResult::Normal(self.div(&s_f.repr, &c_f.repr).map(|v| v.repr))
277    }
278
279    /// Calculate the arcsine of the floating point representation.
280    ///
281    /// # Methodology
282    /// Uses the identity: `asin(x) = atan(x / sqrt(1 - x^2))`
283    /// Returns `NaN` if `|x| > 1`.
284    #[must_use]
285    pub fn asin<const B: Word>(&self, x: &Repr<B>) -> FpResult<B> {
286        if x.is_infinite() {
287            return FpResult::NaN;
288        }
289        assert_limited_precision(self.precision);
290
291        let x_orig = FBig::<R, B>::new(x.clone(), *self);
292        // Domain check: |x| must be <= 1
293        if x_orig.abs_cmp(&FBig::ONE).is_gt() {
294            return FpResult::NaN;
295        }
296
297        let guard_digits = 50;
298        let work_precision = self.precision + guard_digits;
299        let work_context = Self::new(work_precision);
300
301        let x_f = FBig::<R, B>::new(work_context.repr_round(x.clone()).value(), work_context);
302
303        let res = work_context.asin_internal(&x_f);
304        FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr))
305    }
306
307    fn asin_internal<const B: Word>(self, x_f: &FBig<R, B>) -> FBig<R, B> {
308        let one = FBig::<R, B>::ONE.with_precision(self.precision).value();
309        let x2 = x_f.sqr();
310        let d = self.sqrt(&(one - x2).repr).value();
311
312        if d.repr.is_zero() {
313            let pi = self.pi::<B>().value();
314            let half_pi: FBig<R, B> = pi / 2;
315            if x_f.sign() == Sign::Positive {
316                return half_pi;
317            }
318            return -half_pi;
319        }
320
321        self.atan_with_reduction(&(x_f / d))
322    }
323
324    /// Calculate the arccosine of the floating point representation.
325    ///
326    /// # Methodology
327    /// Uses the identity: `acos(x) = pi/2 - asin(x)`.
328    /// Higher precision is used internally to avoid catastrophic cancellation near x ≈ 1.
329    #[must_use]
330    pub fn acos<const B: Word>(&self, x: &Repr<B>) -> FpResult<B> {
331        if x.is_infinite() {
332            return FpResult::NaN;
333        }
334        assert_limited_precision(self.precision);
335
336        let x_orig = FBig::<R, B>::new(x.clone(), *self);
337        // Domain check: |x| must be <= 1
338        if x_orig.abs_cmp(&FBig::ONE).is_gt() {
339            return FpResult::NaN;
340        }
341
342        let guard_digits = 50;
343        let work_precision = self.precision + guard_digits;
344        let work_context = Self::new(work_precision);
345
346        let x_f = FBig::<R, B>::new(work_context.repr_round(x.clone()).value(), work_context);
347
348        let asin_x = work_context.asin_internal(&x_f);
349        let pi = work_context.pi::<B>().value();
350        let half_pi: FBig<R, B> = pi / 2;
351        let res: FBig<R, B> = half_pi - asin_x;
352        FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr))
353    }
354
355    /// Calculate the arctangent of the floating point representation.
356    #[must_use]
357    pub fn atan<const B: Word>(&self, x: &Repr<B>) -> FpResult<B> {
358        if x.is_infinite() {
359            let pi = self.pi::<B>().value();
360            let half_pi: FBig<R, B> = pi / 2;
361            let res: FBig<R, B> = if x.sign() == Sign::Positive {
362                half_pi
363            } else {
364                -half_pi
365            };
366            return FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr));
367        }
368
369        assert_limited_precision(self.precision);
370
371        if x.is_zero() {
372            let res = FBig::<R, B>::ZERO.with_precision(self.precision);
373            return FpResult::Normal(res.map(|v| v.repr));
374        }
375
376        let guard_digits = 50;
377        let work_precision = self.precision + guard_digits;
378        let work_context = Self::new(work_precision);
379
380        let x_f = FBig::<R, B>::new(work_context.repr_round(x.clone()).value(), work_context);
381        let res = work_context.atan_with_reduction(&x_f);
382        FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr))
383    }
384
385    /// Internal arctangent that includes range reduction but no guard digit allocation.
386    fn atan_with_reduction<const B: Word>(self, x_f: &FBig<R, B>) -> FBig<R, B> {
387        let sign = x_f.sign();
388        let mut x_abs = x_f.clone();
389        if sign == Sign::Negative {
390            x_abs = -x_abs;
391        }
392        let mut res = if x_abs >= FBig::<R, B>::ONE.with_precision(self.precision).value() {
393            let pi = self.pi::<B>().value();
394            let inv_x = FBig::<R, B>::ONE.with_precision(self.precision).value() / x_abs;
395            (pi / 2) - self.atan_internal(&inv_x)
396        } else {
397            self.atan_internal(&x_abs)
398        };
399        if sign == Sign::Negative {
400            res = -res;
401        }
402        res
403    }
404
405    /// Internal series for arctangent.
406    /// Evaluates the Euler series for arctangent.
407    fn atan_internal<const B: Word>(self, x: &FBig<R, B>) -> FBig<R, B> {
408        // Euler's series for atan(x)
409        let x2 = x.sqr();
410        let one_plus_x2 = FBig::ONE + &x2;
411        let mut term = x / &one_plus_x2;
412        let mut sum = term.clone();
413        let factor = (2 * &x2) / one_plus_x2;
414        let mut n = 1usize;
415        let threshold = sum.sub_ulp();
416        loop {
417            term *= &factor;
418            term *= n;
419            term /= 2 * n + 1;
420            if term.abs_cmp(&threshold).is_le() {
421                break;
422            }
423            sum += &term;
424            n += 1;
425        }
426        sum
427    }
428
429    /// Calculate the arctangent of y / x.
430    ///
431    /// Handles signed infinities according to IEEE 754 standards.
432    /// Returns `NaN` if both arguments are zero.
433    #[must_use]
434    pub fn atan2<const B: Word>(&self, y: &Repr<B>, x: &Repr<B>) -> FpResult<B> {
435        if y.is_zero() && x.is_zero() {
436            return FpResult::NaN;
437        }
438
439        assert_limited_precision(self.precision);
440
441        let guard_digits = 50;
442        let work_precision = self.precision + guard_digits;
443        let work_context = Self::new(work_precision);
444
445        // Handle Infinities according to IEEE 754
446        if y.is_infinite() || x.is_infinite() {
447            let (sy, sx) = (y.sign() == Sign::Positive, x.sign() == Sign::Positive);
448            let res: FBig<R, B> = match (y.is_infinite(), x.is_infinite(), sy, sx) {
449                (true, true, true, true) => work_context.pi::<B>().value() / 4,
450                (true, true, true, false) => work_context.pi::<B>().value() * 3 / 4,
451                (true, true, false, true) => {
452                    let pi4: FBig<R, B> = work_context.pi::<B>().value() / 4;
453                    -pi4
454                }
455                (true, true, false, false) => {
456                    let pi34: FBig<R, B> = work_context.pi::<B>().value() * 3 / 4;
457                    -pi34
458                }
459                (true, false, true, _) => work_context.pi::<B>().value() / 2,
460                (true, false, false, _) => {
461                    let half_pi: FBig<R, B> = work_context.pi::<B>().value() / 2;
462                    -half_pi
463                }
464                (false, true, _, true) => FBig::<R, B>::ZERO.with_precision(work_precision).value(),
465                (false, true, true, false) => work_context.pi::<B>().value(),
466                (false, true, false, false) => -work_context.pi::<B>().value(),
467                _ => unreachable!(),
468            };
469            // Note: atan2(finite, +inf) returns unsigned ZERO. IEEE 754 requires signed zero,
470            // but `Repr` does not distinguish signed zero.
471            return FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr));
472        }
473
474        let y_f = FBig::<R, B>::new(work_context.repr_round(y.clone()).value(), work_context);
475        let x_f = FBig::<R, B>::new(work_context.repr_round(x.clone()).value(), work_context);
476
477        match x_f.cmp(&FBig::<R, B>::ZERO) {
478            Ordering::Greater => {
479                let res = work_context.atan_with_reduction(&(y_f / x_f));
480                FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr))
481            }
482            Ordering::Less => {
483                let pi = work_context.pi::<B>().value();
484                let y_sign = y_f.sign();
485                let atan_yx = work_context.atan_with_reduction(&(y_f / x_f));
486                let res = if y_sign == Sign::Positive {
487                    atan_yx + pi
488                } else {
489                    atan_yx - pi
490                };
491                FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr))
492            }
493            Ordering::Equal => {
494                // x == 0 case
495                let pi = work_context.pi::<B>().value();
496                let half_pi: FBig<R, B> = pi / 2;
497                if y_f > FBig::<R, B>::ZERO {
498                    FpResult::Normal(half_pi.with_precision(self.precision).map(|v| v.repr))
499                } else {
500                    let res = -half_pi;
501                    FpResult::Normal(res.with_precision(self.precision).map(|v| v.repr))
502                }
503            }
504        }
505    }
506}
507
508impl<R: Round, const B: Word> FBig<R, B> {
509    /// Calculate the sine of the floating point number.
510    ///
511    /// # Panics
512    /// Panics if the input is infinite or the result is not representable as a normal value.
513    #[inline]
514    #[must_use]
515    pub fn sin(&self) -> Self {
516        self.context.sin(&self.repr).value(&self.context)
517    }
518
519    /// Calculate the cosine of the floating point number.
520    ///
521    /// # Panics
522    /// Panics if the input is infinite or the result is not representable as a normal value.
523    #[inline]
524    #[must_use]
525    pub fn cos(&self) -> Self {
526        self.context.cos(&self.repr).value(&self.context)
527    }
528
529    /// Calculate both the sine and cosine of the floating point number.
530    ///
531    /// This is more efficient than calling `sin` and `cos` separately.
532    ///
533    /// # Panics
534    /// Panics if the input is infinite or the results are not representable as normal values.
535    #[inline]
536    #[must_use]
537    pub fn sin_cos(&self) -> (Self, Self) {
538        let (s, c) = self.context.sin_cos(&self.repr);
539        (s.value(&self.context), c.value(&self.context))
540    }
541
542    /// Calculate the tangent of the floating point number.
543    ///
544    /// Returns `FpResult` to safely handle non-finite results (e.g., at singularities).
545    #[inline]
546    #[must_use]
547    pub fn tan(&self) -> FpResult<B> {
548        self.context.tan(&self.repr)
549    }
550
551    /// Calculate the arcsine of the floating point number.
552    ///
553    /// Returns `FpResult` to safely handle domain errors (e.g., |x| > 1).
554    #[inline]
555    #[must_use]
556    pub fn asin(&self) -> FpResult<B> {
557        self.context.asin(&self.repr)
558    }
559
560    /// Calculate the arccosine of the floating point number.
561    ///
562    /// Returns `FpResult` to safely handle domain errors (e.g., |x| > 1).
563    #[inline]
564    #[must_use]
565    pub fn acos(&self) -> FpResult<B> {
566        self.context.acos(&self.repr)
567    }
568
569    /// Calculate the arctangent of the floating point number.
570    ///
571    /// # Panics
572    /// Panics if the result is not representable as a normal value.
573    #[inline]
574    #[must_use]
575    pub fn atan(&self) -> Self {
576        self.context.atan(&self.repr).value(&self.context)
577    }
578
579    /// Calculate the arctangent of y / x.
580    ///
581    /// Returns `FpResult` to safely handle special cases like (0,0) or infinities.
582    #[inline]
583    #[must_use]
584    pub fn atan2(&self, x: &Self) -> FpResult<B> {
585        self.context.atan2(&self.repr, &x.repr)
586    }
587}
dashu_float/math/trig.rs

dashu_float/math/
trig.rs
