00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <config.h>
00024
00025 #include "Polyhedron.defs.hh"
00026 #include "Scalar_Products.defs.hh"
00027
00028 #include <cassert>
00029 #include <iostream>
00030
00031 #ifndef ENSURE_SORTEDNESS
00032 #define ENSURE_SORTEDNESS 0
00033 #endif
00034
00035 namespace PPL = Parma_Polyhedra_Library;
00036
00037 PPL::dimension_type
00038 PPL::Polyhedron::affine_dimension() const {
00039 if (is_empty())
00040 return 0;
00041
00042 const Constraint_System& cs = minimized_constraints();
00043 dimension_type d = space_dim;
00044 for (Constraint_System::const_iterator i = cs.begin(),
00045 cs_end = cs.end(); i != cs_end; ++i)
00046 if (i->is_equality())
00047 --d;
00048 return d;
00049 }
00050
00051 const PPL::Constraint_System&
00052 PPL::Polyhedron::constraints() const {
00053 if (marked_empty()) {
00054
00055
00056 if (con_sys.num_rows() == 0) {
00057
00058
00059 Constraint_System unsat_cs = Constraint_System::zero_dim_empty();
00060 unsat_cs.adjust_topology_and_space_dimension(topology(), space_dim);
00061 const_cast<Constraint_System&>(con_sys).swap(unsat_cs);
00062 }
00063 else {
00064
00065 assert(con_sys.space_dimension() == space_dim);
00066 assert(con_sys.num_rows() == 1);
00067 assert(con_sys[0].is_inconsistent());
00068 }
00069 return con_sys;
00070 }
00071
00072 if (space_dim == 0) {
00073
00074 assert(con_sys.num_columns() == 0 && con_sys.num_rows() == 0);
00075 return con_sys;
00076 }
00077
00078
00079
00080
00081 if (has_pending_generators())
00082 process_pending_generators();
00083 else if (!constraints_are_up_to_date())
00084 update_constraints();
00085
00086
00087 #if ENSURE_SORTEDNESS
00088
00089
00090 if (!has_pending_constraints())
00091 obtain_sorted_constraints();
00092 #endif
00093 return con_sys;
00094 }
00095
00096 const PPL::Constraint_System&
00097 PPL::Polyhedron::minimized_constraints() const {
00098
00099
00100 if (is_necessarily_closed())
00101 minimize();
00102 else
00103 strongly_minimize_constraints();
00104 return constraints();
00105 }
00106
00107 const PPL::Generator_System&
00108 PPL::Polyhedron::generators() const {
00109 if (marked_empty()) {
00110 assert(gen_sys.num_rows() == 0);
00111
00112
00113 if (gen_sys.space_dimension() != space_dim) {
00114 Generator_System gs;
00115 gs.adjust_topology_and_space_dimension(topology(), space_dim);
00116 const_cast<Generator_System&>(gen_sys).swap(gs);
00117 }
00118 return gen_sys;
00119 }
00120
00121 if (space_dim == 0) {
00122 assert(gen_sys.num_columns() == 0 && gen_sys.num_rows() == 0);
00123 return Generator_System::zero_dim_univ();
00124 }
00125
00126
00127
00128
00129 if ((has_pending_constraints() && !process_pending_constraints())
00130 || (!generators_are_up_to_date() && !update_generators())) {
00131
00132 assert(gen_sys.num_rows() == 0);
00133
00134
00135 if (gen_sys.space_dimension() != space_dim) {
00136 Generator_System gs;
00137 gs.adjust_topology_and_space_dimension(topology(), space_dim);
00138 const_cast<Generator_System&>(gen_sys).swap(gs);
00139 }
00140 return gen_sys;
00141 }
00142
00143
00144 #if ENSURE_SORTEDNESS
00145
00146
00147 if (!has_pending_generators())
00148 obtain_sorted_generators();
00149 #else
00150
00151
00152
00153
00154 if (!is_necessarily_closed()
00155 && generators_are_minimized() && !has_pending_generators())
00156 obtain_sorted_generators();
00157 #endif
00158 return gen_sys;
00159 }
00160
00161 const PPL::Generator_System&
00162 PPL::Polyhedron::minimized_generators() const {
00163
00164
00165 if (is_necessarily_closed())
00166 minimize();
00167 else
00168 strongly_minimize_generators();
00169
00170
00171
00172 return generators();
00173 }
00174
00175 PPL::Poly_Con_Relation
00176 PPL::Polyhedron::relation_with(const Constraint& c) const {
00177
00178 if (space_dim < c.space_dimension())
00179 throw_dimension_incompatible("relation_with(c)", "c", c);
00180
00181 if (marked_empty())
00182 return Poly_Con_Relation::saturates()
00183 && Poly_Con_Relation::is_included()
00184 && Poly_Con_Relation::is_disjoint();
00185
00186 if (space_dim == 0)
00187 if (c.is_inconsistent())
00188 if (c.is_strict_inequality() && c.inhomogeneous_term() == 0)
00189
00190
00191 return Poly_Con_Relation::saturates()
00192 && Poly_Con_Relation::is_disjoint();
00193 else
00194 return Poly_Con_Relation::is_disjoint();
00195 else if (c.is_equality() || c.inhomogeneous_term() == 0)
00196 return Poly_Con_Relation::saturates()
00197 && Poly_Con_Relation::is_included();
00198 else
00199
00200
00201
00202 return Poly_Con_Relation::is_included();
00203
00204 if ((has_pending_constraints() && !process_pending_constraints())
00205 || (!generators_are_up_to_date() && !update_generators()))
00206
00207 return Poly_Con_Relation::saturates()
00208 && Poly_Con_Relation::is_included()
00209 && Poly_Con_Relation::is_disjoint();
00210
00211 return gen_sys.relation_with(c);
00212 }
00213
00214 PPL::Poly_Gen_Relation
00215 PPL::Polyhedron::relation_with(const Generator& g) const {
00216
00217 if (space_dim < g.space_dimension())
00218 throw_dimension_incompatible("relation_with(g)", "g", g);
00219
00220
00221 if (marked_empty())
00222 return Poly_Gen_Relation::nothing();
00223
00224
00225
00226 if (space_dim == 0)
00227 return Poly_Gen_Relation::subsumes();
00228
00229 if (has_pending_generators())
00230 process_pending_generators();
00231 else if (!constraints_are_up_to_date())
00232 update_constraints();
00233
00234 return
00235 con_sys.satisfies_all_constraints(g)
00236 ? Poly_Gen_Relation::subsumes()
00237 : Poly_Gen_Relation::nothing();
00238 }
00239
00240 bool
00241 PPL::Polyhedron::is_universe() const {
00242 if (marked_empty())
00243 return false;
00244
00245 if (space_dim == 0)
00246 return true;
00247
00248 if (!has_pending_generators() && constraints_are_up_to_date()) {
00249
00250 for (dimension_type i = con_sys.num_rows(); i-- > 0; )
00251 if (!con_sys[i].is_tautological())
00252 return false;
00253
00254 return true;
00255 }
00256
00257 assert(!has_pending_constraints() && generators_are_up_to_date());
00258
00259
00260 dimension_type num_lines = 0;
00261 dimension_type num_rays = 0;
00262 const dimension_type first_pending = gen_sys.first_pending_row();
00263 for (dimension_type i = first_pending; i-- > 0; )
00264 switch (gen_sys[i].type()) {
00265 case Generator::RAY:
00266 ++num_rays;
00267 break;
00268 case Generator::LINE:
00269 ++num_lines;
00270 break;
00271 default:
00272 break;
00273 }
00274
00275 if (has_pending_generators()) {
00276
00277
00278 assert(generators_are_minimized());
00279 if (num_lines == space_dim) {
00280 assert(num_rays == 0);
00281 return true;
00282 }
00283 assert(num_lines < space_dim);
00284
00285 dimension_type num_pending_lines = 0;
00286 dimension_type num_pending_rays = 0;
00287 const dimension_type gs_num_rows = gen_sys.num_rows();
00288 for (dimension_type i = first_pending; i < gs_num_rows; ++i)
00289 switch (gen_sys[i].type()) {
00290 case Generator::RAY:
00291 ++num_pending_rays;
00292 break;
00293 case Generator::LINE:
00294 ++num_pending_lines;
00295 break;
00296 default:
00297 break;
00298 }
00299
00300
00301 if (num_pending_rays == 0 && num_pending_lines == 0)
00302 return false;
00303
00304
00305 if (num_lines + num_pending_lines < space_dim) {
00306 const dimension_type num_dims_missing
00307 = space_dim - (num_lines + num_pending_lines);
00308
00309
00310 if (num_rays + num_pending_rays <= num_dims_missing)
00311 return false;
00312 }
00313 }
00314 else {
00315
00316 if (generators_are_minimized()) {
00317
00318 assert(num_rays == 0 || num_lines < space_dim);
00319 return num_lines == space_dim;
00320 }
00321 else
00322
00323
00324
00325 if (num_lines < space_dim && num_lines + num_rays <= space_dim)
00326 return false;
00327 }
00328
00329
00330 if (has_pending_generators())
00331 process_pending_generators();
00332 else if (!constraints_are_minimized())
00333 minimize();
00334 if (is_necessarily_closed())
00335 return (con_sys.num_rows() == 1
00336 && con_sys[0].is_inequality()
00337 && con_sys[0].is_tautological());
00338 else {
00339
00340 if (con_sys.num_rows() != 2
00341 || con_sys[0].is_equality()
00342 || con_sys[1].is_equality())
00343 return false;
00344 else {
00345
00346
00347
00348
00349 #ifndef NDEBUG
00350 obtain_sorted_constraints();
00351 const Constraint& eps_leq_one = con_sys[0];
00352 const Constraint& eps_geq_zero = con_sys[1];
00353 const dimension_type eps_index = con_sys.num_columns() - 1;
00354 assert(eps_leq_one[0] > 0 && eps_leq_one[eps_index] < 0
00355 && eps_geq_zero[0] == 0 && eps_geq_zero[eps_index] > 0);
00356 for (dimension_type i = 1; i < eps_index; ++i)
00357 assert(eps_leq_one[i] == 0 && eps_geq_zero[i] == 0);
00358 #endif
00359 return true;
00360 }
00361 }
00362 }
00363
00364 bool
00365 PPL::Polyhedron::is_bounded() const {
00366
00367 if (space_dim == 0
00368 || marked_empty()
00369 || (has_pending_constraints() && !process_pending_constraints())
00370 || (!generators_are_up_to_date() && !update_generators()))
00371 return true;
00372
00373
00374
00375 for (dimension_type i = gen_sys.num_rows(); i-- > 0; )
00376 if (gen_sys[i].is_line_or_ray())
00377 return false;
00378
00379
00380
00381 return true;
00382 }
00383
00384 bool
00385 PPL::Polyhedron::is_topologically_closed() const {
00386
00387 if (is_necessarily_closed())
00388 return true;
00389
00390 if (marked_empty()
00391 || space_dim == 0
00392 || (has_something_pending() && !process_pending()))
00393 return true;
00394
00395
00396 assert(!has_something_pending());
00397
00398 if (generators_are_minimized()) {
00399
00400
00401 const dimension_type n_rows = gen_sys.num_rows();
00402 const dimension_type n_lines = gen_sys.num_lines();
00403 for (dimension_type i = n_rows; i-- > n_lines; ) {
00404 const Generator& gi = gen_sys[i];
00405 if (gi.is_closure_point()) {
00406 bool gi_has_no_matching_point = true;
00407 for (dimension_type j = n_rows; j-- > n_lines; ) {
00408 const Generator& gj = gen_sys[j];
00409 if (i != j
00410 && gj.is_point()
00411 && gi.is_matching_closure_point(gj)) {
00412 gi_has_no_matching_point = false;
00413 break;
00414 }
00415 }
00416 if (gi_has_no_matching_point)
00417 return false;
00418 }
00419 }
00420
00421 return true;
00422 }
00423
00424
00425
00426 strongly_minimize_constraints();
00427 return marked_empty() || !con_sys.has_strict_inequalities();
00428 }
00429
00430 bool
00431 PPL::Polyhedron::OK(bool check_not_empty) const {
00432 #ifndef NDEBUG
00433 using std::endl;
00434 using std::cerr;
00435 #endif
00436
00437
00438
00439 const dimension_type poly_num_columns
00440 = space_dim + (is_necessarily_closed() ? 1 : 2);
00441
00442
00443 if (con_sys.topology() != gen_sys.topology()) {
00444 #ifndef NDEBUG
00445 cerr << "Constraints and generators have different topologies!"
00446 << endl;
00447 #endif
00448 goto bomb;
00449 }
00450
00451
00452 if (!sat_c.OK())
00453 goto bomb;
00454 if (!sat_g.OK())
00455 goto bomb;
00456
00457
00458 if (!status.OK())
00459 goto bomb;
00460
00461 if (marked_empty()) {
00462 if (check_not_empty) {
00463
00464 #ifndef NDEBUG
00465 cerr << "Empty polyhedron!" << endl;
00466 #endif
00467 goto bomb;
00468 }
00469
00470
00471
00472
00473 if (has_something_pending()) {
00474 #ifndef NDEBUG
00475 cerr << "The polyhedron is empty, "
00476 << "but it has something pending" << endl;
00477 #endif
00478 goto bomb;
00479 }
00480 if (con_sys.num_rows() == 0)
00481 return true;
00482 else {
00483 if (con_sys.space_dimension() != space_dim) {
00484 #ifndef NDEBUG
00485 cerr << "The polyhedron is in a space of dimension "
00486 << space_dim
00487 << " while the system of constraints is in a space of dimension "
00488 << con_sys.space_dimension()
00489 << endl;
00490 #endif
00491 goto bomb;
00492 }
00493 if (con_sys.num_rows() != 1) {
00494 #ifndef NDEBUG
00495 cerr << "The system of constraints for an empty polyhedron "
00496 << "has more then one row"
00497 << endl;
00498 #endif
00499 goto bomb;
00500 }
00501 if (!con_sys[0].is_inconsistent()) {
00502 #ifndef NDEBUG
00503 cerr << "Empty polyhedron with a satisfiable system of constraints"
00504 << endl;
00505 #endif
00506 goto bomb;
00507 }
00508
00509 return true;
00510 }
00511 }
00512
00513
00514
00515
00516 if (space_dim == 0) {
00517 if (has_something_pending()) {
00518 #ifndef NDEBUG
00519 cerr << "Zero-dimensional polyhedron with something pending"
00520 << endl;
00521 #endif
00522 goto bomb;
00523 }
00524 if (con_sys.num_rows() != 0 || gen_sys.num_rows() != 0) {
00525 #ifndef NDEBUG
00526 cerr << "Zero-dimensional polyhedron with a non-empty"
00527 << endl
00528 << "system of constraints or generators."
00529 << endl;
00530 #endif
00531 goto bomb;
00532 }
00533 return true;
00534 }
00535
00536
00537
00538 if (!constraints_are_up_to_date() && !generators_are_up_to_date()) {
00539 #ifndef NDEBUG
00540 cerr << "Polyhedron not empty, not zero-dimensional"
00541 << endl
00542 << "and with neither constraints nor generators up-to-date!"
00543 << endl;
00544 #endif
00545 goto bomb;
00546 }
00547
00548
00549
00550
00551
00552
00553
00554 if (constraints_are_up_to_date()) {
00555 if (con_sys.num_columns() != poly_num_columns) {
00556 #ifndef NDEBUG
00557 cerr << "Incompatible size! (con_sys and space_dim)"
00558 << endl;
00559 #endif
00560 goto bomb;
00561 }
00562 if (sat_c_is_up_to_date())
00563 if (con_sys.first_pending_row() != sat_c.num_columns()) {
00564 #ifndef NDEBUG
00565 cerr << "Incompatible size! (con_sys and sat_c)"
00566 << endl;
00567 #endif
00568 goto bomb;
00569 }
00570 if (sat_g_is_up_to_date())
00571 if (con_sys.first_pending_row() != sat_g.num_rows()) {
00572 #ifndef NDEBUG
00573 cerr << "Incompatible size! (con_sys and sat_g)"
00574 << endl;
00575 #endif
00576 goto bomb;
00577 }
00578 if (generators_are_up_to_date())
00579 if (con_sys.num_columns() != gen_sys.num_columns()) {
00580 #ifndef NDEBUG
00581 cerr << "Incompatible size! (con_sys and gen_sys)"
00582 << endl;
00583 #endif
00584 goto bomb;
00585 }
00586 }
00587
00588 if (generators_are_up_to_date()) {
00589 if (gen_sys.num_columns() != poly_num_columns) {
00590 #ifndef NDEBUG
00591 cerr << "Incompatible size! (gen_sys and space_dim)"
00592 << endl;
00593 #endif
00594 goto bomb;
00595 }
00596 if (sat_c_is_up_to_date())
00597 if (gen_sys.first_pending_row() != sat_c.num_rows()) {
00598 #ifndef NDEBUG
00599 cerr << "Incompatible size! (gen_sys and sat_c)"
00600 << endl;
00601 #endif
00602 goto bomb;
00603 }
00604 if (sat_g_is_up_to_date())
00605 if (gen_sys.first_pending_row() != sat_g.num_columns()) {
00606 #ifndef NDEBUG
00607 cerr << "Incompatible size! (gen_sys and sat_g)"
00608 << endl;
00609 #endif
00610 goto bomb;
00611 }
00612
00613
00614 if (!gen_sys.OK())
00615 goto bomb;
00616
00617 if (gen_sys.first_pending_row() == 0) {
00618 #ifndef NDEBUG
00619 cerr << "Up-to-date generator system with all rows pending!"
00620 << endl;
00621 #endif
00622 goto bomb;
00623 }
00624
00625
00626
00627 if (gen_sys.num_rows() > 0 && !gen_sys.has_points()) {
00628 #ifndef NDEBUG
00629 cerr << "Non-empty generator system declared up-to-date "
00630 << "has no points!"
00631 << endl;
00632 #endif
00633 goto bomb;
00634 }
00635
00636 #if 0
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647 if (!is_necessarily_closed()) {
00648 dimension_type num_points = 0;
00649 dimension_type num_closure_points = 0;
00650 dimension_type eps_index = gen_sys.num_columns() - 1;
00651 for (dimension_type i = gen_sys.num_rows(); i-- > 0; )
00652 if (!gen_sys[i].is_line_or_ray())
00653 if (gen_sys[i][eps_index] > 0)
00654 ++num_points;
00655 else
00656 ++num_closure_points;
00657 if (num_points > num_closure_points) {
00658 #ifndef NDEBUG
00659 cerr << "# POINTS > # CLOSURE_POINTS" << endl;
00660 #endif
00661 goto bomb;
00662 }
00663 }
00664
00665 #endif
00666
00667 if (generators_are_minimized()) {
00668
00669
00670
00671
00672 Constraint_System new_con_sys(topology());
00673 Generator_System gs_without_pending = gen_sys;
00674 gs_without_pending.erase_to_end(gen_sys.first_pending_row());
00675 gs_without_pending.unset_pending_rows();
00676 Generator_System copy_of_gen_sys = gs_without_pending;
00677 Saturation_Matrix new_sat_c;
00678 minimize(false, copy_of_gen_sys, new_con_sys, new_sat_c);
00679 const dimension_type copy_num_lines = copy_of_gen_sys.num_lines();
00680 if (gs_without_pending.num_rows() != copy_of_gen_sys.num_rows()
00681 || gs_without_pending.num_lines() != copy_num_lines
00682 || gs_without_pending.num_rays() != copy_of_gen_sys.num_rays()) {
00683 #ifndef NDEBUG
00684 cerr << "Generators are declared minimized, but they are not!\n"
00685 << "Here is the minimized form of the generators:\n";
00686 copy_of_gen_sys.ascii_dump(cerr);
00687 cerr << endl;
00688 #endif
00689 goto bomb;
00690 }
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705 if (copy_num_lines == 0) {
00706 copy_of_gen_sys.strong_normalize();
00707 copy_of_gen_sys.sort_rows();
00708 gs_without_pending.strong_normalize();
00709 gs_without_pending.sort_rows();
00710 if (copy_of_gen_sys != gs_without_pending) {
00711 #ifndef NDEBUG
00712 cerr << "Generators are declared minimized, but they are not!\n"
00713 << "(we are in the case:\n"
00714 << "dimension of lineality space equal to 0)\n"
00715 << "Here is the minimized form of the generators:\n";
00716 copy_of_gen_sys.ascii_dump(cerr);
00717 cerr << endl;
00718 #endif
00719 goto bomb;
00720 }
00721 }
00722 }
00723 }
00724
00725 if (constraints_are_up_to_date()) {
00726
00727 if (!con_sys.OK())
00728 goto bomb;
00729
00730 if (con_sys.first_pending_row() == 0) {
00731 #ifndef NDEBUG
00732 cerr << "Up-to-date constraint system with all rows pending!"
00733 << endl;
00734 #endif
00735 goto bomb;
00736 }
00737
00738
00739
00740
00741
00742
00743
00744 bool no_positivity_constraint = true;
00745 for (dimension_type i = con_sys.num_rows(); i-- > 0; )
00746 if (con_sys[i].inhomogeneous_term() != 0) {
00747 no_positivity_constraint = false;
00748 break;
00749 }
00750 if (no_positivity_constraint) {
00751 #ifndef NDEBUG
00752 cerr << "Non-empty constraint system has no positivity constraint"
00753 << endl;
00754 #endif
00755 goto bomb;
00756 }
00757
00758 if (!is_necessarily_closed()) {
00759
00760
00761
00762 bool no_epsilon_geq_zero = true;
00763 const dimension_type eps_index = con_sys.num_columns() - 1;
00764 for (dimension_type i = con_sys.num_rows(); i-- > 0; )
00765 if (con_sys[i][eps_index] > 0) {
00766 no_epsilon_geq_zero = false;
00767 break;
00768 }
00769 if (no_epsilon_geq_zero) {
00770 #ifndef NDEBUG
00771 cerr << "Non-empty constraint system for NNC polyhedron "
00772 << "has no epsilon >= 0 constraint"
00773 << endl;
00774 #endif
00775 goto bomb;
00776 }
00777 }
00778
00779 Constraint_System cs_without_pending = con_sys;
00780 cs_without_pending.erase_to_end(con_sys.first_pending_row());
00781 cs_without_pending.unset_pending_rows();
00782 Constraint_System copy_of_con_sys = cs_without_pending;
00783 Generator_System new_gen_sys(topology());
00784 Saturation_Matrix new_sat_g;
00785
00786 if (minimize(true, copy_of_con_sys, new_gen_sys, new_sat_g)) {
00787 if (check_not_empty) {
00788
00789 #ifndef NDEBUG
00790 cerr << "Unsatisfiable system of constraints!"
00791 << endl;
00792 #endif
00793 goto bomb;
00794 }
00795
00796 return true;
00797 }
00798
00799 if (constraints_are_minimized()) {
00800
00801
00802
00803
00804 if (cs_without_pending.num_rows() != copy_of_con_sys.num_rows()
00805 || cs_without_pending.num_equalities()
00806 != copy_of_con_sys.num_equalities()) {
00807 #ifndef NDEBUG
00808 cerr << "Constraints are declared minimized, but they are not!\n"
00809 << "Here is the minimized form of the constraints:\n";
00810 copy_of_con_sys.ascii_dump(cerr);
00811 cerr << endl;
00812 #endif
00813 goto bomb;
00814 }
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824 copy_of_con_sys.strong_normalize();
00825 copy_of_con_sys.sort_rows();
00826 cs_without_pending.simplify();
00827 cs_without_pending.strong_normalize();
00828 cs_without_pending.sort_rows();
00829 if (cs_without_pending != copy_of_con_sys) {
00830 #ifndef NDEBUG
00831 cerr << "Constraints are declared minimized, but they are not!\n"
00832 << "Here is the minimized form of the constraints:\n";
00833 copy_of_con_sys.ascii_dump(cerr);
00834 cerr << endl;
00835 #endif
00836 goto bomb;
00837 }
00838 }
00839 }
00840
00841 if (sat_c_is_up_to_date())
00842 for (dimension_type i = sat_c.num_rows(); i-- > 0; ) {
00843 const Generator tmp_gen = gen_sys[i];
00844 const Saturation_Row tmp_sat = sat_c[i];
00845 for (dimension_type j = sat_c.num_columns(); j-- > 0; )
00846 if (Scalar_Products::sign(con_sys[j], tmp_gen) != tmp_sat[j]) {
00847 #ifndef NDEBUG
00848 cerr << "sat_c is declared up-to-date, but it is not!"
00849 << endl;
00850 #endif
00851 goto bomb;
00852 }
00853 }
00854
00855 if (sat_g_is_up_to_date())
00856 for (dimension_type i = sat_g.num_rows(); i-- > 0; ) {
00857 const Constraint tmp_con = con_sys[i];
00858 const Saturation_Row tmp_sat = sat_g[i];
00859 for (dimension_type j = sat_g.num_columns(); j-- > 0; )
00860 if (Scalar_Products::sign(tmp_con, gen_sys[j]) != tmp_sat[j]) {
00861 #ifndef NDEBUG
00862 cerr << "sat_g is declared up-to-date, but it is not!"
00863 << endl;
00864 #endif
00865 goto bomb;
00866 }
00867 }
00868
00869 if (has_pending_constraints()) {
00870 if (con_sys.num_pending_rows() == 0) {
00871 #ifndef NDEBUG
00872 cerr << "The polyhedron is declared to have pending constraints, "
00873 << "but con_sys has no pending rows!"
00874 << endl;
00875 #endif
00876 goto bomb;
00877 }
00878 }
00879
00880 if (has_pending_generators()) {
00881 if (gen_sys.num_pending_rows() == 0) {
00882 #ifndef NDEBUG
00883 cerr << "The polyhedron is declared to have pending generators, "
00884 << "but gen_sys has no pending rows!"
00885 << endl;
00886 #endif
00887 goto bomb;
00888 }
00889 }
00890
00891 return true;
00892
00893 bomb:
00894 #ifndef NDEBUG
00895 cerr << "Here is the guilty polyhedron:"
00896 << endl;
00897 ascii_dump(cerr);
00898 #endif
00899 return false;
00900 }
00901
00902 void
00903 PPL::Polyhedron::add_constraint(const Constraint& c) {
00904
00905 if (c.is_strict_inequality() && is_necessarily_closed())
00906 throw_topology_incompatible("add_constraint(c)", "c", c);
00907
00908
00909 if (space_dim < c.space_dimension())
00910 throw_dimension_incompatible("add_constraint(c)", "c", c);
00911
00912
00913
00914 if (marked_empty())
00915 return;
00916
00917
00918 if (space_dim == 0) {
00919 if (!c.is_tautological())
00920 set_empty();
00921 return;
00922 }
00923
00924
00925 if (has_pending_generators())
00926 process_pending_generators();
00927 else if (!constraints_are_up_to_date())
00928 update_constraints();
00929
00930 const bool adding_pending = can_have_something_pending();
00931
00932
00933 if (c.is_necessarily_closed() || !is_necessarily_closed())
00934
00935
00936 if (adding_pending)
00937 con_sys.insert_pending(c);
00938 else
00939 con_sys.insert(c);
00940 else {
00941
00942
00943
00944
00945
00946 Linear_Expression nc_expr = Linear_Expression(c);
00947 if (c.is_equality())
00948 if (adding_pending)
00949 con_sys.insert_pending(nc_expr == 0);
00950 else
00951 con_sys.insert(nc_expr == 0);
00952 else
00953 if (adding_pending)
00954 con_sys.insert_pending(nc_expr >= 0);
00955 else
00956 con_sys.insert(nc_expr >= 0);
00957 }
00958
00959 if (adding_pending)
00960 set_constraints_pending();
00961 else {
00962
00963 clear_constraints_minimized();
00964 clear_generators_up_to_date();
00965 }
00966
00967
00968 assert(OK());
00969 }
00970
00971 void
00972 PPL::Polyhedron::add_congruence(const Congruence& cg) {
00973
00974
00975 if (space_dim < cg.space_dimension())
00976 throw_dimension_incompatible("add_congruence(cg)", "cg", cg);
00977
00978
00979
00980 if (marked_empty())
00981 return;
00982
00983
00984 if (space_dim == 0) {
00985 if (!cg.is_trivial_true())
00986 set_empty();
00987 return;
00988 }
00989
00990 if (cg.is_equality()) {
00991 Linear_Expression le(cg);
00992 Constraint c(le, Constraint::EQUALITY, NECESSARILY_CLOSED);
00993 add_constraint(c);
00994 }
00995 }
00996
00997 bool
00998 PPL::Polyhedron::add_constraint_and_minimize(const Constraint& c) {
00999
01000 Constraint_System cs(c);
01001 return add_recycled_constraints_and_minimize(cs);
01002 }
01003
01004 void
01005 PPL::Polyhedron::add_generator(const Generator& g) {
01006
01007 if (g.is_closure_point() && is_necessarily_closed())
01008 throw_topology_incompatible("add_generator(g)", "g", g);
01009
01010
01011 const dimension_type g_space_dim = g.space_dimension();
01012 if (space_dim < g_space_dim)
01013 throw_dimension_incompatible("add_generator(g)", "g", g);
01014
01015
01016 if (space_dim == 0) {
01017
01018 assert(g.is_point() || g.is_closure_point());
01019
01020 if (marked_empty())
01021 if (g.type() != Generator::POINT)
01022 throw_invalid_generator("add_generator(g)", "g");
01023 else
01024 status.set_zero_dim_univ();
01025 assert(OK());
01026 return;
01027 }
01028
01029 if (marked_empty()
01030 || (has_pending_constraints() && !process_pending_constraints())
01031 || (!generators_are_up_to_date() && !update_generators())) {
01032
01033
01034 if (!g.is_point())
01035 throw_invalid_generator("add_generator(g)", "g");
01036 if (g.is_necessarily_closed() || !is_necessarily_closed()) {
01037 gen_sys.insert(g);
01038
01039
01040 gen_sys.adjust_topology_and_space_dimension(topology(), space_dim);
01041 if (!is_necessarily_closed()) {
01042
01043
01044
01045
01046 Generator& cp = gen_sys[gen_sys.num_rows() - 1];
01047 cp[space_dim + 1] = 0;
01048 cp.normalize();
01049
01050 gen_sys.insert(g);
01051 }
01052 }
01053 else {
01054
01055
01056
01057
01058
01059 const Linear_Expression nc_expr = Linear_Expression(g);
01060 gen_sys.insert(Generator::point(nc_expr, g.divisor()));
01061
01062
01063 gen_sys.adjust_topology_and_space_dimension(topology(), space_dim);
01064 }
01065
01066 clear_empty();
01067 set_generators_minimized();
01068 }
01069 else {
01070 assert(generators_are_up_to_date());
01071 const bool has_pending = can_have_something_pending();
01072 if (g.is_necessarily_closed() || !is_necessarily_closed()) {
01073
01074
01075 if (has_pending)
01076 gen_sys.insert_pending(g);
01077 else
01078 gen_sys.insert(g);
01079 if (!is_necessarily_closed() && g.is_point()) {
01080
01081
01082
01083
01084 Generator& cp = gen_sys[gen_sys.num_rows() - 1];
01085 cp[space_dim + 1] = 0;
01086 cp.normalize();
01087
01088 if (has_pending)
01089 gen_sys.insert_pending(g);
01090 else
01091 gen_sys.insert(g);
01092 }
01093 }
01094 else {
01095 assert(!g.is_closure_point());
01096
01097
01098
01099
01100
01101 const Linear_Expression nc_expr = Linear_Expression(g);
01102 switch (g.type()) {
01103 case Generator::LINE:
01104 if (has_pending)
01105 gen_sys.insert_pending(Generator::line(nc_expr));
01106 else
01107 gen_sys.insert(Generator::line(nc_expr));
01108 break;
01109 case Generator::RAY:
01110 if (has_pending)
01111 gen_sys.insert_pending(Generator::ray(nc_expr));
01112 else
01113 gen_sys.insert(Generator::ray(nc_expr));
01114 break;
01115 case Generator::POINT:
01116 if (has_pending)
01117 gen_sys.insert_pending(Generator::point(nc_expr, g.divisor()));
01118 else
01119 gen_sys.insert(Generator::point(nc_expr, g.divisor()));
01120 break;
01121 default:
01122 throw_runtime_error("add_generator(const Generator& g)");
01123 }
01124 }
01125
01126 if (has_pending)
01127 set_generators_pending();
01128 else {
01129
01130
01131 clear_generators_minimized();
01132 clear_constraints_up_to_date();
01133 }
01134 }
01135 assert(OK());
01136 }
01137
01138 bool
01139 PPL::Polyhedron::add_generator_and_minimize(const Generator& g) {
01140
01141 Generator_System gs(g);
01142 return add_recycled_generators_and_minimize(gs);
01143 }
01144
01145 void
01146 PPL::Polyhedron::add_recycled_constraints(Constraint_System& cs) {
01147
01148 if (is_necessarily_closed() && cs.has_strict_inequalities())
01149 throw_topology_incompatible("add_constraints(cs)", "cs", cs);
01150
01151
01152 const dimension_type cs_space_dim = cs.space_dimension();
01153 if (space_dim < cs_space_dim)
01154 throw_dimension_incompatible("add_recycled_constraints(cs)", "cs", cs);
01155
01156
01157 if (cs.num_rows() == 0)
01158 return;
01159
01160 if (space_dim == 0) {
01161
01162
01163
01164
01165
01166 if (cs.begin() != cs.end())
01167
01168
01169 status.set_empty();
01170 return;
01171 }
01172
01173 if (marked_empty())
01174 return;
01175
01176
01177 if (has_pending_generators())
01178 process_pending_generators();
01179 else if (!constraints_are_up_to_date())
01180 update_constraints();
01181
01182
01183
01184 cs.adjust_topology_and_space_dimension(topology(), space_dim);
01185
01186 const bool adding_pending = can_have_something_pending();
01187
01188
01189
01190
01191 const dimension_type old_num_rows = con_sys.num_rows();
01192 const dimension_type cs_num_rows = cs.num_rows();
01193 const dimension_type cs_num_columns = cs.num_columns();
01194 con_sys.add_zero_rows(cs_num_rows,
01195 Linear_Row::Flags(topology(),
01196 Linear_Row::RAY_OR_POINT_OR_INEQUALITY));
01197 for (dimension_type i = cs_num_rows; i-- > 0; ) {
01198
01199
01200
01201 Constraint& new_c = con_sys[old_num_rows + i];
01202 Constraint& old_c = cs[i];
01203 if (old_c.is_equality())
01204 new_c.set_is_equality();
01205 for (dimension_type j = cs_num_columns; j-- > 0; )
01206 std::swap(new_c[j], old_c[j]);
01207 }
01208
01209 if (adding_pending)
01210 set_constraints_pending();
01211 else {
01212
01213 con_sys.unset_pending_rows();
01214
01215 con_sys.set_sorted(false);
01216
01217 clear_constraints_minimized();
01218 clear_generators_up_to_date();
01219 }
01220
01221
01222 assert(OK());
01223 }
01224
01225 void
01226 PPL::Polyhedron::add_constraints(const Constraint_System& cs) {
01227
01228 Constraint_System cs_copy = cs;
01229 add_recycled_constraints(cs_copy);
01230 }
01231
01232 bool
01233 PPL::Polyhedron::add_recycled_constraints_and_minimize(Constraint_System& cs) {
01234
01235 if (is_necessarily_closed() && cs.has_strict_inequalities())
01236 throw_topology_incompatible("add_recycled_constraints_and_minimize(cs)",
01237 "cs", cs);
01238
01239
01240 const dimension_type cs_space_dim = cs.space_dimension();
01241 if (space_dim < cs_space_dim)
01242 throw_dimension_incompatible("add_recycled_constraints_and_minimize(cs)",
01243 "cs", cs);
01244
01245
01246 if (cs.num_rows() == 0)
01247 return minimize();
01248
01249
01250 if (space_dim == 0) {
01251
01252
01253
01254
01255
01256 if (cs.begin() == cs.end())
01257 return true;
01258
01259
01260 status.set_empty();
01261 return false;
01262 }
01263
01264
01265
01266 if (!minimize())
01267
01268 return false;
01269
01270 obtain_sorted_constraints_with_sat_c();
01271
01272
01273
01274 if (cs.num_pending_rows() > 0) {
01275 cs.unset_pending_rows();
01276 cs.sort_rows();
01277 }
01278 else if (!cs.is_sorted())
01279 cs.sort_rows();
01280
01281
01282 cs.adjust_topology_and_space_dimension(topology(), space_dim);
01283
01284 const bool empty = add_and_minimize(true, con_sys, gen_sys, sat_c, cs);
01285
01286 if (empty)
01287 set_empty();
01288 else {
01289
01290 set_sat_c_up_to_date();
01291 clear_sat_g_up_to_date();
01292 }
01293 assert(OK());
01294
01295 return !empty;
01296 }
01297
01298 bool
01299 PPL::Polyhedron::add_constraints_and_minimize(const Constraint_System& cs) {
01300
01301 Constraint_System cs_copy = cs;
01302 return add_recycled_constraints_and_minimize(cs_copy);
01303 }
01304
01305 void
01306 PPL::Polyhedron::add_recycled_generators(Generator_System& gs) {
01307
01308 if (is_necessarily_closed() && gs.has_closure_points())
01309 throw_topology_incompatible("add_recycled_generators(gs)", "gs", gs);
01310
01311
01312 const dimension_type gs_space_dim = gs.space_dimension();
01313 if (space_dim < gs_space_dim)
01314 throw_dimension_incompatible("add_recycled_generators(gs)", "gs", gs);
01315
01316
01317 if (gs.num_rows() == 0)
01318 return;
01319
01320
01321
01322 if (space_dim == 0) {
01323 if (marked_empty() && !gs.has_points())
01324 throw_invalid_generators("add_recycled_generators(gs)", "gs");
01325 status.set_zero_dim_univ();
01326 assert(OK(true));
01327 return;
01328 }
01329
01330
01331
01332 gs.adjust_topology_and_space_dimension(topology(), space_dim);
01333
01334
01335 if (!is_necessarily_closed())
01336 gs.add_corresponding_closure_points();
01337
01338
01339 if ((has_pending_constraints() && !process_pending_constraints())
01340 || (!generators_are_up_to_date() && !minimize())) {
01341
01342
01343 if (!gs.has_points())
01344 throw_invalid_generators("add_recycled_generators(gs)", "gs");
01345
01346 std::swap(gen_sys, gs);
01347 if (gen_sys.num_pending_rows() > 0) {
01348
01349
01350
01351
01352 gen_sys.unset_pending_rows();
01353 gen_sys.set_sorted(false);
01354 }
01355 set_generators_up_to_date();
01356 clear_empty();
01357 assert(OK());
01358 return;
01359 }
01360
01361 const bool adding_pending = can_have_something_pending();
01362
01363
01364
01365
01366 const dimension_type old_num_rows = gen_sys.num_rows();
01367 const dimension_type gs_num_rows = gs.num_rows();
01368 const dimension_type gs_num_columns = gs.num_columns();
01369 gen_sys.add_zero_rows(gs_num_rows,
01370 Linear_Row::Flags(topology(),
01371 Linear_Row::RAY_OR_POINT_OR_INEQUALITY));
01372 for (dimension_type i = gs_num_rows; i-- > 0; ) {
01373
01374
01375
01376 Generator& new_g = gen_sys[old_num_rows + i];
01377 Generator& old_g = gs[i];
01378 if (old_g.is_line())
01379 new_g.set_is_line();
01380 for (dimension_type j = gs_num_columns; j-- > 0; )
01381 std::swap(new_g[j], old_g[j]);
01382 }
01383
01384 if (adding_pending)
01385 set_generators_pending();
01386 else {
01387
01388 gen_sys.unset_pending_rows();
01389
01390 gen_sys.set_sorted(false);
01391
01392 clear_constraints_up_to_date();
01393 clear_generators_minimized();
01394 }
01395 assert(OK(true));
01396 }
01397
01398 void
01399 PPL::Polyhedron::add_generators(const Generator_System& gs) {
01400
01401 Generator_System gs_copy = gs;
01402 add_recycled_generators(gs_copy);
01403 }
01404
01405 bool
01406 PPL::Polyhedron::add_recycled_generators_and_minimize(Generator_System& gs) {
01407
01408 if (is_necessarily_closed() && gs.has_closure_points())
01409 throw_topology_incompatible("add_recycled_generators_and_minimize(gs)",
01410 "gs", gs);
01411
01412
01413 const dimension_type gs_space_dim = gs.space_dimension();
01414 if (space_dim < gs_space_dim)
01415 throw_dimension_incompatible("add_recycled_generators_and_minimize(gs)",
01416 "gs", gs);
01417
01418
01419 if (gs.num_rows() == 0)
01420 return minimize();
01421
01422
01423
01424 if (space_dim == 0) {
01425 if (marked_empty() && !gs.has_points())
01426 throw_invalid_generators("add_recycled_generators_and_minimize(gs)",
01427 "gs");
01428 status.set_zero_dim_univ();
01429 assert(OK(true));
01430 return true;
01431 }
01432
01433
01434
01435
01436
01437 gs.adjust_topology_and_space_dimension(topology(), gs_space_dim);
01438
01439
01440
01441 if (!is_necessarily_closed())
01442 gs.add_corresponding_closure_points();
01443
01444
01445 if (gs.num_pending_rows() > 0) {
01446 gs.unset_pending_rows();
01447 gs.sort_rows();
01448 }
01449 else if (!gs.is_sorted())
01450 gs.sort_rows();
01451
01452
01453
01454 gs.adjust_topology_and_space_dimension(topology(), space_dim);
01455
01456 if (minimize()) {
01457 obtain_sorted_generators_with_sat_g();
01458
01459 add_and_minimize(false, gen_sys, con_sys, sat_g, gs);
01460 clear_sat_c_up_to_date();
01461 }
01462 else {
01463
01464 if (!gs.has_points())
01465 throw_invalid_generators("add_recycled_generators_and_minimize(gs)",
01466 "gs");
01467
01468
01469 std::swap(gen_sys, gs);
01470 clear_empty();
01471 set_generators_up_to_date();
01472
01473 minimize();
01474 }
01475 assert(OK(true));
01476 return true;
01477 }
01478
01479 bool
01480 PPL::Polyhedron::add_generators_and_minimize(const Generator_System& gs) {
01481
01482 Generator_System gs_copy = gs;
01483 return add_recycled_generators_and_minimize(gs_copy);
01484 }
01485
01486 void
01487 PPL::Polyhedron::add_congruences(const Congruence_System& cgs) {
01488
01489
01490 if (space_dim < cgs.space_dimension())
01491 throw_dimension_incompatible("add_congruences(cgs)", "cgs", cgs);
01492
01493 Constraint_System cs;
01494 bool inserted = false;
01495 for (Congruence_System::const_iterator i = cgs.begin(),
01496 cgs_end = cgs.end(); i != cgs_end; ++i)
01497 if (i->is_equality()) {
01498 Linear_Expression le(*i);
01499 Constraint c(le, Constraint::EQUALITY, NECESSARILY_CLOSED);
01500
01501 cs.insert(c);
01502 inserted = true;
01503 }
01504
01505
01506 if (inserted)
01507 add_recycled_constraints(cs);
01508 }
01509
01510 void
01511 PPL::Polyhedron::intersection_assign(const Polyhedron& y) {
01512 Polyhedron& x = *this;
01513
01514 if (x.topology() != y.topology())
01515 throw_topology_incompatible("intersection_assign(y)", "y", y);
01516
01517 if (x.space_dim != y.space_dim)
01518 throw_dimension_incompatible("intersection_assign(y)", "y", y);
01519
01520
01521 if (x.marked_empty())
01522 return;
01523 if (y.marked_empty()) {
01524 x.set_empty();
01525 return;
01526 }
01527
01528
01529
01530
01531 if (x.space_dim == 0)
01532 return;
01533
01534
01535
01536 if (x.has_pending_generators())
01537 x.process_pending_generators();
01538 else if (!x.constraints_are_up_to_date())
01539 x.update_constraints();
01540
01541 if (y.has_pending_generators())
01542 y.process_pending_generators();
01543 else if (!y.constraints_are_up_to_date())
01544 y.update_constraints();
01545
01546
01547
01548 assert(!x.has_pending_generators() && x.constraints_are_up_to_date());
01549 assert(!y.has_pending_generators() && y.constraints_are_up_to_date());
01550
01551
01552
01553 if (x.can_have_something_pending()) {
01554 x.con_sys.add_pending_rows(y.con_sys);
01555 x.set_constraints_pending();
01556 }
01557 else {
01558
01559
01560
01561 if (x.con_sys.is_sorted()
01562 && y.con_sys.is_sorted() && !y.has_pending_constraints())
01563 x.con_sys.merge_rows_assign(y.con_sys);
01564 else
01565 x.con_sys.add_rows(y.con_sys);
01566
01567
01568 x.clear_generators_up_to_date();
01569 x.clear_constraints_minimized();
01570 }
01571 assert(x.OK() && y.OK());
01572 }
01573
01574 bool
01575 PPL::Polyhedron::intersection_assign_and_minimize(const Polyhedron& y) {
01576 Polyhedron& x = *this;
01577
01578 if (x.topology() != y.topology())
01579 throw_topology_incompatible("intersection_assign_and_minimize(y)",
01580 "y", y);
01581
01582 if (x.space_dim != y.space_dim)
01583 throw_dimension_incompatible("intersection_assign_and_minimize(y)",
01584 "y", y);
01585
01586
01587 if (x.marked_empty())
01588 return false;
01589 if (y.marked_empty()) {
01590 x.set_empty();
01591 return false;
01592 }
01593
01594
01595
01596
01597 if (x.space_dim == 0)
01598 return true;
01599
01600
01601
01602 if (!x.minimize())
01603
01604 return false;
01605 x.obtain_sorted_constraints_with_sat_c();
01606
01607
01608 if (y.has_pending_generators())
01609 y.process_pending_generators();
01610 else if (!y.constraints_are_up_to_date())
01611 y.update_constraints();
01612
01613 bool empty;
01614 if (y.con_sys.num_pending_rows() > 0) {
01615
01616
01617 x.con_sys.add_pending_rows(y.con_sys);
01618 x.con_sys.sort_pending_and_remove_duplicates();
01619 if (x.con_sys.num_pending_rows() == 0) {
01620
01621 x.clear_pending_constraints();
01622 assert(OK(true));
01623 return true;
01624 }
01625 empty = add_and_minimize(true, x.con_sys, x.gen_sys, x.sat_c);
01626 }
01627 else {
01628 y.obtain_sorted_constraints();
01629 empty = add_and_minimize(true, x.con_sys, x.gen_sys, x.sat_c, y.con_sys);
01630 }
01631
01632 if (empty)
01633 x.set_empty();
01634 else {
01635
01636
01637 x.set_sat_c_up_to_date();
01638 x.clear_sat_g_up_to_date();
01639 }
01640 assert(x.OK(!empty));
01641 return !empty;
01642 }
01643
01644 void
01645 PPL::Polyhedron::poly_hull_assign(const Polyhedron& y) {
01646 Polyhedron& x = *this;
01647
01648 if (x.topology() != y.topology())
01649 throw_topology_incompatible("poly_hull_assign(y)", "y", y);
01650
01651 if (x.space_dim != y.space_dim)
01652 throw_dimension_incompatible("poly_hull_assign(y)", "y", y);
01653
01654
01655 if (y.marked_empty())
01656 return;
01657 if (x.marked_empty()) {
01658 x = y;
01659 return;
01660 }
01661
01662
01663
01664
01665 if (x.space_dim == 0)
01666 return;
01667
01668
01669
01670 if ((x.has_pending_constraints() && !x.process_pending_constraints())
01671 || (!x.generators_are_up_to_date() && !x.update_generators())) {
01672
01673 x = y;
01674 return;
01675 }
01676 if ((y.has_pending_constraints() && !y.process_pending_constraints())
01677 || (!y.generators_are_up_to_date() && !y.update_generators()))
01678
01679 return;
01680
01681
01682
01683 assert(!x.has_pending_constraints() && x.generators_are_up_to_date());
01684 assert(!y.has_pending_constraints() && y.generators_are_up_to_date());
01685
01686
01687
01688 if (x.can_have_something_pending()) {
01689 x.gen_sys.add_pending_rows(y.gen_sys);
01690 x.set_generators_pending();
01691 }
01692 else {
01693
01694
01695
01696 if (x.gen_sys.is_sorted()
01697 && y.gen_sys.is_sorted() && !y.has_pending_generators())
01698 x.gen_sys.merge_rows_assign(y.gen_sys);
01699 else
01700 x.gen_sys.add_rows(y.gen_sys);
01701
01702
01703 x.clear_constraints_up_to_date();
01704 x.clear_generators_minimized();
01705 }
01706
01707 assert(x.OK(true) && y.OK(true));
01708 }
01709
01710 bool
01711 PPL::Polyhedron::poly_hull_assign_and_minimize(const Polyhedron& y) {
01712 Polyhedron& x = *this;
01713
01714 if (x.topology() != y.topology())
01715 throw_topology_incompatible("poly_hull_assign_and_minimize(y)", "y", y);
01716
01717 if (x.space_dim != y.space_dim)
01718 throw_dimension_incompatible("poly_hull_assign_and_minimize(y)", "y", y);
01719
01720
01721 if (y.marked_empty())
01722 return minimize();
01723 if (x.marked_empty()) {
01724 x = y;
01725 return minimize();
01726 }
01727
01728
01729
01730
01731 if (x.space_dim == 0)
01732 return true;
01733
01734
01735
01736 if (!x.minimize()) {
01737
01738 x = y;
01739 return minimize();
01740 }
01741
01742 x.obtain_sorted_generators_with_sat_g();
01743
01744
01745 if ((y.has_pending_constraints() && !y.process_pending_constraints())
01746 || (!y.generators_are_up_to_date() && !y.update_generators()))
01747
01748
01749 return true;
01750
01751 if (y.gen_sys.num_pending_rows() > 0) {
01752
01753
01754 x.gen_sys.add_pending_rows(y.gen_sys);
01755 x.gen_sys.sort_pending_and_remove_duplicates();
01756 if (x.gen_sys.num_pending_rows() == 0) {
01757
01758 x.clear_pending_generators();
01759 assert(OK(true) && y.OK());
01760 return true;
01761 }
01762 add_and_minimize(false, x.gen_sys, x.con_sys, x.sat_g);
01763 }
01764 else {
01765 y.obtain_sorted_generators();
01766 add_and_minimize(false, x.gen_sys, x.con_sys, x.sat_g, y.gen_sys);
01767 }
01768 x.clear_sat_c_up_to_date();
01769
01770 assert(x.OK(true) && y.OK());
01771 return true;
01772 }
01773
01774 void
01775 PPL::Polyhedron::poly_difference_assign(const Polyhedron& y) {
01776 Polyhedron& x = *this;
01777
01778 if (x.topology() != y.topology())
01779 throw_topology_incompatible("poly_difference_assign(y)", "y", y);
01780
01781 if (x.space_dim != y.space_dim)
01782 throw_dimension_incompatible("poly_difference_assign(y)", "y", y);
01783
01784
01785 if (y.marked_empty())
01786 return;
01787
01788 if (x.marked_empty())
01789 return;
01790
01791
01792
01793
01794 if (x.space_dim == 0) {
01795 x.set_empty();
01796 return;
01797 }
01798
01799
01800
01801
01802 if (y.contains(x)) {
01803 x.set_empty();
01804 return;
01805 }
01806
01807 Polyhedron new_polyhedron(topology(), x.space_dim, EMPTY);
01808
01809
01810
01811 x.minimize();
01812 y.minimize();
01813
01814 const Constraint_System& y_cs = y.constraints();
01815 for (Constraint_System::const_iterator i = y_cs.begin(),
01816 y_cs_end = y_cs.end(); i != y_cs_end; ++i) {
01817 const Constraint& c = *i;
01818 assert(!c.is_tautological());
01819 assert(!c.is_inconsistent());
01820
01821
01822
01823
01824
01825
01826 if (x.relation_with(c).implies(Poly_Con_Relation::is_included()))
01827 continue;
01828 Polyhedron z = x;
01829 const Linear_Expression e = Linear_Expression(c);
01830 switch (c.type()) {
01831 case Constraint::NONSTRICT_INEQUALITY:
01832 if (is_necessarily_closed())
01833 z.add_constraint(e <= 0);
01834 else
01835 z.add_constraint(e < 0);
01836 break;
01837 case Constraint::STRICT_INEQUALITY:
01838 z.add_constraint(e <= 0);
01839 break;
01840 case Constraint::EQUALITY:
01841 if (is_necessarily_closed())
01842
01843
01844 return;
01845 else {
01846 Polyhedron w = x;
01847 w.add_constraint(e < 0);
01848 new_polyhedron.poly_hull_assign(w);
01849 z.add_constraint(e > 0);
01850 }
01851 break;
01852 }
01853 new_polyhedron.poly_hull_assign(z);
01854 }
01855 *this = new_polyhedron;
01856
01857 assert(OK());
01858 }
01859
01860 void
01861 PPL::Polyhedron::
01862 affine_image(const Variable var,
01863 const Linear_Expression& expr,
01864 Coefficient_traits::const_reference denominator) {
01865
01866 if (denominator == 0)
01867 throw_invalid_argument("affine_image(v, e, d)", "d == 0");
01868
01869
01870
01871
01872 if (space_dim < expr.space_dimension())
01873 throw_dimension_incompatible("affine_image(v, e, d)", "e", expr);
01874
01875 const dimension_type var_space_dim = var.space_dimension();
01876 if (space_dim < var_space_dim)
01877 throw_dimension_incompatible("affine_image(v, e, d)", "v", var);
01878
01879 if (marked_empty())
01880 return;
01881
01882 if (expr.coefficient(var) != 0) {
01883
01884
01885
01886 if (generators_are_up_to_date()) {
01887
01888
01889 if (denominator > 0)
01890 gen_sys.affine_image(var_space_dim, expr, denominator);
01891 else
01892 gen_sys.affine_image(var_space_dim, -expr, -denominator);
01893 }
01894 if (constraints_are_up_to_date()) {
01895
01896
01897
01898 Linear_Expression inverse;
01899 if (expr[var_space_dim] > 0) {
01900 inverse = -expr;
01901 inverse[var_space_dim] = denominator;
01902 con_sys.affine_preimage(var_space_dim, inverse, expr[var_space_dim]);
01903 }
01904 else {
01905
01906
01907
01908 inverse = expr;
01909 inverse[var_space_dim] = denominator;
01910 neg_assign(inverse[var_space_dim]);
01911 con_sys.affine_preimage(var_space_dim, inverse, -expr[var_space_dim]);
01912 }
01913 }
01914 }
01915 else {
01916
01917
01918 if (has_something_pending())
01919 remove_pending_to_obtain_generators();
01920 else if (!generators_are_up_to_date())
01921 minimize();
01922 if (!marked_empty()) {
01923
01924
01925 if (denominator > 0)
01926 gen_sys.affine_image(var_space_dim, expr, denominator);
01927 else
01928 gen_sys.affine_image(var_space_dim, -expr, -denominator);
01929
01930 clear_constraints_up_to_date();
01931 clear_generators_minimized();
01932 clear_sat_c_up_to_date();
01933 clear_sat_g_up_to_date();
01934 }
01935 }
01936 assert(OK());
01937 }
01938
01939
01940 void
01941 PPL::Polyhedron::
01942 affine_preimage(const Variable var,
01943 const Linear_Expression& expr,
01944 Coefficient_traits::const_reference denominator) {
01945
01946 if (denominator == 0)
01947 throw_invalid_argument("affine_preimage(v, e, d)", "d == 0");
01948
01949
01950
01951
01952 if (space_dim < expr.space_dimension())
01953 throw_dimension_incompatible("affine_preimage(v, e, d)", "e", expr);
01954
01955 const dimension_type var_space_dim = var.space_dimension();
01956 if (space_dim < var_space_dim)
01957 throw_dimension_incompatible("affine_preimage(v, e, d)", "v", var);
01958
01959 if (marked_empty())
01960 return;
01961
01962 if (expr.coefficient(var) != 0) {
01963
01964
01965 if (constraints_are_up_to_date()) {
01966
01967
01968 if (denominator > 0)
01969 con_sys.affine_preimage(var_space_dim, expr, denominator);
01970 else
01971 con_sys.affine_preimage(var_space_dim, -expr, -denominator);
01972 }
01973 if (generators_are_up_to_date()) {
01974
01975
01976
01977 Linear_Expression inverse;
01978 if (expr[var_space_dim] > 0) {
01979 inverse = -expr;
01980 inverse[var_space_dim] = denominator;
01981 gen_sys.affine_image(var_space_dim, inverse, expr[var_space_dim]);
01982 }
01983 else {
01984
01985
01986
01987 inverse = expr;
01988 inverse[var_space_dim] = denominator;
01989 neg_assign(inverse[var_space_dim]);
01990 gen_sys.affine_image(var_space_dim, inverse, -expr[var_space_dim]);
01991 }
01992 }
01993 }
01994 else {
01995
01996
01997 if (has_something_pending())
01998 remove_pending_to_obtain_constraints();
01999 else if (!constraints_are_up_to_date())
02000 minimize();
02001
02002
02003 if (denominator > 0)
02004 con_sys.affine_preimage(var_space_dim, expr, denominator);
02005 else
02006 con_sys.affine_preimage(var_space_dim, -expr, -denominator);
02007
02008 clear_generators_up_to_date();
02009 clear_constraints_minimized();
02010 clear_sat_c_up_to_date();
02011 clear_sat_g_up_to_date();
02012 }
02013 assert(OK());
02014 }
02015
02016 void
02017 PPL::Polyhedron::
02018 bounded_affine_image(const Variable var,
02019 const Linear_Expression& lb_expr,
02020 const Linear_Expression& ub_expr,
02021 Coefficient_traits::const_reference denominator) {
02022
02023 if (denominator == 0)
02024 throw_invalid_argument("bounded_affine_image(v, lb, ub, d)", "d == 0");
02025
02026
02027
02028 const dimension_type var_space_dim = var.space_dimension();
02029 if (space_dim < var_space_dim)
02030 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
02031 "v", var);
02032
02033
02034 const dimension_type lb_space_dim = lb_expr.space_dimension();
02035 if (space_dim < lb_space_dim)
02036 throw_dimension_incompatible("bounded_affine_image(v, lb, ub)",
02037 "lb", lb_expr);
02038 const dimension_type ub_space_dim = ub_expr.space_dimension();
02039 if (space_dim < ub_space_dim)
02040 throw_dimension_incompatible("bounded_affine_image(v, lb, ub)",
02041 "ub", ub_expr);
02042
02043
02044 if (marked_empty())
02045 return;
02046
02047
02048 if (lb_expr.coefficient(var) == 0) {
02049
02050 generalized_affine_image(var,
02051 LESS_THAN_OR_EQUAL,
02052 ub_expr,
02053 denominator);
02054 if (denominator > 0)
02055 add_constraint(lb_expr <= denominator*var);
02056 else
02057 add_constraint(denominator*var <= lb_expr);
02058 }
02059 else if (ub_expr.coefficient(var) == 0) {
02060
02061 generalized_affine_image(var,
02062 GREATER_THAN_OR_EQUAL,
02063 lb_expr,
02064 denominator);
02065 if (denominator > 0)
02066 add_constraint(denominator*var <= ub_expr);
02067 else
02068 add_constraint(ub_expr <= denominator*var);
02069 }
02070 else {
02071
02072
02073 const Variable new_var = Variable(space_dim);
02074 add_space_dimensions_and_embed(1);
02075
02076
02077 add_constraint_and_minimize(denominator*new_var == ub_expr);
02078
02079 generalized_affine_image(var,
02080 GREATER_THAN_OR_EQUAL,
02081 lb_expr,
02082 denominator);
02083
02084
02085 if (denominator > 0)
02086 add_constraint_and_minimize(var <= new_var);
02087 else
02088 add_constraint_and_minimize(new_var <= var);
02089
02090 remove_higher_space_dimensions(space_dim-1);
02091 }
02092 assert(OK());
02093 }
02094
02095 void
02096 PPL::Polyhedron::
02097 bounded_affine_preimage(const Variable var,
02098 const Linear_Expression& lb_expr,
02099 const Linear_Expression& ub_expr,
02100 Coefficient_traits::const_reference denominator) {
02101
02102 if (denominator == 0)
02103 throw_invalid_argument("bounded_affine_preimage(v, lb, ub, d)", "d == 0");
02104
02105
02106
02107 const dimension_type var_space_dim = var.space_dimension();
02108 if (space_dim < var_space_dim)
02109 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
02110 "v", var);
02111
02112
02113 const dimension_type lb_space_dim = lb_expr.space_dimension();
02114 if (space_dim < lb_space_dim)
02115 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub)",
02116 "lb", lb_expr);
02117 const dimension_type ub_space_dim = ub_expr.space_dimension();
02118 if (space_dim < ub_space_dim)
02119 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub)",
02120 "ub", ub_expr);
02121
02122
02123 if (marked_empty())
02124 return;
02125
02126
02127 if (lb_expr.coefficient(var) == 0 && ub_expr.coefficient(var) == 0) {
02128 if (denominator > 0) {
02129 add_constraint(lb_expr <= denominator*var);
02130 add_constraint(denominator*var <= ub_expr);
02131 }
02132 else {
02133 add_constraint(ub_expr <= denominator*var);
02134 add_constraint(denominator*var <= lb_expr);
02135 }
02136
02137
02138 if (is_empty())
02139 return;
02140 add_generator(line(var));
02141 }
02142 else {
02143
02144
02145 const Variable new_var = Variable(space_dim);
02146 add_space_dimensions_and_embed(1);
02147
02148 std::vector<dimension_type> swapping_cycle;
02149 swapping_cycle.push_back(var_space_dim);
02150 swapping_cycle.push_back(space_dim);
02151 swapping_cycle.push_back(0);
02152 if (constraints_are_up_to_date())
02153 con_sys.permute_columns(swapping_cycle);
02154 if (generators_are_up_to_date())
02155 gen_sys.permute_columns(swapping_cycle);
02156
02157
02158 if (denominator > 0) {
02159 add_constraint(lb_expr <= denominator*new_var);
02160 add_constraint(denominator*new_var <= ub_expr);
02161 }
02162 else {
02163 add_constraint(ub_expr <= denominator*new_var);
02164 add_constraint(denominator*new_var <= lb_expr);
02165 }
02166
02167 remove_higher_space_dimensions(space_dim-1);
02168 }
02169 assert(OK());
02170 }
02171
02172 void
02173 PPL::Polyhedron::
02174 generalized_affine_image(const Variable var,
02175 const Relation_Symbol relsym,
02176 const Linear_Expression& expr,
02177 Coefficient_traits::const_reference denominator) {
02178
02179 if (denominator == 0)
02180 throw_invalid_argument("generalized_affine_image(v, r, e, d)", "d == 0");
02181
02182
02183
02184
02185 if (space_dim < expr.space_dimension())
02186 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
02187 "e", expr);
02188
02189 const dimension_type var_space_dim = var.space_dimension();
02190 if (space_dim < var_space_dim)
02191 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
02192 "v", var);
02193
02194
02195 if (is_necessarily_closed()
02196 && (relsym == LESS_THAN || relsym == GREATER_THAN))
02197 throw_invalid_argument("generalized_affine_image(v, r, e, d)",
02198 "r is a strict relation symbol");
02199
02200
02201 affine_image(var, expr, denominator);
02202
02203 if (relsym == EQUAL)
02204
02205 return;
02206
02207
02208
02209 if (is_empty())
02210 return;
02211
02212 switch (relsym) {
02213 case LESS_THAN_OR_EQUAL:
02214 add_generator(ray(-var));
02215 break;
02216 case GREATER_THAN_OR_EQUAL:
02217 add_generator(ray(var));
02218 break;
02219 case LESS_THAN:
02220
02221 case GREATER_THAN:
02222 {
02223
02224 assert(!is_necessarily_closed());
02225
02226
02227 add_generator_and_minimize(ray(relsym == GREATER_THAN ? var : -var));
02228
02229
02230
02231
02232
02233 const dimension_type eps_index = space_dim + 1;
02234 for (dimension_type i = gen_sys.num_rows(); i-- > 0; )
02235 if (gen_sys[i].is_point()) {
02236 Generator& g = gen_sys[i];
02237
02238 gen_sys.add_row(g);
02239 if (relsym == GREATER_THAN)
02240 ++gen_sys[gen_sys.num_rows()-1][var_space_dim];
02241 else
02242 --gen_sys[gen_sys.num_rows()-1][var_space_dim];
02243
02244 g[eps_index] = 0;
02245 }
02246 clear_constraints_up_to_date();
02247 clear_generators_minimized();
02248 gen_sys.set_sorted(false);
02249 clear_sat_c_up_to_date();
02250 clear_sat_g_up_to_date();
02251 }
02252 break;
02253 case EQUAL:
02254
02255 throw std::runtime_error("PPL internal error");
02256 }
02257 assert(OK());
02258 }
02259
02260 void
02261 PPL::Polyhedron::
02262 generalized_affine_preimage(const Variable var,
02263 const Relation_Symbol relsym,
02264 const Linear_Expression& expr,
02265 Coefficient_traits::const_reference denominator) {
02266
02267 if (denominator == 0)
02268 throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
02269 "d == 0");
02270
02271
02272
02273
02274 if (space_dim < expr.space_dimension())
02275 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
02276 "e", expr);
02277
02278 const dimension_type var_space_dim = var.space_dimension();
02279 if (space_dim < var_space_dim)
02280 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
02281 "v", var);
02282
02283
02284 if (is_necessarily_closed()
02285 && (relsym == LESS_THAN || relsym == GREATER_THAN))
02286 throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
02287 "r is a strict relation symbol");
02288
02289
02290 if (relsym == EQUAL) {
02291 affine_preimage(var, expr, denominator);
02292 return;
02293 }
02294
02295
02296 Relation_Symbol reversed_relsym;
02297 switch (relsym) {
02298 case LESS_THAN:
02299 reversed_relsym = GREATER_THAN;
02300 break;
02301 case LESS_THAN_OR_EQUAL:
02302 reversed_relsym = GREATER_THAN_OR_EQUAL;
02303 break;
02304 case GREATER_THAN_OR_EQUAL:
02305 reversed_relsym = LESS_THAN_OR_EQUAL;
02306 break;
02307 case GREATER_THAN:
02308 reversed_relsym = LESS_THAN;
02309 break;
02310 default:
02311
02312 throw std::runtime_error("PPL internal error");
02313 break;
02314 }
02315
02316
02317
02318 const Coefficient& var_coefficient = expr.coefficient(var);
02319 if (var_coefficient != 0) {
02320 Linear_Expression inverse_expr
02321 = expr - (denominator + var_coefficient) * var;
02322 Coefficient inverse_denominator = - var_coefficient;
02323 Relation_Symbol inverse_relsym
02324 = (sgn(denominator) == sgn(inverse_denominator))
02325 ? relsym : reversed_relsym;
02326 generalized_affine_image(var, inverse_relsym, inverse_expr,
02327 inverse_denominator);
02328 return;
02329 }
02330
02331
02332
02333
02334
02335 const Relation_Symbol corrected_relsym
02336 = (denominator > 0) ? relsym : reversed_relsym;
02337 switch (corrected_relsym) {
02338 case LESS_THAN:
02339 add_constraint(denominator*var < expr);
02340 break;
02341 case LESS_THAN_OR_EQUAL:
02342 add_constraint(denominator*var <= expr);
02343 break;
02344 case GREATER_THAN_OR_EQUAL:
02345 add_constraint(denominator*var >= expr);
02346 break;
02347 case GREATER_THAN:
02348 add_constraint(denominator*var > expr);
02349 break;
02350 case EQUAL:
02351
02352 throw std::runtime_error("PPL internal error");
02353 break;
02354 }
02355
02356
02357 if (is_empty())
02358 return;
02359 add_generator(line(var));
02360 assert(OK());
02361 }
02362
02363 void
02364 PPL::Polyhedron::generalized_affine_image(const Linear_Expression& lhs,
02365 const Relation_Symbol relsym,
02366 const Linear_Expression& rhs) {
02367
02368
02369
02370 dimension_type lhs_space_dim = lhs.space_dimension();
02371 if (space_dim < lhs_space_dim)
02372 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
02373 "e1", lhs);
02374
02375
02376 const dimension_type rhs_space_dim = rhs.space_dimension();
02377 if (space_dim < rhs_space_dim)
02378 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
02379 "e2", rhs);
02380
02381
02382 if (is_necessarily_closed()
02383 && (relsym == LESS_THAN || relsym == GREATER_THAN))
02384 throw_invalid_argument("generalized_affine_image(e1, r, e2)",
02385 "r is a strict relation symbol");
02386
02387
02388 if (marked_empty())
02389 return;
02390
02391
02392
02393 for ( ; lhs_space_dim > 0; lhs_space_dim--)
02394 if (lhs.coefficient(Variable(lhs_space_dim - 1)) != 0)
02395 break;
02396
02397
02398 if (lhs_space_dim == 0) {
02399 switch (relsym) {
02400 case LESS_THAN:
02401 add_constraint(lhs < rhs);
02402 break;
02403 case LESS_THAN_OR_EQUAL:
02404 add_constraint(lhs <= rhs);
02405 break;
02406 case EQUAL:
02407 add_constraint(lhs == rhs);
02408 break;
02409 case GREATER_THAN_OR_EQUAL:
02410 add_constraint(lhs >= rhs);
02411 break;
02412 case GREATER_THAN:
02413 add_constraint(lhs > rhs);
02414 break;
02415 }
02416 return;
02417 }
02418
02419
02420
02421
02422
02423 Generator_System new_lines;
02424 bool lhs_vars_intersects_rhs_vars = false;
02425 for (dimension_type i = lhs_space_dim; i-- > 0; )
02426 if (lhs.coefficient(Variable(i)) != 0) {
02427 new_lines.insert(line(Variable(i)));
02428 if (rhs.coefficient(Variable(i)) != 0)
02429 lhs_vars_intersects_rhs_vars = true;
02430 }
02431
02432 if (lhs_vars_intersects_rhs_vars) {
02433
02434
02435 const Variable new_var = Variable(space_dim);
02436 add_space_dimensions_and_embed(1);
02437
02438
02439
02440 if (add_constraint_and_minimize(new_var == rhs)) {
02441
02442
02443 add_recycled_generators_and_minimize(new_lines);
02444
02445
02446
02447
02448 switch (relsym) {
02449 case LESS_THAN:
02450 add_constraint_and_minimize(lhs < new_var);
02451 break;
02452 case LESS_THAN_OR_EQUAL:
02453 add_constraint_and_minimize(lhs <= new_var);
02454 break;
02455 case EQUAL:
02456 add_constraint_and_minimize(lhs == new_var);
02457 break;
02458 case GREATER_THAN_OR_EQUAL:
02459 add_constraint_and_minimize(lhs >= new_var);
02460 break;
02461 case GREATER_THAN:
02462 add_constraint_and_minimize(lhs > new_var);
02463 break;
02464 }
02465 }
02466
02467 remove_higher_space_dimensions(space_dim-1);
02468 }
02469 else {
02470
02471
02472
02473
02474
02475 if (is_empty())
02476 return;
02477
02478
02479
02480 add_recycled_generators_and_minimize(new_lines);
02481
02482
02483
02484 switch (relsym) {
02485 case LESS_THAN:
02486 add_constraint(lhs < rhs);
02487 break;
02488 case LESS_THAN_OR_EQUAL:
02489 add_constraint(lhs <= rhs);
02490 break;
02491 case EQUAL:
02492 add_constraint(lhs == rhs);
02493 break;
02494 case GREATER_THAN_OR_EQUAL:
02495 add_constraint(lhs >= rhs);
02496 break;
02497 case GREATER_THAN:
02498 add_constraint(lhs > rhs);
02499 break;
02500 }
02501 }
02502 assert(OK());
02503 }
02504
02505 void
02506 PPL::Polyhedron::generalized_affine_preimage(const Linear_Expression& lhs,
02507 const Relation_Symbol relsym,
02508 const Linear_Expression& rhs) {
02509
02510
02511
02512 dimension_type lhs_space_dim = lhs.space_dimension();
02513 if (space_dim < lhs_space_dim)
02514 throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
02515 "e1", lhs);
02516
02517
02518 const dimension_type rhs_space_dim = rhs.space_dimension();
02519 if (space_dim < rhs_space_dim)
02520 throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
02521 "e2", rhs);
02522
02523
02524 if (is_necessarily_closed()
02525 && (relsym == LESS_THAN || relsym == GREATER_THAN))
02526 throw_invalid_argument("generalized_affine_preimage(e1, r, e2)",
02527 "r is a strict relation symbol");
02528
02529
02530 if (marked_empty())
02531 return;
02532
02533
02534
02535 for ( ; lhs_space_dim > 0; lhs_space_dim--)
02536 if (lhs.coefficient(Variable(lhs_space_dim - 1)) != 0)
02537 break;
02538
02539
02540
02541 if (lhs_space_dim == 0) {
02542 generalized_affine_image(lhs, relsym, rhs);
02543 return;
02544 }
02545
02546
02547
02548
02549
02550 Generator_System new_lines;
02551 bool lhs_vars_intersects_rhs_vars = false;
02552 for (dimension_type i = lhs_space_dim; i-- > 0; )
02553 if (lhs.coefficient(Variable(i)) != 0) {
02554 new_lines.insert(line(Variable(i)));
02555 if (rhs.coefficient(Variable(i)) != 0)
02556 lhs_vars_intersects_rhs_vars = true;
02557 }
02558
02559 if (lhs_vars_intersects_rhs_vars) {
02560
02561
02562 const Variable new_var = Variable(space_dim);
02563 add_space_dimensions_and_embed(1);
02564
02565
02566
02567 if (add_constraint_and_minimize(new_var == lhs)) {
02568
02569
02570 add_recycled_generators_and_minimize(new_lines);
02571
02572
02573
02574
02575 switch (relsym) {
02576 case LESS_THAN:
02577 add_constraint_and_minimize(new_var < rhs);
02578 break;
02579 case LESS_THAN_OR_EQUAL:
02580 add_constraint_and_minimize(new_var <= rhs);
02581 break;
02582 case EQUAL:
02583 add_constraint_and_minimize(new_var == rhs);
02584 break;
02585 case GREATER_THAN_OR_EQUAL:
02586 add_constraint_and_minimize(new_var >= rhs);
02587 break;
02588 case GREATER_THAN:
02589 add_constraint_and_minimize(new_var > rhs);
02590 break;
02591 }
02592 }
02593
02594 remove_higher_space_dimensions(space_dim-1);
02595 }
02596 else {
02597
02598
02599
02600
02601
02602 switch (relsym) {
02603 case LESS_THAN:
02604 add_constraint(lhs < rhs);
02605 break;
02606 case LESS_THAN_OR_EQUAL:
02607 add_constraint(lhs <= rhs);
02608 break;
02609 case EQUAL:
02610 add_constraint(lhs == rhs);
02611 break;
02612 case GREATER_THAN_OR_EQUAL:
02613 add_constraint(lhs >= rhs);
02614 break;
02615 case GREATER_THAN:
02616 add_constraint(lhs > rhs);
02617 break;
02618 }
02619
02620
02621 if (is_empty())
02622 return;
02623
02624 add_recycled_generators(new_lines);
02625 }
02626 assert(OK());
02627 }
02628
02629 void
02630 PPL::Polyhedron::time_elapse_assign(const Polyhedron& y) {
02631 Polyhedron& x = *this;
02632
02633 if (x.topology() != y.topology())
02634 throw_topology_incompatible("time_elapse_assign(y)", "y", y);
02635
02636 if (x.space_dim != y.space_dim)
02637 throw_dimension_incompatible("time_elapse_assign(y)", "y", y);
02638
02639
02640 if (x.space_dim == 0) {
02641 if (y.marked_empty())
02642 x.set_empty();
02643 return;
02644 }
02645
02646
02647 if (x.marked_empty() || y.marked_empty()
02648 || (x.has_pending_constraints() && !x.process_pending_constraints())
02649 || (!x.generators_are_up_to_date() && !x.update_generators())
02650 || (y.has_pending_constraints() && !y.process_pending_constraints())
02651 || (!y.generators_are_up_to_date() && !y.update_generators())) {
02652 x.set_empty();
02653 return;
02654 }
02655
02656
02657
02658 Generator_System gs = y.gen_sys;
02659 dimension_type gs_num_rows = gs.num_rows();
02660
02661 if (!x.is_necessarily_closed())
02662
02663 for (dimension_type i = gs_num_rows; i-- > 0; )
02664 switch (gs[i].type()) {
02665 case Generator::POINT:
02666
02667
02668 --gs_num_rows;
02669 std::swap(gs[i], gs[gs_num_rows]);
02670 break;
02671 case Generator::CLOSURE_POINT:
02672 {
02673 Generator& cp = gs[i];
02674
02675 if (cp.all_homogeneous_terms_are_zero()) {
02676 --gs_num_rows;
02677 std::swap(cp, gs[gs_num_rows]);
02678 }
02679
02680 else {
02681 cp[0] = 0;
02682
02683 cp.normalize();
02684 }
02685 }
02686 break;
02687 default:
02688
02689 break;
02690 }
02691 else
02692
02693 for (dimension_type i = gs_num_rows; i-- > 0; )
02694 switch (gs[i].type()) {
02695 case Generator::POINT:
02696 {
02697 Generator& p = gs[i];
02698
02699 if (p.all_homogeneous_terms_are_zero()) {
02700 --gs_num_rows;
02701 std::swap(p, gs[gs_num_rows]);
02702 }
02703
02704 else {
02705 p[0] = 0;
02706
02707 p.normalize();
02708 }
02709 }
02710 break;
02711 default:
02712
02713 break;
02714 }
02715
02716
02717
02718
02719
02720 gs.erase_to_end(gs_num_rows);
02721 gs.unset_pending_rows();
02722
02723
02724
02725
02726
02727 if (gs_num_rows == 0)
02728 return;
02729
02730
02731
02732 if (x.can_have_something_pending()) {
02733 x.gen_sys.add_pending_rows(gs);
02734 x.set_generators_pending();
02735 }
02736
02737
02738 else {
02739 if (!x.gen_sys.is_sorted())
02740 x.gen_sys.sort_rows();
02741 gs.sort_rows();
02742 x.gen_sys.merge_rows_assign(gs);
02743
02744 x.clear_constraints_up_to_date();
02745 x.clear_generators_minimized();
02746 }
02747 assert(x.OK(true) && y.OK(true));
02748 }
02749
02750 void
02751 PPL::Polyhedron::topological_closure_assign() {
02752
02753 if (is_necessarily_closed())
02754 return;
02755
02756 if (marked_empty() || space_dim == 0)
02757 return;
02758
02759
02760
02761
02762
02763
02764
02765 if (has_pending_constraints() && !process_pending_constraints())
02766 return;
02767
02768
02769
02770 if (!has_pending_generators() && constraints_are_up_to_date()) {
02771 const dimension_type eps_index = space_dim + 1;
02772 bool changed = false;
02773
02774 for (dimension_type i = con_sys.num_rows(); i-- > 0; ) {
02775 Constraint& c = con_sys[i];
02776 if (c[eps_index] < 0 && !c.is_tautological()) {
02777 c[eps_index] = 0;
02778
02779 c.normalize();
02780 changed = true;
02781 }
02782 }
02783 if (changed) {
02784 con_sys.insert(Constraint::epsilon_leq_one());
02785 con_sys.set_sorted(false);
02786
02787
02788
02789 clear_generators_up_to_date();
02790 clear_constraints_minimized();
02791 }
02792 }
02793 else {
02794
02795 assert(generators_are_up_to_date());
02796
02797 gen_sys.add_corresponding_points();
02798 if (can_have_something_pending())
02799 set_generators_pending();
02800 else {
02801
02802
02803 gen_sys.unset_pending_rows();
02804 gen_sys.set_sorted(false);
02805
02806 clear_constraints_up_to_date();
02807 clear_generators_minimized();
02808 }
02809 }
02810 assert(OK());
02811 }
02812
02814 bool
02815 PPL::operator==(const Polyhedron& x, const Polyhedron& y) {
02816
02817
02818 if (x.topology() != y.topology() || x.space_dim != y.space_dim)
02819 return false;
02820
02821 if (x.marked_empty())
02822 return y.is_empty();
02823 else if (y.marked_empty())
02824 return x.is_empty();
02825 else if (x.space_dim == 0)
02826 return true;
02827
02828 switch (x.quick_equivalence_test(y)) {
02829 case Polyhedron::TVB_TRUE:
02830 return true;
02831
02832 case Polyhedron::TVB_FALSE:
02833 return false;
02834
02835 default:
02836 if (x.is_included_in(y))
02837 if (x.marked_empty())
02838 return y.is_empty();
02839 else
02840 return y.is_included_in(x);
02841 else
02842 return false;
02843 }
02844 }
02845
02846 bool
02847 PPL::Polyhedron::contains(const Polyhedron& y) const {
02848 const Polyhedron& x = *this;
02849
02850
02851 if (x.topology() != y.topology())
02852 throw_topology_incompatible("contains(y)", "y", y);
02853
02854
02855 if (x.space_dim != y.space_dim)
02856 throw_dimension_incompatible("contains(y)", "y", y);
02857
02858 if (y.marked_empty())
02859 return true;
02860 else if (x.marked_empty())
02861 return y.is_empty();
02862 else if (y.space_dim == 0)
02863 return true;
02864 else if (x.quick_equivalence_test(y) == Polyhedron::TVB_TRUE)
02865 return true;
02866 else
02867 return y.is_included_in(x);
02868 }
02869
02870 bool
02871 PPL::Polyhedron::is_disjoint_from(const Polyhedron& y) const {
02872 Polyhedron z = *this;
02873 z.intersection_assign_and_minimize(y);
02874 return z.is_empty();
02875 }
02876
02877 void
02878 PPL::Polyhedron::ascii_dump(std::ostream& s) const {
02879 s << "space_dim " << space_dim << "\n";
02880 status.ascii_dump(s);
02881 s << "\ncon_sys ("
02882 << (constraints_are_up_to_date() ? "" : "not_")
02883 << "up-to-date)"
02884 << "\n";
02885 con_sys.ascii_dump(s);
02886 s << "\ngen_sys ("
02887 << (generators_are_up_to_date() ? "" : "not_")
02888 << "up-to-date)"
02889 << "\n";
02890 gen_sys.ascii_dump(s);
02891 s << "\nsat_c\n";
02892 sat_c.ascii_dump(s);
02893 s << "\nsat_g\n";
02894 sat_g.ascii_dump(s);
02895 s << "\n";
02896 }
02897
02898 PPL_OUTPUT_DEFINITIONS(Polyhedron);
02899
02900 bool
02901 PPL::Polyhedron::ascii_load(std::istream& s) {
02902 std::string str;
02903
02904 if (!(s >> str) || str != "space_dim")
02905 return false;
02906
02907 if (!(s >> space_dim))
02908 return false;
02909
02910 if (!status.ascii_load(s))
02911 return false;
02912
02913 if (!(s >> str) || str != "con_sys")
02914 return false;
02915
02916 if (!(s >> str) || (str != "(not_up-to-date)" && str != "(up-to-date)"))
02917 return false;
02918
02919 if (!con_sys.ascii_load(s))
02920 return false;
02921
02922 if (!(s >> str) || str != "gen_sys")
02923 return false;
02924
02925 if (!(s >> str) || (str != "(not_up-to-date)" && str != "(up-to-date)"))
02926 return false;
02927
02928 if (!gen_sys.ascii_load(s))
02929 return false;
02930
02931 if (!(s >> str) || str != "sat_c")
02932 return false;
02933
02934 if (!sat_c.ascii_load(s))
02935 return false;
02936
02937 if (!(s >> str) || str != "sat_g")
02938 return false;
02939
02940 if (!sat_g.ascii_load(s))
02941 return false;
02942
02943
02944 assert(OK());
02945 return true;
02946 }
02947
02948 PPL::memory_size_type
02949 PPL::Polyhedron::external_memory_in_bytes() const {
02950 return
02951 con_sys.external_memory_in_bytes()
02952 + gen_sys.external_memory_in_bytes()
02953 + sat_c.external_memory_in_bytes()
02954 + sat_g.external_memory_in_bytes();
02955 }
02956
02958 std::ostream&
02959 PPL::IO_Operators::operator<<(std::ostream& s, const Polyhedron& ph) {
02960 if (ph.is_empty())
02961 s << "false";
02962 else
02963 s << ph.minimized_constraints();
02964 return s;
02965 }