rival/interval/
arithmetic.rs

1use super::value::{Endpoint, Ival, IvalClass, classify};
2use crate::{
3    interval::core::endpoint_binary,
4    mpfr::{mpfr_add, mpfr_div, mpfr_hypot, mpfr_mul, mpfr_sub},
5};
6use rug::{Assign, Float, float::Round, ops::AssignRound};
7
8impl Ival {
9    /// Compute the interval sum `a + b`.
10    pub fn add_assign(&mut self, a: &Ival, b: &Ival) {
11        self.lo.immovable =
12            endpoint_binary(mpfr_add, &a.lo, &b.lo, self.lo.as_float_mut(), Round::Down);
13        self.hi.immovable =
14            endpoint_binary(mpfr_add, &a.hi, &b.hi, self.hi.as_float_mut(), Round::Up);
15        self.err = a.err.union(&b.err);
16    }
17
18    /// Compute the interval difference `a - b`.
19    pub fn sub_assign(&mut self, a: &Ival, b: &Ival) {
20        self.lo.immovable =
21            endpoint_binary(mpfr_sub, &a.lo, &b.hi, self.lo.as_float_mut(), Round::Down);
22        self.hi.immovable =
23            endpoint_binary(mpfr_sub, &a.hi, &b.lo, self.hi.as_float_mut(), Round::Up);
24        self.err = a.err.union(&b.err);
25    }
26
27    /// Compute the interval product `a * b`.
28    pub fn mul_assign(&mut self, a: &Ival, b: &Ival) {
29        let class_a = classify(a, false);
30        let class_b = classify(b, false);
31        let err = a.err.union(&b.err);
32
33        let mkmul =
34            |out: &mut Ival, lo_a: &Endpoint, lo_b: &Endpoint, hi_a: &Endpoint, hi_b: &Endpoint| {
35                out.lo.immovable = epmul(
36                    lo_a,
37                    lo_b,
38                    class_a,
39                    class_b,
40                    out.lo.as_float_mut(),
41                    Round::Down,
42                );
43                out.hi.immovable = epmul(
44                    hi_a,
45                    hi_b,
46                    class_a,
47                    class_b,
48                    out.hi.as_float_mut(),
49                    Round::Up,
50                );
51                out.err = err;
52            };
53
54        match (class_a, class_b) {
55            (IvalClass::Pos, IvalClass::Pos) => mkmul(self, &a.lo, &b.lo, &a.hi, &b.hi),
56            (IvalClass::Pos, IvalClass::Neg) => mkmul(self, &a.hi, &b.lo, &a.lo, &b.hi),
57            (IvalClass::Pos, IvalClass::Mix) => mkmul(self, &a.hi, &b.lo, &a.hi, &b.hi),
58            (IvalClass::Neg, IvalClass::Pos) => mkmul(self, &a.lo, &b.hi, &a.hi, &b.lo),
59            (IvalClass::Neg, IvalClass::Neg) => mkmul(self, &a.hi, &b.hi, &a.lo, &b.lo),
60            (IvalClass::Neg, IvalClass::Mix) => mkmul(self, &a.lo, &b.hi, &a.lo, &b.lo),
61            (IvalClass::Mix, IvalClass::Pos) => mkmul(self, &a.lo, &b.hi, &a.hi, &b.hi),
62            (IvalClass::Mix, IvalClass::Neg) => mkmul(self, &a.hi, &b.lo, &a.lo, &b.lo),
63            (IvalClass::Mix, IvalClass::Mix) => {
64                let mut tmp_ival = self.clone();
65                mkmul(self, &a.hi, &b.lo, &a.lo, &b.lo);
66                mkmul(&mut tmp_ival, &a.lo, &b.hi, &a.hi, &b.hi);
67                self.union_assign(tmp_ival);
68            }
69        }
70    }
71
72    /// Compute the interval quotient `a / b`.
73    pub fn div_assign(&mut self, a: &Ival, b: &Ival) {
74        let class_a = classify(a, true);
75        let class_b = classify(b, true);
76
77        let y_lo = b.lo.as_float();
78        let y_hi = b.hi.as_float();
79
80        self.err = a.err.union(&b.err);
81        if (y_lo.is_zero() || y_lo.is_sign_negative())
82            && (y_hi.is_zero() || y_hi.is_sign_positive())
83        {
84            self.err.partial = true;
85        }
86        if y_lo.is_zero() && y_hi.is_zero() {
87            self.err.total = true;
88            self.err.partial = true;
89        }
90
91        if matches!(class_b, IvalClass::Mix) {
92            let immovable = (b.lo.immovable && y_lo.is_zero())
93                || (b.hi.immovable && y_hi.is_zero())
94                || (b.lo.immovable && b.hi.immovable);
95            self.lo.as_float_mut().assign(f64::NEG_INFINITY);
96            self.hi.as_float_mut().assign(f64::INFINITY);
97            self.lo.immovable = immovable;
98            self.hi.immovable = immovable;
99            return;
100        }
101
102        let mut mkdiv = |num_lo: &Endpoint,
103                         den_lo: &Endpoint,
104                         num_hi: &Endpoint,
105                         den_hi: &Endpoint| {
106            self.lo.immovable = epdiv(num_lo, den_lo, class_a, self.lo.as_float_mut(), Round::Down);
107            self.hi.immovable = epdiv(num_hi, den_hi, class_a, self.hi.as_float_mut(), Round::Up);
108        };
109
110        match (class_a, class_b) {
111            (_, IvalClass::Mix) => unreachable!(),
112            (IvalClass::Pos, IvalClass::Pos) => mkdiv(&a.lo, &b.hi, &a.hi, &b.lo),
113            (IvalClass::Pos, IvalClass::Neg) => mkdiv(&a.hi, &b.hi, &a.lo, &b.lo),
114            (IvalClass::Neg, IvalClass::Pos) => mkdiv(&a.lo, &b.lo, &a.hi, &b.hi),
115            (IvalClass::Neg, IvalClass::Neg) => mkdiv(&a.hi, &b.lo, &a.lo, &b.hi),
116            (IvalClass::Mix, IvalClass::Pos) => mkdiv(&a.lo, &b.lo, &a.hi, &b.lo),
117            (IvalClass::Mix, IvalClass::Neg) => mkdiv(&a.hi, &b.hi, &a.lo, &b.hi),
118        }
119    }
120
121    /// Compute the interval fused multiply-add `a * b + c`.
122    pub fn fma_assign(&mut self, a: &Ival, b: &Ival, c: &Ival) {
123        let mut product = Ival::zero(self.prec());
124        product.mul_assign(a, b);
125        self.add_assign(&product, c);
126    }
127
128    /// Compute the interval positive difference `fdim(x, y) = max(x - y, 0)`.
129    pub fn fdim_assign(&mut self, x: &Ival, y: &Ival) {
130        let mut diff = Ival::zero(self.prec());
131        diff.sub_assign(x, y);
132        let zero_ival = Ival::zero(self.prec());
133        self.fmax_assign(&diff, &zero_ival);
134    }
135
136    /// Compute the interval Euclidean distance `hypot(x, y) = sqrt(x² + y²)`.
137    pub fn hypot_assign(&mut self, x: &Ival, y: &Ival) {
138        let mut abs_x = Ival::zero(self.prec());
139        let mut abs_y = Ival::zero(self.prec());
140        abs_x.pre_fabs_assign(x);
141        abs_y.pre_fabs_assign(y);
142
143        self.lo.immovable = endpoint_binary(
144            mpfr_hypot,
145            &abs_x.lo,
146            &abs_y.lo,
147            self.lo.as_float_mut(),
148            Round::Down,
149        );
150        self.hi.immovable = endpoint_binary(
151            mpfr_hypot,
152            &abs_x.hi,
153            &abs_y.hi,
154            self.hi.as_float_mut(),
155            Round::Up,
156        );
157        self.err = x.err.union(&y.err);
158    }
159}
160
161fn epmul(
162    ep1: &Endpoint,
163    ep2: &Endpoint,
164    class1: IvalClass,
165    class2: IvalClass,
166    out: &mut Float,
167    rnd: Round,
168) -> bool {
169    let a = ep1.as_float();
170    let b = ep2.as_float();
171
172    if a.is_zero() || b.is_zero() {
173        out.assign_round(0.0, Round::Nearest);
174        return (ep1.immovable && a.is_zero())
175            || (ep2.immovable && b.is_zero())
176            || (ep1.immovable && ep2.immovable);
177    }
178
179    let exact = mpfr_mul(a, b, out, rnd);
180    (ep1.immovable && ep2.immovable && exact)
181        || (ep1.immovable && a.is_infinite() && !matches!(class2, IvalClass::Mix))
182        || (ep2.immovable && b.is_infinite() && !matches!(class1, IvalClass::Mix))
183}
184
185fn epdiv(ep1: &Endpoint, ep2: &Endpoint, class: IvalClass, out: &mut Float, rnd: Round) -> bool {
186    let a = ep1.as_float();
187    let b = ep2.as_float();
188
189    let exact = mpfr_div(a, b, out, rnd);
190
191    (ep1.immovable && ep2.immovable && exact)
192        || (ep1.immovable && (a.is_zero() || a.is_infinite()))
193        || (ep2.immovable && (b.is_infinite() || (b.is_zero() && !matches!(class, IvalClass::Mix))))
194}
rival/interval/arithmetic.rs

rival/interval/
arithmetic.rs
