zenith_float_num/
ball.rs

1//! Interval enclosures (`Ball`) and Ziv correct-rounding retries.
2
3use crate::common::util::bump_prec_retry;
4use crate::common::util::round_p;
5use crate::defs::WORD_BIT_SIZE;
6use crate::Consts;
7use crate::Error;
8use crate::ExactComplex;
9use crate::ExactNum;
10use crate::RoundingMode;
11
12/// Extra rounding-ulp multiples added to a transcendental Lipschitz radius.
13const BALL_TRANSCENDENTAL_ERROR_TERMS: u32 = 8;
14
15/// Enclosure `mid ± rad` used to certify that a rounded midpoint is unique.
16#[derive(Clone, Debug)]
17pub struct Ball {
18    mid: ExactNum,
19    rad: ExactNum,
20}
21
22impl Ball {
23    /// `mid ± rad`. The radius is taken in absolute value.
24    pub fn new(mid: ExactNum, rad: ExactNum) -> Self {
25        Ball {
26            mid,
27            rad: rad.abs(),
28        }
29    }
30
31    /// Midpoint of the enclosure.
32    pub fn mid(&self) -> &ExactNum {
33        &self.mid
34    }
35
36    /// Non-negative radius.
37    pub fn rad(&self) -> &ExactNum {
38        &self.rad
39    }
40
41    fn rounding_ulp(x: &ExactNum, p: usize) -> ExactNum {
42        if x.is_nan() || x.is_inf() || x.is_zero() {
43            return ExactNum::new(p);
44        }
45        let e = x.exponent().unwrap_or(0);
46        let bits = x.mantissa_max_bit_len().unwrap_or(p) as i32;
47        let mut u = ExactNum::from_word(1, p);
48        u.set_exponent(e.saturating_sub(bits.saturating_sub(2)));
49        u
50    }
51
52    /// Sum of two balls: midpoint add at precision `p`, radius `r1+r2` plus a rounding ulp.
53    pub fn add(&self, other: &Self, p: usize, rm: RoundingMode) -> Self {
54        let mid = self.mid.add(&other.mid, p, rm);
55        let rad = self.rad.add(&other.rad, p, RoundingMode::Up).add(
56            &Self::rounding_ulp(&mid, p),
57            p,
58            RoundingMode::Up,
59        );
60        Ball { mid, rad }
61    }
62
63    /// Product of two balls with a first-order radius bound plus a rounding ulp.
64    pub fn mul(&self, other: &Self, p: usize, rm: RoundingMode) -> Self {
65        let mid = self.mid.mul(&other.mid, p, rm);
66        let a = self.mid.abs().mul(&other.rad, p, RoundingMode::Up);
67        let b = other.mid.abs().mul(&self.rad, p, RoundingMode::Up);
68        let c = self.rad.mul(&other.rad, p, RoundingMode::Up);
69        let rad = a
70            .add(&b, p, RoundingMode::Up)
71            .add(&c, p, RoundingMode::Up)
72            .add(&Self::rounding_ulp(&mid, p), p, RoundingMode::Up);
73        Ball { mid, rad }
74    }
75
76    /// Exponential of a ball. `exp` is increasing; the radius uses
77    /// \(\lvert\exp(m)\rvert(e^{r}-1)\) plus a rounding ulp (same `Up` convention as `add`/`mul`).
78    pub fn exp(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
79        let mid = self.mid.exp(p, rm, cc);
80        let em1 = self.rad.expm1(p, RoundingMode::Up, cc);
81        let rad = mid.abs().mul(&em1, p, RoundingMode::Up).add(
82            &Self::rounding_ulp(&mid, p),
83            p,
84            RoundingMode::Up,
85        );
86        Ball { mid, rad }
87    }
88
89    /// Sine of a ball. \(\lvert\sin'\rvert\le 1\), so the image radius is at most `rad`
90    /// plus a rounding ulp.
91    pub fn sin(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
92        let mid = self.mid.sin(p, rm, cc);
93        let rad = self
94            .rad
95            .add(&Self::rounding_ulp(&mid, p), p, RoundingMode::Up);
96        Ball { mid, rad }
97    }
98
99    /// Cosine of a ball. \(\lvert\cos'\rvert\le 1\).
100    pub fn cos(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
101        let mid = self.mid.cos(p, rm, cc);
102        let rad = self
103            .rad
104            .add(&Self::rounding_ulp(&mid, p), p, RoundingMode::Up);
105        Ball { mid, rad }
106    }
107
108    /// Natural log of a ball. Domain: the ball must lie in \((0,+\infty)\).
109    /// Lipschitz \(\lvert\ln'\rvert=1/x\le 1/(m-r)\).
110    pub fn ln(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
111        if !self.strictly_positive(p) {
112            return Self::nan_ball(p);
113        }
114        let mid = self.mid.ln(p, rm, cc);
115        let den = self.mid.sub(&self.rad, p, RoundingMode::Down);
116        let lip = ExactNum::from_u8(1, p).div(&den, p, RoundingMode::Up);
117        let rad = lip.mul(&self.rad, p, RoundingMode::Up).add(
118            &Self::transcendental_slack(&mid, p),
119            p,
120            RoundingMode::Up,
121        );
122        Ball { mid, rad }
123    }
124
125    /// Square root of a ball. Domain: the ball must lie in \((0,+\infty)\).
126    /// Lipschitz \(1/(2\sqrt{x})\le 1/(2\sqrt{m-r})\).
127    pub fn sqrt(&self, p: usize, rm: RoundingMode) -> Self {
128        if !self.strictly_positive(p) {
129            return Self::nan_ball(p);
130        }
131        let mid = self.mid.sqrt(p, rm);
132        let lo = self
133            .mid
134            .sub(&self.rad, p, RoundingMode::Down)
135            .sqrt(p, RoundingMode::Down);
136        let two = ExactNum::from_u8(2, p);
137        let den = two.mul(&lo, p, RoundingMode::Down);
138        let lip = ExactNum::from_u8(1, p).div(&den, p, RoundingMode::Up);
139        let rad = lip.mul(&self.rad, p, RoundingMode::Up).add(
140            &Self::transcendental_slack(&mid, p),
141            p,
142            RoundingMode::Up,
143        );
144        Ball { mid, rad }
145    }
146
147    /// Error function of a ball. \(\lvert\mathrm{erf}'\rvert\le 2/\sqrt{\pi}\).
148    pub fn erf(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
149        let mid = self.mid.erf(p, rm, cc);
150        let two = ExactNum::from_u8(2, p);
151        let s = cc.pi(p, RoundingMode::Down).sqrt(p, RoundingMode::Down);
152        let lip = two.div(&s, p, RoundingMode::Up);
153        let rad = lip.mul(&self.rad, p, RoundingMode::Up).add(
154            &Self::transcendental_slack(&mid, p),
155            p,
156            RoundingMode::Up,
157        );
158        Ball { mid, rad }
159    }
160
161    /// \(J_0\) of a ball. \(\lvert J_0'\rvert=\lvert J_1\rvert\le 1\).
162    pub fn bessel_j0(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
163        let mid = self.mid.bessel_j(0, p, rm, cc);
164        let rad = self
165            .rad
166            .add(&Self::transcendental_slack(&mid, p), p, RoundingMode::Up);
167        Ball { mid, rad }
168    }
169
170    /// \(J_1\) of a ball. \(\lvert J_1'\rvert=\lvert J_0-J_1/x\rvert\le 1+1/(\lvert m\rvert-r)\)
171    /// when the ball excludes \(0\).
172    pub fn bessel_j1(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
173        if !self.excludes_zero(p) {
174            return Self::nan_ball(p);
175        }
176        let mid = self.mid.bessel_j(1, p, rm, cc);
177        let den = self.mid.abs().sub(&self.rad, p, RoundingMode::Down);
178        let extra = ExactNum::from_u8(1, p).div(&den, p, RoundingMode::Up);
179        let lip = ExactNum::from_u8(1, p).add(&extra, p, RoundingMode::Up);
180        let rad = lip.mul(&self.rad, p, RoundingMode::Up).add(
181            &Self::transcendental_slack(&mid, p),
182            p,
183            RoundingMode::Up,
184        );
185        Ball { mid, rad }
186    }
187
188    fn transcendental_slack(mid: &ExactNum, p: usize) -> ExactNum {
189        let u = Self::rounding_ulp(mid, p);
190        u.mul(
191            &ExactNum::from_u8(BALL_TRANSCENDENTAL_ERROR_TERMS as u8, p),
192            p,
193            RoundingMode::Up,
194        )
195    }
196
197    fn nan_ball(p: usize) -> Self {
198        let n = ExactNum::nan(Some(Error::InvalidArgument));
199        let _ = p;
200        Ball {
201            mid: n.clone(),
202            rad: n,
203        }
204    }
205
206    fn strictly_positive(&self, p: usize) -> bool {
207        if self.mid.is_nan() || self.rad.is_nan() || !self.mid.is_positive() {
208            return false;
209        }
210        matches!(self.mid.cmp(&self.rad), Some(c) if c > 0)
211            && !self.mid.sub(&self.rad, p, RoundingMode::Down).is_negative()
212            && !self.mid.sub(&self.rad, p, RoundingMode::Down).is_zero()
213    }
214
215    fn excludes_zero(&self, p: usize) -> bool {
216        if self.mid.is_nan() || self.rad.is_nan() || self.mid.is_zero() {
217            return false;
218        }
219        matches!(self.mid.abs().cmp(&self.rad), Some(c) if c > 0)
220            && !self
221                .mid
222                .abs()
223                .sub(&self.rad, p, RoundingMode::Down)
224                .is_zero()
225    }
226
227    /// True when `x` lies in `[mid − rad, mid + rad]` (NaN / Inf never contained).
228    pub fn contains(&self, x: &ExactNum, p: usize) -> bool {
229        if x.is_nan() || self.mid.is_nan() || self.rad.is_nan() {
230            return false;
231        }
232        let d = self.mid.sub(x, p, RoundingMode::None).abs();
233        matches!(d.cmp(&self.rad), Some(c) if c <= 0)
234    }
235}
236
237/// Disk enclosure in \(\mathbb{C}\): center `mid`, radius `rad`.
238#[derive(Clone, Debug)]
239pub struct ComplexBall {
240    mid: ExactComplex,
241    rad: ExactNum,
242}
243
244impl ComplexBall {
245    /// Disk `mid` with radius `rad` (taken in absolute value).
246    pub fn new(mid: ExactComplex, rad: ExactNum) -> Self {
247        ComplexBall {
248            mid,
249            rad: rad.abs(),
250        }
251    }
252
253    /// Center.
254    pub fn mid(&self) -> &ExactComplex {
255        &self.mid
256    }
257
258    /// Non-negative radius.
259    pub fn rad(&self) -> &ExactNum {
260        &self.rad
261    }
262
263    fn slack(mid: &ExactComplex, p: usize) -> ExactNum {
264        let u_re = Ball::rounding_ulp(mid.re(), p);
265        let u_im = Ball::rounding_ulp(mid.im(), p);
266        let u = if matches!(u_re.cmp(&u_im), Some(c) if c >= 0) { u_re } else { u_im };
267        u.mul(
268            &ExactNum::from_u8(BALL_TRANSCENDENTAL_ERROR_TERMS as u8, p),
269            p,
270            RoundingMode::Up,
271        )
272    }
273
274    fn nan_disk() -> Self {
275        let n = ExactNum::nan(Some(Error::InvalidArgument));
276        ComplexBall {
277            mid: ExactComplex::new(n.clone(), n.clone()),
278            rad: n,
279        }
280    }
281
282    /// Sum of two disks: radii add, plus rounding slack.
283    pub fn add(&self, other: &Self, p: usize, rm: RoundingMode) -> Self {
284        let mid = self.mid.add(&other.mid, p, rm);
285        let rad = self.rad.add(&other.rad, p, RoundingMode::Up).add(
286            &Self::slack(&mid, p),
287            p,
288            RoundingMode::Up,
289        );
290        ComplexBall { mid, rad }
291    }
292
293    /// Product of two disks: \(\lvert m_1\rvert r_2+\lvert m_2\rvert r_1+r_1 r_2\) plus slack.
294    pub fn mul(&self, other: &Self, p: usize, rm: RoundingMode) -> Self {
295        let mid = self.mid.mul(&other.mid, p, rm);
296        let a = self
297            .mid
298            .abs(p, RoundingMode::Up)
299            .mul(&other.rad, p, RoundingMode::Up);
300        let b = other
301            .mid
302            .abs(p, RoundingMode::Up)
303            .mul(&self.rad, p, RoundingMode::Up);
304        let c = self.rad.mul(&other.rad, p, RoundingMode::Up);
305        let rad = a
306            .add(&b, p, RoundingMode::Up)
307            .add(&c, p, RoundingMode::Up)
308            .add(&Self::slack(&mid, p), p, RoundingMode::Up);
309        ComplexBall { mid, rad }
310    }
311
312    /// Exponential. Lipschitz \(\lvert e^z\rvert\le \exp(\mathrm{Re}\,m+r)\).
313    pub fn exp(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
314        let mid = self.mid.exp(p, rm, cc);
315        let re_hi = self.mid.re().add(&self.rad, p, RoundingMode::Up);
316        let lip = re_hi.exp(p, RoundingMode::Up, cc);
317        let rad =
318            lip.mul(&self.rad, p, RoundingMode::Up)
319                .add(&Self::slack(&mid, p), p, RoundingMode::Up);
320        ComplexBall { mid, rad }
321    }
322
323    /// Principal logarithm. Domain: the disk must exclude \(0\).
324    /// Lipschitz \(1/(\lvert m\rvert-r)\).
325    pub fn ln(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
326        let am = self.mid.abs(p, RoundingMode::Down);
327        if matches!(am.cmp(&self.rad), Some(c) if c <= 0) || am.is_zero() || self.rad.is_nan() {
328            return Self::nan_disk();
329        }
330        let mid = self.mid.ln(p, rm, cc);
331        let den = am.sub(&self.rad, p, RoundingMode::Down);
332        let lip = ExactNum::from_u8(1, p).div(&den, p, RoundingMode::Up);
333        let rad =
334            lip.mul(&self.rad, p, RoundingMode::Up)
335                .add(&Self::slack(&mid, p), p, RoundingMode::Up);
336        ComplexBall { mid, rad }
337    }
338
339    /// Sine. \(\lvert\cos(z)\rvert\le\cosh(\lvert\mathrm{Im}\,m\rvert+r)\).
340    pub fn sin(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
341        let mid = self.mid.sin(p, rm, cc);
342        let im_hi = self.mid.im().abs().add(&self.rad, p, RoundingMode::Up);
343        let lip = im_hi.sinh_cosh(p, RoundingMode::Up, cc).1;
344        let rad =
345            lip.mul(&self.rad, p, RoundingMode::Up)
346                .add(&Self::slack(&mid, p), p, RoundingMode::Up);
347        ComplexBall { mid, rad }
348    }
349
350    /// Cosine. Same Lipschitz bound as `sin`.
351    pub fn cos(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
352        let mid = self.mid.cos(p, rm, cc);
353        let im_hi = self.mid.im().abs().add(&self.rad, p, RoundingMode::Up);
354        let lip = im_hi.sinh_cosh(p, RoundingMode::Up, cc).1;
355        let rad =
356            lip.mul(&self.rad, p, RoundingMode::Up)
357                .add(&Self::slack(&mid, p), p, RoundingMode::Up);
358        ComplexBall { mid, rad }
359    }
360
361    /// True when \(\lvert z-\mathrm{mid}\rvert\le\mathrm{rad}\).
362    pub fn contains(&self, z: &ExactComplex, p: usize) -> bool {
363        if z.is_nan() || self.mid.is_nan() || self.rad.is_nan() {
364            return false;
365        }
366        let d = z
367            .sub(&self.mid, p, RoundingMode::None)
368            .abs(p, RoundingMode::Up);
369        matches!(d.cmp(&self.rad), Some(c) if c <= 0)
370    }
371}
372
373/// Evaluate `compute` at increasing working precision until the result rounds uniquely
374/// to `p` bits (`try_set_precision`). Same retry budget as the transcendental kernel
375/// ([`crate::MAX_PREC_RETRY`]).
376pub fn ziv_round<F>(p: usize, rm: RoundingMode, mut compute: F) -> ExactNum
377where
378    F: FnMut(usize) -> ExactNum,
379{
380    let mut p_inc = WORD_BIT_SIZE;
381    let mut p_wrk = match round_p(p).checked_add(p_inc) {
382        Some(v) => v,
383        None => return ExactNum::nan(Some(Error::InvalidArgument)),
384    };
385    loop {
386        let mut v = compute(p_wrk);
387        if v.try_set_precision(p, rm, p_wrk) {
388            return v;
389        }
390        if bump_prec_retry(&mut p_wrk, &mut p_inc, p).is_err() {
391            return ExactNum::nan(Some(Error::PrecisionRetryExhausted));
392        }
393    }
394}
395
396/// Same Ziv loop as [`ziv_round`], with every input lifted to `p_wrk` before `compute`.
397///
398/// `MAX_PREC_RETRY` still bounds the number of bumps via [`bump_prec_retry`].
399pub fn ziv_round_vec<F>(p: usize, rm: RoundingMode, inputs: &[ExactNum], mut compute: F) -> ExactNum
400where
401    F: FnMut(usize, &[ExactNum]) -> ExactNum,
402{
403    let mut p_inc = WORD_BIT_SIZE;
404    let mut p_wrk = match round_p(p).checked_add(p_inc) {
405        Some(v) => v,
406        None => return ExactNum::nan(Some(Error::InvalidArgument)),
407    };
408    loop {
409        let mut xs = alloc::vec::Vec::with_capacity(inputs.len());
410        for x in inputs {
411            let mut y = x.clone();
412            if y.set_precision(p_wrk, RoundingMode::None).is_err() {
413                y = x.clone();
414            }
415            xs.push(y);
416        }
417        let mut v = compute(p_wrk, &xs);
418        if v.try_set_precision(p, rm, p_wrk) {
419            return v;
420        }
421        if bump_prec_retry(&mut p_wrk, &mut p_inc, p).is_err() {
422            return ExactNum::nan(Some(Error::PrecisionRetryExhausted));
423        }
424    }
425}
426
427#[cfg(test)]
428mod tests {
429    use super::*;
430
431    #[test]
432    fn ball_add_contains_true_sum() {
433        let p = 128;
434        let rm = RoundingMode::ToEven;
435        let a = ExactNum::from(3);
436        let b = ExactNum::from(4);
437        let u = ExactNum::from_word(1, p);
438        let ba = Ball::new(a.clone(), u.clone());
439        let bb = Ball::new(b.clone(), u.clone());
440        let sum = ba.add(&bb, p, rm);
441        let true_sum = a.add(&b, p, rm);
442        assert!(sum.contains(&true_sum, p));
443    }
444
445    #[test]
446    fn ziv_round_vec_hypot_atan2() {
447        let rm = RoundingMode::ToEven;
448        for p in [64usize, 128, 256] {
449            let a = ExactNum::from_u8(3, p);
450            let b = ExactNum::from_u8(4, p);
451            let got = ziv_round_vec(p, rm, &[a, b], |pw, xs| {
452                xs[0].hypot(&xs[1], pw, RoundingMode::None)
453            });
454            assert_eq!(got.cmp(&ExactNum::from_u8(5, p)), Some(0));
455        }
456
457        let p = 256;
458        let mut cc = Consts::new().unwrap();
459        let one = ExactNum::from_u8(1, p);
460        let got = ziv_round_vec(p, rm, &[one.clone(), one], |pw, xs| {
461            xs[0].atan2(&xs[1], pw, RoundingMode::None, &mut cc)
462        });
463        let quarter = cc.pi(p, rm).div(&ExactNum::from_u8(4, p), p, rm);
464        assert_eq!(got.cmp(&quarter), Some(0));
465    }
466
467    #[test]
468    fn ziv_round_sqrt_matches_direct() {
469        let p = 128;
470        let rm = RoundingMode::ToEven;
471        let two = ExactNum::from(2);
472        let via_ziv = ziv_round(p, rm, |pw| two.sqrt(pw, RoundingMode::None));
473        let direct = two.sqrt(p, rm);
474        assert_eq!(via_ziv.cmp(&direct), Some(0));
475    }
476
477    #[test]
478    fn ball_exp_contains_one_and_two() {
479        let p = 128;
480        let rm = RoundingMode::ToEven;
481        let mut cc = Consts::new().unwrap();
482        let two = ExactNum::from_u8(2, p);
483        let rad = two.powsi(-20, p, rm);
484        let z = Ball::new(ExactNum::from_u8(0, p), rad.clone());
485        let ez = z.exp(p, rm, &mut cc);
486        assert!(ez.contains(&ExactNum::from_u8(1, p), p));
487
488        let ln2 = cc.ln_2(p, rm);
489        let bln = Ball::new(ln2, rad);
490        let e2 = bln.exp(p, rm, &mut cc);
491        assert!(e2.contains(&two, p));
492    }
493
494    #[test]
495    fn ball_sin_contains_zero_at_origin_and_pi() {
496        let p = 128;
497        let rm = RoundingMode::ToEven;
498        let mut cc = Consts::new().unwrap();
499        let two = ExactNum::from_u8(2, p);
500        let rad = two.powsi(-20, p, rm);
501        let z = Ball::new(ExactNum::from_u8(0, p), rad.clone());
502        let sz = z.sin(p, rm, &mut cc);
503        assert!(sz.contains(&ExactNum::from_u8(0, p), p));
504
505        let pi = cc.pi(p, rm);
506        let bpi = Ball::new(pi, rad);
507        let sp = bpi.sin(p, rm, &mut cc);
508        assert!(sp.contains(&ExactNum::from_u8(0, p), p));
509    }
510
511    #[test]
512    fn ball_exp_sin_contain_scalar_at_a_point() {
513        let p = 128;
514        let rm = RoundingMode::ToEven;
515        let mut cc = Consts::new().unwrap();
516        let x = ExactNum::from_u8(1, p);
517        let rad = ExactNum::from_u8(2, p).powsi(-12, p, rm);
518        let b = Ball::new(x.clone(), rad);
519        let hi = x.exp(256, rm, &mut cc);
520        assert!(b.exp(p, rm, &mut cc).contains(&hi, p));
521        let hs = x.sin(256, rm, &mut cc);
522        assert!(b.sin(p, rm, &mut cc).contains(&hs, p));
523    }
524
525    #[test]
526    fn ball_plan_transcendental_golds() {
527        let p = 256;
528        let rm = RoundingMode::ToEven;
529        let mut cc = Consts::new().unwrap();
530        let rad = ExactNum::from_u8(1, p).ldexp(-(p as i32), p, RoundingMode::None);
531        let six = ExactNum::from_u8(6, p);
532        let pi6 = cc.pi(p, rm).div(&six, p, rm);
533        let bsin = Ball::new(pi6, rad.clone());
534        let half = ExactNum::from_u8(1, p).div(&ExactNum::from_u8(2, p), p, rm);
535        assert!(bsin.sin(p, rm, &mut cc).contains(&half, p));
536
537        let one = ExactNum::from_u8(1, p);
538        let be = Ball::new(one.clone(), rad.clone());
539        let e = cc.e(p, rm);
540        assert!(be.exp(p, rm, &mut cc).contains(&e, p));
541
542        let erf1 = one.erf(p, rm, &mut cc);
543        assert!(be.erf(p, rm, &mut cc).contains(&erf1, p));
544
545        let ln1 = one.ln(256, rm, &mut cc);
546        assert!(be.ln(p, rm, &mut cc).contains(&ln1, p));
547        let sq = one.sqrt(256, rm);
548        assert!(be.sqrt(p, rm).contains(&sq, p));
549        assert!(be
550            .cos(p, rm, &mut cc)
551            .contains(&one.cos(256, rm, &mut cc), p));
552        assert!(be
553            .bessel_j0(p, rm, &mut cc)
554            .contains(&one.bessel_j(0, 256, rm, &mut cc), p));
555        assert!(be
556            .bessel_j1(p, rm, &mut cc)
557            .contains(&one.bessel_j(1, 256, rm, &mut cc), p));
558
559        let composed = be.sin(p, rm, &mut cc).exp(p, rm, &mut cc);
560        let true_c = one.sin(256, rm, &mut cc).exp(256, rm, &mut cc);
561        assert!(composed.contains(&true_c, p));
562    }
563
564    #[test]
565    fn complex_ball_exp_and_pythagoras() {
566        let p = 256;
567        let rm = RoundingMode::ToEven;
568        let mut cc = Consts::new().unwrap();
569        let mid = ExactComplex::new(
570            ExactNum::from_u8(3, p).div(&ExactNum::from_u8(10, p), p, rm),
571            ExactNum::from_u8(2, p).div(&ExactNum::from_u8(10, p), p, rm),
572        );
573        let rad = ExactNum::from_u8(1, p);
574        let unit = ComplexBall::new(ExactComplex::zero(p), rad);
575        let e_mid = mid.exp(p, rm, &mut cc);
576        assert!(unit.exp(p, rm, &mut cc).contains(&e_mid, p));
577
578        let small = ExactNum::from_u8(1, p).ldexp(-40, p, RoundingMode::None);
579        let d = ComplexBall::new(mid.clone(), small);
580        let s = d.sin(p, rm, &mut cc);
581        let c = d.cos(p, rm, &mut cc);
582        let ss = s.mul(&s, p, rm);
583        let cc2 = c.mul(&c, p, rm);
584        let py = ss.add(&cc2, p, rm);
585        assert!(py.contains(&ExactComplex::one(p), p));
586    }
587}
zenith_float_num/ball.rs

zenith_float_num/
ball.rs
