From 8eb87fbeae276f0d25f8a5f1bdb6ebfc59228f9c Mon Sep 17 00:00:00 2001
From: Soonho Kong <soonhok@cs.cmu.edu>
Date: Tue, 13 Aug 2013 20:05:54 -0700
Subject: [PATCH] Fix interval<T>::operator- and interval<T>::operator/

---
 src/interval/interval.h     |  4 ++--
 src/interval/interval_def.h | 32 ++++++++++++++++++++++++++++++++
 2 files changed, 34 insertions(+), 2 deletions(-)

diff --git a/src/interval/interval.h b/src/interval/interval.h
index dd5f327ecb..58f60c1411 100644
--- a/src/interval/interval.h
+++ b/src/interval/interval.h
@@ -214,9 +214,9 @@ public:
     friend interval<T> operator/(interval<T> a, T const & b) { return a /= b; }
 
     friend interval<T> operator+(T const & a, interval<T> b) { return b += a; }
-    friend interval<T> operator-(T const & a, interval<T> b) { return b += -a; }
+    friend interval<T> operator-(T const & a, interval<T> b) { b.neg(); return b += a; }
     friend interval<T> operator*(T const & a, interval<T> b) { return b *= a; }
-    friend interval<T> operator/(T const & a, interval<T> b) { b = b / a; return b; }
+    friend interval<T> operator/(T const & a, interval<T> b) { b.inv(); return b *= a; }
 
 
     bool check_invariant() const;
diff --git a/src/interval/interval_def.h b/src/interval/interval_def.h
index 047a6e4489..ebfb4e4d05 100644
--- a/src/interval/interval_def.h
+++ b/src/interval/interval_def.h
@@ -625,6 +625,38 @@ interval<T> & interval<T>::operator*=(T const & o) {
 
 template<typename T>
 interval<T> & interval<T>::operator/=(T const & o) {
+    xnumeral_kind new_l_kind, new_u_kind;
+    static thread_local T tmp1;
+    if (this->is_zero()) {
+        return *this;
+    }
+    if (numeric_traits<T>::is_zero(o)) {
+        numeric_traits<T>::reset(m_lower);
+        numeric_traits<T>::reset(m_upper);
+        m_lower_open = m_upper_open = true;
+        m_lower_inf  = m_upper_inf  = true;
+        return *this;
+    }
+
+    if(numeric_traits<T>::is_pos(o)) {
+        // [a, b] / c = [a/c, b/c] when c > 0
+        round_to_minus_inf();
+        div(m_lower, new_l_kind, m_lower, lower_kind(), o, XN_NUMERAL);
+        round_to_plus_inf();
+        div(m_upper, new_u_kind, m_upper, upper_kind(), o, XN_NUMERAL);
+        m_lower_inf = new_l_kind == XN_MINUS_INFINITY;
+        m_upper_inf = new_u_kind == XN_PLUS_INFINITY;
+    }
+    else {
+        // [a, b] / c = [b/c, a/c] when c < 0
+        round_to_minus_inf();
+        div(tmp1, new_l_kind, m_upper, upper_kind(), o, XN_NUMERAL);
+        round_to_plus_inf();
+        div(m_upper, new_u_kind, m_lower, lower_kind(), o, XN_NUMERAL);
+        m_lower = tmp1;
+        m_lower_inf = new_l_kind == XN_MINUS_INFINITY;
+        m_upper_inf = new_u_kind == XN_PLUS_INFINITY;
+    }
     return *this;
 }