00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifndef PPL_Box_templates_hh
00024 #define PPL_Box_templates_hh 1
00025
00026 #include "Variables_Set.defs.hh"
00027 #include "Constraint_System.defs.hh"
00028 #include "Constraint_System.inlines.hh"
00029 #include "Generator_System.defs.hh"
00030 #include "Generator_System.inlines.hh"
00031 #include "Poly_Con_Relation.defs.hh"
00032 #include "Poly_Gen_Relation.defs.hh"
00033 #include "Polyhedron.defs.hh"
00034 #include "Grid.defs.hh"
00035 #include "BD_Shape.defs.hh"
00036 #include "Octagonal_Shape.defs.hh"
00037 #include "MIP_Problem.defs.hh"
00038 #include "Rational_Interval.hh"
00039 #include <iostream>
00040
00041 namespace Parma_Polyhedra_Library {
00042
00043 template <typename ITV>
00044 inline
00045 Box<ITV>::Box(dimension_type num_dimensions, Degenerate_Element kind)
00046 : seq(num_dimensions <= max_space_dimension()
00047 ? num_dimensions
00048 : (throw_space_dimension_overflow("Box(n, k)",
00049 "n exceeds the maximum "
00050 "allowed space dimension"),
00051 num_dimensions)),
00052 status() {
00053
00054
00055 if (kind == UNIVERSE) {
00056 for (dimension_type i = num_dimensions; i-- > 0; )
00057 seq[i].assign(UNIVERSE);
00058 set_empty_up_to_date();
00059 }
00060 else
00061 set_empty();
00062 assert(OK());
00063 }
00064
00065 template <typename ITV>
00066 inline
00067 Box<ITV>::Box(const Constraint_System& cs)
00068 : seq(cs.space_dimension() <= max_space_dimension()
00069 ? cs.space_dimension()
00070 : (throw_space_dimension_overflow("Box(cs)",
00071 "cs exceeds the maximum "
00072 "allowed space dimension"),
00073 cs.space_dimension())),
00074 status() {
00075
00076 for (dimension_type i = cs.space_dimension(); i-- > 0; )
00077 seq[i].assign(UNIVERSE);
00078 add_constraints_no_check(cs);
00079 }
00080
00081 template <typename ITV>
00082 inline
00083 Box<ITV>::Box(const Congruence_System& cgs)
00084 : seq(cgs.space_dimension() <= max_space_dimension()
00085 ? cgs.space_dimension()
00086 : (throw_space_dimension_overflow("Box(cgs)",
00087 "cgs exceeds the maximum "
00088 "allowed space dimension"),
00089 cgs.space_dimension())),
00090 status() {
00091
00092 for (dimension_type i = cgs.space_dimension(); i-- > 0; )
00093 seq[i].assign(UNIVERSE);
00094 add_congruences_no_check(cgs);
00095 }
00096
00097 template <typename ITV>
00098 template <typename Other_ITV>
00099 inline
00100 Box<ITV>::Box(const Box<Other_ITV>& y, Complexity_Class)
00101 : seq(y.space_dimension()),
00102
00103
00104 status() {
00105
00106 if (y.marked_empty())
00107 set_empty();
00108
00109 if (!y.marked_empty())
00110 for (dimension_type k = y.space_dimension(); k-- > 0; )
00111 seq[k].assign(y.seq[k]);
00112 assert(OK());
00113 }
00114
00115 template <typename ITV>
00116 Box<ITV>::Box(const Generator_System& gs)
00117 : seq(gs.space_dimension() <= max_space_dimension()
00118 ? gs.space_dimension()
00119 : (throw_space_dimension_overflow("Box(gs)",
00120 "gs exceeds the maximum "
00121 "allowed space dimension"),
00122 gs.space_dimension())),
00123 status() {
00124 const Generator_System::const_iterator gs_begin = gs.begin();
00125 const Generator_System::const_iterator gs_end = gs.end();
00126 if (gs_begin == gs_end) {
00127
00128 set_empty();
00129 return;
00130 }
00131
00132
00133 set_empty_up_to_date();
00134
00135 const dimension_type space_dim = space_dimension();
00136 DIRTY_TEMP0(mpq_class, q);
00137 bool point_seen = false;
00138
00139 for (Generator_System::const_iterator
00140 gs_i = gs_begin; gs_i != gs_end; ++gs_i) {
00141 const Generator& g = *gs_i;
00142 if (g.is_point()) {
00143 const Coefficient& d = g.divisor();
00144 if (point_seen) {
00145
00146 for (dimension_type i = space_dim; i-- > 0; ) {
00147 assign_r(q.get_num(), g.coefficient(Variable(i)), ROUND_NOT_NEEDED);
00148 assign_r(q.get_den(), d, ROUND_NOT_NEEDED);
00149 q.canonicalize();
00150 seq[i].join_assign(q);
00151 }
00152 }
00153 else {
00154
00155 point_seen = true;
00156 for (dimension_type i = space_dim; i-- > 0; ) {
00157 assign_r(q.get_num(), g.coefficient(Variable(i)), ROUND_NOT_NEEDED);
00158 assign_r(q.get_den(), d, ROUND_NOT_NEEDED);
00159 q.canonicalize();
00160 seq[i].assign(q);
00161 }
00162 }
00163 }
00164 }
00165
00166 if (!point_seen)
00167
00168 throw std::invalid_argument("PPL::Box<ITV>::Box(gs):\n"
00169 "the non-empty generator system gs "
00170 "contains no points.");
00171
00172
00173 ITV q_interval;
00174 for (Generator_System::const_iterator gs_i = gs_begin;
00175 gs_i != gs_end; ++gs_i) {
00176 const Generator& g = *gs_i;
00177 switch (g.type()) {
00178 case Generator::LINE:
00179 for (dimension_type i = space_dim; i-- > 0; )
00180 if (g.coefficient(Variable(i)) != 0)
00181 seq[i].assign(UNIVERSE);
00182 break;
00183 case Generator::RAY:
00184 for (dimension_type i = space_dim; i-- > 0; )
00185 switch (sgn(g.coefficient(Variable(i)))) {
00186 case 1:
00187 seq[i].upper_set(UNBOUNDED);
00188 break;
00189 case -1:
00190 seq[i].lower_set(UNBOUNDED);
00191 break;
00192 default:
00193 break;
00194 }
00195 break;
00196 case Generator::CLOSURE_POINT:
00197 {
00198 const Coefficient& d = g.divisor();
00199 for (dimension_type i = space_dim; i-- > 0; ) {
00200 assign_r(q.get_num(), g.coefficient(Variable(i)), ROUND_NOT_NEEDED);
00201 assign_r(q.get_den(), d, ROUND_NOT_NEEDED);
00202 q.canonicalize();
00203 ITV& seq_i = seq[i];
00204 seq_i.lower_widen(q, true);
00205 seq_i.upper_widen(q, true);
00206 }
00207 }
00208 break;
00209 default:
00210
00211 break;
00212 }
00213 }
00214 assert(OK());
00215 }
00216
00217 template <typename ITV>
00218 template <typename T>
00219 Box<ITV>::Box(const BD_Shape<T>& bds, Complexity_Class)
00220 : seq(bds.space_dimension() <= max_space_dimension()
00221 ? bds.space_dimension()
00222 : (throw_space_dimension_overflow("Box(bds)",
00223 "bds exceeds the maximum "
00224 "allowed space dimension"),
00225 bds.space_dimension())),
00226 status() {
00227
00228 bds.shortest_path_closure_assign();
00229 if (bds.marked_empty()) {
00230 set_empty();
00231 assert(OK());
00232 return;
00233 }
00234
00235
00236 set_empty_up_to_date();
00237
00238 const dimension_type space_dim = space_dimension();
00239 if (space_dim == 0) {
00240 assert(OK());
00241 return;
00242 }
00243
00244 DIRTY_TEMP(typename BD_Shape<T>::coefficient_type, tmp);
00245 const DB_Row<typename BD_Shape<T>::coefficient_type>& dbm_0 = bds.dbm[0];
00246 for (dimension_type i = space_dim; i-- > 0; ) {
00247 ITV& seq_i = seq[i];
00248
00249 const typename BD_Shape<T>::coefficient_type& u = dbm_0[i+1];
00250 if (is_plus_infinity(u))
00251 seq_i.upper_set_uninit(UNBOUNDED);
00252 else
00253 seq_i.upper_set_uninit(u);
00254
00255
00256 const typename BD_Shape<T>::coefficient_type& negated_l = bds.dbm[i+1][0];
00257 if (is_plus_infinity(negated_l))
00258 seq_i.lower_set_uninit(UNBOUNDED);
00259 else {
00260 neg_assign_r(tmp, negated_l, ROUND_DOWN);
00261 seq_i.lower_set_uninit(tmp);
00262 }
00263
00264
00265 seq_i.complete_init();
00266 }
00267 assert(OK());
00268 }
00269
00270 template <typename ITV>
00271 template <typename T>
00272 Box<ITV>::Box(const Octagonal_Shape<T>& oct, Complexity_Class)
00273 : seq(oct.space_dimension() <= max_space_dimension()
00274 ? oct.space_dimension()
00275 : (throw_space_dimension_overflow("Box(oct)",
00276 "oct exceeds the maximum "
00277 "allowed space dimension"),
00278 oct.space_dimension())),
00279 status() {
00280
00281 oct.strong_closure_assign();
00282 if (oct.marked_empty()) {
00283 set_empty();
00284 return;
00285 }
00286
00287
00288 set_empty_up_to_date();
00289
00290 const dimension_type space_dim = space_dimension();
00291 if (space_dim == 0)
00292 return;
00293
00294 DIRTY_TEMP0(mpq_class, bound);
00295 for (dimension_type i = space_dim; i-- > 0; ) {
00296 ITV& seq_i = seq[i];
00297 const dimension_type ii = 2*i;
00298 const dimension_type cii = ii + 1;
00299
00300
00301 const typename Octagonal_Shape<T>::coefficient_type& twice_ub
00302 = oct.matrix[cii][ii];
00303 if (!is_plus_infinity(twice_ub)) {
00304 assign_r(bound, twice_ub, ROUND_NOT_NEEDED);
00305 div2exp_assign_r(bound, bound, 1, ROUND_NOT_NEEDED);
00306 seq_i.upper_set_uninit(bound);
00307 }
00308 else
00309 seq_i.upper_set_uninit(UNBOUNDED);
00310
00311
00312 const typename Octagonal_Shape<T>::coefficient_type& twice_lb
00313 = oct.matrix[ii][cii];
00314 if (!is_plus_infinity(twice_lb)) {
00315 assign_r(bound, twice_lb, ROUND_NOT_NEEDED);
00316 neg_assign_r(bound, bound, ROUND_NOT_NEEDED);
00317 div2exp_assign_r(bound, bound, 1, ROUND_NOT_NEEDED);
00318 seq_i.lower_set_uninit(bound);
00319 }
00320 else
00321 seq_i.lower_set_uninit(UNBOUNDED);
00322 seq_i.complete_init();
00323 }
00324 }
00325
00326 template <typename ITV>
00327 Box<ITV>::Box(const Polyhedron& ph, Complexity_Class complexity)
00328 : seq(ph.space_dimension() <= max_space_dimension()
00329 ? ph.space_dimension()
00330 : (throw_space_dimension_overflow("Box(ph)",
00331 "ph exceeds the maximum "
00332 "allowed space dimension"),
00333 ph.space_dimension())),
00334 status() {
00335
00336 set_empty_up_to_date();
00337
00338
00339
00340 if (ph.marked_empty()) {
00341 set_empty();
00342 return;
00343 }
00344
00345
00346 const dimension_type space_dim = ph.space_dimension();
00347 if (space_dim == 0)
00348 return;
00349
00350
00351 if (ph.generators_are_up_to_date() && !ph.has_pending_constraints()) {
00352 Box tmp(ph.generators());
00353 swap(tmp);
00354 return;
00355 }
00356
00357
00358 assert(ph.constraints_are_up_to_date());
00359
00360 if (complexity == POLYNOMIAL_COMPLEXITY) {
00361
00362 for (dimension_type i = space_dimension(); i-- > 0; )
00363 seq[i].assign(UNIVERSE);
00364
00365 refine_with_constraints(ph.simplified_constraints());
00366 }
00367 else if (complexity == SIMPLEX_COMPLEXITY) {
00368 MIP_Problem lp(space_dim);
00369 const Constraint_System& ph_cs = ph.constraints();
00370 if (!ph_cs.has_strict_inequalities())
00371 lp.add_constraints(ph_cs);
00372 else
00373
00374 for (Constraint_System::const_iterator i = ph_cs.begin(),
00375 ph_cs_end = ph_cs.end(); i != ph_cs_end; ++i) {
00376 const Constraint& c = *i;
00377 if (c.is_strict_inequality())
00378 lp.add_constraint(Linear_Expression(c) >= 0);
00379 else
00380 lp.add_constraint(c);
00381 }
00382
00383 if (!lp.is_satisfiable()) {
00384 set_empty();
00385 return;
00386 }
00387
00388 Generator g(point());
00389 DIRTY_TEMP0(mpq_class, bound);
00390 DIRTY_TEMP(Coefficient, bound_num);
00391 DIRTY_TEMP(Coefficient, bound_den);
00392 for (dimension_type i = space_dim; i-- > 0; ) {
00393 ITV& seq_i = seq[i];
00394 lp.set_objective_function(Variable(i));
00395
00396 lp.set_optimization_mode(MAXIMIZATION);
00397 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00398 g = lp.optimizing_point();
00399 lp.evaluate_objective_function(g, bound_num, bound_den);
00400 assign_r(bound.get_num(), bound_num, ROUND_NOT_NEEDED);
00401 assign_r(bound.get_den(), bound_den, ROUND_NOT_NEEDED);
00402 assert(is_canonical(bound));
00403 seq_i.upper_set_uninit(bound);
00404 }
00405 else
00406 seq_i.upper_set_uninit(UNBOUNDED);
00407
00408 lp.set_optimization_mode(MINIMIZATION);
00409 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00410 g = lp.optimizing_point();
00411 lp.evaluate_objective_function(g, bound_num, bound_den);
00412 assign_r(bound.get_num(), bound_num, ROUND_NOT_NEEDED);
00413 assign_r(bound.get_den(), bound_den, ROUND_NOT_NEEDED);
00414 assert(is_canonical(bound));
00415 seq_i.lower_set_uninit(bound);
00416 }
00417 else
00418 seq_i.lower_set_uninit(UNBOUNDED);
00419 seq_i.complete_init();
00420 }
00421 }
00422 else {
00423 assert(complexity == ANY_COMPLEXITY);
00424 if (ph.is_empty())
00425 set_empty();
00426 else {
00427 Box tmp(ph.generators());
00428 swap(tmp);
00429 }
00430 }
00431 }
00432
00433 template <typename ITV>
00434 Box<ITV>::Box(const Grid& gr, Complexity_Class)
00435 : seq(gr.space_dimension() <= max_space_dimension()
00436 ? gr.space_dimension()
00437 : (throw_space_dimension_overflow("Box(gr)",
00438 "gr exceeds the maximum "
00439 "allowed space dimension"),
00440 gr.space_dimension())),
00441 status() {
00442
00443
00444
00445 if (gr.marked_empty()) {
00446 set_empty();
00447 return;
00448 }
00449
00450
00451 set_empty_up_to_date();
00452
00453 const dimension_type space_dim = gr.space_dimension();
00454
00455 if (space_dim == 0)
00456 return;
00457
00458 if (!gr.generators_are_up_to_date() && !gr.update_generators()) {
00459
00460 set_empty();
00461 return;
00462 }
00463
00464 assert(!gr.gen_sys.empty());
00465
00466
00467 std::vector<bool> bounded_interval(space_dim, true);
00468
00469 const Grid_Generator *first_point = 0;
00470
00471
00472
00473
00474
00475
00476 for (Grid_Generator_System::const_iterator gs_i = gr.gen_sys.begin(),
00477 gs_end = gr.gen_sys.end(); gs_i != gs_end; ++gs_i) {
00478 Grid_Generator& g = const_cast<Grid_Generator&>(*gs_i);
00479 if (g.is_point()) {
00480 if (first_point == 0) {
00481 first_point = &g;
00482 continue;
00483 }
00484 const Grid_Generator& point = *first_point;
00485
00486 for (dimension_type dim = space_dim; dim-- > 0; )
00487 g[dim] -= point[dim];
00488 g.set_divisor(point.divisor());
00489 }
00490 for (dimension_type col = space_dim; col > 0; )
00491 if (g[col--] != 0)
00492 bounded_interval[col] = false;
00493 }
00494
00495
00496
00497
00498 assert(first_point != 0);
00499 const Grid_Generator& point = *first_point;
00500 DIRTY_TEMP0(mpq_class, bound);
00501 const Coefficient& divisor = point.divisor();
00502 for (dimension_type i = space_dim; i-- > 0; ) {
00503 ITV& seq_i = seq[i];
00504 if (bounded_interval[i]) {
00505 assign_r(bound.get_num(), point[i+1], ROUND_NOT_NEEDED);
00506 assign_r(bound.get_den(), divisor, ROUND_NOT_NEEDED);
00507 bound.canonicalize();
00508 seq_i.assign(bound);
00509 }
00510 else
00511 seq_i.assign(UNIVERSE);
00512 }
00513 }
00514
00515 template <typename ITV>
00516 template <typename D1, typename D2, typename R>
00517 Box<ITV>::Box(const Partially_Reduced_Product<D1, D2, R>& dp,
00518 Complexity_Class complexity)
00519 : seq(), status() {
00520 if (dp.space_dimension() > max_space_dimension())
00521 throw_space_dimension_overflow("Box(dp)",
00522 "dp exceeds the maximum "
00523 "allowed space dimension");
00524 Box tmp1(dp.domain1(), complexity);
00525 Box tmp2(dp.domain2(), complexity);
00526 tmp1.intersection_assign(tmp2);
00527 swap(tmp1);
00528 }
00529
00530 template <typename ITV>
00531 inline void
00532 Box<ITV>::add_space_dimensions_and_embed(const dimension_type m) {
00533
00534 if (m == 0)
00535 return;
00536
00537
00538
00539 seq.insert(seq.end(), m, ITV());
00540 for (dimension_type sz = seq.size(), i = sz - m; i < sz; ++i)
00541 seq[i].assign(UNIVERSE);
00542 assert(OK());
00543 }
00544
00545 template <typename ITV>
00546 inline void
00547 Box<ITV>::add_space_dimensions_and_project(const dimension_type m) {
00548
00549 if (m == 0)
00550 return;
00551
00552
00553 seq.insert(seq.end(), m, ITV());
00554 for (dimension_type sz = seq.size(), i = sz - m; i < sz; ++i)
00555 seq[i].assign(0);
00556
00557 assert(OK());
00558 }
00559
00560 template <typename ITV>
00561 bool
00562 operator==(const Box<ITV>& x, const Box<ITV>& y) {
00563 const dimension_type x_space_dim = x.space_dimension();
00564 if (x_space_dim != y.space_dimension())
00565 return false;
00566
00567 if (x.is_empty())
00568 return y.is_empty();
00569
00570 if (y.is_empty())
00571 return x.is_empty();
00572
00573 for (dimension_type k = x_space_dim; k-- > 0; )
00574 if (x.seq[k] != y.seq[k])
00575 return false;
00576 return true;
00577 }
00578
00579 template <typename ITV>
00580 bool
00581 Box<ITV>::bounds(const Linear_Expression& expr, const bool from_above) const {
00582
00583 const dimension_type expr_space_dim = expr.space_dimension();
00584 const dimension_type space_dim = space_dimension();
00585 if (space_dim < expr_space_dim)
00586 throw_dimension_incompatible((from_above
00587 ? "bounds_from_above(e)"
00588 : "bounds_from_below(e)"), "e", expr);
00589
00590 if (space_dim == 0 || is_empty())
00591 return true;
00592
00593 const int from_above_sign = from_above ? 1 : -1;
00594 for (dimension_type i = expr_space_dim; i-- > 0; )
00595 switch (sgn(expr.coefficient(Variable(i))) * from_above_sign) {
00596 case 1:
00597 if (seq[i].upper_is_unbounded())
00598 return false;
00599 break;
00600 case 0:
00601
00602 break;
00603 case -1:
00604 if (seq[i].lower_is_unbounded())
00605 return false;
00606 break;
00607 }
00608 return true;
00609 }
00610
00611 template <typename ITV>
00612 Poly_Con_Relation
00613 interval_relation(const ITV& i,
00614 const Constraint::Type constraint_type,
00615 Coefficient_traits::const_reference num,
00616 Coefficient_traits::const_reference den) {
00617
00618 if (i.is_universe())
00619 return Poly_Con_Relation::strictly_intersects();
00620
00621 DIRTY_TEMP0(mpq_class, bound);
00622 assign_r(bound.get_num(), num, ROUND_NOT_NEEDED);
00623 assign_r(bound.get_den(), den, ROUND_NOT_NEEDED);
00624 bound.canonicalize();
00625 neg_assign_r(bound, bound, ROUND_NOT_NEEDED);
00626 const bool is_lower_bound = (den > 0);
00627
00628 DIRTY_TEMP0(mpq_class, bound_diff);
00629 if (constraint_type == Constraint::EQUALITY) {
00630 if (i.lower_is_unbounded()) {
00631 assert(!i.upper_is_unbounded());
00632 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
00633 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
00634 switch (sgn(bound_diff)) {
00635 case 1:
00636 return Poly_Con_Relation::strictly_intersects();
00637 case 0:
00638 return i.upper_is_open()
00639 ? Poly_Con_Relation::is_disjoint()
00640 : Poly_Con_Relation::strictly_intersects();
00641 case -1:
00642 return Poly_Con_Relation::is_disjoint();
00643 }
00644 }
00645 else {
00646 assign_r(bound_diff, i.lower(), ROUND_NOT_NEEDED);
00647 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
00648 switch (sgn(bound_diff)) {
00649 case 1:
00650 return Poly_Con_Relation::is_disjoint();
00651 case 0:
00652 if (i.lower_is_open())
00653 return Poly_Con_Relation::is_disjoint();
00654 if (i.is_singleton())
00655 return Poly_Con_Relation::is_included()
00656 && Poly_Con_Relation::saturates();
00657 return Poly_Con_Relation::strictly_intersects();
00658 case -1:
00659 if (i.upper_is_unbounded())
00660 return Poly_Con_Relation::strictly_intersects();
00661 else {
00662 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
00663 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
00664 switch (sgn(bound_diff)) {
00665 case 1:
00666 return Poly_Con_Relation::strictly_intersects();
00667 case 0:
00668 if (i.upper_is_open())
00669 return Poly_Con_Relation::is_disjoint();
00670 else
00671 return Poly_Con_Relation::strictly_intersects();
00672 case -1:
00673 return Poly_Con_Relation::is_disjoint();
00674 }
00675 }
00676 }
00677 }
00678 }
00679
00680 assert(constraint_type != Constraint::EQUALITY);
00681 if (is_lower_bound) {
00682 if (i.lower_is_unbounded()) {
00683 assert(!i.upper_is_unbounded());
00684 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
00685 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
00686 switch (sgn(bound_diff)) {
00687 case 1:
00688 return Poly_Con_Relation::strictly_intersects();
00689 case 0:
00690 if (constraint_type == Constraint::STRICT_INEQUALITY
00691 || i.upper_is_open())
00692 return Poly_Con_Relation::is_disjoint();
00693 else
00694 return Poly_Con_Relation::strictly_intersects();
00695 case -1:
00696 return Poly_Con_Relation::is_disjoint();
00697 }
00698 }
00699 else {
00700 assign_r(bound_diff, i.lower(), ROUND_NOT_NEEDED);
00701 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
00702 switch (sgn(bound_diff)) {
00703 case 1:
00704 return Poly_Con_Relation::is_included();
00705 case 0:
00706 if (constraint_type == Constraint::NONSTRICT_INEQUALITY
00707 || i.lower_is_open()) {
00708 Poly_Con_Relation result = Poly_Con_Relation::is_included();
00709 if (i.is_singleton())
00710 result = result && Poly_Con_Relation::saturates();
00711 return result;
00712 }
00713 else {
00714 assert(constraint_type == Constraint::STRICT_INEQUALITY
00715 && !i.lower_is_open());
00716 if (i.is_singleton())
00717 return Poly_Con_Relation::is_disjoint()
00718 && Poly_Con_Relation::saturates();
00719 else
00720 return Poly_Con_Relation::strictly_intersects();
00721 }
00722 case -1:
00723 if (i.upper_is_unbounded())
00724 return Poly_Con_Relation::strictly_intersects();
00725 else {
00726 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
00727 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
00728 switch (sgn(bound_diff)) {
00729 case 1:
00730 return Poly_Con_Relation::strictly_intersects();
00731 case 0:
00732 if (constraint_type == Constraint::STRICT_INEQUALITY
00733 || i.upper_is_open())
00734 return Poly_Con_Relation::is_disjoint();
00735 else
00736 return Poly_Con_Relation::strictly_intersects();
00737 case -1:
00738 return Poly_Con_Relation::is_disjoint();
00739 }
00740 }
00741 }
00742 }
00743 }
00744 else {
00745
00746 if (i.upper_is_unbounded())
00747 return Poly_Con_Relation::strictly_intersects();
00748 else {
00749 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
00750 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
00751 switch (sgn(bound_diff)) {
00752 case -1:
00753 return Poly_Con_Relation::is_included();
00754 case 0:
00755 if (constraint_type == Constraint::NONSTRICT_INEQUALITY
00756 || i.upper_is_open()) {
00757 Poly_Con_Relation result = Poly_Con_Relation::is_included();
00758 if (i.is_singleton())
00759 result = result && Poly_Con_Relation::saturates();
00760 return result;
00761 }
00762 else {
00763 assert(constraint_type == Constraint::STRICT_INEQUALITY
00764 && !i.upper_is_open());
00765 if (i.is_singleton())
00766 return Poly_Con_Relation::is_disjoint()
00767 && Poly_Con_Relation::saturates();
00768 else
00769 return Poly_Con_Relation::strictly_intersects();
00770 }
00771 case 1:
00772 if (i.lower_is_unbounded())
00773 return Poly_Con_Relation::strictly_intersects();
00774 else {
00775 assign_r(bound_diff, i.lower(), ROUND_NOT_NEEDED);
00776 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
00777 switch (sgn(bound_diff)) {
00778 case -1:
00779 return Poly_Con_Relation::strictly_intersects();
00780 case 0:
00781 if (constraint_type == Constraint::STRICT_INEQUALITY
00782 || i.lower_is_open())
00783 return Poly_Con_Relation::is_disjoint();
00784 else
00785 return Poly_Con_Relation::strictly_intersects();
00786 case 1:
00787 return Poly_Con_Relation::is_disjoint();
00788 }
00789 }
00790 }
00791 }
00792 }
00793
00794
00795 throw std::runtime_error("PPL internal error");
00796 }
00797
00798 template <typename ITV>
00799 Poly_Con_Relation
00800 Box<ITV>::relation_with(const Congruence& cg) const {
00801 const dimension_type cg_space_dim = cg.space_dimension();
00802 const dimension_type space_dim = space_dimension();
00803
00804
00805 if (cg_space_dim > space_dim)
00806 throw_dimension_incompatible("relation_with(cg)", cg);
00807
00808 if (is_empty())
00809 return Poly_Con_Relation::saturates()
00810 && Poly_Con_Relation::is_included()
00811 && Poly_Con_Relation::is_disjoint();
00812
00813 if (space_dim == 0) {
00814 if (cg.is_inconsistent())
00815 return Poly_Con_Relation::is_disjoint();
00816 else
00817 return Poly_Con_Relation::saturates()
00818 && Poly_Con_Relation::is_included();
00819 }
00820
00821 if (cg.is_equality()) {
00822 const Constraint c(cg);
00823 return relation_with(c);
00824 }
00825
00826 DIRTY_TEMP0(Rational_Interval, r);
00827 DIRTY_TEMP0(Rational_Interval, t);
00828 DIRTY_TEMP0(mpq_class, m);
00829 r = 0;
00830 for (dimension_type i = cg.space_dimension(); i-- > 0; ) {
00831 const Coefficient& cg_i = cg.coefficient(Variable(i));
00832 if (sgn(cg_i) != 0) {
00833 assign_r(m, cg_i, ROUND_NOT_NEEDED);
00834
00835 t = seq[i];
00836 t *= m;
00837 r += t;
00838 }
00839 }
00840
00841 if (r.lower_is_unbounded() || r.upper_is_unbounded())
00842 return Poly_Con_Relation::strictly_intersects();
00843
00844
00845
00846
00847
00848 TEMP_INTEGER(lower);
00849 TEMP_INTEGER(mod);
00850 TEMP_INTEGER(v);
00851 mod = cg.modulus();
00852 v = cg.inhomogeneous_term() % mod;
00853 assign_r(lower, r.lower(), ROUND_DOWN);
00854 v -= ((lower / mod) * mod);
00855 if (v + lower > 0)
00856 v -= mod;
00857 return interval_relation(r, Constraint::EQUALITY, v);
00858 }
00859
00860 template <typename ITV>
00861 Poly_Con_Relation
00862 Box<ITV>::relation_with(const Constraint& c) const {
00863 const dimension_type c_space_dim = c.space_dimension();
00864 const dimension_type space_dim = space_dimension();
00865
00866
00867 if (c_space_dim > space_dim)
00868 throw_dimension_incompatible("relation_with(c)", c);
00869
00870 if (is_empty())
00871 return Poly_Con_Relation::saturates()
00872 && Poly_Con_Relation::is_included()
00873 && Poly_Con_Relation::is_disjoint();
00874
00875 if (space_dim == 0) {
00876 if ((c.is_equality() && c.inhomogeneous_term() != 0)
00877 || (c.is_inequality() && c.inhomogeneous_term() < 0))
00878 return Poly_Con_Relation::is_disjoint();
00879 else if (c.is_strict_inequality() && c.inhomogeneous_term() == 0)
00880
00881
00882 return Poly_Con_Relation::saturates()
00883 && Poly_Con_Relation::is_disjoint();
00884 else if (c.is_equality() || c.inhomogeneous_term() == 0)
00885 return Poly_Con_Relation::saturates()
00886 && Poly_Con_Relation::is_included();
00887 else
00888
00889
00890
00891 return Poly_Con_Relation::is_included();
00892 }
00893
00894 dimension_type c_num_vars = 0;
00895 dimension_type c_only_var = 0;
00896
00897 if (extract_interval_constraint(c, c_space_dim, c_num_vars, c_only_var))
00898 if (c_num_vars == 0)
00899
00900 switch (sgn(c.inhomogeneous_term())) {
00901 case -1:
00902 return Poly_Con_Relation::is_disjoint();
00903 case 0:
00904 if (c.is_strict_inequality())
00905 return Poly_Con_Relation::saturates()
00906 && Poly_Con_Relation::is_disjoint();
00907 else
00908 return Poly_Con_Relation::saturates()
00909 && Poly_Con_Relation::is_included();
00910 case 1:
00911 return Poly_Con_Relation::is_included();
00912 }
00913 else {
00914
00915 return interval_relation(seq[c_only_var],
00916 c.type(),
00917 c.inhomogeneous_term(),
00918 c.coefficient(Variable(c_only_var)));
00919 }
00920 else {
00921
00922 DIRTY_TEMP0(Rational_Interval, r);
00923 DIRTY_TEMP0(Rational_Interval, t);
00924 DIRTY_TEMP0(mpq_class, m);
00925 r = 0;
00926 for (dimension_type i = c.space_dimension(); i-- > 0; ) {
00927 const Coefficient& c_i = c.coefficient(Variable(i));
00928 if (sgn(c_i) != 0) {
00929 assign_r(m, c_i, ROUND_NOT_NEEDED);
00930
00931 t = seq[i];
00932 t *= m;
00933 r += t;
00934 }
00935 }
00936 return interval_relation(r,
00937 c.type(),
00938 c.inhomogeneous_term());
00939 }
00940
00941
00942 throw std::runtime_error("PPL internal error");
00943 }
00944
00945 template <typename ITV>
00946 Poly_Gen_Relation
00947 Box<ITV>::relation_with(const Generator& g) const {
00948 const dimension_type space_dim = space_dimension();
00949 const dimension_type g_space_dim = g.space_dimension();
00950
00951
00952 if (space_dim < g_space_dim)
00953 throw_dimension_incompatible("relation_with(g)", g);
00954
00955
00956 if (is_empty())
00957 return Poly_Gen_Relation::nothing();
00958
00959
00960
00961 if (space_dim == 0)
00962 return Poly_Gen_Relation::subsumes();
00963
00964 if (g.is_line_or_ray()) {
00965 if (g.is_line()) {
00966 for (dimension_type i = g_space_dim; i-- > 0; )
00967 if (g.coefficient(Variable(i)) != 0 && !seq[i].is_universe())
00968 return Poly_Gen_Relation::nothing();
00969 return Poly_Gen_Relation::subsumes();
00970 }
00971 else {
00972 assert(g.is_ray());
00973 for (dimension_type i = g_space_dim; i-- > 0; )
00974 switch (sgn(g.coefficient(Variable(i)))) {
00975 case 1:
00976 if (!seq[i].upper_is_unbounded())
00977 return Poly_Gen_Relation::nothing();
00978 break;
00979 case 0:
00980 break;
00981 case -1:
00982 if (!seq[i].lower_is_unbounded())
00983 return Poly_Gen_Relation::nothing();
00984 break;
00985 }
00986 return Poly_Gen_Relation::subsumes();
00987 }
00988 }
00989
00990
00991 const Coefficient& g_divisor = g.divisor();
00992 DIRTY_TEMP0(mpq_class, g_coord);
00993 DIRTY_TEMP0(mpq_class, bound);
00994 for (dimension_type i = g_space_dim; i-- > 0; ) {
00995 const ITV& seq_i = seq[i];
00996 if (seq_i.is_universe())
00997 continue;
00998 assign_r(g_coord.get_num(), g.coefficient(Variable(i)), ROUND_NOT_NEEDED);
00999 assign_r(g_coord.get_den(), g_divisor, ROUND_NOT_NEEDED);
01000 g_coord.canonicalize();
01001
01002 if (!seq_i.lower_is_unbounded()) {
01003 assign_r(bound, seq_i.lower(), ROUND_NOT_NEEDED);
01004 if (g_coord <= bound) {
01005 if (seq_i.lower_is_open()) {
01006 if (g.is_point() || g_coord != bound)
01007 return Poly_Gen_Relation::nothing();
01008 }
01009 else if (g_coord != bound)
01010 return Poly_Gen_Relation::nothing();
01011 }
01012 }
01013
01014 if (!seq_i.upper_is_unbounded()) {
01015 assign_r(bound, seq_i.upper(), ROUND_NOT_NEEDED);
01016 if (g_coord >= bound) {
01017 if (seq_i.upper_is_open()) {
01018 if (g.is_point() || g_coord != bound)
01019 return Poly_Gen_Relation::nothing();
01020 }
01021 else if (g_coord != bound)
01022 return Poly_Gen_Relation::nothing();
01023 }
01024 }
01025 }
01026 return Poly_Gen_Relation::subsumes();
01027 }
01028
01029
01030 template <typename ITV>
01031 bool
01032 Box<ITV>::max_min(const Linear_Expression& expr,
01033 const bool maximize,
01034 Coefficient& ext_n, Coefficient& ext_d,
01035 bool& included) const {
01036
01037 const dimension_type space_dim = space_dimension();
01038 const dimension_type expr_space_dim = expr.space_dimension();
01039 if (space_dim < expr_space_dim)
01040 throw_dimension_incompatible((maximize
01041 ? "maximize(e, ...)"
01042 : "minimize(e, ...)"), "e", expr);
01043
01044 if (space_dim == 0) {
01045 if (marked_empty())
01046 return false;
01047 else {
01048 ext_n = expr.inhomogeneous_term();
01049 ext_d = 1;
01050 included = true;
01051 return true;
01052 }
01053 }
01054
01055
01056 if (is_empty())
01057 return false;
01058
01059 DIRTY_TEMP0(mpq_class, result);
01060 assign_r(result, expr.inhomogeneous_term(), ROUND_NOT_NEEDED);
01061 bool is_included = true;
01062 const int maximize_sign = maximize ? 1 : -1;
01063 DIRTY_TEMP0(mpq_class, bound_i);
01064 DIRTY_TEMP0(mpq_class, expr_i);
01065 for (dimension_type i = expr_space_dim; i-- > 0; ) {
01066 const ITV& seq_i = seq[i];
01067 assign_r(expr_i, expr.coefficient(Variable(i)), ROUND_NOT_NEEDED);
01068 switch (sgn(expr_i) * maximize_sign) {
01069 case 1:
01070 if (seq_i.upper_is_unbounded())
01071 return false;
01072 assign_r(bound_i, seq_i.upper(), ROUND_NOT_NEEDED);
01073 add_mul_assign_r(result, bound_i, expr_i, ROUND_NOT_NEEDED);
01074 if (seq_i.upper_is_open())
01075 is_included = false;
01076 break;
01077 case 0:
01078
01079 break;
01080 case -1:
01081 if (seq_i.lower_is_unbounded())
01082 return false;
01083 assign_r(bound_i, seq_i.lower(), ROUND_NOT_NEEDED);
01084 add_mul_assign_r(result, bound_i, expr_i, ROUND_NOT_NEEDED);
01085 if (seq_i.lower_is_open())
01086 is_included = false;
01087 break;
01088 }
01089 }
01090
01091 assert(is_canonical(result));
01092 ext_n = result.get_num();
01093 ext_d = result.get_den();
01094 included = is_included;
01095 return true;
01096 }
01097
01098 template <typename ITV>
01099 bool
01100 Box<ITV>::max_min(const Linear_Expression& expr,
01101 const bool maximize,
01102 Coefficient& ext_n, Coefficient& ext_d,
01103 bool& included,
01104 Generator& g) const {
01105 if (!max_min(expr, maximize, ext_n, ext_d, included))
01106 return false;
01107
01108
01109 Linear_Expression g_expr;
01110 DIRTY_TEMP(Coefficient, g_divisor);
01111 g_divisor = 1;
01112 const int maximize_sign = maximize ? 1 : -1;
01113 DIRTY_TEMP0(mpq_class, g_coord);
01114 DIRTY_TEMP(Coefficient, num);
01115 DIRTY_TEMP(Coefficient, den);
01116 DIRTY_TEMP(Coefficient, lcm);
01117 DIRTY_TEMP(Coefficient, factor);
01118 for (dimension_type i = space_dimension(); i-- > 0; ) {
01119 const ITV& seq_i = seq[i];
01120 switch (sgn(expr.coefficient(Variable(i))) * maximize_sign) {
01121 case 1:
01122 assign_r(g_coord, seq_i.upper(), ROUND_NOT_NEEDED);
01123 break;
01124 case 0:
01125
01126
01127
01128 if (seq_i.contains(0))
01129 continue;
01130 if (!seq_i.lower_is_unbounded())
01131 if (seq_i.lower_is_open())
01132 if (!seq_i.upper_is_unbounded())
01133 if (seq_i.upper_is_open()) {
01134
01135 assign_r(g_coord, seq_i.lower(), ROUND_NOT_NEEDED);
01136 DIRTY_TEMP0(mpq_class, q_seq_i_upper);
01137 assign_r(q_seq_i_upper, seq_i.upper(), ROUND_NOT_NEEDED);
01138 g_coord += q_seq_i_upper;
01139 g_coord /= 2;
01140 }
01141 else
01142
01143 assign_r(g_coord, seq_i.upper(), ROUND_NOT_NEEDED);
01144 else {
01145
01146 assign_r(g_coord, seq_i.lower(), ROUND_NOT_NEEDED);
01147 ++g_coord;
01148 }
01149 else
01150
01151 assign_r(g_coord, seq_i.lower(), ROUND_NOT_NEEDED);
01152 else {
01153
01154
01155 assert(!seq_i.upper_is_unbounded());
01156 assign_r(g_coord, seq_i.upper(), ROUND_NOT_NEEDED);
01157 if (seq_i.upper_is_open())
01158 --g_coord;
01159 }
01160 break;
01161 case -1:
01162 assign_r(g_coord, seq_i.lower(), ROUND_NOT_NEEDED);
01163 break;
01164 }
01165
01166 assign_r(den, g_coord.get_den(), ROUND_NOT_NEEDED);
01167 lcm_assign(lcm, g_divisor, den);
01168 exact_div_assign(factor, lcm, g_divisor);
01169 g_expr *= factor;
01170 exact_div_assign(factor, lcm, den);
01171 assign_r(num, g_coord.get_num(), ROUND_NOT_NEEDED);
01172 num *= factor;
01173 g_expr += num * Variable(i);
01174 g_divisor = lcm;
01175 }
01176 g = Generator::point(g_expr, g_divisor);
01177 return true;
01178 }
01179
01180 template <typename ITV>
01181 bool
01182 Box<ITV>::contains(const Box& y) const {
01183 const Box& x = *this;
01184
01185 if (x.space_dimension() != y.space_dimension())
01186 x.throw_dimension_incompatible("contains(y)", y);
01187
01188
01189 if (y.is_empty())
01190 return true;
01191
01192
01193 if (x.is_empty())
01194 return false;
01195
01196 for (dimension_type k = x.seq.size(); k-- > 0; )
01197
01198 if (!x.seq[k].contains(y.seq[k]))
01199 return false;
01200 return true;
01201 }
01202
01203 template <typename ITV>
01204 bool
01205 Box<ITV>::is_disjoint_from(const Box& y) const {
01206 const Box& x = *this;
01207
01208 if (x.space_dimension() != y.space_dimension())
01209 x.throw_dimension_incompatible("is_disjoint_from(y)", y);
01210
01211
01212
01213 if (x.marked_empty() || y.marked_empty())
01214 return true;
01215
01216 for (dimension_type k = x.seq.size(); k-- > 0; )
01217
01218 if (x.seq[k].is_disjoint_from(y.seq[k]))
01219 return true;
01220 return false;
01221 }
01222
01223 template <typename ITV>
01224 bool
01225 Box<ITV>::OK() const {
01226 if (status.test_empty_up_to_date() && !status.test_empty()) {
01227 Box tmp = *this;
01228 tmp.reset_empty_up_to_date();
01229 if (tmp.check_empty()) {
01230 #ifndef NDEBUG
01231 std::cerr << "The box is empty, but it is marked as non-empty."
01232 << std::endl;
01233 #endif // NDEBUG
01234 return false;
01235 }
01236 }
01237
01238
01239 if (!marked_empty()) {
01240 for (dimension_type k = seq.size(); k-- > 0; )
01241 if (!seq[k].OK())
01242 return false;
01243 }
01244
01245 return true;
01246 }
01247
01248 template <typename ITV>
01249 dimension_type
01250 Box<ITV>::affine_dimension() const {
01251 dimension_type d = space_dimension();
01252
01253 if (d == 0)
01254 return 0;
01255
01256
01257 if (is_empty())
01258 return 0;
01259
01260 for (dimension_type k = d; k-- > 0; )
01261 if (seq[k].is_singleton())
01262 --d;
01263
01264 return d;
01265 }
01266
01267 template <typename ITV>
01268 bool
01269 Box<ITV>::check_empty() const {
01270 assert(!marked_empty());
01271 Box<ITV>& x = const_cast<Box<ITV>&>(*this);
01272 for (dimension_type k = seq.size(); k-- > 0; )
01273 if (seq[k].is_empty()) {
01274 x.set_empty();
01275 return true;
01276 }
01277 x.set_nonempty();;
01278 return false;
01279 }
01280
01281 template <typename ITV>
01282 bool
01283 Box<ITV>::is_universe() const {
01284 if (marked_empty())
01285 return false;
01286 for (dimension_type k = seq.size(); k-- > 0; )
01287 if (!seq[k].is_universe())
01288 return false;
01289 return true;
01290 }
01291
01292 template <typename ITV>
01293 bool
01294 Box<ITV>::is_topologically_closed() const {
01295 if (!ITV::info_type::store_open || is_empty())
01296 return true;
01297
01298 for (dimension_type k = seq.size(); k-- > 0; )
01299 if (!seq[k].is_topologically_closed())
01300 return false;
01301 return true;
01302 }
01303
01304 template <typename ITV>
01305 bool
01306 Box<ITV>::is_discrete() const {
01307 if (is_empty())
01308 return true;
01309 for (dimension_type k = seq.size(); k-- > 0; )
01310 if (!seq[k].is_singleton())
01311 return false;
01312 return true;
01313 }
01314
01315 template <typename ITV>
01316 bool
01317 Box<ITV>::is_bounded() const {
01318 if (is_empty())
01319 return true;
01320 for (dimension_type k = seq.size(); k-- > 0; )
01321 if (seq[k].is_unbounded())
01322 return false;
01323 return true;
01324 }
01325
01326 template <typename ITV>
01327 bool
01328 Box<ITV>::contains_integer_point() const {
01329 if (marked_empty())
01330 return false;
01331 for (dimension_type k = seq.size(); k-- > 0; )
01332 if (!seq[k].contains_integer_point())
01333 return false;
01334 return true;
01335 }
01336
01337 template <typename ITV>
01338 bool
01339 Box<ITV>::constrains(Variable var) const {
01340
01341 const dimension_type var_space_dim = var.space_dimension();
01342 if (space_dimension() < var_space_dim)
01343 throw_dimension_incompatible("constrains(v)", "v", var);
01344
01345 if (marked_empty() || !seq[var_space_dim-1].is_universe())
01346 return true;
01347
01348 return is_empty();
01349 }
01350
01351 template <typename ITV>
01352 void
01353 Box<ITV>::unconstrain(const Variables_Set& to_be_unconstrained) {
01354
01355
01356
01357 if (to_be_unconstrained.empty())
01358 return;
01359
01360
01361 const dimension_type min_space_dim = to_be_unconstrained.space_dimension();
01362 if (space_dimension() < min_space_dim)
01363 throw_dimension_incompatible("unconstrain(vs)", min_space_dim);
01364
01365
01366 if (marked_empty())
01367 return;
01368
01369
01370
01371
01372 for (Variables_Set::const_iterator tbu = to_be_unconstrained.begin(),
01373 tbu_end = to_be_unconstrained.end(); tbu != tbu_end; ++tbu) {
01374 ITV& seq_tbu = seq[*tbu];
01375 if (!seq_tbu.is_empty())
01376 seq_tbu.assign(UNIVERSE);
01377 else {
01378 set_empty();
01379 break;
01380 }
01381 }
01382 assert(OK());
01383 }
01384
01385 template <typename ITV>
01386 void
01387 Box<ITV>::topological_closure_assign() {
01388 if (!ITV::info_type::store_open || is_empty())
01389 return;
01390
01391 for (dimension_type k = seq.size(); k-- > 0; )
01392 seq[k].topological_closure_assign();
01393 }
01394
01395 template <typename ITV>
01396 void
01397 Box<ITV>::intersection_assign(const Box& y) {
01398 Box& x = *this;
01399 const dimension_type space_dim = space_dimension();
01400
01401
01402 if (space_dim != y.space_dimension())
01403 x.throw_dimension_incompatible("intersection_assign(y)", y);
01404
01405
01406 if (x.marked_empty())
01407 return;
01408 if (y.marked_empty()) {
01409 x.set_empty();
01410 return;
01411 }
01412
01413
01414
01415 if (space_dim == 0)
01416 return;
01417
01418
01419
01420 reset_empty_up_to_date();
01421
01422 for (dimension_type k = space_dim; k-- > 0; )
01423 x.seq[k].intersect_assign(y.seq[k]);
01424
01425 assert(x.OK());
01426 }
01427
01428 template <typename ITV>
01429 void
01430 Box<ITV>::upper_bound_assign(const Box& y) {
01431 Box& x = *this;
01432
01433
01434 if (x.space_dimension() != y.space_dimension())
01435 x.throw_dimension_incompatible("upper_bound_assign(y)", y);
01436
01437
01438 if (y.marked_empty())
01439 return;
01440 if (x.marked_empty()) {
01441 x = y;
01442 return;
01443 }
01444
01445 for (dimension_type k = x.seq.size(); k-- > 0; )
01446 x.seq[k].join_assign(y.seq[k]);
01447
01448 assert(x.OK());
01449 }
01450
01451 template <typename ITV>
01452 void
01453 Box<ITV>::concatenate_assign(const Box& y) {
01454 Box& x = *this;
01455 const dimension_type x_space_dim = x.space_dimension();
01456 const dimension_type y_space_dim = y.space_dimension();
01457
01458
01459 if (y.marked_empty())
01460 x.set_empty();
01461
01462
01463 if (y_space_dim == 0)
01464 return;
01465
01466
01467
01468 x.seq.reserve(x_space_dim + y_space_dim);
01469
01470
01471
01472 if (x.marked_empty()) {
01473 x.seq.insert(x.seq.end(), y_space_dim, ITV());
01474 assert(x.OK());
01475 return;
01476 }
01477
01478
01479 std::copy(y.seq.begin(), y.seq.end(),
01480 std::back_insert_iterator<Sequence>(x.seq));
01481
01482 if (!y.status.test_empty_up_to_date())
01483 reset_empty_up_to_date();
01484
01485 assert(x.OK());
01486 }
01487
01488 template <typename ITV>
01489 void
01490 Box<ITV>::difference_assign(const Box& y) {
01491 const dimension_type space_dim = space_dimension();
01492
01493
01494 if (space_dim != y.space_dimension())
01495 throw_dimension_incompatible("difference_assign(y)", y);
01496
01497 Box& x = *this;
01498 if (x.is_empty() || y.is_empty())
01499 return;
01500
01501 switch (space_dim) {
01502 case 0:
01503
01504
01505 x.set_empty();
01506 break;
01507
01508 case 1:
01509 x.seq[0].difference_assign(y.seq[0]);
01510 if (x.seq[0].is_empty())
01511 x.set_empty();
01512 break;
01513
01514 default:
01515 {
01516 dimension_type index_non_contained = space_dim;
01517 dimension_type number_non_contained = 0;
01518 for (dimension_type i = space_dim; i-- > 0; )
01519 if (!y.seq[i].contains(x.seq[i])) {
01520 if (++number_non_contained == 1)
01521 index_non_contained = i;
01522 else
01523 break;
01524 }
01525
01526 switch (number_non_contained) {
01527 case 0:
01528
01529 x.set_empty();
01530 break;
01531 case 1:
01532 x.seq[index_non_contained]
01533 .difference_assign(y.seq[index_non_contained]);
01534 if (x.seq[index_non_contained].is_empty())
01535 x.set_empty();
01536 break;
01537 default:
01538
01539 break;
01540 }
01541 }
01542 break;
01543 }
01544 assert(OK());
01545 }
01546
01547 template <typename ITV>
01548 bool
01549 Box<ITV>::simplify_using_context_assign(const Box& y) {
01550
01551 used(y);
01552 return true;
01553 }
01554
01555 template <typename ITV>
01556 void
01557 Box<ITV>::time_elapse_assign(const Box& y) {
01558 Box& x = *this;
01559 const dimension_type x_space_dim = x.space_dimension();
01560
01561
01562 if (x_space_dim != y.space_dimension())
01563 x.throw_dimension_incompatible("time_elapse_assign(y)", y);
01564
01565
01566 if (x_space_dim == 0) {
01567 if (y.marked_empty())
01568 x.set_empty();
01569 return;
01570 }
01571
01572
01573
01574 if (x.marked_empty() || y.marked_empty()
01575 || x.is_empty() || y.is_empty()) {
01576 x.set_empty();
01577 return;
01578 }
01579
01580 for (dimension_type i = x_space_dim; i-- > 0; ) {
01581 ITV& x_seq_i = x.seq[i];
01582 const ITV& y_seq_i = y.seq[i];
01583 if (!x_seq_i.lower_is_unbounded())
01584 if (y_seq_i.lower_is_unbounded() || y_seq_i.lower() < 0)
01585 x_seq_i.lower_set(UNBOUNDED);
01586 if (!x_seq_i.upper_is_unbounded())
01587 if (y_seq_i.upper_is_unbounded() || y_seq_i.upper() > 0)
01588 x_seq_i.upper_set(UNBOUNDED);
01589 }
01590 assert(x.OK());
01591 }
01592
01593 template <typename ITV>
01594 inline void
01595 Box<ITV>::remove_space_dimensions(const Variables_Set& to_be_removed) {
01596
01597
01598
01599 if (to_be_removed.empty()) {
01600 assert(OK());
01601 return;
01602 }
01603
01604 const dimension_type old_space_dim = space_dimension();
01605
01606
01607 const dimension_type tbr_space_dim = to_be_removed.space_dimension();
01608 if (old_space_dim < tbr_space_dim)
01609 throw_dimension_incompatible("remove_space_dimensions(vs)",
01610 tbr_space_dim);
01611
01612 const dimension_type new_space_dim = old_space_dim - to_be_removed.size();
01613
01614
01615
01616
01617 if (is_empty() || new_space_dim == 0) {
01618 seq.resize(new_space_dim);
01619 assert(OK());
01620 return;
01621 }
01622
01623
01624
01625 Variables_Set::const_iterator tbr = to_be_removed.begin();
01626 Variables_Set::const_iterator tbr_end = to_be_removed.end();
01627 dimension_type dst = *tbr;
01628 dimension_type src = dst + 1;
01629 for (++tbr; tbr != tbr_end; ++tbr) {
01630 const dimension_type tbr_next = *tbr;
01631
01632 while (src < tbr_next)
01633 seq[dst++].swap(seq[src++]);
01634 ++src;
01635 }
01636
01637 while (src < old_space_dim)
01638 seq[dst++].swap(seq[src++]);
01639
01640 assert(dst == new_space_dim);
01641 seq.resize(new_space_dim);
01642
01643 assert(OK());
01644 }
01645
01646 template <typename ITV>
01647 void
01648 Box<ITV>::remove_higher_space_dimensions(const dimension_type new_dim) {
01649
01650
01651 const dimension_type old_dim = space_dimension();
01652 if (new_dim > old_dim)
01653 throw_dimension_incompatible("remove_higher_space_dimensions(nd)",
01654 new_dim);
01655
01656
01657
01658
01659 if (new_dim == old_dim) {
01660 assert(OK());
01661 return;
01662 }
01663
01664 seq.erase(seq.begin() + new_dim, seq.end());
01665 assert(OK());
01666 }
01667
01668 template <typename ITV>
01669 template <typename Partial_Function>
01670 void
01671 Box<ITV>::map_space_dimensions(const Partial_Function& pfunc) {
01672 const dimension_type space_dim = space_dimension();
01673 if (space_dim == 0)
01674 return;
01675
01676 if (pfunc.has_empty_codomain()) {
01677
01678 remove_higher_space_dimensions(0);
01679 return;
01680 }
01681
01682 const dimension_type new_space_dim = pfunc.max_in_codomain() + 1;
01683
01684 if (is_empty()) {
01685 remove_higher_space_dimensions(new_space_dim);
01686 return;
01687 }
01688
01689
01690 Box<ITV> tmp(new_space_dim);
01691
01692 for (dimension_type i = 0; i < space_dim; ++i) {
01693 dimension_type new_i;
01694 if (pfunc.maps(i, new_i))
01695 seq[i].swap(tmp.seq[new_i]);
01696 }
01697 swap(tmp);
01698 assert(OK());
01699 }
01700
01701 template <typename ITV>
01702 void
01703 Box<ITV>::fold_space_dimensions(const Variables_Set& to_be_folded,
01704 const Variable var) {
01705 const dimension_type space_dim = space_dimension();
01706
01707 if (var.space_dimension() > space_dim)
01708 throw_dimension_incompatible("fold_space_dimensions(tbf, v)", "v", var);
01709
01710
01711 if (to_be_folded.empty())
01712 return;
01713
01714
01715 if (to_be_folded.space_dimension() > space_dim)
01716 throw_dimension_incompatible("fold_space_dimensions(tbf, ...)",
01717 to_be_folded.space_dimension());
01718
01719
01720 if (to_be_folded.find(var.id()) != to_be_folded.end())
01721 throw_generic("fold_space_dimensions(tbf, v)",
01722 "v should not occur in tbf");
01723
01724
01725 if (!is_empty()) {
01726
01727
01728 ITV& seq_v = seq[var.id()];
01729 for (Variables_Set::const_iterator i = to_be_folded.begin(),
01730 tbf_end = to_be_folded.end(); i != tbf_end; ++i)
01731 seq_v.join_assign(seq[*i]);
01732 }
01733 remove_space_dimensions(to_be_folded);
01734 }
01735
01736 template <typename ITV>
01737 void
01738 Box<ITV>::add_constraint_no_check(const Constraint& c) {
01739 const dimension_type c_space_dim = c.space_dimension();
01740 assert(c_space_dim <= space_dimension());
01741
01742 dimension_type c_num_vars = 0;
01743 dimension_type c_only_var = 0;
01744
01745 if (!extract_interval_constraint(c, c_space_dim, c_num_vars, c_only_var))
01746 throw_generic("add_constraint(c)", "c is not an interval constraint");
01747
01748
01749
01750 if (c.is_strict_inequality() && c_num_vars != 0
01751 && !Box::interval_type::info_type::store_open)
01752 throw_generic("add_constraint(c)", "c is a nontrivial strict constraint");
01753
01754
01755 if (marked_empty())
01756 return;
01757
01758 const Coefficient& n = c.inhomogeneous_term();
01759 if (c_num_vars == 0) {
01760
01761 if (n < 0
01762 || (c.is_equality() && n != 0)
01763 || (c.is_strict_inequality() && n == 0))
01764 set_empty();
01765 return;
01766 }
01767
01768 assert(c_num_vars == 1);
01769 const Coefficient& d = c.coefficient(Variable(c_only_var));
01770 add_interval_constraint_no_check(c_only_var, c.type(), n, d);
01771 }
01772
01773 template <typename ITV>
01774 void
01775 Box<ITV>::add_constraints_no_check(const Constraint_System& cs) {
01776 assert(cs.space_dimension() <= space_dimension());
01777
01778
01779
01780 for (Constraint_System::const_iterator i = cs.begin(),
01781 cs_end = cs.end(); i != cs_end; ++i)
01782 add_constraint_no_check(*i);
01783 assert(OK());
01784 }
01785
01786 template <typename ITV>
01787 void
01788 Box<ITV>::add_congruence_no_check(const Congruence& cg) {
01789 const dimension_type cg_space_dim = cg.space_dimension();
01790 assert(cg_space_dim <= space_dimension());
01791
01792
01793 if (cg.is_proper_congruence()) {
01794 if (cg.is_inconsistent()) {
01795 set_empty();
01796 return;
01797 }
01798 else if (cg.is_tautological())
01799 return;
01800 else
01801
01802 throw_generic("add_congruence(cg)",
01803 "cg is a nontrivial proper congruence");
01804 }
01805
01806 assert(cg.is_equality());
01807 dimension_type cg_num_vars = 0;
01808 dimension_type cg_only_var = 0;
01809
01810 if (!extract_interval_congruence(cg, cg_space_dim, cg_num_vars, cg_only_var))
01811 throw_generic("add_congruence(cg)", "cg is not an interval congruence");
01812
01813
01814 if (marked_empty())
01815 return;
01816
01817 const Coefficient& n = cg.inhomogeneous_term();
01818 if (cg_num_vars == 0) {
01819
01820 if (n != 0)
01821 set_empty();
01822 return;
01823 }
01824
01825 assert(cg_num_vars == 1);
01826 const Coefficient& d = cg.coefficient(Variable(cg_only_var));
01827 add_interval_constraint_no_check(cg_only_var, Constraint::EQUALITY, n, d);
01828 }
01829
01830 template <typename ITV>
01831 void
01832 Box<ITV>::add_congruences_no_check(const Congruence_System& cgs) {
01833 assert(cgs.space_dimension() <= space_dimension());
01834
01835
01836
01837 for (Congruence_System::const_iterator i = cgs.begin(),
01838 cgs_end = cgs.end(); i != cgs_end; ++i)
01839 add_congruence_no_check(*i);
01840 assert(OK());
01841 }
01842
01843 template <typename ITV>
01844 void
01845 Box<ITV>::refine_no_check(const Constraint& c) {
01846 const dimension_type c_space_dim = c.space_dimension();
01847 assert(c_space_dim <= space_dimension());
01848 assert(!marked_empty());
01849
01850 dimension_type c_num_vars = 0;
01851 dimension_type c_only_var = 0;
01852
01853 if (!extract_interval_constraint(c, c_space_dim, c_num_vars, c_only_var))
01854 return;
01855
01856 const Coefficient& n = c.inhomogeneous_term();
01857 if (c_num_vars == 0) {
01858
01859 if (n < 0
01860 || (c.is_equality() && n != 0)
01861 || (c.is_strict_inequality() && n == 0))
01862 set_empty();
01863 return;
01864 }
01865
01866 assert(c_num_vars == 1);
01867 const Coefficient& d = c.coefficient(Variable(c_only_var));
01868 add_interval_constraint_no_check(c_only_var, c.type(), n, d);
01869 }
01870
01871 template <typename ITV>
01872 void
01873 Box<ITV>::refine_no_check(const Constraint_System& cs) {
01874 assert(cs.space_dimension() <= space_dimension());
01875 for (Constraint_System::const_iterator i = cs.begin(),
01876 cs_end = cs.end(); !marked_empty() && i != cs_end; ++i)
01877 refine_no_check(*i);
01878 assert(OK());
01879 }
01880
01881 template <typename ITV>
01882 void
01883 Box<ITV>::refine_no_check(const Congruence& cg) {
01884 assert(!marked_empty());
01885
01886 const dimension_type cg_space_dim = cg.space_dimension();
01887 assert(cg_space_dim <= space_dimension());
01888
01889 if (cg.is_proper_congruence()) {
01890
01891
01892
01893 if (cg.is_inconsistent())
01894 set_empty();
01895 return;
01896 }
01897
01898 assert(cg.is_equality());
01899 dimension_type cg_num_vars = 0;
01900 dimension_type cg_only_var = 0;
01901
01902 if (!extract_interval_congruence(cg, cg_space_dim, cg_num_vars, cg_only_var))
01903 return;
01904
01905 if (cg_num_vars == 0) {
01906
01907 if (cg.inhomogeneous_term() != 0)
01908 set_empty();
01909 return;
01910 }
01911
01912 assert(cg_num_vars == 1);
01913 const Coefficient& n = cg.inhomogeneous_term();
01914 const Coefficient& d = cg.coefficient(Variable(cg_only_var));
01915 add_interval_constraint_no_check(cg_only_var, Constraint::EQUALITY, n, d);
01916 }
01917
01918 template <typename ITV>
01919 void
01920 Box<ITV>::refine_no_check(const Congruence_System& cgs) {
01921 assert(cgs.space_dimension() <= space_dimension());
01922 for (Congruence_System::const_iterator i = cgs.begin(),
01923 cgs_end = cgs.end(); !marked_empty() && i != cgs_end; ++i)
01924 refine_no_check(*i);
01925 assert(OK());
01926 }
01927
01928 #if 1 // Alternative implementations for propagate_constraint_no_check.
01929 namespace {
01930
01931 inline bool
01932 propagate_constraint_check_result(Result r, Ternary& open) {
01933 switch (r) {
01934 case V_NEG_OVERFLOW:
01935 case V_POS_OVERFLOW:
01936 case V_UNKNOWN_NEG_OVERFLOW:
01937 case V_UNKNOWN_POS_OVERFLOW:
01938 return true;
01939 case V_LT:
01940 case V_GT:
01941 open = T_YES;
01942 return false;
01943 case V_LE:
01944 case V_GE:
01945 if (open == T_NO)
01946 open = T_MAYBE;
01947 return false;
01948 case V_EQ:
01949 return false;
01950 default:
01951 assert(false);
01952 return true;
01953 }
01954 }
01955
01956 }
01957
01958 template <typename ITV>
01959 void
01960 Box<ITV>::propagate_constraint_no_check(const Constraint& c) {
01961 assert(c.space_dimension() <= space_dimension());
01962
01963 typedef
01964 typename Select_Temp_Boundary_Type<typename ITV::boundary_type>::type
01965 Temp_Boundary_Type;
01966
01967 const dimension_type c_space_dim = c.space_dimension();
01968 const Constraint::Type c_type = c.type();
01969 const Coefficient& c_inhomogeneous_term = c.inhomogeneous_term();
01970
01971
01972 dimension_type last_k = c_space_dim;
01973 for (dimension_type k = c_space_dim; k-- > 0; ) {
01974 if (c.coefficient(Variable(k)) != 0) {
01975 last_k = k;
01976 break;
01977 }
01978 }
01979 if (last_k == c_space_dim) {
01980
01981 if (c_inhomogeneous_term < 0
01982 || (c_inhomogeneous_term == 0
01983 && c_type != Constraint::NONSTRICT_INEQUALITY))
01984 set_empty();
01985 return;
01986 }
01987
01988
01989 assert(last_k < c_space_dim);
01990 Result r;
01991 Temp_Boundary_Type t_bound;
01992 Temp_Boundary_Type t_a;
01993 Temp_Boundary_Type t_x;
01994 Ternary open;
01995 for (dimension_type k = last_k+1; k-- > 0; ) {
01996 const Coefficient& a_k = c.coefficient(Variable(k));
01997 int sgn_a_k = sgn(a_k);
01998 if (sgn_a_k == 0)
01999 continue;
02000 if (sgn_a_k > 0) {
02001 open = (c_type == Constraint::STRICT_INEQUALITY) ? T_YES : T_NO;
02002 if (open == T_NO)
02003 maybe_reset_fpu_inexact<Temp_Boundary_Type>();
02004 r = assign_r(t_bound, c_inhomogeneous_term, ROUND_UP);
02005 if (propagate_constraint_check_result(r, open))
02006 goto maybe_refine_upper_1;
02007 r = neg_assign_r(t_bound, t_bound, ROUND_DOWN);
02008 if (propagate_constraint_check_result(r, open))
02009 goto maybe_refine_upper_1;
02010 for (dimension_type i = last_k+1; i-- > 0; ) {
02011 if (i == k)
02012 continue;
02013 const Coefficient& a_i = c.coefficient(Variable(i));
02014 int sgn_a_i = sgn(a_i);
02015 if (sgn_a_i == 0)
02016 continue;
02017 ITV& x_i = seq[i];
02018 if (sgn_a_i < 0) {
02019 if (x_i.lower_is_unbounded())
02020 goto maybe_refine_upper_1;
02021 r = assign_r(t_a, a_i, ROUND_DOWN);
02022 if (propagate_constraint_check_result(r, open))
02023 goto maybe_refine_upper_1;
02024 r = assign_r(t_x, x_i.lower(), ROUND_DOWN);
02025 if (propagate_constraint_check_result(r, open))
02026 goto maybe_refine_upper_1;
02027 if (x_i.lower_is_open())
02028 open = T_YES;
02029 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_DOWN);
02030 if (propagate_constraint_check_result(r, open))
02031 goto maybe_refine_upper_1;
02032 }
02033 else {
02034 assert(sgn_a_i > 0);
02035 if (x_i.upper_is_unbounded())
02036 goto maybe_refine_upper_1;
02037 r = assign_r(t_a, a_i, ROUND_UP);
02038 if (propagate_constraint_check_result(r, open))
02039 goto maybe_refine_upper_1;
02040 r = assign_r(t_x, x_i.upper(), ROUND_UP);
02041 if (propagate_constraint_check_result(r, open))
02042 goto maybe_refine_upper_1;
02043 if (x_i.upper_is_open())
02044 open = T_YES;
02045 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_DOWN);
02046 if (propagate_constraint_check_result(r, open))
02047 goto maybe_refine_upper_1;
02048 }
02049 }
02050 r = assign_r(t_a, a_k, ROUND_UP);
02051 if (propagate_constraint_check_result(r, open))
02052 goto maybe_refine_upper_1;
02053 r = div_assign_r(t_bound, t_bound, t_a, ROUND_DOWN);
02054 if (propagate_constraint_check_result(r, open))
02055 goto maybe_refine_upper_1;
02056
02057
02058 if (open == T_MAYBE
02059 && maybe_check_fpu_inexact<Temp_Boundary_Type>() == 1)
02060 open = T_YES;
02061 seq[k].lower_narrow(t_bound, open == T_YES);
02062 reset_empty_up_to_date();
02063 maybe_refine_upper_1:
02064 if (c_type != Constraint::EQUALITY)
02065 continue;
02066 open = T_NO;
02067 maybe_reset_fpu_inexact<Temp_Boundary_Type>();
02068 r = assign_r(t_bound, c_inhomogeneous_term, ROUND_DOWN);
02069 if (propagate_constraint_check_result(r, open))
02070 goto next_k;
02071 r = neg_assign_r(t_bound, t_bound, ROUND_UP);
02072 if (propagate_constraint_check_result(r, open))
02073 goto next_k;
02074 for (dimension_type i = c_space_dim; i-- > 0; ) {
02075 if (i == k)
02076 continue;
02077 const Coefficient& a_i = c.coefficient(Variable(i));
02078 int sgn_a_i = sgn(a_i);
02079 if (sgn_a_i == 0)
02080 continue;
02081 ITV& x_i = seq[i];
02082 if (sgn_a_i < 0) {
02083 if (x_i.upper_is_unbounded())
02084 goto next_k;
02085 r = assign_r(t_a, a_i, ROUND_UP);
02086 if (propagate_constraint_check_result(r, open))
02087 goto next_k;
02088 r = assign_r(t_x, x_i.upper(), ROUND_UP);
02089 if (propagate_constraint_check_result(r, open))
02090 goto next_k;
02091 if (x_i.upper_is_open())
02092 open = T_YES;
02093 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
02094 if (propagate_constraint_check_result(r, open))
02095 goto next_k;
02096 }
02097 else {
02098 assert(sgn_a_i > 0);
02099 if (x_i.lower_is_unbounded())
02100 goto next_k;
02101 r = assign_r(t_a, a_i, ROUND_DOWN);
02102 if (propagate_constraint_check_result(r, open))
02103 goto next_k;
02104 r = assign_r(t_x, x_i.lower(), ROUND_DOWN);
02105 if (propagate_constraint_check_result(r, open))
02106 goto next_k;
02107 if (x_i.lower_is_open())
02108 open = T_YES;
02109 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
02110 if (propagate_constraint_check_result(r, open))
02111 goto next_k;
02112 }
02113 }
02114 r = assign_r(t_a, a_k, ROUND_DOWN);
02115 if (propagate_constraint_check_result(r, open))
02116 goto next_k;
02117 r = div_assign_r(t_bound, t_bound, t_a, ROUND_UP);
02118 if (propagate_constraint_check_result(r, open))
02119 goto next_k;
02120
02121
02122 if (open == T_MAYBE
02123 && maybe_check_fpu_inexact<Temp_Boundary_Type>() == 1)
02124 open = T_YES;
02125 seq[k].upper_narrow(t_bound, open == T_YES);
02126 reset_empty_up_to_date();
02127 }
02128 else {
02129 assert(sgn_a_k < 0);
02130 open = (c_type == Constraint::STRICT_INEQUALITY) ? T_YES : T_NO;
02131 if (open == T_NO)
02132 maybe_reset_fpu_inexact<Temp_Boundary_Type>();
02133 r = assign_r(t_bound, c_inhomogeneous_term, ROUND_UP);
02134 if (propagate_constraint_check_result(r, open))
02135 goto maybe_refine_upper_2;
02136 r = neg_assign_r(t_bound, t_bound, ROUND_DOWN);
02137 if (propagate_constraint_check_result(r, open))
02138 goto maybe_refine_upper_2;
02139 for (dimension_type i = c_space_dim; i-- > 0; ) {
02140 if (i == k)
02141 continue;
02142 const Coefficient& a_i = c.coefficient(Variable(i));
02143 int sgn_a_i = sgn(a_i);
02144 if (sgn_a_i == 0)
02145 continue;
02146 ITV& x_i = seq[i];
02147 if (sgn_a_i < 0) {
02148 if (x_i.lower_is_unbounded())
02149 goto maybe_refine_upper_2;
02150 r = assign_r(t_a, a_i, ROUND_DOWN);
02151 if (propagate_constraint_check_result(r, open))
02152 goto maybe_refine_upper_2;
02153 r = assign_r(t_x, x_i.lower(), ROUND_DOWN);
02154 if (propagate_constraint_check_result(r, open))
02155 goto maybe_refine_upper_2;
02156 if (x_i.lower_is_open())
02157 open = T_YES;
02158 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
02159 if (propagate_constraint_check_result(r, open))
02160 goto maybe_refine_upper_2;
02161 }
02162 else {
02163 assert(sgn_a_i > 0);
02164 if (x_i.upper_is_unbounded())
02165 goto maybe_refine_upper_2;
02166 r = assign_r(t_a, a_i, ROUND_UP);
02167 if (propagate_constraint_check_result(r, open))
02168 goto maybe_refine_upper_2;
02169 r = assign_r(t_x, x_i.upper(), ROUND_UP);
02170 if (propagate_constraint_check_result(r, open))
02171 goto maybe_refine_upper_2;
02172 if (x_i.upper_is_open())
02173 open = T_YES;
02174 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
02175 if (propagate_constraint_check_result(r, open))
02176 goto maybe_refine_upper_2;
02177 }
02178 }
02179 r = assign_r(t_a, a_k, ROUND_UP);
02180 if (propagate_constraint_check_result(r, open))
02181 goto maybe_refine_upper_2;
02182 r = div_assign_r(t_bound, t_bound, t_a, ROUND_UP);
02183 if (propagate_constraint_check_result(r, open))
02184 goto maybe_refine_upper_2;
02185
02186
02187 if (open == T_MAYBE
02188 && maybe_check_fpu_inexact<Temp_Boundary_Type>() == 1)
02189 open = T_YES;
02190 seq[k].upper_narrow(t_bound, open == T_YES);
02191 reset_empty_up_to_date();
02192 maybe_refine_upper_2:
02193 if (c_type != Constraint::EQUALITY)
02194 continue;
02195 open = T_NO;
02196 maybe_reset_fpu_inexact<Temp_Boundary_Type>();
02197 r = assign_r(t_bound, c_inhomogeneous_term, ROUND_DOWN);
02198 if (propagate_constraint_check_result(r, open))
02199 goto next_k;
02200 r = neg_assign_r(t_bound, t_bound, ROUND_UP);
02201 if (propagate_constraint_check_result(r, open))
02202 goto next_k;
02203 for (dimension_type i = c_space_dim; i-- > 0; ) {
02204 if (i == k)
02205 continue;
02206 const Coefficient& a_i = c.coefficient(Variable(i));
02207 int sgn_a_i = sgn(a_i);
02208 if (sgn_a_i == 0)
02209 continue;
02210 ITV& x_i = seq[i];
02211 if (sgn_a_i < 0) {
02212 if (x_i.upper_is_unbounded())
02213 goto next_k;
02214 r = assign_r(t_a, a_i, ROUND_UP);
02215 if (propagate_constraint_check_result(r, open))
02216 goto next_k;
02217 r = assign_r(t_x, x_i.upper(), ROUND_UP);
02218 if (propagate_constraint_check_result(r, open))
02219 goto next_k;
02220 if (x_i.upper_is_open())
02221 open = T_YES;
02222 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
02223 if (propagate_constraint_check_result(r, open))
02224 goto next_k;
02225 }
02226 else {
02227 assert(sgn_a_i > 0);
02228 if (x_i.lower_is_unbounded())
02229 goto next_k;
02230 r = assign_r(t_a, a_i, ROUND_DOWN);
02231 if (propagate_constraint_check_result(r, open))
02232 goto next_k;
02233 r = assign_r(t_x, x_i.lower(), ROUND_DOWN);
02234 if (propagate_constraint_check_result(r, open))
02235 goto next_k;
02236 if (x_i.lower_is_open())
02237 open = T_YES;
02238 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
02239 if (propagate_constraint_check_result(r, open))
02240 goto next_k;
02241 }
02242 }
02243 r = assign_r(t_a, a_k, ROUND_DOWN);
02244 if (propagate_constraint_check_result(r, open))
02245 goto next_k;
02246 r = div_assign_r(t_bound, t_bound, t_a, ROUND_DOWN);
02247 if (propagate_constraint_check_result(r, open))
02248 goto next_k;
02249
02250
02251 if (open == T_MAYBE
02252 && maybe_check_fpu_inexact<Temp_Boundary_Type>() == 1)
02253 open = T_YES;
02254 seq[k].lower_narrow(t_bound, open == T_YES);
02255 reset_empty_up_to_date();
02256 }
02257 next_k:
02258 ;
02259 }
02260 }
02261
02262 #else // Alternative implementations for propagate_constraint_no_check.
02263
02264 template <typename ITV>
02265 void
02266 Box<ITV>::propagate_constraint_no_check(const Constraint& c) {
02267 assert(c.space_dimension() <= space_dimension());
02268
02269 dimension_type c_space_dim = c.space_dimension();
02270 ITV k[c_space_dim];
02271 ITV p[c_space_dim];
02272 for (dimension_type i = c_space_dim; i-- > 0; ) {
02273 k[i] = c.coefficient(Variable(i));
02274 ITV& p_i = p[i];
02275 p_i = seq[i];
02276 p_i.mul_assign(p_i, k[i]);
02277 }
02278 const Coefficient& inhomogeneous_term = c.inhomogeneous_term();
02279 for (dimension_type i = c_space_dim; i-- > 0; ) {
02280 int sgn_coefficient_i = sgn(c.coefficient(Variable(i)));
02281 if (sgn_coefficient_i == 0)
02282 continue;
02283 ITV q(inhomogeneous_term);
02284 for (dimension_type j = c_space_dim; j-- > 0; ) {
02285 if (i == j)
02286 continue;
02287 q.add_assign(q, p[j]);
02288 }
02289 q.div_assign(q, k[i]);
02290 q.neg_assign(q);
02291 Relation_Symbol rel;
02292 switch (c.type()) {
02293 case Constraint::EQUALITY:
02294 rel = EQUAL;
02295 break;
02296 case Constraint::NONSTRICT_INEQUALITY:
02297 rel = (sgn_coefficient_i > 0) ? GREATER_OR_EQUAL : LESS_OR_EQUAL;
02298 break;
02299 case Constraint::STRICT_INEQUALITY:
02300 rel = (sgn_coefficient_i > 0) ? GREATER_THAN : LESS_THAN;
02301 break;
02302 }
02303 seq[i].refine_existential(rel, q);
02304
02305
02306
02307
02308
02309 if (seq[i].is_empty()) {
02310 set_empty();
02311 break;
02312 }
02313 }
02314
02315 assert(OK());
02316 }
02317
02318 #endif // Alternative implementations for propagate_constraint_no_check.
02319
02320 template <typename ITV>
02321 void
02322 Box<ITV>::propagate_constraints_no_check(const Constraint_System& cs) {
02323 assert(cs.space_dimension() <= space_dimension());
02324
02325 bool changed;
02326 do {
02327 Sequence copy(seq);
02328 for (Constraint_System::const_iterator i = cs.begin(),
02329 cs_end = cs.end(); i != cs_end; ++i)
02330 propagate_constraint_no_check(*i);
02331
02332
02333
02334
02335 maybe_abandon();
02336
02337 changed = (copy != seq);
02338 } while (changed);
02339 }
02340
02341 template <typename ITV>
02342 void
02343 Box<ITV>::affine_image(const Variable var,
02344 const Linear_Expression& expr,
02345 Coefficient_traits::const_reference denominator) {
02346
02347 if (denominator == 0)
02348 throw_generic("affine_image(v, e, d)", "d == 0");
02349
02350
02351 const dimension_type space_dim = space_dimension();
02352 const dimension_type expr_space_dim = expr.space_dimension();
02353 if (space_dim < expr_space_dim)
02354 throw_dimension_incompatible("affine_image(v, e, d)", "e", expr);
02355
02356 const dimension_type var_space_dim = var.space_dimension();
02357 if (space_dim < var_space_dim)
02358 throw_dimension_incompatible("affine_image(v, e, d)", "v", var);
02359
02360 if (is_empty())
02361 return;
02362
02363 Tmp_Interval_Type expr_value, temp0, temp1;
02364 expr_value.assign(expr.inhomogeneous_term());
02365 for (dimension_type i = expr_space_dim; i-- > 0; ) {
02366 const Coefficient& coeff = expr.coefficient(Variable(i));
02367 if (coeff != 0) {
02368 temp0.assign(coeff);
02369 temp1.assign(seq[i]);
02370 temp0.mul_assign(temp0, temp1);
02371 expr_value.add_assign(expr_value, temp0);
02372 }
02373 }
02374 if (denominator != 1) {
02375 temp0.assign(denominator);
02376 expr_value.div_assign(expr_value, temp0);
02377 }
02378 seq[var.id()].assign(expr_value);
02379
02380 assert(OK());
02381 }
02382
02383 template <typename ITV>
02384 void
02385 Box<ITV>::affine_preimage(const Variable var,
02386 const Linear_Expression& expr,
02387 Coefficient_traits::const_reference
02388 denominator) {
02389
02390 if (denominator == 0)
02391 throw_generic("affine_preimage(v, e, d)", "d == 0");
02392
02393
02394 const dimension_type x_space_dim = space_dimension();
02395 const dimension_type expr_space_dim = expr.space_dimension();
02396 if (x_space_dim < expr_space_dim)
02397 throw_dimension_incompatible("affine_preimage(v, e, d)", "e", expr);
02398
02399 const dimension_type var_space_dim = var.space_dimension();
02400 if (x_space_dim < var_space_dim)
02401 throw_dimension_incompatible("affine_preimage(v, e, d)", "v", var);
02402
02403 if (is_empty())
02404 return;
02405
02406 const Coefficient& expr_v = expr.coefficient(var);
02407 const bool invertible = (expr_v != 0);
02408 if (!invertible) {
02409 Tmp_Interval_Type expr_value, temp0, temp1;
02410 expr_value.assign(expr.inhomogeneous_term());
02411 for (dimension_type i = expr_space_dim; i-- > 0; ) {
02412 const Coefficient& coeff = expr.coefficient(Variable(i));
02413 if (coeff != 0) {
02414 temp0.assign(coeff);
02415 temp1.assign(seq[i]);
02416 temp0.mul_assign(temp0, temp1);
02417 expr_value.add_assign(expr_value, temp0);
02418 }
02419 }
02420 if (denominator != 1) {
02421 temp0.assign(denominator);
02422 expr_value.div_assign(expr_value, temp0);
02423 }
02424 ITV& x_seq_v = seq[var.id()];
02425 expr_value.intersect_assign(x_seq_v);
02426 if (expr_value.is_empty())
02427 set_empty();
02428 else
02429 x_seq_v.assign(UNIVERSE);
02430 }
02431 else {
02432
02433
02434
02435
02436 Linear_Expression inverse;
02437 inverse -= expr;
02438 inverse += (expr_v + denominator) * var;
02439 affine_image(var, inverse, expr_v);
02440 }
02441 assert(OK());
02442 }
02443
02444 template <typename ITV>
02445 void
02446 Box<ITV>
02447 ::bounded_affine_image(const Variable var,
02448 const Linear_Expression& lb_expr,
02449 const Linear_Expression& ub_expr,
02450 Coefficient_traits::const_reference denominator) {
02451
02452 if (denominator == 0)
02453 throw_generic("bounded_affine_image(v, lb, ub, d)", "d == 0");
02454
02455
02456 const dimension_type space_dim = space_dimension();
02457
02458
02459 const dimension_type lb_space_dim = lb_expr.space_dimension();
02460 if (space_dim < lb_space_dim)
02461 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
02462 "lb", lb_expr);
02463 const dimension_type ub_space_dim = ub_expr.space_dimension();
02464 if (space_dim < ub_space_dim)
02465 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
02466 "ub", ub_expr);
02467
02468 const dimension_type var_space_dim = var.space_dimension();
02469 if (space_dim < var_space_dim)
02470 throw_dimension_incompatible("affine_image(v, e, d)", "v", var);
02471
02472
02473 if (is_empty())
02474 return;
02475
02476
02477 if (denominator > 0)
02478 refine_with_constraint(lb_expr <= ub_expr);
02479 else
02480 refine_with_constraint(lb_expr >= ub_expr);
02481
02482
02483 if (lb_expr.coefficient(var) == 0) {
02484
02485 generalized_affine_image(var,
02486 LESS_OR_EQUAL,
02487 ub_expr,
02488 denominator);
02489 if (denominator > 0)
02490 refine_with_constraint(lb_expr <= denominator*var);
02491 else
02492 refine_with_constraint(denominator*var <= lb_expr);
02493 }
02494 else if (ub_expr.coefficient(var) == 0) {
02495
02496 generalized_affine_image(var,
02497 GREATER_OR_EQUAL,
02498 lb_expr,
02499 denominator);
02500 if (denominator > 0)
02501 refine_with_constraint(denominator*var <= ub_expr);
02502 else
02503 refine_with_constraint(ub_expr <= denominator*var);
02504 }
02505 else {
02506
02507
02508
02509
02510 DIRTY_TEMP(Coefficient, max_num);
02511 DIRTY_TEMP(Coefficient, max_den);
02512 bool max_included;
02513 DIRTY_TEMP(Coefficient, min_num);
02514 DIRTY_TEMP(Coefficient, min_den);
02515 bool min_included;
02516 ITV& seq_v = seq[var.id()];
02517 if (maximize(ub_expr, max_num, max_den, max_included)) {
02518 if (minimize(lb_expr, min_num, min_den, min_included)) {
02519
02520
02521
02522 min_den *= denominator;
02523 DIRTY_TEMP0(mpq_class, q);
02524 assign_r(q.get_num(), min_num, ROUND_NOT_NEEDED);
02525 assign_r(q.get_den(), min_den, ROUND_NOT_NEEDED);
02526 q.canonicalize();
02527 (denominator > 0)
02528 ? seq_v.lower_set(q, !min_included)
02529 : seq_v.upper_set(q, !min_included);
02530
02531
02532 max_den *= denominator;
02533 assign_r(q.get_num(), max_num, ROUND_NOT_NEEDED);
02534 assign_r(q.get_den(), max_den, ROUND_NOT_NEEDED);
02535 q.canonicalize();
02536 (denominator > 0)
02537 ? seq_v.upper_set(q, !max_included)
02538 : seq_v.lower_set(q, !max_included);
02539 }
02540 else {
02541
02542
02543
02544 DIRTY_TEMP0(mpq_class, q);
02545 max_den *= denominator;
02546 assign_r(q.get_num(), max_num, ROUND_NOT_NEEDED);
02547 assign_r(q.get_den(), max_den, ROUND_NOT_NEEDED);
02548 q.canonicalize();
02549 if (denominator > 0) {
02550 seq_v.lower_set(UNBOUNDED);
02551 seq_v.upper_set(q, !max_included);
02552 }
02553 else {
02554 seq_v.upper_set(UNBOUNDED);
02555 seq_v.lower_set(q, !max_included);
02556 }
02557 }
02558 }
02559 else if (minimize(lb_expr, min_num, min_den, min_included)) {
02560
02561
02562
02563 min_den *= denominator;
02564 DIRTY_TEMP0(mpq_class, q);
02565 assign_r(q.get_num(), min_num, ROUND_NOT_NEEDED);
02566 assign_r(q.get_den(), min_den, ROUND_NOT_NEEDED);
02567 q.canonicalize();
02568 if (denominator > 0) {
02569 seq_v.upper_set(UNBOUNDED);
02570 seq_v.lower_set(q, !min_included);
02571 }
02572 else {
02573 seq_v.lower_set(UNBOUNDED);
02574 seq_v.upper_set(q, !min_included);
02575 }
02576 }
02577 else {
02578
02579
02580
02581 seq_v.upper_set(UNBOUNDED);
02582 seq_v.lower_set(UNBOUNDED);
02583 }
02584 }
02585 assert(OK());
02586 }
02587
02588 template <typename ITV>
02589 void
02590 Box<ITV>
02591 ::bounded_affine_preimage(const Variable var,
02592 const Linear_Expression& lb_expr,
02593 const Linear_Expression& ub_expr,
02594 Coefficient_traits::const_reference denominator) {
02595
02596 const dimension_type space_dim = space_dimension();
02597 if (denominator == 0)
02598 throw_generic("bounded_affine_preimage(v, lb, ub, d)", "d == 0");
02599
02600
02601
02602 const dimension_type var_space_dim = var.space_dimension();
02603 if (space_dim < var_space_dim)
02604 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
02605 "v", var);
02606
02607
02608 const dimension_type lb_space_dim = lb_expr.space_dimension();
02609 if (space_dim < lb_space_dim)
02610 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub)",
02611 "lb", lb_expr);
02612 const dimension_type ub_space_dim = ub_expr.space_dimension();
02613 if (space_dim < ub_space_dim)
02614 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub)",
02615 "ub", ub_expr);
02616
02617
02618 if (marked_empty())
02619 return;
02620
02621 const bool negative_denom = (denominator < 0);
02622 const Coefficient& lb_var_coeff = lb_expr.coefficient(var);
02623 const Coefficient& ub_var_coeff = ub_expr.coefficient(var);
02624
02625
02626
02627 if (lb_var_coeff == ub_var_coeff) {
02628 if (negative_denom)
02629 refine_with_constraint(lb_expr >= ub_expr);
02630 else
02631 refine_with_constraint(lb_expr <= ub_expr);
02632 }
02633
02634 ITV& seq_var = seq[var.id()];
02635 if (!seq_var.is_universe()) {
02636
02637
02638 TEMP_INTEGER(pos_denominator);
02639 pos_denominator = denominator;
02640 if (negative_denom)
02641 neg_assign(pos_denominator, pos_denominator);
02642
02643
02644 bool open_lower = seq_var.lower_is_open();
02645 bool unbounded_lower = seq_var.lower_is_unbounded();
02646 DIRTY_TEMP0(mpq_class, q_seq_var_lower);
02647 DIRTY_TEMP(Coefficient, num_lower);
02648 DIRTY_TEMP(Coefficient, den_lower);
02649 if (!unbounded_lower) {
02650 assign_r(q_seq_var_lower, seq_var.lower(), ROUND_NOT_NEEDED);
02651 assign_r(num_lower, q_seq_var_lower.get_num(), ROUND_NOT_NEEDED);
02652 assign_r(den_lower, q_seq_var_lower.get_den(), ROUND_NOT_NEEDED);
02653 if (negative_denom)
02654 neg_assign(den_lower, den_lower);
02655 num_lower *= pos_denominator;
02656 seq_var.lower_set(UNBOUNDED);
02657 }
02658 bool open_upper = seq_var.upper_is_open();
02659 bool unbounded_upper = seq_var.upper_is_unbounded();
02660 DIRTY_TEMP0(mpq_class, q_seq_var_upper);
02661 DIRTY_TEMP(Coefficient, num_upper);
02662 DIRTY_TEMP(Coefficient, den_upper);
02663 if (!unbounded_upper) {
02664 assign_r(q_seq_var_upper, seq_var.upper(), ROUND_NOT_NEEDED);
02665 assign_r(num_upper, q_seq_var_upper.get_num(), ROUND_NOT_NEEDED);
02666 assign_r(den_upper, q_seq_var_upper.get_den(), ROUND_NOT_NEEDED);
02667 if (negative_denom)
02668 neg_assign(den_upper, den_upper);
02669 num_upper *= pos_denominator;
02670 seq_var.upper_set(UNBOUNDED);
02671 }
02672
02673 if (!unbounded_lower) {
02674
02675
02676
02677 Linear_Expression revised_lb_expr(ub_expr);
02678 revised_lb_expr -= ub_var_coeff * var;
02679 DIRTY_TEMP(Coefficient, d);
02680 neg_assign(d, den_lower);
02681 revised_lb_expr *= d;
02682 revised_lb_expr += num_lower;
02683
02684
02685
02686 bool included;
02687 DIRTY_TEMP(Coefficient, den);
02688 if (minimize(revised_lb_expr, num_lower, den, included)) {
02689 den_lower *= (den * ub_var_coeff);
02690 DIRTY_TEMP0(mpq_class, q);
02691 assign_r(q.get_num(), num_lower, ROUND_NOT_NEEDED);
02692 assign_r(q.get_den(), den_lower, ROUND_NOT_NEEDED);
02693 q.canonicalize();
02694 open_lower |= !included;
02695 if ((ub_var_coeff >= 0) ? !negative_denom : negative_denom)
02696 seq_var.lower_narrow(q, open_lower);
02697 else
02698 seq_var.upper_narrow(q, open_lower);
02699 if (seq_var.is_empty()) {
02700 set_empty();
02701 return;
02702 }
02703 }
02704 }
02705
02706 if (!unbounded_upper) {
02707
02708
02709
02710 Linear_Expression revised_ub_expr(lb_expr);
02711 revised_ub_expr -= lb_var_coeff * var;
02712 DIRTY_TEMP(Coefficient, d);
02713 neg_assign(d, den_upper);
02714 revised_ub_expr *= d;
02715 revised_ub_expr += num_upper;
02716
02717
02718
02719 bool included;
02720 DIRTY_TEMP(Coefficient, den);
02721 if (maximize(revised_ub_expr, num_upper, den, included)) {
02722 den_upper *= (den * lb_var_coeff);
02723 DIRTY_TEMP0(mpq_class, q);
02724 assign_r(q.get_num(), num_upper, ROUND_NOT_NEEDED);
02725 assign_r(q.get_den(), den_upper, ROUND_NOT_NEEDED);
02726 q.canonicalize();
02727 open_upper |= !included;
02728 if ((lb_var_coeff >= 0) ? !negative_denom : negative_denom)
02729 seq_var.upper_narrow(q, open_upper);
02730 else
02731 seq_var.lower_narrow(q, open_upper);
02732 if (seq_var.is_empty()) {
02733 set_empty();
02734 return;
02735 }
02736 }
02737 }
02738 }
02739
02740
02741
02742 if (lb_var_coeff != ub_var_coeff) {
02743 if (denominator > 0)
02744 refine_with_constraint(lb_expr <= ub_expr);
02745 else
02746 refine_with_constraint(lb_expr >= ub_expr);
02747 }
02748
02749 assert(OK());
02750 }
02751
02752 template <typename ITV>
02753 void
02754 Box<ITV>
02755 ::generalized_affine_image(const Variable var,
02756 const Relation_Symbol relsym,
02757 const Linear_Expression& expr,
02758 Coefficient_traits::const_reference denominator) {
02759
02760 if (denominator == 0)
02761 throw_generic("generalized_affine_image(v, r, e, d)", "d == 0");
02762
02763
02764 const dimension_type space_dim = space_dimension();
02765
02766
02767 if (space_dim < expr.space_dimension())
02768 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
02769 "e", expr);
02770
02771 const dimension_type var_space_dim = var.space_dimension();
02772 if (space_dim < var_space_dim)
02773 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
02774 "v", var);
02775
02776
02777 if (relsym == NOT_EQUAL)
02778 throw_generic("generalized_affine_image(v, r, e, d)",
02779 "r is the disequality relation symbol");
02780
02781
02782 affine_image(var, expr, denominator);
02783
02784 if (relsym == EQUAL)
02785
02786 return;
02787
02788
02789 if (is_empty())
02790 return;
02791
02792 ITV& seq_var = seq[var.id()];
02793 switch (relsym) {
02794 case LESS_OR_EQUAL:
02795 seq_var.lower_set(UNBOUNDED);
02796 break;
02797 case LESS_THAN:
02798 seq_var.lower_set(UNBOUNDED);
02799 if (!seq_var.upper_is_unbounded())
02800 seq_var.refine_existential(LESS_THAN, seq_var.upper());
02801 break;
02802 case GREATER_OR_EQUAL:
02803 seq_var.upper_set(UNBOUNDED);
02804 break;
02805 case GREATER_THAN:
02806 seq_var.upper_set(UNBOUNDED);
02807 if (!seq_var.lower_is_unbounded())
02808 seq_var.refine_existential(GREATER_THAN, seq_var.lower());
02809 break;
02810 default:
02811
02812 throw std::runtime_error("PPL internal error");
02813 }
02814 assert(OK());
02815 }
02816
02817 template <typename ITV>
02818 void
02819 Box<ITV>
02820 ::generalized_affine_preimage(const Variable var,
02821 const Relation_Symbol relsym,
02822 const Linear_Expression& expr,
02823 Coefficient_traits::const_reference denominator)
02824 {
02825
02826 if (denominator == 0)
02827 throw_generic("generalized_affine_preimage(v, r, e, d)",
02828 "d == 0");
02829
02830
02831 const dimension_type space_dim = space_dimension();
02832
02833
02834 if (space_dim < expr.space_dimension())
02835 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
02836 "e", expr);
02837
02838 const dimension_type var_space_dim = var.space_dimension();
02839 if (space_dim < var_space_dim)
02840 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
02841 "v", var);
02842
02843 if (relsym == NOT_EQUAL)
02844 throw_generic("generalized_affine_preimage(v, r, e, d)",
02845 "r is the disequality relation symbol");
02846
02847
02848 if (relsym == EQUAL) {
02849 affine_preimage(var, expr, denominator);
02850 return;
02851 }
02852
02853
02854 Relation_Symbol reversed_relsym;
02855 switch (relsym) {
02856 case LESS_THAN:
02857 reversed_relsym = GREATER_THAN;
02858 break;
02859 case LESS_OR_EQUAL:
02860 reversed_relsym = GREATER_OR_EQUAL;
02861 break;
02862 case GREATER_OR_EQUAL:
02863 reversed_relsym = LESS_OR_EQUAL;
02864 break;
02865 case GREATER_THAN:
02866 reversed_relsym = LESS_THAN;
02867 break;
02868 default:
02869
02870 throw std::runtime_error("PPL internal error");
02871 }
02872
02873
02874
02875 const Coefficient& var_coefficient = expr.coefficient(var);
02876 if (var_coefficient != 0) {
02877 Linear_Expression inverse_expr
02878 = expr - (denominator + var_coefficient) * var;
02879 TEMP_INTEGER(inverse_denominator);
02880 neg_assign(inverse_denominator, var_coefficient);
02881 Relation_Symbol inverse_relsym
02882 = (sgn(denominator) == sgn(inverse_denominator))
02883 ? relsym : reversed_relsym;
02884 generalized_affine_image(var, inverse_relsym, inverse_expr,
02885 inverse_denominator);
02886 return;
02887 }
02888
02889
02890
02891
02892
02893
02894
02895
02896 DIRTY_TEMP(Coefficient, max_num);
02897 DIRTY_TEMP(Coefficient, max_den);
02898 bool max_included;
02899 bool bound_above = maximize(denominator*var, max_num, max_den, max_included);
02900 DIRTY_TEMP(Coefficient, min_num);
02901 DIRTY_TEMP(Coefficient, min_den);
02902 bool min_included;
02903 bool bound_below = minimize(denominator*var, min_num, min_den, min_included);
02904
02905 const Relation_Symbol corrected_relsym
02906 = (denominator > 0) ? relsym : reversed_relsym;
02907
02908
02909 DIRTY_TEMP(Linear_Expression, revised_expr);
02910 dimension_type dim = space_dim;
02911 TEMP_INTEGER(d);
02912 if (corrected_relsym == LESS_THAN || corrected_relsym == LESS_OR_EQUAL) {
02913 if (bound_below) {
02914 for ( ; dim > 0; dim--) {
02915 d = min_den * expr.coefficient(Variable(dim - 1));
02916 revised_expr += d * Variable(dim - 1);
02917 }
02918 }
02919 }
02920 else {
02921 if (bound_above) {
02922 for ( ; dim > 0; dim--) {
02923 d = max_den * expr.coefficient(Variable(dim - 1));
02924 revised_expr += d * Variable(dim - 1);
02925 }
02926 }
02927 }
02928
02929 switch (corrected_relsym) {
02930 case LESS_THAN:
02931 if (bound_below)
02932 refine_with_constraint(min_num < revised_expr);
02933 break;
02934 case LESS_OR_EQUAL:
02935 if (bound_below)
02936 (min_included)
02937 ? refine_with_constraint(min_num <= revised_expr)
02938 : refine_with_constraint(min_num < revised_expr);
02939 break;
02940 case GREATER_OR_EQUAL:
02941 if (bound_above)
02942 (max_included)
02943 ? refine_with_constraint(max_num >= revised_expr)
02944 : refine_with_constraint(max_num > revised_expr);
02945 break;
02946 case GREATER_THAN:
02947 if (bound_above)
02948 refine_with_constraint(max_num > revised_expr);
02949 break;
02950 default:
02951
02952 throw std::runtime_error("PPL internal error");
02953 }
02954
02955 if (is_empty())
02956 return;
02957 ITV& seq_v = seq[var.id()];
02958 seq_v.lower_set(UNBOUNDED);
02959 seq_v.upper_set(UNBOUNDED);
02960 assert(OK());
02961 }
02962
02963 template <typename ITV>
02964 void
02965 Box<ITV>
02966 ::generalized_affine_image(const Linear_Expression& lhs,
02967 const Relation_Symbol relsym,
02968 const Linear_Expression& rhs) {
02969
02970
02971
02972 dimension_type lhs_space_dim = lhs.space_dimension();
02973 const dimension_type space_dim = space_dimension();
02974 if (space_dim < lhs_space_dim)
02975 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
02976 "e1", lhs);
02977
02978
02979 const dimension_type rhs_space_dim = rhs.space_dimension();
02980 if (space_dim < rhs_space_dim)
02981 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
02982 "e2", rhs);
02983
02984
02985 if (relsym == NOT_EQUAL)
02986 throw_generic("generalized_affine_image(e1, r, e2)",
02987 "r is the disequality relation symbol");
02988
02989
02990 if (marked_empty())
02991 return;
02992
02993
02994 DIRTY_TEMP(Coefficient, max_num);
02995 DIRTY_TEMP(Coefficient, max_den);
02996 bool max_included;
02997 bool max_rhs = maximize(rhs, max_num, max_den, max_included);
02998 DIRTY_TEMP(Coefficient, min_num);
02999 DIRTY_TEMP(Coefficient, min_den);
03000 bool min_included;
03001 bool min_rhs = minimize(rhs, min_num, min_den, min_included);
03002
03003
03004
03005
03006
03007 bool has_var = false;
03008 bool has_more_than_one_var = false;
03009
03010 dimension_type has_var_id = 0;
03011 for ( ; lhs_space_dim > 0; --lhs_space_dim)
03012 if (lhs.coefficient(Variable(lhs_space_dim - 1)) != 0) {
03013 if (has_var) {
03014 ITV& seq_i = seq[lhs_space_dim - 1];
03015 seq_i.lower_set(UNBOUNDED);
03016 seq_i.upper_set(UNBOUNDED);
03017 has_more_than_one_var = true;
03018 }
03019 else {
03020 has_var = true;
03021 has_var_id = lhs_space_dim - 1;
03022 }
03023 }
03024
03025 if (has_more_than_one_var) {
03026
03027
03028
03029
03030
03031 ITV& seq_var = seq[has_var_id];
03032 seq_var.lower_set(UNBOUNDED);
03033 seq_var.upper_set(UNBOUNDED);
03034 assert(OK());
03035 return;
03036 }
03037
03038 if (has_var) {
03039
03040 ITV& seq_var = seq[has_var_id];
03041
03042
03043
03044 const Coefficient& inhomo = lhs.inhomogeneous_term();
03045 const Coefficient& coeff = lhs.coefficient(Variable(has_var_id));
03046 DIRTY_TEMP0(mpq_class, q_max);
03047 DIRTY_TEMP0(mpq_class, q_min);
03048 if (max_rhs) {
03049 max_num -= inhomo * max_den;
03050 max_den *= coeff;
03051 assign_r(q_max.get_num(), max_num, ROUND_NOT_NEEDED);
03052 assign_r(q_max.get_den(), max_den, ROUND_NOT_NEEDED);
03053 q_max.canonicalize();
03054 }
03055 if (min_rhs) {
03056 min_num -= inhomo * min_den;
03057 min_den *= coeff;
03058 assign_r(q_min.get_num(), min_num, ROUND_NOT_NEEDED);
03059 assign_r(q_min.get_den(), min_den, ROUND_NOT_NEEDED);
03060 q_min.canonicalize();
03061 }
03062
03063
03064
03065 if (coeff > 0)
03066
03067 switch (relsym) {
03068 case LESS_OR_EQUAL:
03069 seq_var.lower_set(UNBOUNDED);
03070 max_rhs
03071 ? seq_var.upper_set(q_max, !max_included)
03072 : seq_var.upper_set(UNBOUNDED);
03073 break;
03074 case LESS_THAN:
03075 seq_var.lower_set(UNBOUNDED);
03076 max_rhs
03077 ? seq_var.upper_set(q_max, true)
03078 : seq_var.upper_set(UNBOUNDED);
03079 break;
03080 case EQUAL:
03081 max_rhs
03082 ? seq_var.upper_set(q_max, !max_included)
03083 : seq_var.upper_set(UNBOUNDED);
03084 min_rhs
03085 ? seq_var.lower_set(q_min, !min_included)
03086 : seq_var.lower_set(UNBOUNDED);
03087 break;
03088 case GREATER_OR_EQUAL:
03089 seq_var.upper_set(UNBOUNDED);
03090 min_rhs
03091 ? seq_var.lower_set(q_min, !min_included)
03092 : seq_var.lower_set(UNBOUNDED);
03093 break;
03094 case GREATER_THAN:
03095 seq_var.upper_set(UNBOUNDED);
03096 min_rhs
03097 ? seq_var.lower_set(q_min, true)
03098 : seq_var.lower_set(UNBOUNDED);
03099 break;
03100 default:
03101
03102 throw std::runtime_error("PPL internal error");
03103 }
03104 else
03105
03106 switch (relsym) {
03107 case GREATER_OR_EQUAL:
03108 seq_var.lower_set(UNBOUNDED);
03109 min_rhs
03110 ? seq_var.upper_set(q_min, !min_included)
03111 : seq_var.upper_set(UNBOUNDED);
03112 break;
03113 case GREATER_THAN:
03114 seq_var.lower_set(UNBOUNDED);
03115 min_rhs
03116 ? seq_var.upper_set(q_min, true)
03117 : seq_var.upper_set(UNBOUNDED);
03118 break;
03119 case EQUAL:
03120 max_rhs
03121 ? seq_var.lower_set(q_max, !max_included)
03122 : seq_var.lower_set(UNBOUNDED);
03123 min_rhs
03124 ? seq_var.upper_set(q_min, !min_included)
03125 : seq_var.upper_set(UNBOUNDED);
03126 break;
03127 case LESS_OR_EQUAL:
03128 seq_var.upper_set(UNBOUNDED);
03129 max_rhs
03130 ? seq_var.lower_set(q_max, !max_included)
03131 : seq_var.lower_set(UNBOUNDED);
03132 break;
03133 case LESS_THAN:
03134 seq_var.upper_set(UNBOUNDED);
03135 max_rhs
03136 ? seq_var.lower_set(q_max, true)
03137 : seq_var.lower_set(UNBOUNDED);
03138 break;
03139 default:
03140
03141 throw std::runtime_error("PPL internal error");
03142 }
03143 }
03144
03145 else {
03146
03147
03148 const Coefficient& inhomo = lhs.inhomogeneous_term();
03149 switch (relsym) {
03150 case LESS_THAN:
03151 refine_with_constraint(inhomo < rhs);
03152 break;
03153 case LESS_OR_EQUAL:
03154 refine_with_constraint(inhomo <= rhs);
03155 break;
03156 case EQUAL:
03157 refine_with_constraint(inhomo == rhs);
03158 break;
03159 case GREATER_OR_EQUAL:
03160 refine_with_constraint(inhomo >= rhs);
03161 break;
03162 case GREATER_THAN:
03163 refine_with_constraint(inhomo > rhs);
03164 break;
03165 default:
03166
03167 throw std::runtime_error("PPL internal error");
03168 }
03169 }
03170 assert(OK());
03171 }
03172
03173 template <typename ITV>
03174 void
03175 Box<ITV>::generalized_affine_preimage(const Linear_Expression& lhs,
03176 const Relation_Symbol relsym,
03177 const Linear_Expression& rhs) {
03178
03179
03180
03181 dimension_type lhs_space_dim = lhs.space_dimension();
03182 const dimension_type space_dim = space_dimension();
03183 if (space_dim < lhs_space_dim)
03184 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
03185 "e1", lhs);
03186
03187
03188 const dimension_type rhs_space_dim = rhs.space_dimension();
03189 if (space_dim < rhs_space_dim)
03190 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
03191 "e2", rhs);
03192
03193
03194 if (relsym == NOT_EQUAL)
03195 throw_generic("generalized_affine_image(e1, r, e2)",
03196 "r is the disequality relation symbol");
03197
03198
03199 if (marked_empty())
03200 return;
03201
03202
03203
03204
03205 Linear_Expression revised_lhs = lhs;
03206 Linear_Expression revised_rhs = rhs;
03207 for (dimension_type d = lhs_space_dim; d-- > 0; ) {
03208 const Variable& var = Variable(d);
03209 if (lhs.coefficient(var) != 0) {
03210 DIRTY_TEMP(Coefficient, temp);
03211 temp = rhs.coefficient(var) + lhs.coefficient(var);
03212 revised_rhs -= temp * var;
03213 revised_lhs -= temp * var;
03214 }
03215 }
03216 generalized_affine_image(revised_lhs, relsym, revised_rhs);
03217 assert(OK());
03218 }
03219
03220 template <typename ITV>
03221 template <typename Iterator>
03222 void
03223 Box<ITV>::CC76_widening_assign(const Box& y, Iterator first, Iterator last) {
03224 if (y.is_empty())
03225 return;
03226
03227 for (dimension_type i = seq.size(); i-- > 0; )
03228 seq[i].CC76_widening_assign(y.seq[i], first, last);
03229
03230 assert(OK());
03231 }
03232
03233 template <typename ITV>
03234 void
03235 Box<ITV>::CC76_widening_assign(const Box& y, unsigned* tp) {
03236 static typename ITV::boundary_type stop_points[] = {
03237 typename ITV::boundary_type(-2),
03238 typename ITV::boundary_type(-1),
03239 typename ITV::boundary_type(0),
03240 typename ITV::boundary_type(1),
03241 typename ITV::boundary_type(2)
03242 };
03243
03244 Box& x = *this;
03245
03246 if (tp != 0 && *tp > 0) {
03247 Box<ITV> x_tmp(x);
03248 x_tmp.CC76_widening_assign(y, 0);
03249
03250 if (!x.contains(x_tmp))
03251 --(*tp);
03252 return;
03253 }
03254 x.CC76_widening_assign(y,
03255 stop_points,
03256 stop_points
03257 + sizeof(stop_points)/sizeof(stop_points[0]));
03258 }
03259
03260 template <typename ITV>
03261 void
03262 Box<ITV>::limited_CC76_extrapolation_assign(const Box& y,
03263 const Constraint_System& cs,
03264 unsigned* tp) {
03265
03266 used(cs);
03267 Box& x = *this;
03268 x.CC76_widening_assign(y, tp);
03269 }
03270
03271 template <typename ITV>
03272 void
03273 Box<ITV>::CC76_narrowing_assign(const Box& y) {
03274 const dimension_type space_dim = space_dimension();
03275
03276
03277 if (space_dim != y.space_dimension())
03278 throw_dimension_incompatible("CC76_narrowing_assign(y)", y);
03279
03280 #ifndef NDEBUG
03281 {
03282
03283 const Box x_copy = *this;
03284 const Box y_copy = y;
03285 assert(y_copy.contains(x_copy));
03286 }
03287 #endif
03288
03289
03290
03291 if (space_dim == 0)
03292 return;
03293
03294
03295 if (y.is_empty())
03296 return;
03297
03298 if (is_empty())
03299 return;
03300
03301
03302
03303 for (dimension_type i = space_dim; i-- > 0; ) {
03304 ITV& x_i = seq[i];
03305 const ITV& y_i = y.seq[i];
03306 if (!x_i.lower_is_unbounded()
03307 && !y_i.lower_is_unbounded()
03308 && x_i.lower() != y_i.lower())
03309 x_i.lower() = y_i.lower();
03310 if (!x_i.upper_is_unbounded()
03311 && !y_i.upper_is_unbounded()
03312 && x_i.upper() != y_i.upper())
03313 x_i.upper() = y_i.upper();
03314 }
03315 assert(OK());
03316 }
03317
03318 template <typename ITV>
03319 Constraint_System
03320 Box<ITV>::constraints() const {
03321 Constraint_System cs;
03322 const dimension_type space_dim = space_dimension();
03323 if (space_dim == 0) {
03324 if (marked_empty())
03325 cs = Constraint_System::zero_dim_empty();
03326 }
03327 else if (marked_empty())
03328 cs.insert(0*Variable(space_dim-1) <= -1);
03329 else {
03330
03331
03332 cs.insert(0*Variable(space_dim-1) <= 0);
03333
03334 for (dimension_type k = 0; k < space_dim; ++k) {
03335 bool closed = false;
03336 DIRTY_TEMP(Coefficient, n);
03337 DIRTY_TEMP(Coefficient, d);
03338 if (get_lower_bound(k, closed, n, d)) {
03339 if (closed)
03340 cs.insert(d*Variable(k) >= n);
03341 else
03342 cs.insert(d*Variable(k) > n);
03343 }
03344 if (get_upper_bound(k, closed, n, d)) {
03345 if (closed)
03346 cs.insert(d*Variable(k) <= n);
03347 else
03348 cs.insert(d*Variable(k) < n);
03349 }
03350 }
03351 }
03352 return cs;
03353 }
03354
03355 template <typename ITV>
03356 Constraint_System
03357 Box<ITV>::minimized_constraints() const {
03358 Constraint_System cs;
03359 const dimension_type space_dim = space_dimension();
03360 if (space_dim == 0) {
03361 if (marked_empty())
03362 cs = Constraint_System::zero_dim_empty();
03363 }
03364
03365 else if (is_empty())
03366 cs.insert(0*Variable(space_dim-1) <= -1);
03367 else {
03368
03369
03370 cs.insert(0*Variable(space_dim-1) <= 0);
03371
03372 for (dimension_type k = 0; k < space_dim; ++k) {
03373 bool closed = false;
03374 DIRTY_TEMP(Coefficient, n);
03375 DIRTY_TEMP(Coefficient, d);
03376 if (get_lower_bound(k, closed, n, d)) {
03377 if (closed)
03378
03379 if (seq[k].is_singleton()) {
03380 cs.insert(d*Variable(k) == n);
03381 continue;
03382 }
03383 else
03384 cs.insert(d*Variable(k) >= n);
03385 else
03386 cs.insert(d*Variable(k) > n);
03387 }
03388 if (get_upper_bound(k, closed, n, d)) {
03389 if (closed)
03390 cs.insert(d*Variable(k) <= n);
03391 else
03392 cs.insert(d*Variable(k) < n);
03393 }
03394 }
03395 }
03396 return cs;
03397 }
03398
03399 template <typename ITV>
03400 Congruence_System
03401 Box<ITV>::congruences() const {
03402 Congruence_System cgs;
03403 const dimension_type space_dim = space_dimension();
03404 if (space_dim == 0) {
03405 if (marked_empty())
03406 cgs = Congruence_System::zero_dim_empty();
03407 }
03408
03409 else if (is_empty())
03410 cgs.insert((0*Variable(space_dim-1) %= -1) / 0);
03411 else {
03412
03413
03414 cgs.insert(0*Variable(space_dim-1) %= 0);
03415
03416 for (dimension_type k = 0; k < space_dim; ++k) {
03417 bool closed = false;
03418 DIRTY_TEMP(Coefficient, n);
03419 DIRTY_TEMP(Coefficient, d);
03420 if (get_lower_bound(k, closed, n, d) && closed)
03421
03422 if (seq[k].is_singleton())
03423 cgs.insert((d*Variable(k) %= n) / 0);
03424 }
03425 }
03426 return cgs;
03427 }
03428
03429 template <typename ITV>
03430 memory_size_type
03431 Box<ITV>::external_memory_in_bytes() const {
03432 memory_size_type n = seq.capacity() * sizeof(ITV);
03433 for (dimension_type k = seq.size(); k-- > 0; )
03434 n += seq[k].external_memory_in_bytes();
03435 return n;
03436 }
03437
03439 template <typename ITV>
03440 std::ostream&
03441 IO_Operators::operator<<(std::ostream& s, const Box<ITV>& box) {
03442 if (box.is_empty())
03443 s << "false";
03444 else if (box.is_universe())
03445 s << "true";
03446 else
03447 for (dimension_type k = 0,
03448 space_dim = box.space_dimension(); k < space_dim; ) {
03449 s << Variable(k) << " in " << box[k];
03450 ++k;
03451 if (k < space_dim)
03452 s << ", ";
03453 else
03454 break;
03455 }
03456 return s;
03457 }
03458
03459 template <typename ITV>
03460 void
03461 Box<ITV>::ascii_dump(std::ostream& s) const {
03462 const char separator = ' ';
03463 status.ascii_dump(s);
03464 const dimension_type space_dim = space_dimension();
03465 s << "space_dim" << separator << space_dim;
03466 s << "\n";
03467 for (dimension_type i = 0; i < space_dim; ++i)
03468 seq[i].ascii_dump(s);
03469 }
03470
03471 PPL_OUTPUT_TEMPLATE_DEFINITIONS(ITV, Box<ITV>)
03472
03473 template <typename ITV>
03474 bool
03475 Box<ITV>::ascii_load(std::istream& s) {
03476 if (!status.ascii_load(s))
03477 return false;
03478
03479 std::string str;
03480 dimension_type space_dim;
03481 if (!(s >> str) || str != "space_dim")
03482 return false;
03483 if (!(s >> space_dim))
03484 return false;
03485
03486 seq.clear();
03487 ITV seq_i;
03488 for (dimension_type i = 0; i < space_dim; ++i) {
03489 if (seq_i.ascii_load(s))
03490 seq.push_back(seq_i);
03491 else
03492 return false;
03493 }
03494
03495
03496 assert(OK());
03497 return true;
03498 }
03499
03500 template <typename ITV>
03501 void
03502 Box<ITV>::throw_dimension_incompatible(const char* method,
03503 const Box& y) const {
03504 std::ostringstream s;
03505 s << "PPL::Box::" << method << ":" << std::endl
03506 << "this->space_dimension() == " << this->space_dimension()
03507 << ", y->space_dimension() == " << y.space_dimension() << ".";
03508 throw std::invalid_argument(s.str());
03509 }
03510
03511 template <typename ITV>
03512 void
03513 Box<ITV>
03514 ::throw_dimension_incompatible(const char* method,
03515 dimension_type required_dim) const {
03516 std::ostringstream s;
03517 s << "PPL::Box::" << method << ":" << std::endl
03518 << "this->space_dimension() == " << space_dimension()
03519 << ", required dimension == " << required_dim << ".";
03520 throw std::invalid_argument(s.str());
03521 }
03522
03523 template <typename ITV>
03524 void
03525 Box<ITV>::throw_dimension_incompatible(const char* method,
03526 const Constraint& c) const {
03527 std::ostringstream s;
03528 s << "PPL::Box::" << method << ":" << std::endl
03529 << "this->space_dimension() == " << space_dimension()
03530 << ", c->space_dimension == " << c.space_dimension() << ".";
03531 throw std::invalid_argument(s.str());
03532 }
03533
03534 template <typename ITV>
03535 void
03536 Box<ITV>::throw_dimension_incompatible(const char* method,
03537 const Congruence& cg) const {
03538 std::ostringstream s;
03539 s << "PPL::Box::" << method << ":" << std::endl
03540 << "this->space_dimension() == " << space_dimension()
03541 << ", cg->space_dimension == " << cg.space_dimension() << ".";
03542 throw std::invalid_argument(s.str());
03543 }
03544
03545 template <typename ITV>
03546 void
03547 Box<ITV>::throw_dimension_incompatible(const char* method,
03548 const Constraint_System& cs) const {
03549 std::ostringstream s;
03550 s << "PPL::Box::" << method << ":" << std::endl
03551 << "this->space_dimension() == " << space_dimension()
03552 << ", cs->space_dimension == " << cs.space_dimension() << ".";
03553 throw std::invalid_argument(s.str());
03554 }
03555
03556 template <typename ITV>
03557 void
03558 Box<ITV>::throw_dimension_incompatible(const char* method,
03559 const Congruence_System& cgs) const {
03560 std::ostringstream s;
03561 s << "PPL::Box::" << method << ":" << std::endl
03562 << "this->space_dimension() == " << space_dimension()
03563 << ", cgs->space_dimension == " << cgs.space_dimension() << ".";
03564 throw std::invalid_argument(s.str());
03565 }
03566
03567 template <typename ITV>
03568 void
03569 Box<ITV>::throw_dimension_incompatible(const char* method,
03570 const Generator& g) const {
03571 std::ostringstream s;
03572 s << "PPL::Box::" << method << ":" << std::endl
03573 << "this->space_dimension() == " << space_dimension()
03574 << ", g->space_dimension == " << g.space_dimension() << ".";
03575 throw std::invalid_argument(s.str());
03576 }
03577
03578 template <typename ITV>
03579 void
03580 Box<ITV>::throw_constraint_incompatible(const char* method) {
03581 std::ostringstream s;
03582 s << "PPL::Box::" << method << ":" << std::endl
03583 << "the constraint is incompatible.";
03584 throw std::invalid_argument(s.str());
03585 }
03586
03587 template <typename ITV>
03588 void
03589 Box<ITV>::throw_expression_too_complex(const char* method,
03590 const Linear_Expression& e) {
03591 using namespace IO_Operators;
03592 std::ostringstream s;
03593 s << "PPL::Box::" << method << ":" << std::endl
03594 << e << " is too complex.";
03595 throw std::invalid_argument(s.str());
03596 }
03597
03598 template <typename ITV>
03599 void
03600 Box<ITV>::throw_dimension_incompatible(const char* method,
03601 const char* name_row,
03602 const Linear_Expression& e) const {
03603 std::ostringstream s;
03604 s << "PPL::Box::" << method << ":" << std::endl
03605 << "this->space_dimension() == " << space_dimension()
03606 << ", " << name_row << "->space_dimension() == "
03607 << e.space_dimension() << ".";
03608 throw std::invalid_argument(s.str());
03609 }
03610
03611 template <typename ITV>
03612 void
03613 Box<ITV>::throw_generic(const char* method, const char* reason) {
03614 std::ostringstream s;
03615 s << "PPL::Box::" << method << ":" << std::endl
03616 << reason;
03617 throw std::invalid_argument(s.str());
03618 }
03619
03620 template <typename ITV>
03621 void
03622 Box<ITV>::throw_space_dimension_overflow(const char* method,
03623 const char* reason) {
03624 std::ostringstream s;
03625 s << "PPL::Box::" << method << ":" << std::endl
03626 << reason;
03627 throw std::length_error(s.str());
03628 }
03629
03630 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
03631
03632 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
03633 template <typename Specialization,
03634 typename Temp, typename To, typename ITV>
03635 bool
03636 l_m_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
03637 const Box<ITV>& x, const Box<ITV>& y,
03638 const Rounding_Dir dir,
03639 Temp& tmp0, Temp& tmp1, Temp& tmp2) {
03640 const dimension_type x_space_dim = x.space_dimension();
03641
03642 if (x_space_dim != y.space_dimension())
03643 return false;
03644
03645
03646 if (x_space_dim == 0) {
03647 if (x.marked_empty() == y.marked_empty())
03648 assign_r(r, 0, ROUND_NOT_NEEDED);
03649 else
03650 assign_r(r, PLUS_INFINITY, ROUND_NOT_NEEDED);
03651 return true;
03652 }
03653
03654
03655 (void) x.is_empty();
03656 (void) y.is_empty();
03657
03658
03659 if (x.marked_empty() || y.marked_empty()) {
03660 if (x.marked_empty() == y.marked_empty()) {
03661 assign_r(r, 0, ROUND_NOT_NEEDED);
03662 return true;
03663 }
03664 else
03665 goto pinf;
03666 }
03667
03668 assign_r(tmp0, 0, ROUND_NOT_NEEDED);
03669 for (dimension_type i = x_space_dim; i-- > 0; ) {
03670 const ITV& x_i = x.seq[i];
03671 const ITV& y_i = y.seq[i];
03672
03673 if (x_i.lower_is_unbounded()) {
03674 if (!y_i.lower_is_unbounded())
03675 goto pinf;
03676 }
03677 else if (y_i.lower_is_unbounded())
03678 goto pinf;
03679 else {
03680 const Temp* tmp1p;
03681 const Temp* tmp2p;
03682 if (x_i.lower() > y_i.lower()) {
03683 maybe_assign(tmp1p, tmp1, x_i.lower(), dir);
03684 maybe_assign(tmp2p, tmp2, y_i.lower(), inverse(dir));
03685 }
03686 else {
03687 maybe_assign(tmp1p, tmp1, y_i.lower(), dir);
03688 maybe_assign(tmp2p, tmp2, x_i.lower(), inverse(dir));
03689 }
03690 sub_assign_r(tmp1, *tmp1p, *tmp2p, dir);
03691 assert(sgn(tmp1) >= 0);
03692 Specialization::combine(tmp0, tmp1, dir);
03693 }
03694
03695 if (x_i.upper_is_unbounded())
03696 if (y_i.upper_is_unbounded())
03697 continue;
03698 else
03699 goto pinf;
03700 else if (y_i.upper_is_unbounded())
03701 goto pinf;
03702 else {
03703 const Temp* tmp1p;
03704 const Temp* tmp2p;
03705 if (x_i.upper() > y_i.upper()) {
03706 maybe_assign(tmp1p, tmp1, x_i.upper(), dir);
03707 maybe_assign(tmp2p, tmp2, y_i.upper(), inverse(dir));
03708 }
03709 else {
03710 maybe_assign(tmp1p, tmp1, y_i.upper(), dir);
03711 maybe_assign(tmp2p, tmp2, x_i.upper(), inverse(dir));
03712 }
03713 sub_assign_r(tmp1, *tmp1p, *tmp2p, dir);
03714 assert(sgn(tmp1) >= 0);
03715 Specialization::combine(tmp0, tmp1, dir);
03716 }
03717 }
03718 Specialization::finalize(tmp0, dir);
03719 assign_r(r, tmp0, dir);
03720 return true;
03721
03722 pinf:
03723 assign_r(r, PLUS_INFINITY, ROUND_NOT_NEEDED);
03724 return true;
03725 }
03726
03727 }
03728
03729 #endif // !defined(PPL_Box_templates_hh)