ocas_domain/
real_ball.rs

1//! Real ball domain implementation.
2//!
3//! A real ball represents a real number together with a rigorous error bound.
4//! The default build uses lightweight conservative `f64` balls. When the `mpfr`
5//! feature is enabled, the implementation is backed by `rug::Float` and uses
6//! directed rounding to produce rigorous arbitrary-precision enclosures.
7//!
8//! The default `f64` version is suitable for templates and demonstration; it is
9//! not a rigorous arbitrary-precision interval library.
10
11#[cfg(feature = "mpfr")]
12use rug::ops::{AddAssignRound, DivAssignRound, Pow, SubAssignRound};
13
14use crate::domain::Domain;
15
16/// Precision used for the MPFR-backed real ball in bits.
17#[cfg(feature = "mpfr")]
18const BALL_PRECISION: u32 = 53;
19
20/// A real ball: midpoint with radius.
21///
22/// The default `f64` build implements `Copy`. The `mpfr` build stores an
23/// owned `rug::Float` and therefore only implements `Clone`.
24#[derive(Debug, Clone, PartialEq)]
25pub struct RealBall {
26    #[cfg(not(feature = "mpfr"))]
27    mid: f64,
28    #[cfg(not(feature = "mpfr"))]
29    rad: f64,
30    #[cfg(feature = "mpfr")]
31    mid: rug::Float,
32    #[cfg(feature = "mpfr")]
33    rad: rug::Float,
34}
35
36impl Eq for RealBall {}
37
38#[cfg(not(feature = "mpfr"))]
39impl Copy for RealBall {}
40
41impl RealBall {
42    /// Create a new ball from midpoint and radius.
43    ///
44    /// The radius is clamped to be non-negative.
45    #[cfg(not(feature = "mpfr"))]
46    pub fn new(mid: f64, rad: f64) -> Self {
47        Self {
48            mid,
49            rad: rad.max(0.0),
50        }
51    }
52
53    /// Create a new ball from midpoint and radius.
54    ///
55    /// The radius is clamped to be non-negative.
56    #[cfg(feature = "mpfr")]
57    pub fn new(mid: rug::Float, rad: rug::Float) -> Self {
58        let zero = mpfr_zero();
59        if rad < zero {
60            Self { mid, rad: zero }
61        } else {
62            Self { mid, rad }
63        }
64    }
65
66    /// Create a ball from an `f64` value with zero radius.
67    #[cfg(not(feature = "mpfr"))]
68    pub fn from_f64(value: f64) -> Self {
69        Self::new(value, 0.0)
70    }
71
72    /// Create a ball from an `f64` value with zero radius.
73    #[cfg(feature = "mpfr")]
74    pub fn from_f64(value: f64) -> Self {
75        Self::new(rug::Float::with_val(BALL_PRECISION, value), mpfr_zero())
76    }
77
78    /// Access the midpoint.
79    #[cfg(not(feature = "mpfr"))]
80    pub fn mid(&self) -> f64 {
81        self.mid
82    }
83
84    /// Access the midpoint.
85    #[cfg(feature = "mpfr")]
86    pub fn mid(&self) -> &rug::Float {
87        &self.mid
88    }
89
90    /// Access the radius.
91    #[cfg(not(feature = "mpfr"))]
92    pub fn rad(&self) -> f64 {
93        self.rad
94    }
95
96    /// Access the radius.
97    #[cfg(feature = "mpfr")]
98    pub fn rad(&self) -> &rug::Float {
99        &self.rad
100    }
101
102    /// Return a conservative lower bound.
103    #[cfg(not(feature = "mpfr"))]
104    pub fn lower(&self) -> f64 {
105        self.mid - self.rad
106    }
107
108    /// Return a conservative lower bound.
109    #[cfg(feature = "mpfr")]
110    pub fn lower(&self) -> rug::Float {
111        let mut lo = self.mid.clone();
112        lo.sub_assign_round(&self.rad, rug::float::Round::Down);
113        lo
114    }
115
116    /// Return a conservative upper bound.
117    #[cfg(not(feature = "mpfr"))]
118    pub fn upper(&self) -> f64 {
119        self.mid + self.rad
120    }
121
122    /// Return a conservative upper bound.
123    #[cfg(feature = "mpfr")]
124    pub fn upper(&self) -> rug::Float {
125        let mut hi = self.mid.clone();
126        hi.add_assign_round(&self.rad, rug::float::Round::Up);
127        hi
128    }
129
130    /// Return the precision of the MPFR backing store in bits.
131    #[cfg(feature = "mpfr")]
132    pub fn precision(&self) -> u32 {
133        self.mid.prec()
134    }
135}
136
137/// The real ball domain.
138///
139/// In the default build this wraps `f64`-based balls. When the `mpfr` feature
140/// is enabled, this domain uses `rug::Float` with rigorous rounding.
141#[derive(Debug, Clone, Copy, PartialEq, Eq)]
142pub struct RealBallDomain;
143
144impl Domain for RealBallDomain {
145    type Element = RealBall;
146
147    fn zero(&self) -> Self::Element {
148        RealBall::from_f64(0.0)
149    }
150
151    fn one(&self) -> Self::Element {
152        RealBall::from_f64(1.0)
153    }
154
155    fn add(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
156        #[cfg(not(feature = "mpfr"))]
157        {
158            RealBall::new(a.mid + b.mid, a.rad + b.rad)
159        }
160        #[cfg(feature = "mpfr")]
161        {
162            let mut mid = a.mid.clone();
163            let mut rad = a.rad.clone();
164            mid.add_assign_round(&b.mid, rug::float::Round::Down);
165            rad.add_assign_round(&b.rad, rug::float::Round::Up);
166            rad.add_assign_round(&rounding_unit(&mid), rug::float::Round::Up);
167            RealBall::new(mid, rad)
168        }
169    }
170
171    fn sub(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
172        #[cfg(not(feature = "mpfr"))]
173        {
174            RealBall::new(a.mid - b.mid, a.rad + b.rad)
175        }
176        #[cfg(feature = "mpfr")]
177        {
178            let mut mid = a.mid.clone();
179            let mut rad = a.rad.clone();
180            mid.sub_assign_round(&b.mid, rug::float::Round::Up);
181            rad.add_assign_round(&b.rad, rug::float::Round::Up);
182            rad.add_assign_round(&rounding_unit(&mid), rug::float::Round::Up);
183            RealBall::new(mid, rad)
184        }
185    }
186
187    fn neg(&self, a: &Self::Element) -> Self::Element {
188        #[cfg(not(feature = "mpfr"))]
189        {
190            RealBall::new(-a.mid, a.rad)
191        }
192        #[cfg(feature = "mpfr")]
193        {
194            RealBall::new(-a.mid.clone(), a.rad.clone())
195        }
196    }
197
198    fn mul(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
199        #[cfg(not(feature = "mpfr"))]
200        {
201            let a_lo = a.lower();
202            let a_hi = a.upper();
203            let b_lo = b.lower();
204            let b_hi = b.upper();
205            let p1 = a_lo * b_lo;
206            let p2 = a_lo * b_hi;
207            let p3 = a_hi * b_lo;
208            let p4 = a_hi * b_hi;
209            let lo = p1.min(p2).min(p3).min(p4);
210            let hi = p1.max(p2).max(p3).max(p4);
211            let mid = (lo + hi) / 2.0;
212            let rad = (hi - lo) / 2.0;
213            RealBall::new(mid, rad)
214        }
215        #[cfg(feature = "mpfr")]
216        {
217            let a_lo = a.lower();
218            let a_hi = a.upper();
219            let b_lo = b.lower();
220            let b_hi = b.upper();
221            let mut products = [
222                a_lo.clone() * b_lo.clone(),
223                a_lo * b_hi.clone(),
224                a_hi.clone() * b_lo,
225                a_hi * b_hi,
226            ];
227            for p in &mut products {
228                p.set_prec(BALL_PRECISION);
229            }
230            let lo = products.iter().cloned().reduce(mpfr_min).unwrap();
231            let hi = products.iter().cloned().reduce(mpfr_max).unwrap();
232            let mut mid = (lo.clone() + &hi) / 2u32;
233            mid.set_prec(BALL_PRECISION);
234            let mut rad = (hi - &lo).abs() / 2u32;
235            rad.set_prec(BALL_PRECISION);
236            rad.add_assign_round(&rounding_unit(&mid), rug::float::Round::Up);
237            RealBall::new(mid, rad)
238        }
239    }
240
241    fn div(&self, a: &Self::Element, b: &Self::Element) -> Option<Self::Element> {
242        let inv = self.inv(b)?;
243        Some(self.mul(a, &inv))
244    }
245
246    fn inv(&self, a: &Self::Element) -> Option<Self::Element> {
247        #[cfg(not(feature = "mpfr"))]
248        {
249            if a.lower() <= 0.0 && a.upper() >= 0.0 {
250                return None;
251            }
252            let lo = a.lower();
253            let hi = a.upper();
254            let inv_lo = 1.0 / hi;
255            let inv_hi = 1.0 / lo;
256            let mid = (inv_lo + inv_hi) / 2.0;
257            let rad = (inv_hi - inv_lo).abs() / 2.0;
258            Some(RealBall::new(mid, rad))
259        }
260        #[cfg(feature = "mpfr")]
261        {
262            let lo = a.lower();
263            let hi = a.upper();
264            if lo <= mpfr_zero() && hi >= mpfr_zero() {
265                return None;
266            }
267            let mut inv_lo = rug::Float::with_val(BALL_PRECISION, 1.0);
268            inv_lo.div_assign_round(&hi, rug::float::Round::Down);
269            let mut inv_hi = rug::Float::with_val(BALL_PRECISION, 1.0);
270            inv_hi.div_assign_round(&lo, rug::float::Round::Up);
271            let mut mid = (inv_lo.clone() + &inv_hi) / 2u32;
272            mid.set_prec(BALL_PRECISION);
273            let mut rad = (inv_hi - &inv_lo).abs() / 2u32;
274            rad.set_prec(BALL_PRECISION);
275            rad.add_assign_round(&rounding_unit(&mid), rug::float::Round::Up);
276            Some(RealBall::new(mid, rad))
277        }
278    }
279}
280
281#[cfg(feature = "mpfr")]
282fn mpfr_zero() -> rug::Float {
283    rug::Float::with_val(BALL_PRECISION, 0.0)
284}
285
286/// Rounding-unit helper for `rug::Float`.
287#[cfg(feature = "mpfr")]
288fn rounding_unit(x: &rug::Float) -> rug::Float {
289    let p = x.prec();
290    let mut u = rug::Float::with_val(p, 1.0);
291    let denom = rug::Float::with_val(p, 2.0).pow(p - 1);
292    u.div_assign_round(&denom, rug::float::Round::Up);
293    u
294}
295
296/// Minimum of two `rug::Float` values.
297#[cfg(feature = "mpfr")]
298fn mpfr_min(a: rug::Float, b: rug::Float) -> rug::Float {
299    if a < b { a } else { b }
300}
301
302/// Maximum of two `rug::Float` values.
303#[cfg(feature = "mpfr")]
304fn mpfr_max(a: rug::Float, b: rug::Float) -> rug::Float {
305    if a > b { a } else { b }
306}
307
308#[cfg(test)]
309mod tests {
310    use super::*;
311
312    #[cfg(not(feature = "mpfr"))]
313    #[test]
314    fn ball_default_addition_widens_radius() {
315        let domain = RealBallDomain;
316        let a = RealBall::from_f64(1.0);
317        let b = RealBall::from_f64(2.0);
318        let sum = domain.add(&a, &b);
319        assert!((sum.mid() - 3.0).abs() < 1e-12);
320        assert!((sum.rad() - 0.0).abs() < 1e-12);
321    }
322
323    #[cfg(not(feature = "mpfr"))]
324    #[test]
325    fn ball_default_multiplication_contains_true_product() {
326        let domain = RealBallDomain;
327        let a = RealBall::from_f64(2.0);
328        let b = RealBall::from_f64(3.0);
329        let prod = domain.mul(&a, &b);
330        assert!(prod.lower() <= 6.0 && 6.0 <= prod.upper());
331    }
332
333    #[cfg(not(feature = "mpfr"))]
334    #[test]
335    fn ball_default_inverse_excludes_zero() {
336        let domain = RealBallDomain;
337        let a = RealBall::from_f64(2.0);
338        let inv = domain.inv(&a).expect("ball does not contain zero");
339        let prod = domain.mul(&a, &inv);
340        assert!(prod.lower() <= 1.0 && 1.0 <= prod.upper());
341    }
342
343    #[cfg(not(feature = "mpfr"))]
344    #[test]
345    fn ball_default_inverse_of_zero_ball_is_none() {
346        let domain = RealBallDomain;
347        let a = RealBall::from_f64(0.0);
348        assert!(domain.inv(&a).is_none());
349    }
350
351    #[cfg(not(feature = "mpfr"))]
352    #[test]
353    fn ball_default_contains_zero_leads_to_none() {
354        let domain = RealBallDomain;
355        let a = RealBall::new(1.0, 0.0);
356        let b = RealBall::new(-0.05, 0.1);
357        assert!(domain.div(&a, &b).is_none());
358    }
359
360    #[cfg(feature = "mpfr")]
361    mod mpfr_tests {
362        use super::*;
363
364        #[test]
365        fn ball_mpfr_contains_exact_value() {
366            let domain = RealBallDomain;
367            let a = RealBall::from_f64(1.0);
368            let b = RealBall::from_f64(2.0);
369            let sum = domain.add(&a, &b);
370            let lo = sum.lower();
371            let hi = sum.upper();
372            assert!(lo <= rug::Float::with_val(BALL_PRECISION, 3.0));
373            assert!(hi >= rug::Float::with_val(BALL_PRECISION, 3.0));
374        }
375
376        #[test]
377        fn ball_mpfr_inverse_excludes_zero() {
378            let domain = RealBallDomain;
379            let a = RealBall::from_f64(2.0);
380            let inv = domain.inv(&a).expect("invertible");
381            let prod = domain.mul(&a, &inv);
382            let lo = prod.lower();
383            let hi = prod.upper();
384            assert!(lo <= rug::Float::with_val(BALL_PRECISION, 1.0));
385            assert!(hi >= rug::Float::with_val(BALL_PRECISION, 1.0));
386        }
387    }
388}
ocas_domain/real_ball.rs

ocas_domain/
real_ball.rs
