1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

use crate::classify::*;
use crate::interval::*;
use crate::simd::*;
use core::arch::x86_64::*;

impl Interval {
    pub fn abs(self) -> Self {
        match self.classify() {
            C_E | C_P0 | C_P1 | C_Z => self,
            C_M => {
                // [0, max(-a, b)] = [max(-a, b); 0]
                let r0 = self.rep; // [b; -a]
                let r1 = unsafe { _mm_max_sd(r0, swap(r0)) }; // [_; max(-a, b)]
                let r = unsafe { _mm_unpacklo_pd(_mm_setzero_pd(), r1) };
                Self { rep: r }
            }
            C_N0 | C_N1 => {
                // [-b, -a] = [-a; b]
                Self {
                    rep: swap(self.rep),
                }
            }
            _ => unreachable!(),
        }
    }

    pub fn max(self, rhs: Self) -> Self {
        // IEEE 754's min/max do not propagate nan.
        // https://www.agner.org/optimize/blog/read.php?i=1012
        if self.is_empty() || rhs.is_empty() {
            return Self::empty();
        }

        unsafe {
            let max = _mm_max_pd(self.rep, rhs.rep); // [max(b, d); max(-a, -c)]
            let min = _mm_min_pd(self.rep, rhs.rep); // [min(b, d); min(-a, -c)]
            let r = _mm_move_sd(max, min); // [min(-a, -c); max(b, d)]
            Self { rep: r }
        }
    }

    pub fn min(self, rhs: Self) -> Self {
        if self.is_empty() || rhs.is_empty() {
            return Self::empty();
        }

        unsafe {
            let min = _mm_min_pd(self.rep, rhs.rep); // [min(b, d); min(-a, -c)]
            let max = _mm_max_pd(self.rep, rhs.rep); // [max(b, d); max(-a, -c)]
            let r = _mm_move_sd(min, max); // [max(-a, -c); min(b, d)]
            Self { rep: r }
        }
    }
}

impl DecoratedInterval {
    pub fn abs(self) -> Self {
        if self.is_nai() {
            return self;
        }

        Self::set_dec(self.x.abs(), self.d)
    }

    pub fn max(self, rhs: Self) -> Self {
        if self.is_nai() {
            return self;
        }

        Self::set_dec(self.x.max(rhs.x), self.d.min(rhs.d))
    }

    pub fn min(self, rhs: Self) -> Self {
        if self.is_nai() {
            return self;
        }

        Self::set_dec(self.x.min(rhs.x), self.d.min(rhs.d))
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    type DI = DecoratedInterval;
    type I = Interval;

    #[test]
    fn empty() {
        assert!(I::empty().abs().is_empty());

        assert!(I::empty().max(I::PI).is_empty());
        assert!(I::PI.max(I::empty()).is_empty());

        assert!(I::empty().min(I::PI).is_empty());
        assert!(I::PI.min(I::empty()).is_empty());

        assert!(DI::empty().abs().is_empty());

        assert!(DI::empty().max(DI::PI).is_empty());
        assert!(DI::PI.max(DI::empty()).is_empty());

        assert!(DI::empty().min(DI::PI).is_empty());
        assert!(DI::PI.min(DI::empty()).is_empty());
    }

    #[test]
    fn nai() {
        assert!(DI::nai().abs().is_nai());

        assert!(DI::nai().max(DI::PI).is_nai());
        assert!(DI::PI.max(DI::nai()).is_nai());

        assert!(DI::nai().min(DI::PI).is_nai());
        assert!(DI::PI.min(DI::nai()).is_nai());
    }
}