diff --git a/src/math/interval/interval.h b/src/math/interval/interval.h index 273d9a626..a4d6def83 100644 --- a/src/math/interval/interval.h +++ b/src/math/interval/interval.h @@ -80,14 +80,14 @@ public: #define DEP_IN_LOWER2 4 #define DEP_IN_UPPER2 8 -typedef short bound_deps; -inline bool dep_in_lower1(bound_deps d) { return (d & DEP_IN_LOWER1) != 0; } -inline bool dep_in_lower2(bound_deps d) { return (d & DEP_IN_LOWER2) != 0; } -inline bool dep_in_upper1(bound_deps d) { return (d & DEP_IN_UPPER1) != 0; } -inline bool dep_in_upper2(bound_deps d) { return (d & DEP_IN_UPPER2) != 0; } -inline bound_deps dep1_to_dep2(bound_deps d) { +typedef short deps_combine_rule; +inline bool dep_in_lower1(deps_combine_rule d) { return (d & DEP_IN_LOWER1) != 0; } +inline bool dep_in_lower2(deps_combine_rule d) { return (d & DEP_IN_LOWER2) != 0; } +inline bool dep_in_upper1(deps_combine_rule d) { return (d & DEP_IN_UPPER1) != 0; } +inline bool dep_in_upper2(deps_combine_rule d) { return (d & DEP_IN_UPPER2) != 0; } +inline deps_combine_rule dep1_to_dep2(deps_combine_rule d) { SASSERT(!dep_in_lower2(d) && !dep_in_upper2(d)); - bound_deps r = d << 2; + deps_combine_rule r = d << 2; SASSERT(dep_in_lower1(d) == dep_in_lower2(r)); SASSERT(dep_in_upper1(d) == dep_in_upper2(r)); SASSERT(!dep_in_lower1(r) && !dep_in_upper1(r)); @@ -99,9 +99,9 @@ inline bound_deps dep1_to_dep2(bound_deps d) { It contains the dependencies for the output lower and upper bounds for the resultant interval. */ -struct interval_deps { - bound_deps m_lower_deps; - bound_deps m_upper_deps; +struct interval_deps_combine_rule { + deps_combine_rule m_lower_combine; + deps_combine_rule m_upper_combine; }; template @@ -239,30 +239,30 @@ public: /** \brief b <- -a */ - void neg(interval const & a, interval & b, interval_deps & b_deps); + void neg(interval const & a, interval & b, interval_deps_combine_rule & b_deps); void neg(interval const & a, interval & b); - void neg_jst(interval const & a, interval_deps & b_deps); + void neg_jst(interval const & a, interval_deps_combine_rule & b_deps); /** \brief c <- a + b */ - void add(interval const & a, interval const & b, interval & c, interval_deps & c_deps); + void add(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps); void add(interval const & a, interval const & b, interval & c); - void add_jst(interval const & a, interval const & b, interval_deps & c_deps); + void add_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps); /** \brief c <- a - b */ - void sub(interval const & a, interval const & b, interval & c, interval_deps & c_deps); + void sub(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps); void sub(interval const & a, interval const & b, interval & c); - void sub_jst(interval const & a, interval const & b, interval_deps & c_deps); + void sub_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps); /** \brief b <- k * a */ - void mul(numeral const & k, interval const & a, interval & b, interval_deps & b_deps); + void mul(numeral const & k, interval const & a, interval & b, interval_deps_combine_rule & b_deps); void mul(numeral const & k, interval const & a, interval & b) { div_mul(k, a, b, false); } - void mul_jst(numeral const & k, interval const & a, interval_deps & b_deps); + void mul_jst(numeral const & k, interval const & a, interval_deps_combine_rule & b_deps); /** \brief b <- (n/d) * a */ @@ -278,58 +278,58 @@ public: That is, we must invert k rounding towards +oo or -oo depending whether we are computing a lower or upper bound. */ - void div(interval const & a, numeral const & k, interval & b, interval_deps & b_deps); + void div(interval const & a, numeral const & k, interval & b, interval_deps_combine_rule & b_deps); void div(interval const & a, numeral const & k, interval & b) { div_mul(k, a, b, true); } - void div_jst(interval const & a, numeral const & k, interval_deps & b_deps) { mul_jst(k, a, b_deps); } + void div_jst(interval const & a, numeral const & k, interval_deps_combine_rule & b_deps) { mul_jst(k, a, b_deps); } /** \brief c <- a * b */ - void mul(interval const & a, interval const & b, interval & c, interval_deps & c_deps); + void mul(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps); void mul(interval const & a, interval const & b, interval & c); - void mul_jst(interval const & a, interval const & b, interval_deps & c_deps); + void mul_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps); /** \brief b <- a^n */ - void power(interval const & a, unsigned n, interval & b, interval_deps & b_deps); + void power(interval const & a, unsigned n, interval & b, interval_deps_combine_rule & b_deps); void power(interval const & a, unsigned n, interval & b); - void power_jst(interval const & a, unsigned n, interval_deps & b_deps); + void power_jst(interval const & a, unsigned n, interval_deps_combine_rule & b_deps); /** \brief b <- a^(1/n) with precision p. \pre if n is even, then a must not contain negative numbers. */ - void nth_root(interval const & a, unsigned n, numeral const & p, interval & b, interval_deps & b_deps); + void nth_root(interval const & a, unsigned n, numeral const & p, interval & b, interval_deps_combine_rule & b_deps); void nth_root(interval const & a, unsigned n, numeral const & p, interval & b); - void nth_root_jst(interval const & a, unsigned n, numeral const & p, interval_deps & b_deps); + void nth_root_jst(interval const & a, unsigned n, numeral const & p, interval_deps_combine_rule & b_deps); /** \brief Given an equation x^n = y and an interval for y, compute the solution set for x with precision p. \pre if n is even, then !lower_is_neg(y) */ - void xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x, interval_deps & x_deps); + void xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x, interval_deps_combine_rule & x_deps); void xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x); - void xn_eq_y_jst(interval const & y, unsigned n, numeral const & p, interval_deps & x_deps); + void xn_eq_y_jst(interval const & y, unsigned n, numeral const & p, interval_deps_combine_rule & x_deps); /** \brief b <- 1/a \pre !contains_zero(a) */ - void inv(interval const & a, interval & b, interval_deps & b_deps); + void inv(interval const & a, interval & b, interval_deps_combine_rule & b_deps); void inv(interval const & a, interval & b); - void inv_jst(interval const & a, interval_deps & b_deps); + void inv_jst(interval const & a, interval_deps_combine_rule & b_deps); /** \brief c <- a/b \pre !contains_zero(b) \pre &a == &c (that is, c should not be an alias for a) */ - void div(interval const & a, interval const & b, interval & c, interval_deps & c_deps); + void div(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps); void div(interval const & a, interval const & b, interval & c); - void div_jst(interval const & a, interval const & b, interval_deps & c_deps); + void div_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps); /** \brief Store in r an interval that contains the number pi. diff --git a/src/math/interval/interval_def.h b/src/math/interval/interval_def.h index 1b9bb473e..07fd242cd 100644 --- a/src/math/interval/interval_def.h +++ b/src/math/interval/interval_def.h @@ -713,31 +713,31 @@ bool interval_manager::before(interval const & a, interval const & b) const { } template -void interval_manager::neg_jst(interval const & a, interval_deps & b_deps) { +void interval_manager::neg_jst(interval const & a, interval_deps_combine_rule & b_deps) { if (lower_is_inf(a)) { if (upper_is_inf(a)) { - b_deps.m_lower_deps = 0; - b_deps.m_upper_deps = 0; + b_deps.m_lower_combine = 0; + b_deps.m_upper_combine = 0; } else { - b_deps.m_lower_deps = DEP_IN_UPPER1; - b_deps.m_upper_deps = 0; + b_deps.m_lower_combine = DEP_IN_UPPER1; + b_deps.m_upper_combine = 0; } } else { if (upper_is_inf(a)) { - b_deps.m_lower_deps = 0; - b_deps.m_upper_deps = DEP_IN_LOWER1; + b_deps.m_lower_combine = 0; + b_deps.m_upper_combine = DEP_IN_LOWER1; } else { - b_deps.m_lower_deps = DEP_IN_UPPER1; - b_deps.m_upper_deps = DEP_IN_LOWER1; + b_deps.m_lower_combine = DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_LOWER1; } } } template -void interval_manager::neg(interval const & a, interval & b, interval_deps & b_deps) { +void interval_manager::neg(interval const & a, interval & b, interval_deps_combine_rule & b_deps) { neg_jst(a, b_deps); neg(a, b); } @@ -792,13 +792,13 @@ void interval_manager::neg(interval const & a, interval & b) { } template -void interval_manager::add_jst(interval const & a, interval const & b, interval_deps & c_deps) { - c_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2; - c_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2; +void interval_manager::add_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps) { + c_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2; + c_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2; } template -void interval_manager::add(interval const & a, interval const & b, interval & c, interval_deps & c_deps) { +void interval_manager::add(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps) { add_jst(a, b, c_deps); add(a, b, c); } @@ -818,13 +818,13 @@ void interval_manager::add(interval const & a, interval const & b, interval & } template -void interval_manager::sub_jst(interval const & a, interval const & b, interval_deps & c_deps) { - c_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2; - c_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2; +void interval_manager::sub_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps) { + c_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2; + c_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2; } template -void interval_manager::sub(interval const & a, interval const & b, interval & c, interval_deps & c_deps) { +void interval_manager::sub(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps) { sub_jst(a, b, c_deps); sub(a, b, c); } @@ -844,18 +844,18 @@ void interval_manager::sub(interval const & a, interval const & b, interval & } template -void interval_manager::mul_jst(numeral const & k, interval const & a, interval_deps & b_deps) { +void interval_manager::mul_jst(numeral const & k, interval const & a, interval_deps_combine_rule & b_deps) { if (m().is_zero(k)) { - b_deps.m_lower_deps = 0; - b_deps.m_upper_deps = 0; + b_deps.m_lower_combine = 0; + b_deps.m_upper_combine = 0; } else if (m().is_neg(k)) { - b_deps.m_lower_deps = DEP_IN_UPPER1; - b_deps.m_upper_deps = DEP_IN_LOWER1; + b_deps.m_lower_combine = DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_LOWER1; } else { - b_deps.m_lower_deps = DEP_IN_LOWER1; - b_deps.m_upper_deps = DEP_IN_UPPER1; + b_deps.m_lower_combine = DEP_IN_LOWER1; + b_deps.m_upper_combine = DEP_IN_UPPER1; } } @@ -918,7 +918,7 @@ void interval_manager::div_mul(numeral const & k, interval const & a, interva } template -void interval_manager::mul(numeral const & k, interval const & a, interval & b, interval_deps & b_deps) { +void interval_manager::mul(numeral const & k, interval const & a, interval & b, interval_deps_combine_rule & b_deps) { mul_jst(k, a, b_deps); mul(k, a, b); } @@ -931,57 +931,57 @@ void interval_manager::mul(int n, int d, interval const & a, interval & b) { } template -void interval_manager::div(interval const & a, numeral const & k, interval & b, interval_deps & b_deps) { +void interval_manager::div(interval const & a, numeral const & k, interval & b, interval_deps_combine_rule & b_deps) { div_jst(a, k, b_deps); div(a, k, b); } template -void interval_manager::mul_jst(interval const & i1, interval const & i2, interval_deps & r_deps) { +void interval_manager::mul_jst(interval const & i1, interval const & i2, interval_deps_combine_rule & r_deps) { if (is_zero(i1)) { - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; } else if (is_zero(i2)) { - r_deps.m_lower_deps = DEP_IN_LOWER2 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_LOWER2 | DEP_IN_UPPER2; } else if (is_N(i1)) { if (is_N(i2)) { // x <= b <= 0, y <= d <= 0 --> b*d <= x*y // a <= x <= b <= 0, c <= y <= d <= 0 --> x*y <= a*c (we can use the fact that x or y is always negative (i.e., b is neg or d is neg)) - r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER1; // we can replace DEP_IN_UPPER1 with DEP_IN_UPPER2 + r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER1; // we can replace DEP_IN_UPPER1 with DEP_IN_UPPER2 } else if (is_M(i2)) { // a <= x <= b <= 0, y <= d, d > 0 --> a*d <= x*y (uses the fact that b is not positive) // a <= x <= b <= 0, c <= y, c < 0 --> x*y <= a*c (uses the fact that b is not positive) - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_UPPER1; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER1; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_UPPER1; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER1; } else { // a <= x <= b <= 0, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that x is neg (b is not positive) or y is pos (c is not negative)) // x <= b <= 0, 0 <= c <= y --> x*y <= b*c - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_UPPER1; // we can replace DEP_IN_UPPER1 with DEP_IN_UPPER2 - r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_UPPER1; // we can replace DEP_IN_UPPER1 with DEP_IN_UPPER2 + r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2; } } else if (is_M(i1)) { if (is_N(i2)) { // b > 0, x <= b, c <= y <= d <= 0 --> b*c <= x*y (uses the fact that d is not positive) // a < 0, a <= x, c <= y <= d <= 0 --> x*y <= a*c (uses the fact that d is not positive) - r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; } else if (is_M(i2)) { - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; } else { // a < 0, a <= x, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that c is not negative) // b > 0, x <= b, 0 <= c <= y <= d --> x*y <= b*d (uses the fact that c is not negative) - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_LOWER2; - r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER2; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_LOWER2; + r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER2; } } else { @@ -989,27 +989,27 @@ void interval_manager::mul_jst(interval const & i1, interval const & i2, inte if (is_N(i2)) { // 0 <= a <= x <= b, c <= y <= d <= 0 --> x*y <= b*c (uses the fact that x is pos (a is not neg) or y is neg (d is not pos)) // 0 <= a <= x, y <= d <= 0 --> a*d <= x*y - r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_UPPER2 - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_UPPER2 + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2; } else if (is_M(i2)) { // 0 <= a <= x <= b, c <= y --> b*c <= x*y (uses the fact that a is not negative) // 0 <= a <= x <= b, y <= d --> x*y <= b*d (uses the fact that a is not negative) - r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_LOWER1; - r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1; + r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_LOWER1; + r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1; } else { SASSERT(is_P(i2)); // 0 <= a <= x, 0 <= c <= y --> a*c <= x*y // x <= b, y <= d --> x*y <= b*d (uses the fact that x is pos (a is not negative) or y is pos (c is not negative)) - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2; - r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_LOWER2 + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2; + r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_LOWER2 } } } template -void interval_manager::mul(interval const & i1, interval const & i2, interval & r, interval_deps & r_deps) { +void interval_manager::mul(interval const & i1, interval const & i2, interval & r, interval_deps_combine_rule & r_deps) { mul_jst(i1, i2, r_deps); mul(i1, i2, r); } @@ -1223,55 +1223,55 @@ void interval_manager::mul(interval const & i1, interval const & i2, interval } template -void interval_manager::power_jst(interval const & a, unsigned n, interval_deps & b_deps) { +void interval_manager::power_jst(interval const & a, unsigned n, interval_deps_combine_rule & b_deps) { if (n == 1) { - b_deps.m_lower_deps = DEP_IN_LOWER1; - b_deps.m_upper_deps = DEP_IN_UPPER1; + b_deps.m_lower_combine = DEP_IN_LOWER1; + b_deps.m_upper_combine = DEP_IN_UPPER1; } else if (n % 2 == 0) { if (lower_is_pos(a)) { // [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) - b_deps.m_lower_deps = DEP_IN_LOWER1; + b_deps.m_lower_combine = DEP_IN_LOWER1; if (upper_is_inf(a)) - b_deps.m_upper_deps = 0; + b_deps.m_upper_combine = 0; else - b_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; } else if (upper_is_neg(a)) { // [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 - b_deps.m_lower_deps = DEP_IN_UPPER1; + b_deps.m_lower_combine = DEP_IN_UPPER1; if (lower_is_inf(a)) - b_deps.m_upper_deps = 0; + b_deps.m_upper_combine = 0; else - b_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; } else { // [l, u]^n = [0, max{l^n, u^n}] otherwise // we need both bounds to justify upper bound - b_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; - b_deps.m_lower_deps = 0; + b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; + b_deps.m_lower_combine = 0; } } else { // Remark: when n is odd x^n is monotonic. if (lower_is_inf(a)) - b_deps.m_lower_deps = 0; + b_deps.m_lower_combine = 0; else - b_deps.m_lower_deps = DEP_IN_LOWER1; + b_deps.m_lower_combine = DEP_IN_LOWER1; if (upper_is_inf(a)) - b_deps.m_upper_deps = 0; + b_deps.m_upper_combine = 0; else - b_deps.m_upper_deps = DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_UPPER1; } } template -void interval_manager::power(interval const & a, unsigned n, interval & b, interval_deps & b_deps) { +void interval_manager::power(interval const & a, unsigned n, interval & b, interval_deps_combine_rule & b_deps) { power_jst(a, n, b_deps); power(a, n, b); } @@ -1392,7 +1392,7 @@ void interval_manager::power(interval const & a, unsigned n, interval & b) { template -void interval_manager::nth_root(interval const & a, unsigned n, numeral const & p, interval & b, interval_deps & b_deps) { +void interval_manager::nth_root(interval const & a, unsigned n, numeral const & p, interval & b, interval_deps_combine_rule & b_deps) { nth_root_jst(a, n, p, b_deps); nth_root(a, n, p, b); } @@ -1437,16 +1437,16 @@ void interval_manager::nth_root(interval const & a, unsigned n, numeral const } template -void interval_manager::nth_root_jst(interval const & a, unsigned n, numeral const & p, interval_deps & b_deps) { - b_deps.m_lower_deps = DEP_IN_LOWER1; +void interval_manager::nth_root_jst(interval const & a, unsigned n, numeral const & p, interval_deps_combine_rule & b_deps) { + b_deps.m_lower_combine = DEP_IN_LOWER1; if (n % 2 == 0) - b_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; else - b_deps.m_upper_deps = DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_UPPER1; } template -void interval_manager::xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x, interval_deps & x_deps) { +void interval_manager::xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x, interval_deps_combine_rule & x_deps) { xn_eq_y_jst(y, n, p, x_deps); xn_eq_y(y, n, p, x); } @@ -1486,29 +1486,29 @@ void interval_manager::xn_eq_y(interval const & y, unsigned n, numeral const } template -void interval_manager::xn_eq_y_jst(interval const & y, unsigned n, numeral const & p, interval_deps & x_deps) { +void interval_manager::xn_eq_y_jst(interval const & y, unsigned n, numeral const & p, interval_deps_combine_rule & x_deps) { if (n % 2 == 0) { - x_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; - x_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; + x_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; + x_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; } else { - x_deps.m_lower_deps = DEP_IN_LOWER1; - x_deps.m_upper_deps = DEP_IN_UPPER1; + x_deps.m_lower_combine = DEP_IN_LOWER1; + x_deps.m_upper_combine = DEP_IN_UPPER1; } } template -void interval_manager::inv_jst(interval const & a, interval_deps & b_deps) { +void interval_manager::inv_jst(interval const & a, interval_deps_combine_rule & b_deps) { SASSERT(!contains_zero(a)); if (is_P1(a)) { - b_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; - b_deps.m_upper_deps = DEP_IN_LOWER1; + b_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_LOWER1; } else if (is_N1(a)) { // x <= u < 0 --> 1/u <= 1/x // l <= x <= u < 0 --> 1/l <= 1/x (use lower and upper bounds) - b_deps.m_lower_deps = DEP_IN_UPPER1; - b_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER1; + b_deps.m_lower_combine = DEP_IN_UPPER1; + b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1; } else { UNREACHABLE(); @@ -1516,7 +1516,7 @@ void interval_manager::inv_jst(interval const & a, interval_deps & b_deps) { } template -void interval_manager::inv(interval const & a, interval & b, interval_deps & b_deps) { +void interval_manager::inv(interval const & a, interval & b, interval_deps_combine_rule & b_deps) { inv_jst(a, b_deps); inv(a, b); } @@ -1602,16 +1602,16 @@ void interval_manager::inv(interval const & a, interval & b) { } template -void interval_manager::div_jst(interval const & i1, interval const & i2, interval_deps & r_deps) { +void interval_manager::div_jst(interval const & i1, interval const & i2, interval_deps_combine_rule & r_deps) { SASSERT(!contains_zero(i2)); if (is_zero(i1)) { if (is_P1(i2)) { - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2; - r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2; + r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2; } else { - r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2; } } else { @@ -1619,28 +1619,28 @@ void interval_manager::div_jst(interval const & i1, interval const & i2, inte if (is_N1(i2)) { // x <= b <= 0, c <= y <= d < 0 --> b/c <= x/y // a <= x <= b <= 0, y <= d < 0 --> x/y <= a/d - r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2; } else { // a <= x, a < 0, 0 < c <= y --> a/c <= x/y // x <= b <= 0, 0 < c <= y <= d --> x/y <= b/d - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2; - r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2; + r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; } } else if (is_M(i1)) { if (is_N1(i2)) { // 0 < a <= x <= b < 0, y <= d < 0 --> b/d <= x/y // 0 < a <= x <= b < 0, y <= d < 0 --> x/y <= a/d - r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2; } 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 - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2; - r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2; + r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2; } } else { @@ -1648,22 +1648,22 @@ void interval_manager::div_jst(interval const & i1, interval const & i2, inte if (is_N1(i2)) { // b > 0, x <= b, c <= y <= d < 0 --> b/d <= x/y // 0 <= a <= x, c <= y <= d < 0 --> x/y <= a/c - r_deps.m_lower_deps = DEP_IN_UPPER1 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; } else { SASSERT(is_P1(i2)); // 0 <= a <= x, 0 < c <= y <= d --> a/d <= x/y // b > 0 x <= b, 0 < c <= y --> x/y <= b/c - r_deps.m_lower_deps = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; - r_deps.m_upper_deps = DEP_IN_UPPER1 | DEP_IN_LOWER2; + r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2; + r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2; } } } } template -void interval_manager::div(interval const & i1, interval const & i2, interval & r, interval_deps & r_deps) { +void interval_manager::div(interval const & i1, interval const & i2, interval & r, interval_deps_combine_rule & r_deps) { div_jst(i1, i2, r_deps); div(i1, i2, r); } diff --git a/src/math/lp/horner.cpp b/src/math/lp/horner.cpp index e07ddcde3..fce8086a0 100644 --- a/src/math/lp/horner.cpp +++ b/src/math/lp/horner.cpp @@ -145,29 +145,31 @@ interv horner::interval_of_expr(const nex& e) { interv horner::interval_of_mul(const nex& e) { SASSERT(e.is_mul()); auto & es = e.children(); - SASSERT(es.size()); - TRACE("nla_horner_details", tout << "e=" << e << "\n";); + interval_deps_combine_rule comb_rule; interv a = interval_of_expr(es[0]); + if (m_intervals.is_zero(a)) { + m_intervals.set_zero_interval_deps_for_mult(a); + TRACE("nla_horner_details", tout << "es[0]= "<< es[0] << std::endl << "a = "; m_intervals.display(tout, a); ); + return a; + } TRACE("nla_horner_details", tout << "es[0]= "<< es[0] << std::endl << "a = "; m_intervals.display(tout, a); ); for (unsigned k = 1; k < es.size(); k++) { - interval_deps deps; - if (m_intervals.is_zero(a)) { - interv t; - set_interval_for_scalar(t, rational(1)); - m_intervals.mul(a, t, a, deps); - m_intervals.add_deps(a, t, deps, a); - TRACE("nla_horner_details", tout << "got zero\n"; ); - break; - } interv b = interval_of_expr(es[k]); + if (m_intervals.is_zero(b)) { + m_intervals.set_zero_interval_deps_for_mult(b); + TRACE("nla_horner_details", tout << "es[k]= "<< es[k] << std::endl << ", "; m_intervals.display(tout, b); ); + TRACE("nla_horner_details", tout << "got zero\n"; ); + return b; + } + TRACE("nla_horner_details", tout << "es[k]= "<< es[k] << std::endl << ", "; m_intervals.display(tout, b); ); interv c; - m_intervals.mul(a, b, c, deps); - m_intervals.add_deps(a, b, deps, a); + m_intervals.mul(a, b, c, comb_rule); + m_intervals.combine_deps(a, b, comb_rule, a); m_intervals.set(a, c); - TRACE("nla_horner_details", tout << "es[" << k << "]=" << es[k] << ", "; m_intervals.display(tout, a) << "\n";); + TRACE("nla_horner_details", tout << "part mult "; m_intervals.display(tout, a) << "\n";); } TRACE("nla_horner_details", tout << "e=" << e << "\n"; - tout << " interv = "; m_intervals.display(tout, a);); + tout << " return "; m_intervals.display(tout, a);); return a; } @@ -190,10 +192,10 @@ interv horner::interval_of_sum(const nex& e) { return b; } interv c; - interval_deps deps; + interval_deps_combine_rule combine_rule; TRACE("nla_horner_details_sum", tout << "a = "; m_intervals.display(tout, a) << "\nb = "; m_intervals.display(tout, b) << "\n";); - m_intervals.add(a, b, c, deps); - m_intervals.add_deps(a, b, deps, a); + m_intervals.add(a, b, c, combine_rule); + m_intervals.combine_deps(a, b, combine_rule, a); m_intervals.set(a, c); TRACE("nla_horner_details_sum", tout << es[k] << ", "; m_intervals.display(tout, a); tout << "\n";); @@ -209,7 +211,7 @@ interv horner::interval_of_sum(const nex& e) { // sets the dependencies also void horner::set_var_interval(lpvar v, interv& b) { m_intervals.set_var_interval_with_deps(v, b); - TRACE("nla_horner_details", tout << "v = "; print_var(v, tout) << "\n"; m_intervals.display(tout, b);); + TRACE("nla_horner_details_var", tout << "v = "; print_var(v, tout) << "\n"; m_intervals.display(tout, b);); } } diff --git a/src/math/lp/nla_intervals.cpp b/src/math/lp/nla_intervals.cpp index e2baa8f44..e63b00bbe 100644 --- a/src/math/lp/nla_intervals.cpp +++ b/src/math/lp/nla_intervals.cpp @@ -32,6 +32,11 @@ void intervals::set_var_interval_with_deps(lpvar v, interval& b) { } } +void intervals::set_zero_interval_deps_for_mult(interval& a) { + a.m_lower_dep = m_dep_manager.mk_join(a.m_lower_dep, a.m_upper_dep); + a.m_upper_dep = a.m_lower_dep; +} + bool intervals::check_interval_for_conflict_on_zero(const interval & i) { return check_interval_for_conflict_on_zero_lower(i) || check_interval_for_conflict_on_zero_upper(i); } diff --git a/src/math/lp/nla_intervals.h b/src/math/lp/nla_intervals.h index 2c2865bd7..348fc8abc 100644 --- a/src/math/lp/nla_intervals.h +++ b/src/math/lp/nla_intervals.h @@ -74,9 +74,9 @@ class intervals : common { }; - void add_deps(interval const& a, interval const& b, interval_deps const& deps, interval& i) const { - i.m_lower_dep = mk_dependency(a, b, deps.m_lower_deps); - i.m_upper_dep = mk_dependency(a, b, deps.m_upper_deps); + void add_deps(interval const& a, interval const& b, interval_deps_combine_rule const& deps, interval& i) const { + i.m_lower_dep = mk_dependency(a, b, deps.m_lower_combine); + i.m_upper_dep = mk_dependency(a, b, deps.m_upper_combine); } @@ -118,7 +118,7 @@ class intervals : common { im_config(numeral_manager & m, ci_dependency_manager& d):m_manager(m), m_dep_manager(d) {} private: - ci_dependency* mk_dependency(interval const& a, interval const& b, bound_deps bd) const { + ci_dependency* mk_dependency(interval const& a, interval const& b, deps_combine_rule bd) const { ci_dependency* dep = nullptr; if (dep_in_lower1(bd)) { dep = m_dep_manager.mk_join(dep, a.m_lower_dep); @@ -173,16 +173,19 @@ public: bool is_zero(const interval& a) const { return m_config.is_zero(a); } void mul(const interval& a, const interval& b, interval& c) { m_imanager.mul(a, b, c); } void add(const interval& a, const interval& b, interval& c) { m_imanager.add(a, b, c); } - void add(const interval& a, const interval& b, interval& c, interval_deps& deps) { m_imanager.add(a, b, c, deps); } + void add(const interval& a, const interval& b, interval& c, interval_deps_combine_rule& deps) { m_imanager.add(a, b, c, deps); } void set(interval& a, const interval& b) { m_imanager.set(a, b); } - void mul(const interval& a, const interval& b, interval& c, interval_deps& deps) { m_imanager.mul(a, b, c, deps); } - void add_deps(interval const& a, interval const& b, interval_deps const& deps, interval& i) const { + void mul(const interval& a, const interval& b, interval& c, interval_deps_combine_rule& deps) { m_imanager.mul(a, b, c, deps); } + void combine_deps(interval const& a, interval const& b, interval_deps_combine_rule const& deps, interval& i) const { m_config.add_deps(a, b, deps, i); } + + bool upper_is_inf(const interval& a) const { return m_config.upper_is_inf(a); } bool lower_is_inf(const interval& a) const { return m_config.lower_is_inf(a); } void set_var_interval_with_deps(lpvar, interval &); + void set_zero_interval_deps_for_mult(interval&); bool is_inf(const interval& i) const { return m_config.is_inf(i); } bool check_interval_for_conflict_on_zero(const interval & i); bool check_interval_for_conflict_on_zero_lower(const interval & i);