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