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