diff --git a/src/interval/CMakeLists.txt b/src/interval/CMakeLists.txt index 9cc261e8f4..94f78a29e7 100644 --- a/src/interval/CMakeLists.txt +++ b/src/interval/CMakeLists.txt @@ -1 +1,2 @@ add_library(interval interval.cpp) +target_link_libraries(interval ${EXTRA_LIBS}) diff --git a/src/interval/interval.cpp b/src/interval/interval.cpp index 51dde8db35..fdaf31a59b 100644 --- a/src/interval/interval.cpp +++ b/src/interval/interval.cpp @@ -8,10 +8,8 @@ Author: Leonardo de Moura #include "mpq.h" namespace lean { - template class interval; template class interval; - } void pp(lean::interval const & i) { std::cout << i << std::endl; } diff --git a/src/interval/interval.h b/src/interval/interval.h index 331719b40c..844bd60519 100644 --- a/src/interval/interval.h +++ b/src/interval/interval.h @@ -27,6 +27,9 @@ class interval { static void round_to_minus_inf() { numeric_traits::set_rounding(false); } static void reset(T & v) { numeric_traits::reset(v); } static void inv(T & v) { numeric_traits::inv(v); } + static void neg(T & v) { numeric_traits::neg(v); } + static bool is_zero(T const & v) { return numeric_traits::is_zero(v); } + static void power(T & v, unsigned k) { return numeric_traits::power(v, k); } void _swap(interval & b); bool _eq(interval const & b) const; @@ -36,19 +39,36 @@ public: interval & operator=(interval const & n); interval & operator=(interval && n); + // (-oo, oo) interval(); + // [n,n] template interval(T2 const & n):m_lower(n), m_upper(n) { set_closed_endpoints();} + // copy constructor interval(interval const & n); - interval(interval && n); + // move constructor + interval(interval && src); + // [l,u], (l,u], [l,u), (l,u) template interval(T2 const & l, T2 const & u, bool l_open = false, bool u_open = false):m_lower(l), m_upper(u) { m_lower_open = l_open; m_upper_open = u_open; m_lower_inf = false; m_upper_inf = false; } + // [l,u], (l,u], [l,u), (l,u) template interval(bool l_open, T2 const & l, T2 const & u, bool u_open):interval(l, u, l_open, u_open) {} + // (-oo, u], (-oo, u] + template interval(T2 const & u, bool u_open):m_upper(u) { + m_lower_open = true; m_upper_open = u_open; m_lower_inf = true; m_upper_inf = false; + } + // [u, +oo), (u, +oo) + template interval(bool l_open, T2 const & l):m_lower(l) { + m_lower_open = l_open; m_upper_open = true; m_lower_inf = false; m_upper_inf = true; + } ~interval(); friend void swap(interval & a, interval & b) { a._swap(b); } + T const & lower() const { return m_lower; } + T const & upper() const { return m_upper; } + bool is_lower_neg() const { return ::lean::is_neg(m_lower, lower_kind()); } bool is_lower_pos() const { return ::lean::is_pos(m_lower, lower_kind()); } bool is_lower_zero() const { return ::lean::is_zero(m_lower, lower_kind()); } @@ -62,16 +82,40 @@ public: bool is_lower_inf() const { return m_lower_inf; } bool is_upper_inf() const { return m_upper_inf; } + // all values in the interval are non-negative bool is_P() const { return is_lower_pos() || is_lower_zero(); } + // interval of the form [0, ... bool is_P0() const { return is_lower_zero() && !is_lower_open(); } + // all values in the interval are positive bool is_P1() const { return is_lower_pos() || (is_lower_zero() && is_lower_open()); } + // all values in the interval are non-positive bool is_N() const { return is_upper_neg() || is_upper_zero(); } + // interval of the form ..., 0] bool is_N0() const { return is_upper_zero() && !is_upper_open(); } + // all values in the interval are negative bool is_N1() const { return is_upper_neg() || (is_upper_zero() && is_upper_open()); } + // lower is negative and upper is positive bool is_M() const { return is_lower_neg() && is_upper_pos(); } + // [0,0] bool is_zero() const { return is_lower_zero() && is_upper_zero(); } + // Interval contains the value zero bool contains_zero() const; + /** + \brief Return true iff this contains b. + That is, every value in b is a value of this. + */ + bool contains(interval & b) const; + + /** + \brief Return true is the interval contains only one value. + */ + bool is_singleton() const; + /** + \brief Return true is the interval contains only v. + */ + template bool is_singleton(T2 const & v) const { return is_singleton() && m_lower == v; } + friend bool operator==(interval const & a, interval const & b) { return a._eq(b); } friend bool operator!=(interval const & a, interval const & b) { return !operator==(a, b); } @@ -88,6 +132,9 @@ public: void set_is_lower_inf(bool v) { m_lower_inf = v; } void set_is_upper_inf(bool v) { m_upper_inf = v; } + void neg(); + friend interval neg(interval o) { o.neg(); return o; } + interval & operator+=(interval const & o); interval & operator-=(interval const & o); interval & operator*=(interval const & o); @@ -96,6 +143,9 @@ public: void inv(); friend interval inv(interval o) { o.inv(); return o; } + void power(unsigned n); + friend interval power(interval o, unsigned k) { o.power(k); return o; } + friend interval operator+(interval a, interval const & b) { return a += b; } friend interval operator-(interval a, interval const & b) { return a -= b; } friend interval operator*(interval a, interval const & b) { return a *= b; } diff --git a/src/interval/interval_def.h b/src/interval/interval_def.h index 2d22cf2d6f..7b4c9de225 100644 --- a/src/interval/interval_def.h +++ b/src/interval/interval_def.h @@ -7,6 +7,7 @@ Author: Leonardo de Moura #pragma once #include #include "interval.h" +#include "trace.h" namespace lean { @@ -44,7 +45,11 @@ template interval::interval(): m_lower(0), m_upper(0) { - set_closed_endpoints(); + m_lower_inf = true; + m_lower_open = true; + m_upper_inf = true; + m_upper_open = true; + lean_assert(check_invariant()); } template @@ -55,6 +60,7 @@ interval::interval(interval const & n): m_upper_open(n.m_upper_open), m_lower_inf(n.m_lower_inf), m_upper_inf(n.m_upper_inf) { + lean_assert(check_invariant()); } template @@ -65,6 +71,7 @@ interval::interval(interval && n): m_upper_open(n.m_upper_open), m_lower_inf(n.m_lower_inf), m_upper_inf(n.m_upper_inf) { + lean_assert(check_invariant()); } template @@ -90,6 +97,32 @@ bool interval::contains_zero() const { (is_upper_pos() || (is_upper_zero() && !is_upper_open())); } +template +bool interval::contains(interval & b) const { + if (!m_lower_inf) { + if (b.m_lower_inf) + return false; + if (m_lower > b.m_lower) + return false; + if (m_lower == b.m_lower && m_lower_open && !b.m_lower_open) + return false; + } + if (!m_upper_inf) { + if (b.m_upper_inf) + return false; + if (m_upper < b.m_upper) + return false; + if (m_upper == b.m_upper && m_upper_open && !b.m_upper_open) + return false; + } + return true; +} + +template +bool interval::is_singleton() const { + return !m_lower_inf && !m_upper_inf && !m_lower_open && !m_upper_open && m_lower == m_upper; +} + template bool interval::_eq(interval const & b) const { return @@ -101,13 +134,61 @@ bool interval::_eq(interval const & b) const { template bool interval::before(interval const & b) const { - // TODO - //if (is_upper_inf() || b.lower_is_inf()) - // return false; - // return m_upper < b.m_lower || (is_upper_open() && - return true; + if (is_upper_inf() || b.is_lower_inf()) + return false; + return m_upper < b.m_lower || (is_upper_open() && m_upper == b.m_lower); } +template +void interval::neg() { + using std::swap; + if (is_lower_inf()) { + if (is_upper_inf()) { + // (-oo, oo) case + // do nothing + } + else { + // (-oo, a| --> |-a, oo) + swap(m_lower, m_upper); + neg(m_lower); + m_lower_inf = false; + m_lower_open = m_upper_open; + + reset(m_upper); + m_upper_inf = true; + m_upper_open = true; + } + } + else { + if (is_upper_inf()) { + // |a, oo) --> (-oo, -a| + swap(m_upper, m_lower); + neg(m_upper); + m_upper_inf = false; + m_upper_open = m_lower_open; + + reset(m_lower); + m_lower_inf = true; + m_lower_open = true; + } + else { + // |a, b| --> |-b, -a| + swap(m_lower, m_upper); + neg(m_lower); + neg(m_upper); + + unsigned lo = m_lower_open; + unsigned li = m_lower_inf; + m_lower_open = m_upper_open; + m_lower_inf = m_upper_inf; + m_upper_open = lo; + m_upper_inf = li; + } + } + lean_assert(check_invariant()); +} + + template interval & interval::operator+=(interval const & o) { xnumeral_kind new_l_kind, new_u_kind; @@ -125,15 +206,23 @@ interval & interval::operator+=(interval const & o) { template interval & interval::operator-=(interval const & o) { + using std::swap; + static thread_local T new_l_val; + static thread_local T new_u_val; xnumeral_kind new_l_kind, new_u_kind; + lean_trace("numerics", tout << "this: " << *this << " o: " << o << "\n";); round_to_minus_inf(); - sub(m_lower, new_l_kind, m_lower, lower_kind(), o.m_upper, o.upper_kind()); + sub(new_l_val, new_l_kind, m_lower, lower_kind(), o.m_upper, o.upper_kind()); round_to_plus_inf(); - add(m_upper, new_u_kind, m_upper, upper_kind(), o.m_lower, o.lower_kind()); + sub(new_u_val, new_u_kind, m_upper, upper_kind(), o.m_lower, o.lower_kind()); + lean_trace("numerics", tout << "new: " << new_l_val << " " << new_u_val << "\n";); + swap(new_l_val, m_lower); + swap(new_u_val, m_upper); m_lower_inf = new_l_kind == XN_MINUS_INFINITY; m_upper_inf = new_u_kind == XN_PLUS_INFINITY; + bool o_l = o.m_lower_open; m_lower_open = m_lower_open || o.m_upper_open; - m_upper_open = m_upper_open || o.m_lower_open; + m_upper_open = m_upper_open || o_l; lean_assert(check_invariant()); return *this; } @@ -321,12 +410,165 @@ interval & interval::operator*=(interval const & o) { set_is_lower_inf(new_l_kind == XN_MINUS_INFINITY); set_is_upper_inf(new_u_kind == XN_PLUS_INFINITY); lean_assert(!(i1_contains_zero || i2_contains_zero) || contains_zero()); + lean_assert(check_invariant()); return *this; } template interval & interval::operator/=(interval const & o) { - // TODO + using std::swap; + interval const & i1 = *this; + interval const & i2 = o; + lean_assert(!i2.contains_zero()); + + if (i1.is_zero()) { + // 0/other = 0 if other != 0 + // do nothing + } + else { + T const & a = i1.m_lower; xnumeral_kind a_k = i1.lower_kind(); + T const & b = i1.m_upper; xnumeral_kind b_k = i1.upper_kind(); + T const & c = i2.m_lower; xnumeral_kind c_k = i2.lower_kind(); + T const & d = i2.m_upper; xnumeral_kind d_k = i2.upper_kind(); + + bool a_o = i1.m_lower_open; + bool b_o = i1.m_upper_open; + bool c_o = i2.m_lower_open; + bool d_o = i2.m_upper_open; + + static thread_local T new_l_val; + static thread_local T new_u_val; + xnumeral_kind new_l_kind, new_u_kind; + + if (i1.is_N()) { + if (i2.is_N1()) { + // x <= b <= 0, c <= y <= d < 0 --> b/c <= x/y + // a <= x <= b <= 0, y <= d < 0 --> x/y <= a/d + set_is_lower_open(i1.is_N0() ? false : b_o || c_o); + set_is_upper_open(a_o || d_o); + + round_to_minus_inf(); + div(new_l_val, new_l_kind, b, b_k, c, c_k); + if (is_zero(d)) { + lean_assert(d_o); + reset(new_u_val); + new_u_kind = XN_PLUS_INFINITY; + } + else { + round_to_plus_inf(); + div(new_u_val, new_u_kind, a, a_k, d, d_k); + } + } + else { + // a <= x, a < 0, 0 < c <= y --> a/c <= x/y + // x <= b <= 0, 0 < c <= y <= d --> x/y <= b/d + lean_assert(i2.is_P1()); + + set_is_upper_open(i1.is_N0() ? false : (b_o || d_o)); + set_is_lower_open(a_o || c_o); + + if (is_zero(c)) { + lean_assert(c_o); + reset(new_l_val); + new_l_kind = XN_MINUS_INFINITY; + } + else { + round_to_minus_inf(); + div(new_l_val, new_l_kind, a, a_k, c, c_k); + } + round_to_plus_inf(); + div(new_u_val, new_u_kind, b, b_k, d, d_k); + } + } + else if (i1.is_M()) { + if (i2.is_N1()) { + // 0 < a <= x <= b < 0, y <= d < 0 --> b/d <= x/y + // 0 < a <= x <= b < 0, y <= d < 0 --> x/y <= a/d + set_is_lower_open(b_o || d_o); + set_is_upper_open(a_o || d_o); + + if (is_zero(d)) { + lean_assert(d_o); + reset(new_l_val); + reset(new_u_val); + new_l_kind = XN_MINUS_INFINITY; + new_u_kind = XN_PLUS_INFINITY; + } + else { + round_to_minus_inf(); + div(new_l_val, new_l_kind, b, b_k, d, d_k); + round_to_plus_inf(); + div(new_u_val, new_u_kind, a, a_k, d, d_k); + } + } + else { + // 0 < a <= x <= b < 0, 0 < c <= y --> a/c <= x/y + // 0 < a <= x <= b < 0, 0 < c <= y --> x/y <= b/c + lean_assert(i1.is_P1()); + + set_is_lower_open(a_o || c_o); + set_is_upper_open(b_o || c_o); + + if (is_zero(c)) { + lean_assert(c_o); + reset(new_l_val); + reset(new_u_val); + new_l_kind = XN_MINUS_INFINITY; + new_u_kind = XN_PLUS_INFINITY; + } + else { + round_to_minus_inf(); + div(new_l_val, new_l_kind, a, a_k, c, c_k); + round_to_plus_inf(); + div(new_u_val, new_u_kind, b, b_k, c, c_k); + } + } + } + else { + lean_assert(i1.is_P()); + if (i2.is_N1()) { + // b > 0, x <= b, c <= y <= d < 0 --> b/d <= x/y + // 0 <= a <= x, c <= y <= d < 0 --> x/y <= a/c + set_is_upper_open(i1.is_P0() ? false : a_o || c_o); + set_is_lower_open(b_o || d_o); + + if (is_zero(d)) { + lean_assert(d_o); + reset(new_l_val); + new_l_kind = XN_MINUS_INFINITY; + } + else { + round_to_minus_inf(); + div(new_l_val, new_l_kind, b, b_k, d, d_k); + } + round_to_plus_inf(); + div(new_u_val, new_u_kind, a, a_k, c, c_k); + } + else { + lean_assert(i2.is_P1()); + // 0 <= a <= x, 0 < c <= y <= d --> a/d <= x/y + // b > 0 x <= b, 0 < c <= y --> x/y <= b/c + set_is_lower_open(i1.is_P0() ? false : a_o || d_o); + set_is_upper_open(b_o || c_o); + + round_to_minus_inf(); + div(new_l_val, new_l_kind, a, a_k, d, d_k); + if (is_zero(c)) { + lean_assert(c_o); + reset(new_u_val); + new_u_kind = XN_PLUS_INFINITY; + } + else { + round_to_plus_inf(); + div(new_u_val, new_u_kind, b, b_k, c, c_k); + } + } + } + swap(m_lower, new_l_val); + swap(m_upper, new_u_val); + m_lower_inf = (new_l_kind == XN_MINUS_INFINITY); + m_upper_inf = (new_u_kind == XN_PLUS_INFINITY); + } return *this; } @@ -404,6 +646,115 @@ void interval::inv() { else { lean_unreachable(); } + lean_assert(check_invariant()); +} + +template +void interval::power(unsigned n) { + using std::swap; + lean_assert(n > 0); + if (n == 1) { + // a^1 = a + // nothing to be done + return; + } + else if (n % 2 == 0) { + if (is_lower_pos()) { + // [l, u]^n = [l^n, u^n] if l > 0 + // 0 < l <= x --> l^n <= x^n (lower bound guarantees that is positive) + // 0 < l <= x <= u --> x^n <= u^n (use lower and upper bound -- need the fact that x is positive) + lean_assert(!is_lower_inf()); + round_to_minus_inf(); + power(m_lower, n); + + if (!m_upper_inf) { + round_to_plus_inf(); + power(m_upper, n); + } + } + else if (is_upper_neg()) { + // [l, u]^n = [u^n, l^n] if u < 0 + // l <= x <= u < 0 --> x^n <= l^n (use lower and upper bound -- need the fact that x is negative) + // x <= u < 0 --> u^n <= x^n + lean_assert(!is_upper_inf()); + bool lo = m_lower_open; + bool li = m_lower_inf; + + swap(m_lower, m_upper); + round_to_minus_inf(); + power(m_lower, n); + m_lower_open = m_upper_open; + m_lower_inf = false; + + if (li) { + reset(m_upper); + } + else { + round_to_plus_inf(); + power(m_upper, n); + } + m_upper_inf = li; + m_upper_open = lo; + } + else { + // [l, u]^n = [0, max{l^n, u^n}] otherwise + // we need both bounds to justify upper bound + xnumeral_kind un1_kind = lower_kind(); + xnumeral_kind un2_kind = upper_kind(); + static thread_local T un1; + static thread_local T un2; + swap(un1, m_lower); + swap(un2, m_upper); + round_to_plus_inf(); + ::lean::power(un1, un1_kind, n); + ::lean::power(un2, un2_kind, n); + + if (gt(un1, un1_kind, un2, un2_kind) || (eq(un1, un1_kind, un2, un2_kind) && !m_lower_open && m_upper_open)) { + swap(m_upper, un1); + m_upper_inf = (un1_kind == XN_PLUS_INFINITY); + m_upper_open = m_lower_open; + } + else { + swap(m_upper, un2); + m_upper_inf = (un2_kind == XN_PLUS_INFINITY); + } + + reset(m_lower); + m_lower_inf = false; + m_lower_open = false; + } + } + else { + // Remark: when n is odd x^n is monotonic. + if (!m_lower_inf) + power(m_lower, n); + if (!m_upper_inf) + power(m_upper, n); + } +} + +/** + return a/(x^n) + + If to_plus_inf, then result >= a/(x^n) + If not to_plus_inf, then result <= a/(x^n) + +*/ +template +T a_div_x_n(T a, T const & x, unsigned n, bool to_plus_inf) { + if (n == 1) { + numeric_traits::set_rounding(to_plus_inf); + a /= x; + } + else { + static thread_local T tmp; + numeric_traits::set_rounding(!to_plus_inf); + tmp = x; + numeric_traits::power(tmp, n); + numeric_traits::set_rounding(to_plus_inf); + a /= x; + } + return a; } template diff --git a/src/tests/interval/interval.cpp b/src/tests/interval/interval.cpp index f9bd3da4e7..d85a21f4b5 100644 --- a/src/tests/interval/interval.cpp +++ b/src/tests/interval/interval.cpp @@ -4,19 +4,91 @@ Released under Apache 2.0 license as described in the file LICENSE. Author: Leonardo de Moura */ -#include "interval_def.h" +#include +#include "interval.h" #include "test.h" +#include "trace.h" #include "mpq.h" using namespace lean; -void tst1() { - interval t(1, 3); - std::cout << t + interval(2, 4, false, true) << "\n"; +typedef interval qi; +typedef std::vector qiv; + +qiv mk_some_intervals(int low, int hi) { + qiv r; + for (unsigned lo = 0; lo < 2; lo++) + for (unsigned uo = 0; uo < 2; uo++) + for (int l = low; l <= hi; l++) + for (int u = l; u <= hi; u++) { + if ((lo || uo) && l == u) + continue; + r.push_back(qi(lo, l, u, uo)); + } + return r; +} + +template bool closed_endpoints(interval const & i) { return !i.is_lower_open() && !i.is_upper_open(); } + +static void tst1() { + qi t(1, 3); + std::cout << t + qi(2, 4, false, true) << "\n"; std::cout << t << " --> " << inv(t) << "\n"; + lean_assert(neg(t) == qi(-3, -1)); + lean_assert(neg(neg(t)) == t); + lean_assert(qi(1, 2) + qi(2, 3) == qi(3, 5)); + lean_assert(qi(1, 5) + qi(-2, -3) == qi(-1, 2)); + for (auto i1 : mk_some_intervals(-2, 2)) + for (auto i2 : mk_some_intervals(-2, 2)) { + auto r = i1 + i2; + auto s = i1; + s += i2; + lean_assert(r == s); + lean_assert(r.lower() == i1.lower() + i2.lower()); + lean_assert(r.upper() == i1.upper() + i2.upper()); + lean_assert(r.is_lower_open() == i1.is_lower_open() || i2.is_lower_open()); + lean_assert(r.is_upper_open() == i1.is_upper_open() || i2.is_upper_open()); + r -= i2; + lean_assert(r.contains(i1)); + r = i1 - i2; + s = i1; + s -= i2; + lean_assert(r == s); + lean_assert(r.lower() == i1.lower() - i2.upper()); + lean_assert(r.upper() == i1.upper() - i2.lower()); + lean_assert(r.is_lower_open() == i1.is_lower_open() || i2.is_upper_open()); + lean_assert(r.is_upper_open() == i1.is_upper_open() || i2.is_lower_open()); + r -= r; + lean_assert(r.contains_zero()); + r = i1 * i2; + s = i1; + s *= i2; + lean_assert(r == s); + lean_assert(r.lower() == std::min(i1.lower()*i2.lower(), std::min(i1.lower()*i2.upper(), std::min(i1.upper()*i2.lower(), i1.upper()*i2.upper())))); + lean_assert(r.upper() == std::max(i1.lower()*i2.lower(), std::max(i1.lower()*i2.upper(), std::max(i1.upper()*i2.lower(), i1.upper()*i2.upper())))); + } + lean_assert(qi(1, 3).before(qi(4, 6))); + lean_assert(!qi(1, 3).before(qi(3, 6))); + lean_assert(qi(1, 3, true, true).before(qi(3, 6))); +} + +static void tst2() { + lean_assert(power(qi(2, 3), 2) == qi(4, 9)); + lean_assert(power(qi(-2, 3), 2) == qi(0, 9)); + lean_assert(power(qi(true, -2, 3, true), 2) == qi(false, 0, 9, true)); + lean_assert(power(qi(true, -4, 3, true), 2) == qi(false, 0, 16, true)); + lean_assert(power(qi(-3, -2), 2) == qi(4, 9)); + std::cout << power(qi(false, -3, -2, true), 2) << " --> " << qi(true, 4, 9, false) << "\n"; + lean_assert(power(qi(false, -3, -2, true), 2) == qi(true, 4, 9, false)); + lean_assert(power(qi(-3,-1), 3) == qi(-27, -1)); + lean_assert(power(qi(-3, 4), 3) == qi(-27, 64)); + lean_assert(power(qi(),3) == qi()); + lean_assert(power(qi(),2) == qi(false, 0)); // (-oo, oo)^2 == [0, oo) } int main() { continue_on_violation(true); + enable_trace("numerics"); tst1(); + tst2(); return 0; }