inari/
arith.rs

1use crate::{classify::*, interval::*, simd::*};
2use forward_ref::*;
3use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
4
5impl Neg for Interval {
6    type Output = Self;
7
8    fn neg(self) -> Self {
9        // [-b, -a] = [b; -a]
10        Self {
11            rep: swap(self.rep),
12        }
13    }
14}
15
16forward_ref_unop!(impl Neg, neg for Interval);
17
18impl Add for Interval {
19    type Output = Self;
20
21    fn add(self, rhs: Self) -> Self {
22        // [a + c, b + d] = [-a - c; b + d] = [-a; b] .+ [-c; d]
23        let x = self.rep; // [-a; b]
24        let y = rhs.rep; // [-c; d]
25        Self { rep: add_ru(x, y) }
26    }
27}
28
29forward_ref_binop!(impl Add, add for Interval, Interval);
30
31impl Sub for Interval {
32    type Output = Self;
33
34    fn sub(self, rhs: Self) -> Self {
35        // [a - d, b - c] = [-a + d; b - c] = [-a; b] .+ [d; -c]
36        let x = self.rep; // [-a; b]
37        let y = swap(rhs.rep); // [d; -c]
38        Self { rep: add_ru(x, y) }
39    }
40}
41
42forward_ref_binop!(impl Sub, sub for Interval, Interval);
43
44impl Mul for Interval {
45    type Output = Self;
46
47    #[allow(clippy::many_single_char_names)]
48    fn mul(self, rhs: Self) -> Self {
49        // [a, b] * [c, d] =
50        //
51        //    |      M     |      N     |      P     |  Z
52        // ---+------------+------------+------------+-----
53        //  M |     *1     | [b*c, a*c] | [a*d, b*d] | {0}
54        //  N | [a*d, a*c] | [b*d, a*c] | [a*d, b*c] | {0}
55        //  P | [b*c, b*d] | [b*c, a*d] | [a*c, b*d] | {0}
56        //  Z |     {0}    |     {0}    |     {0}    | {0}
57        // *1 [min{a*d, b*c}, max{a*c, b*d}]
58
59        use IntervalClass2::*;
60        match self.classify2(rhs) {
61            E_E | E_M | E_N0 | E_N1 | E_P0 | E_P1 | E_Z | M_E | N0_E | N1_E | P0_E | P1_E | Z_E => {
62                Self::EMPTY
63            }
64            M_Z | N0_Z | N1_Z | P0_Z | P1_Z | Z_M | Z_N0 | Z_N1 | Z_P0 | Z_P1 | Z_Z => Self::zero(),
65            M_M => {
66                // M * M => [min(a*d, b*c), max(a*c, b*d)]
67                //   = [-min(a*d, b*c); max(a*c, b*d)]
68                //   = [max(-a*d, -b*c); max(a*c, b*d)]
69                //   = [max(-a*d, b*-c); max(-a*-c, b*d)]
70                //   = .max([-a*d; -a*-c], [b*-c; b*d])
71                //   = .max([-a; -a] .* [d; -c], [b; b] .* [-c; d])
72                let x = shuffle02(self.rep, self.rep); // [-a; -a]
73                let y = swap(rhs.rep); // [d; -c]
74                let xy = mul_ru(x, y);
75                let z = shuffle13(self.rep, self.rep); // [b; b]
76                let w = rhs.rep; // [-c; d]
77                let zw = mul_ru(z, w);
78                let r = max(xy, zw);
79                Self { rep: r }
80            }
81            M_N0 | M_N1 => {
82                // M * N => [b*c, a*c] = [-b*c; a*c] = [b*-c; -a*-c] = [b; -a] .* [-c; -c]
83                let x = swap(self.rep); // [b; -a]
84                let y = shuffle02(rhs.rep, rhs.rep); // [-c; -c]
85                Self { rep: mul_ru(x, y) }
86            }
87            M_P0 | M_P1 => {
88                // M * P => [a*d, b*d] = [-a*d; b*d] = [-a; b] .* [d; d]
89                let x = self.rep; // [-a; b]
90                let y = shuffle13(rhs.rep, rhs.rep); // [d; d]
91                Self { rep: mul_ru(x, y) }
92            }
93            N0_M | N1_M => {
94                // N * M => [a*d, a*c] = [-a*d; a*c] = [-a*d; -a*-c] = [-a; -a] .* [d; -c]
95                let x = shuffle02(self.rep, self.rep); // [-a; -a]
96                let y = swap(rhs.rep); // [d; -c]
97                Self { rep: mul_ru(x, y) }
98            }
99            N0_N0 | N0_N1 | N1_N0 | N1_N1 => {
100                // N * N => [b*d, a*c] = [-b*d; a*c] = [-b*d; -a*-c] = [-b; -a] .* [d; -c]
101                let x = swap(self.rep); // [b; -a]
102                let x = neg0(x); // [-b; -a]
103                let y = swap(rhs.rep); // [d; -c]
104                Self { rep: mul_ru(x, y) }
105            }
106            N0_P0 | N0_P1 | N1_P0 | N1_P1 => {
107                // N * P => [a*d, b*c] = [-a*d; b*c] = [-a; b] .* [d; c]
108                let x = self.rep; // [-a; b]
109                let y = neg0(rhs.rep); // [c; d]
110                let y = swap(y); // [d; c]
111                Self { rep: mul_ru(x, y) }
112            }
113            P0_M | P1_M => {
114                // P * M => [b*c, b*d] = [-b*c; b*d] = [b*-c; b*d] = [b; b] .* [-c; d]
115                let x = shuffle13(self.rep, self.rep); // [b; b]
116                let y = rhs.rep; // [-c; d]
117                Self { rep: mul_ru(x, y) }
118            }
119            P0_N0 | P0_N1 | P1_N0 | P1_N1 => {
120                // P * N => [b*c, a*d] = [-b*c; a*d] = [b*-c; a*d] = [b; a] .* [-c; d]
121                let x = neg0(self.rep); // [a; b]
122                let x = swap(x); // [b; a]
123                let y = rhs.rep; // [-c; d]
124                Self { rep: mul_ru(x, y) }
125            }
126            P0_P0 | P0_P1 | P1_P0 | P1_P1 => {
127                // P * P => [a*c, b*d] = [-a*c; b*d] = [-a; b] .* [c; d]
128                let x = self.rep; // [-a; b]
129                let y = neg0(rhs.rep); // [c; d]
130                Self { rep: mul_ru(x, y) }
131            }
132        }
133    }
134}
135
136forward_ref_binop!(impl Mul, mul for Interval, Interval);
137
138impl Div for Interval {
139    type Output = Self;
140
141    fn div(self, rhs: Self) -> Self {
142        // [a, b] / [c, d] =
143        //
144        //    |  M  |     N0    |     N1     |     P0    |     P1     | Z
145        // ---+-----+-----------+------------+-----------+------------+---
146        //  M |  ℝ  |     ℝ     | [b/d, a/d] |     ℝ     | [a/c, b/c] | ∅
147        //  N |  ℝ  | [b/c, +∞] | [b/c, a/d] | [-∞, b/d] | [a/c, b/d] | ∅
148        //  P |  ℝ  | [-∞, a/c] | [b/d, a/c] | [a/d, +∞] | [a/d, b/c] | ∅
149        //  Z | {0} |    {0}    |     {0}    |    {0}    |     {0}    | ∅
150
151        use IntervalClass2::*;
152        match self.classify2(rhs) {
153            E_E | E_M | E_N0 | E_N1 | E_P0 | E_P1 | E_Z | M_E | M_Z | N0_E | N0_Z | N1_E | N1_Z
154            | P0_E | P0_Z | P1_E | P1_Z | Z_E | Z_Z => Self::EMPTY,
155            M_M | M_N0 | M_P0 | N0_M | N1_M | P0_M | P1_M => Self::ENTIRE,
156            Z_M | Z_N0 | Z_N1 | Z_P0 | Z_P1 => Self::zero(),
157            M_N1 => self.div_m_n1(rhs),
158            M_P1 => self.div_m_p1(rhs),
159            N0_N0 | N1_N0 => self.div_n_n0(rhs),
160            N0_N1 | N1_N1 => self.div_n_n1(rhs),
161            N0_P0 | N1_P0 => self.div_n_p0(rhs),
162            N0_P1 | N1_P1 => self.div_n_p1(rhs),
163            P0_N0 | P1_N0 => self.div_p_n0(rhs),
164            P0_N1 | P1_N1 => self.div_p_n1(rhs),
165            P0_P0 | P1_P0 => self.div_p_p0(rhs),
166            P0_P1 | P1_P1 => self.div_p_p1(rhs),
167        }
168    }
169}
170
171forward_ref_binop!(impl Div, div for Interval, Interval);
172
173impl Neg for DecInterval {
174    type Output = Self;
175
176    fn neg(self) -> Self {
177        if self.is_nai() {
178            return self;
179        }
180
181        Self::set_dec(-self.x, self.d)
182    }
183}
184
185forward_ref_unop!(impl Neg, neg for DecInterval);
186
187impl Add for DecInterval {
188    type Output = Self;
189
190    fn add(self, rhs: Self) -> Self {
191        if self.is_nai() || rhs.is_nai() {
192            return Self::NAI;
193        }
194
195        Self::set_dec(self.x + rhs.x, self.d.min(rhs.d))
196    }
197}
198
199forward_ref_binop!(impl Add, add for DecInterval, DecInterval);
200
201impl Sub for DecInterval {
202    type Output = Self;
203
204    fn sub(self, rhs: Self) -> Self {
205        if self.is_nai() || rhs.is_nai() {
206            return Self::NAI;
207        }
208
209        Self::set_dec(self.x - rhs.x, self.d.min(rhs.d))
210    }
211}
212
213forward_ref_binop!(impl Sub, sub for DecInterval, DecInterval);
214
215impl Mul for DecInterval {
216    type Output = Self;
217
218    fn mul(self, rhs: Self) -> Self {
219        if self.is_nai() || rhs.is_nai() {
220            return Self::NAI;
221        }
222
223        Self::set_dec(self.x * rhs.x, self.d.min(rhs.d))
224    }
225}
226
227forward_ref_binop!(impl Mul, mul for DecInterval, DecInterval);
228
229impl Div for DecInterval {
230    type Output = Self;
231
232    fn div(self, rhs: Self) -> Self {
233        if self.is_nai() || rhs.is_nai() {
234            return Self::NAI;
235        }
236
237        let d = if rhs.x.contains(0.0) {
238            Decoration::Trv
239        } else {
240            self.d.min(rhs.d)
241        };
242        Self::set_dec(self.x / rhs.x, d)
243    }
244}
245
246forward_ref_binop!(impl Div, div for DecInterval, DecInterval);
247
248macro_rules! impl_op_assign {
249    ($OpAssign:ident, $op_assign:ident, $op:ident) => {
250        impl $OpAssign for Interval {
251            fn $op_assign(&mut self, rhs: Self) {
252                *self = self.$op(rhs);
253            }
254        }
255
256        forward_ref_op_assign!(impl $OpAssign, $op_assign for Interval,
257                               Interval);
258
259        impl $OpAssign for DecInterval {
260            fn $op_assign(&mut self, rhs: Self) {
261                *self = self.$op(rhs);
262            }
263        }
264
265        forward_ref_op_assign!(impl $OpAssign, $op_assign for DecInterval,
266                               DecInterval);
267    };
268}
269
270impl_op_assign!(AddAssign, add_assign, add);
271impl_op_assign!(SubAssign, sub_assign, sub);
272impl_op_assign!(MulAssign, mul_assign, mul);
273impl_op_assign!(DivAssign, div_assign, div);
274
275#[cfg(test)]
276mod tests {
277    use crate::*;
278    use DecInterval as DI;
279    use Interval as I;
280
281    #[test]
282    fn add_assign() {
283        let mut i = const_interval!(3.0, 4.0);
284        i += const_interval!(1.0, 2.0);
285        assert_eq!(i, const_interval!(4.0, 6.0));
286
287        let mut i = const_dec_interval!(3.0, 4.0);
288        i += const_dec_interval!(1.0, 2.0);
289        assert_eq!(i, const_dec_interval!(4.0, 6.0));
290    }
291
292    #[test]
293    fn sub_assign() {
294        let mut i = const_interval!(3.0, 4.0);
295        i -= const_interval!(1.0, 2.0);
296        assert_eq!(i, const_interval!(1.0, 3.0));
297
298        let mut i = const_dec_interval!(3.0, 4.0);
299        i -= const_dec_interval!(1.0, 2.0);
300        assert_eq!(i, const_dec_interval!(1.0, 3.0));
301    }
302
303    #[test]
304    fn mul_assign() {
305        let mut i = const_interval!(3.0, 4.0);
306        i *= const_interval!(1.0, 2.0);
307        assert_eq!(i, const_interval!(3.0, 8.0));
308
309        let mut i = const_dec_interval!(3.0, 4.0);
310        i *= const_dec_interval!(1.0, 2.0);
311        assert_eq!(i, const_dec_interval!(3.0, 8.0));
312    }
313
314    #[test]
315    fn div_assign() {
316        let mut i = const_interval!(3.0, 4.0);
317        i /= const_interval!(1.0, 2.0);
318        assert_eq!(i, const_interval!(1.5, 4.0));
319
320        let mut i = const_dec_interval!(3.0, 4.0);
321        i /= const_dec_interval!(1.0, 2.0);
322        assert_eq!(i, const_dec_interval!(1.5, 4.0));
323    }
324
325    #[test]
326    fn empty() {
327        assert!((-I::EMPTY).is_empty());
328
329        assert!((I::EMPTY + I::PI).is_empty());
330        assert!((I::PI + I::EMPTY).is_empty());
331
332        assert!((I::EMPTY - I::PI).is_empty());
333        assert!((I::PI - I::EMPTY).is_empty());
334
335        assert!((I::EMPTY * I::PI).is_empty());
336        assert!((I::PI * I::EMPTY).is_empty());
337
338        assert!((I::EMPTY / I::PI).is_empty());
339        assert!((I::PI / I::EMPTY).is_empty());
340
341        assert!((-DI::EMPTY).is_empty());
342
343        assert!((DI::EMPTY + DI::PI).is_empty());
344        assert!((DI::PI + DI::EMPTY).is_empty());
345
346        assert!((DI::EMPTY - DI::PI).is_empty());
347        assert!((DI::PI - DI::EMPTY).is_empty());
348
349        assert!((DI::EMPTY * DI::PI).is_empty());
350        assert!((DI::PI * DI::EMPTY).is_empty());
351
352        assert!((DI::EMPTY / DI::PI).is_empty());
353        assert!((DI::PI / DI::EMPTY).is_empty());
354    }
355
356    #[test]
357    fn nai() {
358        assert!((-DI::NAI).is_nai());
359
360        assert!((DI::NAI + DI::PI).is_nai());
361        assert!((DI::PI + DI::NAI).is_nai());
362
363        assert!((DI::NAI - DI::PI).is_nai());
364        assert!((DI::PI - DI::NAI).is_nai());
365
366        assert!((DI::NAI * DI::PI).is_nai());
367        assert!((DI::PI * DI::NAI).is_nai());
368
369        assert!((DI::NAI / DI::PI).is_nai());
370        assert!((DI::PI / DI::NAI).is_nai());
371    }
372
373    #[allow(clippy::op_ref)]
374    #[test]
375    fn ref_type_args() {
376        const E: I = I::EMPTY;
377        const DE: DI = DI::EMPTY;
378
379        let _ = -&E;
380
381        let _ = &E + E;
382        let _ = E + &E;
383        let _ = &E + &E;
384
385        let _ = &E - E;
386        let _ = E - &E;
387        let _ = &E - &E;
388
389        let _ = &E * E;
390        let _ = E * &E;
391        let _ = &E * &E;
392
393        let _ = &E / E;
394        let _ = E / &E;
395        let _ = &E / &E;
396
397        let _ = -&DE;
398
399        let _ = &DE + DE;
400        let _ = DE + &DE;
401        let _ = &DE + &DE;
402
403        let _ = &DE - DE;
404        let _ = DE - &DE;
405        let _ = &DE - &DE;
406
407        let _ = &DE * DE;
408        let _ = DE * &DE;
409        let _ = &DE * &DE;
410
411        let _ = &DE / DE;
412        let _ = DE / &DE;
413        let _ = &DE / &DE;
414
415        let mut e = I::EMPTY;
416        e += &E;
417        e -= &E;
418        e *= &E;
419        e /= &E;
420
421        let mut de = DI::EMPTY;
422        de += &DE;
423        de -= &DE;
424        de *= &DE;
425        de /= &DE;
426    }
427}
inari/arith.rs

inari/
arith.rs
