00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include <ppl-config.h>
00025
00026 #include "Polyhedron.defs.hh"
00027 #include "Scalar_Products.defs.hh"
00028 #include <cassert>
00029 #include <string>
00030 #include <iostream>
00031 #include <sstream>
00032 #include <stdexcept>
00033
00034 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00035
00044 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
00045 #define BE_LAZY 1
00046
00047 namespace PPL = Parma_Polyhedra_Library;
00048
00049 PPL::Polyhedron::Polyhedron(const Topology topol,
00050 const dimension_type num_dimensions,
00051 const Degenerate_Element kind)
00052 : con_sys(topol),
00053 gen_sys(topol),
00054 sat_c(),
00055 sat_g() {
00056
00057 assert(num_dimensions <= max_space_dimension());
00058
00059 if (kind == EMPTY)
00060 status.set_empty();
00061 else if (num_dimensions > 0) {
00062 con_sys.add_low_level_constraints();
00063 con_sys.adjust_topology_and_space_dimension(topol, num_dimensions);
00064 set_constraints_minimized();
00065 }
00066 space_dim = num_dimensions;
00067 assert(OK());
00068 }
00069
00070 PPL::Polyhedron::Polyhedron(const Polyhedron& y, Complexity_Class)
00071 : con_sys(y.topology()),
00072 gen_sys(y.topology()),
00073 status(y.status),
00074 space_dim(y.space_dim) {
00075
00076 assert(topology() == y.topology());
00077 if (y.constraints_are_up_to_date())
00078 con_sys.assign_with_pending(y.con_sys);
00079 if (y.generators_are_up_to_date())
00080 gen_sys.assign_with_pending(y.gen_sys);
00081 if (y.sat_c_is_up_to_date())
00082 sat_c = y.sat_c;
00083 if (y.sat_g_is_up_to_date())
00084 sat_g = y.sat_g;
00085 }
00086
00087 PPL::Polyhedron::Polyhedron(const Topology topol, const Constraint_System& ccs)
00088 : con_sys(topol),
00089 gen_sys(topol),
00090 sat_c(),
00091 sat_g() {
00092
00093 assert(ccs.space_dimension() <= max_space_dimension());
00094
00095
00096 Constraint_System cs = ccs;
00097
00098
00099 const dimension_type cs_space_dim = cs.space_dimension();
00100 if (!cs.adjust_topology_and_space_dimension(topol, cs_space_dim))
00101 throw_topology_incompatible((topol == NECESSARILY_CLOSED)
00102 ? "C_Polyhedron(cs)"
00103 : "NNC_Polyhedron(cs)", "cs", cs);
00104
00105
00106 space_dim = cs_space_dim;
00107
00108 if (space_dim > 0) {
00109
00110 std::swap(con_sys, cs);
00111 if (con_sys.num_pending_rows() > 0) {
00112
00113
00114
00115
00116 con_sys.unset_pending_rows();
00117 con_sys.set_sorted(false);
00118 }
00119 con_sys.add_low_level_constraints();
00120 set_constraints_up_to_date();
00121 }
00122 else {
00123
00124 if (cs.num_columns() > 0)
00125
00126 for (dimension_type i = cs.num_rows(); i-- > 0; )
00127 if (cs[i].is_inconsistent()) {
00128
00129 set_empty();
00130 break;
00131 }
00132 }
00133 assert(OK());
00134 }
00135
00136 PPL::Polyhedron::Polyhedron(const Topology topol,
00137 Constraint_System& cs,
00138 Recycle_Input)
00139 : con_sys(topol),
00140 gen_sys(topol),
00141 sat_c(),
00142 sat_g() {
00143
00144 assert(cs.space_dimension() <= max_space_dimension());
00145
00146
00147 const dimension_type cs_space_dim = cs.space_dimension();
00148 if (!cs.adjust_topology_and_space_dimension(topol, cs_space_dim))
00149 throw_topology_incompatible((topol == NECESSARILY_CLOSED)
00150 ? "C_Polyhedron(cs, recycle)"
00151 : "NNC_Polyhedron(cs, recycle)", "cs", cs);
00152
00153
00154 space_dim = cs_space_dim;
00155
00156 if (space_dim > 0) {
00157
00158 std::swap(con_sys, cs);
00159 if (con_sys.num_pending_rows() > 0) {
00160
00161
00162
00163
00164 con_sys.unset_pending_rows();
00165 con_sys.set_sorted(false);
00166 }
00167 con_sys.add_low_level_constraints();
00168 set_constraints_up_to_date();
00169 }
00170 else {
00171
00172 if (cs.num_columns() > 0)
00173
00174 for (dimension_type i = cs.num_rows(); i-- > 0; )
00175 if (cs[i].is_inconsistent()) {
00176
00177 set_empty();
00178 break;
00179 }
00180 }
00181 assert(OK());
00182 }
00183
00184 PPL::Polyhedron::Polyhedron(const Topology topol, const Generator_System& cgs)
00185 : con_sys(topol),
00186 gen_sys(topol),
00187 sat_c(),
00188 sat_g() {
00189
00190 assert(cgs.space_dimension() <= max_space_dimension());
00191
00192
00193 Generator_System gs = cgs;
00194
00195
00196 if (gs.has_no_rows()) {
00197 space_dim = gs.space_dimension();
00198 status.set_empty();
00199 assert(OK());
00200 return;
00201 }
00202
00203
00204 if (!gs.has_points())
00205 throw_invalid_generators((topol == NECESSARILY_CLOSED)
00206 ? "C_Polyhedron(gs)"
00207 : "NNC_Polyhedron(gs)", "gs");
00208
00209 const dimension_type gs_space_dim = gs.space_dimension();
00210
00211 if (!gs.adjust_topology_and_space_dimension(topol, gs_space_dim))
00212 throw_topology_incompatible((topol == NECESSARILY_CLOSED)
00213 ? "C_Polyhedron(gs)"
00214 : "NNC_Polyhedron(gs)", "gs", gs);
00215
00216 if (gs_space_dim > 0) {
00217
00218 std::swap(gen_sys, gs);
00219
00220
00221 if (topol == NOT_NECESSARILY_CLOSED)
00222 gen_sys.add_corresponding_closure_points();
00223 if (gen_sys.num_pending_rows() > 0) {
00224
00225
00226
00227
00228 gen_sys.unset_pending_rows();
00229 gen_sys.set_sorted(false);
00230 }
00231
00232 set_generators_up_to_date();
00233
00234
00235 space_dim = gs_space_dim;
00236 assert(OK());
00237 return;
00238 }
00239
00240
00241
00242
00243 space_dim = 0;
00244 assert(OK());
00245 }
00246
00247 PPL::Polyhedron::Polyhedron(const Topology topol,
00248 Generator_System& gs,
00249 Recycle_Input)
00250 : con_sys(topol),
00251 gen_sys(topol),
00252 sat_c(),
00253 sat_g() {
00254
00255 assert(gs.space_dimension() <= max_space_dimension());
00256
00257
00258 if (gs.has_no_rows()) {
00259 space_dim = gs.space_dimension();
00260 status.set_empty();
00261 assert(OK());
00262 return;
00263 }
00264
00265
00266 if (!gs.has_points())
00267 throw_invalid_generators((topol == NECESSARILY_CLOSED)
00268 ? "C_Polyhedron(gs, recycle)"
00269 : "NNC_Polyhedron(gs, recycle)", "gs");
00270
00271 const dimension_type gs_space_dim = gs.space_dimension();
00272
00273 if (!gs.adjust_topology_and_space_dimension(topol, gs_space_dim))
00274 throw_topology_incompatible((topol == NECESSARILY_CLOSED)
00275 ? "C_Polyhedron(gs, recycle)"
00276 : "NNC_Polyhedron(gs, recycle)", "gs", gs);
00277
00278 if (gs_space_dim > 0) {
00279
00280 std::swap(gen_sys, gs);
00281
00282
00283 if (topol == NOT_NECESSARILY_CLOSED)
00284 gen_sys.add_corresponding_closure_points();
00285 if (gen_sys.num_pending_rows() > 0) {
00286
00287
00288
00289
00290 gen_sys.unset_pending_rows();
00291 gen_sys.set_sorted(false);
00292 }
00293
00294 set_generators_up_to_date();
00295
00296
00297 space_dim = gs_space_dim;
00298 assert(OK());
00299 return;
00300 }
00301
00302
00303
00304
00305 space_dim = 0;
00306 assert(OK());
00307 }
00308
00309 PPL::Polyhedron&
00310 PPL::Polyhedron::operator=(const Polyhedron& y) {
00311
00312 assert(topology() == y.topology());
00313 space_dim = y.space_dim;
00314 if (y.marked_empty())
00315 set_empty();
00316 else if (space_dim == 0)
00317 set_zero_dim_univ();
00318 else {
00319 status = y.status;
00320 if (y.constraints_are_up_to_date())
00321 con_sys.assign_with_pending(y.con_sys);
00322 if (y.generators_are_up_to_date())
00323 gen_sys.assign_with_pending(y.gen_sys);
00324 if (y.sat_c_is_up_to_date())
00325 sat_c = y.sat_c;
00326 if (y.sat_g_is_up_to_date())
00327 sat_g = y.sat_g;
00328 }
00329 return *this;
00330 }
00331
00332 PPL::Polyhedron::Three_Valued_Boolean
00333 PPL::Polyhedron::quick_equivalence_test(const Polyhedron& y) const {
00334
00335 assert(topology() == y.topology());
00336 assert(space_dim == y.space_dim);
00337 assert(!marked_empty() && !y.marked_empty() && space_dim > 0);
00338
00339 const Polyhedron& x = *this;
00340
00341 if (x.is_necessarily_closed()) {
00342 if (!x.has_something_pending() && !y.has_something_pending()) {
00343 bool css_normalized = false;
00344 if (x.constraints_are_minimized() && y.constraints_are_minimized()) {
00345
00346
00347 if (x.con_sys.num_rows() != y.con_sys.num_rows())
00348 return Polyhedron::TVB_FALSE;
00349
00350 dimension_type x_num_equalities = x.con_sys.num_equalities();
00351 if (x_num_equalities != y.con_sys.num_equalities())
00352 return Polyhedron::TVB_FALSE;
00353
00354
00355 css_normalized = (x_num_equalities == 0);
00356 }
00357
00358 if (x.generators_are_minimized() && y.generators_are_minimized()) {
00359
00360
00361 if (x.gen_sys.num_rows() != y.gen_sys.num_rows())
00362 return Polyhedron::TVB_FALSE;
00363
00364 const dimension_type x_num_lines = x.gen_sys.num_lines();
00365 if (x_num_lines != y.gen_sys.num_lines())
00366 return Polyhedron::TVB_FALSE;
00367
00368 if (x_num_lines == 0) {
00369
00370 x.obtain_sorted_generators();
00371 y.obtain_sorted_generators();
00372 if (x.gen_sys == y.gen_sys)
00373 return Polyhedron::TVB_TRUE;
00374 else
00375 return Polyhedron::TVB_FALSE;
00376 }
00377 }
00378
00379 if (css_normalized) {
00380
00381 x.obtain_sorted_constraints();
00382 y.obtain_sorted_constraints();
00383 if (x.con_sys == y.con_sys)
00384 return Polyhedron::TVB_TRUE;
00385 else
00386 return Polyhedron::TVB_FALSE;
00387 }
00388 }
00389 }
00390 return Polyhedron::TVB_DONT_KNOW;
00391 }
00392
00393 bool
00394 PPL::Polyhedron::is_included_in(const Polyhedron& y) const {
00395
00396 assert(topology() == y.topology());
00397 assert(space_dim == y.space_dim);
00398 assert(!marked_empty() && !y.marked_empty() && space_dim > 0);
00399
00400 const Polyhedron& x = *this;
00401
00402
00403 if (x.has_pending_constraints() && !x.process_pending_constraints())
00404 return true;
00405
00406 if (y.has_pending_generators())
00407 y.process_pending_generators();
00408
00409 #if BE_LAZY
00410 if (!x.generators_are_up_to_date() && !x.update_generators())
00411 return true;
00412 if (!y.constraints_are_up_to_date())
00413 y.update_constraints();
00414 #else
00415 if (!x.generators_are_minimized())
00416 x.minimize();
00417 if (!y.constraints_are_minimized())
00418 y.minimize();
00419 #endif
00420
00421 assert(x.OK());
00422 assert(y.OK());
00423
00424 const Generator_System& gs = x.gen_sys;
00425 const Constraint_System& cs = y.con_sys;
00426
00427 if (x.is_necessarily_closed())
00428
00429
00430
00431
00432
00433
00434 for (dimension_type i = cs.num_rows(); i-- > 0; ) {
00435 const Constraint& c = cs[i];
00436 if (c.is_inequality()) {
00437 for (dimension_type j = gs.num_rows(); j-- > 0; ) {
00438 const Generator& g = gs[j];
00439 const int sp_sign = Scalar_Products::sign(c, g);
00440 if (g.is_line()) {
00441 if (sp_sign != 0)
00442 return false;
00443 }
00444 else
00445
00446 if (sp_sign < 0)
00447 return false;
00448 }
00449 }
00450 else {
00451
00452 for (dimension_type j = gs.num_rows(); j-- > 0; )
00453 if (Scalar_Products::sign(c, gs[j]) != 0)
00454 return false;
00455 }
00456 }
00457 else {
00458
00459
00460 for (dimension_type i = cs.num_rows(); i-- > 0; ) {
00461 const Constraint& c = cs[i];
00462 switch (c.type()) {
00463 case Constraint::NONSTRICT_INEQUALITY:
00464 for (dimension_type j = gs.num_rows(); j-- > 0; ) {
00465 const Generator& g = gs[j];
00466 const int sp_sign = Scalar_Products::reduced_sign(c, g);
00467 if (g.is_line()) {
00468 if (sp_sign != 0)
00469 return false;
00470 }
00471 else
00472
00473 if (sp_sign < 0)
00474 return false;
00475 }
00476 break;
00477 case Constraint::EQUALITY:
00478 for (dimension_type j = gs.num_rows(); j-- > 0; )
00479 if (Scalar_Products::reduced_sign(c, gs[j]) != 0)
00480 return false;
00481 break;
00482 case Constraint::STRICT_INEQUALITY:
00483 for (dimension_type j = gs.num_rows(); j-- > 0; ) {
00484 const Generator& g = gs[j];
00485 const int sp_sign = Scalar_Products::reduced_sign(c, g);
00486 switch (g.type()) {
00487 case Generator::POINT:
00488
00489
00490
00491 if (sp_sign <= 0)
00492 return false;
00493 break;
00494 case Generator::LINE:
00495
00496 if (sp_sign != 0)
00497 return false;
00498 break;
00499 case Generator::RAY:
00500
00501 case Generator::CLOSURE_POINT:
00502
00503 if (sp_sign < 0)
00504 return false;
00505 break;
00506 }
00507 }
00508 break;
00509 }
00510 }
00511 }
00512
00513
00514 return true;
00515 }
00516
00517 bool
00518 PPL::Polyhedron::bounds(const Linear_Expression& expr,
00519 const bool from_above) const {
00520
00521
00522 const dimension_type expr_space_dim = expr.space_dimension();
00523 if (space_dim < expr_space_dim)
00524 throw_dimension_incompatible((from_above
00525 ? "bounds_from_above(e)"
00526 : "bounds_from_below(e)"), "e", expr);
00527
00528
00529 if (space_dim == 0
00530 || marked_empty()
00531 || (has_pending_constraints() && !process_pending_constraints())
00532 || (!generators_are_up_to_date() && !update_generators()))
00533 return true;
00534
00535
00536 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
00537 const Generator& g = gen_sys[i];
00538
00539 if (g.is_line_or_ray()) {
00540 const int sp_sign = Scalar_Products::homogeneous_sign(expr, g);
00541 if (sp_sign != 0
00542 && (g.is_line()
00543 || (from_above && sp_sign > 0)
00544 || (!from_above && sp_sign < 0)))
00545
00546 return false;
00547 }
00548 }
00549
00550
00551 return true;
00552 }
00553
00554 bool
00555 PPL::Polyhedron::max_min(const Linear_Expression& expr,
00556 const bool maximize,
00557 Coefficient& ext_n, Coefficient& ext_d,
00558 bool& included,
00559 Generator& g) const {
00560
00561
00562 const dimension_type expr_space_dim = expr.space_dimension();
00563 if (space_dim < expr_space_dim)
00564 throw_dimension_incompatible((maximize
00565 ? "maximize(e, ...)"
00566 : "minimize(e, ...)"), "e", expr);
00567
00568
00569 if (space_dim == 0) {
00570 if (marked_empty())
00571 return false;
00572 else {
00573 ext_n = expr.inhomogeneous_term();
00574 ext_d = 1;
00575 included = true;
00576 g = point();
00577 return true;
00578 }
00579 }
00580
00581
00582 if (marked_empty()
00583 || (has_pending_constraints() && !process_pending_constraints())
00584 || (!generators_are_up_to_date() && !update_generators()))
00585 return false;
00586
00587
00588
00589
00590 DIRTY_TEMP0(mpq_class, extremum);
00591
00592
00593 bool first_candidate = true;
00594
00595
00596 PPL_UNINITIALIZED(dimension_type, ext_position);
00597
00598
00599 PPL_UNINITIALIZED(bool, ext_included);
00600
00601 TEMP_INTEGER(sp);
00602 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
00603 const Generator& gen_sys_i = gen_sys[i];
00604 Scalar_Products::homogeneous_assign(sp, expr, gen_sys_i);
00605
00606 if (gen_sys_i.is_line_or_ray()) {
00607 const int sp_sign = sgn(sp);
00608 if (sp_sign != 0
00609 && (gen_sys_i.is_line()
00610 || (maximize && sp_sign > 0)
00611 || (!maximize && sp_sign < 0)))
00612
00613 return false;
00614 }
00615 else {
00616
00617 assert(gen_sys_i.is_point() || gen_sys_i.is_closure_point());
00618
00619
00620 DIRTY_TEMP0(mpq_class, candidate);
00621 assign_r(candidate.get_num(), sp, ROUND_NOT_NEEDED);
00622 assign_r(candidate.get_den(), gen_sys_i[0], ROUND_NOT_NEEDED);
00623 candidate.canonicalize();
00624 const bool g_is_point = gen_sys_i.is_point();
00625 if (first_candidate
00626 || (maximize
00627 && (candidate > extremum
00628 || (g_is_point
00629 && !ext_included
00630 && candidate == extremum)))
00631 || (!maximize
00632 && (candidate < extremum
00633 || (g_is_point
00634 && !ext_included
00635 && candidate == extremum)))) {
00636
00637 first_candidate = false;
00638 extremum = candidate;
00639 ext_position = i;
00640 ext_included = g_is_point;
00641 }
00642 }
00643 }
00644
00645
00646 DIRTY_TEMP0(mpz_class, n);
00647 assign_r(n, expr.inhomogeneous_term(), ROUND_NOT_NEEDED);
00648 extremum += n;
00649
00650
00651
00652 assert(!first_candidate);
00653 ext_n = Coefficient(extremum.get_num());
00654 ext_d = Coefficient(extremum.get_den());
00655 included = ext_included;
00656 g = gen_sys[ext_position];
00657
00658 return true;
00659 }
00660
00661 void
00662 PPL::Polyhedron::set_zero_dim_univ() {
00663 status.set_zero_dim_univ();
00664 space_dim = 0;
00665 con_sys.clear();
00666 gen_sys.clear();
00667 }
00668
00669 void
00670 PPL::Polyhedron::set_empty() {
00671 status.set_empty();
00672
00673 con_sys.clear();
00674 gen_sys.clear();
00675 sat_c.clear();
00676 sat_g.clear();
00677 }
00678
00679 bool
00680 PPL::Polyhedron::process_pending_constraints() const {
00681 assert(space_dim > 0 && !marked_empty());
00682 assert(has_pending_constraints() && !has_pending_generators());
00683
00684 Polyhedron& x = const_cast<Polyhedron&>(*this);
00685
00686
00687
00688 if (!x.sat_c_is_up_to_date())
00689 x.sat_c.transpose_assign(x.sat_g);
00690 if (!x.con_sys.is_sorted())
00691 x.obtain_sorted_constraints_with_sat_c();
00692
00693
00694 x.con_sys.sort_pending_and_remove_duplicates();
00695 if (x.con_sys.num_pending_rows() == 0) {
00696
00697 x.clear_pending_constraints();
00698 assert(OK(true));
00699 return true;
00700 }
00701
00702 const bool empty = add_and_minimize(true, x.con_sys, x.gen_sys, x.sat_c);
00703 assert(x.con_sys.num_pending_rows() == 0);
00704
00705 if (empty)
00706 x.set_empty();
00707 else {
00708 x.clear_pending_constraints();
00709 x.clear_sat_g_up_to_date();
00710 x.set_sat_c_up_to_date();
00711 }
00712 assert(OK(!empty));
00713 return !empty;
00714 }
00715
00716 void
00717 PPL::Polyhedron::process_pending_generators() const {
00718 assert(space_dim > 0 && !marked_empty());
00719 assert(has_pending_generators() && !has_pending_constraints());
00720
00721 Polyhedron& x = const_cast<Polyhedron&>(*this);
00722
00723
00724
00725 if (!x.sat_g_is_up_to_date())
00726 x.sat_g.transpose_assign(x.sat_c);
00727 if (!x.gen_sys.is_sorted())
00728 x.obtain_sorted_generators_with_sat_g();
00729
00730
00731 x.gen_sys.sort_pending_and_remove_duplicates();
00732 if (x.gen_sys.num_pending_rows() == 0) {
00733
00734 x.clear_pending_generators();
00735 assert(OK(true));
00736 return;
00737 }
00738
00739 add_and_minimize(false, x.gen_sys, x.con_sys, x.sat_g);
00740 assert(x.gen_sys.num_pending_rows() == 0);
00741
00742 x.clear_pending_generators();
00743 x.clear_sat_c_up_to_date();
00744 x.set_sat_g_up_to_date();
00745 }
00746
00747 void
00748 PPL::Polyhedron::remove_pending_to_obtain_constraints() const {
00749 assert(has_something_pending());
00750
00751 Polyhedron& x = const_cast<Polyhedron&>(*this);
00752
00753
00754 if (x.has_pending_constraints()) {
00755
00756 x.con_sys.unset_pending_rows();
00757 x.con_sys.set_sorted(false);
00758 x.clear_pending_constraints();
00759 x.clear_constraints_minimized();
00760 x.clear_generators_up_to_date();
00761 }
00762 else {
00763 assert(x.has_pending_generators());
00764
00765
00766 x.process_pending_generators();
00767 }
00768 assert(OK(true));
00769 }
00770
00771 bool
00772 PPL::Polyhedron::remove_pending_to_obtain_generators() const {
00773 assert(has_something_pending());
00774
00775 Polyhedron& x = const_cast<Polyhedron&>(*this);
00776
00777
00778 if (x.has_pending_generators()) {
00779
00780 x.gen_sys.unset_pending_rows();
00781 x.gen_sys.set_sorted(false);
00782 x.clear_pending_generators();
00783 x.clear_generators_minimized();
00784 x.clear_constraints_up_to_date();
00785 assert(OK(true));
00786 return true;
00787 }
00788 else {
00789 assert(x.has_pending_constraints());
00790
00791
00792 return x.process_pending_constraints();
00793 }
00794 }
00795
00796 void
00797 PPL::Polyhedron::update_constraints() const {
00798 assert(space_dim > 0);
00799 assert(!marked_empty());
00800 assert(generators_are_up_to_date());
00801
00802 assert(!has_something_pending());
00803
00804 Polyhedron& x = const_cast<Polyhedron&>(*this);
00805 minimize(false, x.gen_sys, x.con_sys, x.sat_c);
00806
00807 x.set_sat_c_up_to_date();
00808 x.clear_sat_g_up_to_date();
00809
00810
00811 x.set_constraints_minimized();
00812 x.set_generators_minimized();
00813 }
00814
00815 bool
00816 PPL::Polyhedron::update_generators() const {
00817 assert(space_dim > 0);
00818 assert(!marked_empty());
00819 assert(constraints_are_up_to_date());
00820
00821 assert(!has_something_pending());
00822
00823 Polyhedron& x = const_cast<Polyhedron&>(*this);
00824
00825
00826 const bool empty = minimize(true, x.con_sys, x.gen_sys, x.sat_g);
00827 if (empty)
00828 x.set_empty();
00829 else {
00830
00831 x.set_sat_g_up_to_date();
00832 x.clear_sat_c_up_to_date();
00833
00834
00835 x.set_constraints_minimized();
00836 x.set_generators_minimized();
00837 }
00838 return !empty;
00839 }
00840
00841 void
00842 PPL::Polyhedron::update_sat_c() const {
00843 assert(constraints_are_minimized());
00844 assert(generators_are_minimized());
00845 assert(!sat_c_is_up_to_date());
00846
00847
00848 const dimension_type csr = con_sys.first_pending_row();
00849 const dimension_type gsr = gen_sys.first_pending_row();
00850 Polyhedron& x = const_cast<Polyhedron&>(*this);
00851
00852
00853
00854 x.sat_c.resize(gsr, csr);
00855 for (dimension_type i = gsr; i-- > 0; )
00856 for (dimension_type j = csr; j-- > 0; ) {
00857 const int sp_sign = Scalar_Products::sign(con_sys[j], gen_sys[i]);
00858
00859
00860
00861
00862 assert(sp_sign >= 0);
00863 if (sp_sign > 0)
00864
00865 x.sat_c[i].set(j);
00866 else
00867
00868 x.sat_c[i].clear(j);
00869 }
00870 x.set_sat_c_up_to_date();
00871 }
00872
00873 void
00874 PPL::Polyhedron::update_sat_g() const {
00875 assert(constraints_are_minimized());
00876 assert(generators_are_minimized());
00877 assert(!sat_g_is_up_to_date());
00878
00879
00880 const dimension_type csr = con_sys.first_pending_row();
00881 const dimension_type gsr = gen_sys.first_pending_row();
00882 Polyhedron& x = const_cast<Polyhedron&>(*this);
00883
00884
00885
00886 x.sat_g.resize(csr, gsr);
00887 for (dimension_type i = csr; i-- > 0; )
00888 for (dimension_type j = gsr; j-- > 0; ) {
00889 const int sp_sign = Scalar_Products::sign(con_sys[i], gen_sys[j]);
00890
00891
00892
00893
00894 assert(sp_sign >= 0);
00895 if (sp_sign > 0)
00896
00897 x.sat_g[i].set(j);
00898 else
00899
00900 x.sat_g[i].clear(j);
00901 }
00902 x.set_sat_g_up_to_date();
00903 }
00904
00905 void
00906 PPL::Polyhedron::obtain_sorted_constraints() const {
00907 assert(constraints_are_up_to_date());
00908
00909 Polyhedron& x = const_cast<Polyhedron&>(*this);
00910 if (!x.con_sys.is_sorted()) {
00911 if (x.sat_g_is_up_to_date()) {
00912
00913 x.con_sys.sort_and_remove_with_sat(x.sat_g);
00914
00915 x.clear_sat_c_up_to_date();
00916 }
00917 else if (x.sat_c_is_up_to_date()) {
00918
00919 x.sat_g.transpose_assign(x.sat_c);
00920 x.con_sys.sort_and_remove_with_sat(x.sat_g);
00921 x.set_sat_g_up_to_date();
00922 x.clear_sat_c_up_to_date();
00923 }
00924 else
00925
00926
00927 x.con_sys.sort_rows();
00928 }
00929
00930 assert(con_sys.check_sorted());
00931 }
00932
00933 void
00934 PPL::Polyhedron::obtain_sorted_generators() const {
00935 assert(generators_are_up_to_date());
00936
00937 Polyhedron& x = const_cast<Polyhedron&>(*this);
00938 if (!x.gen_sys.is_sorted()) {
00939 if (x.sat_c_is_up_to_date()) {
00940
00941 x.gen_sys.sort_and_remove_with_sat(x.sat_c);
00942
00943 x.clear_sat_g_up_to_date();
00944 }
00945 else if (x.sat_g_is_up_to_date()) {
00946
00947 x.sat_c.transpose_assign(x.sat_g);
00948 x.gen_sys.sort_and_remove_with_sat(x.sat_c);
00949 x.set_sat_c_up_to_date();
00950 x.clear_sat_g_up_to_date();
00951 }
00952 else
00953
00954
00955 x.gen_sys.sort_rows();
00956 }
00957
00958 assert(gen_sys.check_sorted());
00959 }
00960
00961 void
00962 PPL::Polyhedron::obtain_sorted_constraints_with_sat_c() const {
00963 assert(constraints_are_up_to_date());
00964 assert(constraints_are_minimized());
00965
00966 Polyhedron& x = const_cast<Polyhedron&>(*this);
00967
00968 if (!x.sat_c_is_up_to_date() && !x.sat_g_is_up_to_date())
00969 x.update_sat_c();
00970
00971 if (x.con_sys.is_sorted()) {
00972 if (x.sat_c_is_up_to_date())
00973
00974
00975 return;
00976 }
00977 else {
00978 if (!x.sat_g_is_up_to_date()) {
00979
00980
00981 x.sat_g.transpose_assign(x.sat_c);
00982 x.set_sat_g_up_to_date();
00983 }
00984
00985 x.con_sys.sort_and_remove_with_sat(x.sat_g);
00986 }
00987
00988 x.sat_c.transpose_assign(x.sat_g);
00989 x.set_sat_c_up_to_date();
00990
00991 x.con_sys.set_sorted(true);
00992
00993 assert(con_sys.check_sorted());
00994 }
00995
00996 void
00997 PPL::Polyhedron::obtain_sorted_generators_with_sat_g() const {
00998 assert(generators_are_up_to_date());
00999
01000 Polyhedron& x = const_cast<Polyhedron&>(*this);
01001
01002 if (!x.sat_c_is_up_to_date() && !x.sat_g_is_up_to_date())
01003 x.update_sat_g();
01004
01005 if (x.gen_sys.is_sorted()) {
01006 if (x.sat_g_is_up_to_date())
01007
01008
01009 return;
01010 }
01011 else {
01012 if (!x.sat_c_is_up_to_date()) {
01013
01014
01015 x.sat_c.transpose_assign(x.sat_g);
01016 x.set_sat_c_up_to_date();
01017 }
01018
01019 x.gen_sys.sort_and_remove_with_sat(x.sat_c);
01020 }
01021
01022 x.sat_g.transpose_assign(sat_c);
01023 x.set_sat_g_up_to_date();
01024
01025 x.gen_sys.set_sorted(true);
01026
01027 assert(gen_sys.check_sorted());
01028 }
01029
01030 bool
01031 PPL::Polyhedron::minimize() const {
01032
01033 if (marked_empty())
01034 return false;
01035 if (space_dim == 0)
01036 return true;
01037
01038
01039 if (has_something_pending()) {
01040 const bool not_empty = process_pending();
01041 assert(OK());
01042 return not_empty;
01043 }
01044
01045
01046
01047 if (constraints_are_minimized() && generators_are_minimized())
01048 return true;
01049
01050
01051
01052
01053
01054
01055
01056 if (constraints_are_up_to_date()) {
01057
01058 const bool ret = update_generators();
01059 assert(OK());
01060 return ret;
01061 }
01062 else {
01063 assert(generators_are_up_to_date());
01064 update_constraints();
01065 assert(OK());
01066 return true;
01067 }
01068 }
01069
01070 bool
01071 PPL::Polyhedron::strongly_minimize_constraints() const {
01072 assert(!is_necessarily_closed());
01073
01074
01075 Polyhedron& x = const_cast<Polyhedron&>(*this);
01076
01077
01078
01079 if (!minimize())
01080 return false;
01081
01082
01083
01084 if (x.space_dim == 0)
01085 return true;
01086
01087
01088 if (!sat_g_is_up_to_date()) {
01089 assert(sat_c_is_up_to_date());
01090 x.sat_g.transpose_assign(sat_c);
01091 }
01092
01093
01094
01095
01096 Bit_Row sat_all_but_rays;
01097 Bit_Row sat_all_but_points;
01098 Bit_Row sat_all_but_closure_points;
01099
01100 const dimension_type gs_rows = gen_sys.num_rows();
01101 const dimension_type n_lines = gen_sys.num_lines();
01102 for (dimension_type i = gs_rows; i-- > n_lines; )
01103 switch (gen_sys[i].type()) {
01104 case Generator::RAY:
01105 sat_all_but_rays.set(i);
01106 break;
01107 case Generator::POINT:
01108 sat_all_but_points.set(i);
01109 break;
01110 case Generator::CLOSURE_POINT:
01111 sat_all_but_closure_points.set(i);
01112 break;
01113 default:
01114
01115 throw std::runtime_error("PPL internal error: "
01116 "strongly_minimize_constraints.");
01117 }
01118 Bit_Row sat_lines_and_rays;
01119 set_union(sat_all_but_points, sat_all_but_closure_points,
01120 sat_lines_and_rays);
01121 Bit_Row sat_lines_and_closure_points;
01122 set_union(sat_all_but_rays, sat_all_but_points,
01123 sat_lines_and_closure_points);
01124 Bit_Row sat_lines;
01125 set_union(sat_lines_and_rays, sat_lines_and_closure_points,
01126 sat_lines);
01127
01128
01129
01130
01131 bool changed = false;
01132 bool found_eps_leq_one = false;
01133
01134
01135
01136
01137 Constraint_System& cs = x.con_sys;
01138 Bit_Matrix& sat = x.sat_g;
01139 dimension_type cs_rows = cs.num_rows();
01140 const dimension_type eps_index = cs.num_columns() - 1;
01141 for (dimension_type i = 0; i < cs_rows; )
01142 if (cs[i].is_strict_inequality()) {
01143
01144 Bit_Row sat_ci;
01145 set_union(sat[i], sat_lines_and_closure_points, sat_ci);
01146 if (sat_ci == sat_lines) {
01147
01148 if (!found_eps_leq_one) {
01149
01150 const Constraint& c = cs[i];
01151 bool all_zeroes = true;
01152 for (dimension_type k = eps_index; k-- > 1; )
01153 if (c[k] != 0) {
01154 all_zeroes = false;
01155 break;
01156 }
01157 if (all_zeroes && (c[0] + c[eps_index] == 0)) {
01158
01159 found_eps_leq_one = true;
01160
01161 ++i;
01162 continue;
01163 }
01164 }
01165
01166
01167
01168
01169 --cs_rows;
01170 std::swap(cs[i], cs[cs_rows]);
01171 std::swap(sat[i], sat[cs_rows]);
01172
01173 changed = true;
01174
01175
01176 continue;
01177 }
01178
01179
01180
01181 sat_ci.clear();
01182 set_union(sat[i], sat_all_but_points, sat_ci);
01183 bool eps_redundant = false;
01184 for (dimension_type j = 0; j < cs_rows; ++j)
01185 if (i != j && cs[j].is_strict_inequality()
01186 && subset_or_equal(sat[j], sat_ci)) {
01187
01188
01189
01190 --cs_rows;
01191 std::swap(cs[i], cs[cs_rows]);
01192 std::swap(sat[i], sat[cs_rows]);
01193 eps_redundant = true;
01194
01195 changed = true;
01196 break;
01197 }
01198
01199
01200 if (!eps_redundant)
01201 ++i;
01202 }
01203 else
01204
01205 ++i;
01206
01207 if (changed) {
01208
01209
01210 assert(cs_rows < cs.num_rows());
01211 cs.erase_to_end(cs_rows);
01212
01213 cs.unset_pending_rows();
01214
01215 cs.set_sorted(false);
01216
01217 x.clear_generators_up_to_date();
01218
01219
01220
01221
01222 if (!found_eps_leq_one) {
01223 MIP_Problem lp;
01224
01225
01226
01227
01228 cs.set_necessarily_closed();
01229 try {
01230 lp.add_space_dimensions_and_embed(cs.space_dimension());
01231 lp.add_constraints(cs);
01232 cs.set_not_necessarily_closed();
01233 }
01234 catch (...) {
01235 cs.set_not_necessarily_closed();
01236 throw;
01237 }
01238
01239 lp.set_objective_function(Variable(x.space_dim));
01240 lp.set_optimization_mode(MAXIMIZATION);
01241 MIP_Problem_Status status = lp.solve();
01242 assert(status != UNFEASIBLE_MIP_PROBLEM);
01243
01244
01245 if (status == UNBOUNDED_MIP_PROBLEM)
01246 cs.insert(Constraint::epsilon_leq_one());
01247 }
01248 }
01249
01250 assert(OK());
01251 return true;
01252 }
01253
01254 bool
01255 PPL::Polyhedron::strongly_minimize_generators() const {
01256 assert(!is_necessarily_closed());
01257
01258
01259 Polyhedron& x = const_cast<Polyhedron&>(*this);
01260
01261
01262
01263 if (!minimize())
01264 return false;
01265
01266
01267
01268 if (x.space_dim == 0)
01269 return true;
01270
01271
01272 if (!sat_c_is_up_to_date()) {
01273 assert(sat_g_is_up_to_date());
01274 x.sat_c.transpose_assign(sat_g);
01275 }
01276
01277
01278
01279 Bit_Row sat_all_but_strict_ineq;
01280 const dimension_type cs_rows = con_sys.num_rows();
01281 const dimension_type n_equals = con_sys.num_equalities();
01282 for (dimension_type i = cs_rows; i-- > n_equals; )
01283 if (con_sys[i].is_strict_inequality())
01284 sat_all_but_strict_ineq.set(i);
01285
01286
01287 bool changed = false;
01288
01289
01290
01291 Generator_System& gs = const_cast<Generator_System&>(gen_sys);
01292 Bit_Matrix& sat = const_cast<Bit_Matrix&>(sat_c);
01293 dimension_type gs_rows = gs.num_rows();
01294 const dimension_type n_lines = gs.num_lines();
01295 const dimension_type eps_index = gs.num_columns() - 1;
01296 for (dimension_type i = n_lines; i < gs_rows; )
01297 if (gs[i].is_point()) {
01298
01299
01300 Bit_Row sat_gi;
01301 set_union(sat[i], sat_all_but_strict_ineq, sat_gi);
01302
01303
01304
01305 bool eps_redundant = false;
01306 for (dimension_type j = n_lines; j < gs_rows; ++j)
01307 if (i != j && gs[j].is_point() && subset_or_equal(sat[j], sat_gi)) {
01308
01309
01310
01311 --gs_rows;
01312 std::swap(gs[i], gs[gs_rows]);
01313 std::swap(sat[i], sat[gs_rows]);
01314 eps_redundant = true;
01315 changed = true;
01316 break;
01317 }
01318 if (!eps_redundant) {
01319
01320 Generator& gi = gs[i];
01321 if (gi[eps_index] != gi[0]) {
01322 gi[eps_index] = gi[0];
01323
01324 gi.normalize();
01325 changed = true;
01326 }
01327
01328 ++i;
01329 }
01330 }
01331 else
01332
01333 ++i;
01334
01335
01336
01337 if (gs_rows < gs.num_rows()) {
01338 gs.erase_to_end(gs_rows);
01339 gs.unset_pending_rows();
01340 }
01341
01342 if (changed) {
01343
01344 x.gen_sys.set_sorted(false);
01345
01346 x.clear_constraints_up_to_date();
01347 }
01348
01349 assert(OK());
01350 return true;
01351 }
01352
01353 void
01354 PPL::Polyhedron::refine_no_check(const Constraint& c) {
01355 assert(!marked_empty());
01356 assert(space_dim >= c.space_dimension());
01357
01358
01359 if (space_dim == 0) {
01360 if (c.is_inconsistent())
01361 set_empty();
01362 return;
01363 }
01364
01365
01366 if (has_pending_generators())
01367 process_pending_generators();
01368 else if (!constraints_are_up_to_date())
01369 update_constraints();
01370
01371 const bool adding_pending = can_have_something_pending();
01372
01373 if (c.is_necessarily_closed() || !is_necessarily_closed())
01374
01375
01376 if (adding_pending)
01377 con_sys.insert_pending(c);
01378 else
01379 con_sys.insert(c);
01380 else {
01381
01382
01383
01384
01385 Linear_Expression nc_expr = Linear_Expression(c);
01386 if (c.is_equality())
01387 if (adding_pending)
01388 con_sys.insert_pending(nc_expr == 0);
01389 else
01390 con_sys.insert(nc_expr == 0);
01391 else
01392 if (adding_pending)
01393 con_sys.insert_pending(nc_expr >= 0);
01394 else
01395 con_sys.insert(nc_expr >= 0);
01396 }
01397
01398 if (adding_pending)
01399 set_constraints_pending();
01400 else {
01401
01402 clear_constraints_minimized();
01403 clear_generators_up_to_date();
01404 }
01405
01406
01407
01408 assert(OK());
01409 }
01410
01411 void
01412 PPL::Polyhedron::throw_runtime_error(const char* method) const {
01413 std::ostringstream s;
01414 s << "PPL::";
01415 if (is_necessarily_closed())
01416 s << "C_";
01417 else
01418 s << "NNC_";
01419 s << "Polyhedron::" << method << "." << std::endl;
01420 throw std::runtime_error(s.str());
01421 }
01422
01423 void
01424 PPL::Polyhedron::throw_invalid_argument(const char* method,
01425 const char* reason) const {
01426 std::ostringstream s;
01427 s << "PPL::";
01428 if (is_necessarily_closed())
01429 s << "C_";
01430 else
01431 s << "NNC_";
01432 s << "Polyhedron::" << method << ":" << std::endl
01433 << reason << ".";
01434 throw std::invalid_argument(s.str());
01435 }
01436
01437 void
01438 PPL::Polyhedron::throw_topology_incompatible(const char* method,
01439 const char* ph_name,
01440 const Polyhedron& ph) const {
01441 std::ostringstream s;
01442 s << "PPL::";
01443 if (is_necessarily_closed())
01444 s << "C_";
01445 else
01446 s << "NNC_";
01447 s << "Polyhedron::" << method << ":" << std::endl
01448 << ph_name << " is a ";
01449 if (ph.is_necessarily_closed())
01450 s << "C_";
01451 else
01452 s << "NNC_";
01453 s << "Polyhedron." << std::endl;
01454 throw std::invalid_argument(s.str());
01455 }
01456
01457 void
01458 PPL::Polyhedron::throw_topology_incompatible(const char* method,
01459 const char* c_name,
01460 const Constraint&) const {
01461 assert(is_necessarily_closed());
01462 std::ostringstream s;
01463 s << "PPL::C_Polyhedron::" << method << ":" << std::endl
01464 << c_name << " is a strict inequality.";
01465 throw std::invalid_argument(s.str());
01466 }
01467
01468 void
01469 PPL::Polyhedron::throw_topology_incompatible(const char* method,
01470 const char* g_name,
01471 const Generator&) const {
01472 assert(is_necessarily_closed());
01473 std::ostringstream s;
01474 s << "PPL::C_Polyhedron::" << method << ":" << std::endl
01475 << g_name << " is a closure point.";
01476 throw std::invalid_argument(s.str());
01477 }
01478
01479 void
01480 PPL::Polyhedron::throw_topology_incompatible(const char* method,
01481 const char* cs_name,
01482 const Constraint_System&) const {
01483 assert(is_necessarily_closed());
01484 std::ostringstream s;
01485 s << "PPL::C_Polyhedron::" << method << ":" << std::endl
01486 << cs_name << " contains strict inequalities.";
01487 throw std::invalid_argument(s.str());
01488 }
01489
01490 void
01491 PPL::Polyhedron::throw_topology_incompatible(const char* method,
01492 const char* gs_name,
01493 const Generator_System&) const {
01494 std::ostringstream s;
01495 s << "PPL::C_Polyhedron::" << method << ":" << std::endl
01496 << gs_name << " contains closure points.";
01497 throw std::invalid_argument(s.str());
01498 }
01499
01500 void
01501 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01502 const char* other_name,
01503 dimension_type other_dim) const {
01504 std::ostringstream s;
01505 s << "PPL::"
01506 << (is_necessarily_closed() ? "C_" : "NNC_")
01507 << "Polyhedron::" << method << ":\n"
01508 << "this->space_dimension() == " << space_dimension() << ", "
01509 << other_name << ".space_dimension() == " << other_dim << ".";
01510 throw std::invalid_argument(s.str());
01511 }
01512
01513 void
01514 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01515 const char* ph_name,
01516 const Polyhedron& ph) const {
01517 throw_dimension_incompatible(method, ph_name, ph.space_dimension());
01518 }
01519
01520 void
01521 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01522 const char* e_name,
01523 const Linear_Expression& e) const {
01524 throw_dimension_incompatible(method, e_name, e.space_dimension());
01525 }
01526
01527 void
01528 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01529 const char* c_name,
01530 const Constraint& c) const {
01531 throw_dimension_incompatible(method, c_name, c.space_dimension());
01532 }
01533
01534 void
01535 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01536 const char* g_name,
01537 const Generator& g) const {
01538 throw_dimension_incompatible(method, g_name, g.space_dimension());
01539 }
01540
01541 void
01542 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01543 const char* cg_name,
01544 const Congruence& cg) const {
01545 throw_dimension_incompatible(method, cg_name, cg.space_dimension());
01546 }
01547
01548 void
01549 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01550 const char* cs_name,
01551 const Constraint_System& cs) const {
01552 throw_dimension_incompatible(method, cs_name, cs.space_dimension());
01553 }
01554
01555 void
01556 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01557 const char* gs_name,
01558 const Generator_System& gs) const {
01559 throw_dimension_incompatible(method, gs_name, gs.space_dimension());
01560 }
01561
01562 void
01563 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01564 const char* cgs_name,
01565 const Congruence_System& cgs) const {
01566 throw_dimension_incompatible(method, cgs_name, cgs.space_dimension());
01567 }
01568
01569 void
01570 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
01571 const char* var_name,
01572 const Variable var) const {
01573 std::ostringstream s;
01574 s << "PPL::";
01575 if (is_necessarily_closed())
01576 s << "C_";
01577 else
01578 s << "NNC_";
01579 s << "Polyhedron::" << method << ":" << std::endl
01580 << "this->space_dimension() == " << space_dimension() << ", "
01581 << var_name << ".space_dimension() == " << var.space_dimension() << ".";
01582 throw std::invalid_argument(s.str());
01583 }
01584
01585 void
01586 PPL::Polyhedron::
01587 throw_dimension_incompatible(const char* method,
01588 dimension_type required_space_dim) const {
01589 std::ostringstream s;
01590 s << "PPL::";
01591 if (is_necessarily_closed())
01592 s << "C_";
01593 else
01594 s << "NNC_";
01595 s << "Polyhedron::" << method << ":" << std::endl
01596 << "this->space_dimension() == " << space_dimension()
01597 << ", required space dimension == " << required_space_dim << ".";
01598 throw std::invalid_argument(s.str());
01599 }
01600
01601 void
01602 PPL::Polyhedron::throw_space_dimension_overflow(const Topology topol,
01603 const char* method,
01604 const char* reason) {
01605 std::ostringstream s;
01606 s << "PPL::";
01607 if (topol == NECESSARILY_CLOSED)
01608 s << "C_";
01609 else
01610 s << "NNC_";
01611 s << "Polyhedron::" << method << ":" << std::endl
01612 << reason << ".";
01613 throw std::length_error(s.str());
01614 }
01615
01616 void
01617 PPL::Polyhedron::throw_invalid_generator(const char* method,
01618 const char* g_name) const {
01619 std::ostringstream s;
01620 s << "PPL::";
01621 if (is_necessarily_closed())
01622 s << "C_";
01623 else
01624 s << "NNC_";
01625 s << "Polyhedron::" << method << ":" << std::endl
01626 << "*this is an empty polyhedron and "
01627 << g_name << " is not a point.";
01628 throw std::invalid_argument(s.str());
01629 }
01630
01631 void
01632 PPL::Polyhedron::throw_invalid_generators(const char* method,
01633 const char* gs_name) const {
01634 std::ostringstream s;
01635 s << "PPL::";
01636 if (is_necessarily_closed())
01637 s << "C_";
01638 else
01639 s << "NNC_";
01640 s << "Polyhedron::" << method << ":" << std::endl
01641 << "*this is an empty polyhedron and" << std::endl
01642 << "the non-empty generator system " << gs_name << " contains no points.";
01643 throw std::invalid_argument(s.str());
01644 }