00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifndef PPL_Polyhedron_templates_hh
00024 #define PPL_Polyhedron_templates_hh 1
00025
00026 #include "Interval.defs.hh"
00027 #include "Generator.defs.hh"
00028 #include "LP_Problem.defs.hh"
00029 #include <algorithm>
00030 #include <deque>
00031
00032 namespace Parma_Polyhedra_Library {
00033
00034 template <typename Box>
00035 Polyhedron::Polyhedron(Topology topol, const Box& box)
00036 : con_sys(topol),
00037 gen_sys(topol),
00038 sat_c(),
00039 sat_g() {
00040
00041 space_dim = box.space_dimension();
00042
00043
00044 if (box.is_empty()) {
00045 set_empty();
00046 return;
00047 }
00048
00049
00050 if (space_dim == 0) {
00051 set_zero_dim_univ();
00052 return;
00053 }
00054
00055
00056
00057
00058 con_sys.insert(Variable(space_dim - 1) >= 0);
00059
00060 for (dimension_type k = space_dim; k-- > 0; ) {
00061
00062 bool l_closed = false;
00063 Coefficient l_n, l_d;
00064 bool l_bounded = box.get_lower_bound(k, l_closed, l_n, l_d);
00065 if (l_bounded && topol == NECESSARILY_CLOSED && !l_closed)
00066 throw_invalid_argument("C_Polyhedron(const Box& box):",
00067 " box has an open lower bound");
00068
00069 bool u_closed = false;
00070 Coefficient u_n, u_d;
00071 bool u_bounded = box.get_upper_bound(k, u_closed, u_n, u_d);
00072 if (u_bounded && topol == NECESSARILY_CLOSED && !u_closed)
00073 throw_invalid_argument("C_Polyhedron(const Box& box):",
00074 " box has an open upper bound");
00075
00076
00077 if (l_bounded && u_bounded
00078 && l_closed && u_closed
00079 && l_n == u_n && l_d == u_d) {
00080
00081 con_sys.insert(l_d * Variable(k) == l_n);
00082 }
00083 else {
00084
00085 if (l_bounded) {
00086 if (l_closed)
00087
00088 con_sys.insert(l_d * Variable(k) >= l_n);
00089 else
00090
00091 con_sys.insert(l_d * Variable(k) > l_n);
00092 }
00093
00094 if (u_bounded) {
00095 if (u_closed)
00096
00097 con_sys.insert(u_d * Variable(k) <= u_n);
00098 else
00099
00100 con_sys.insert(u_d * Variable(k) < u_n);
00101 }
00102 }
00103 }
00104
00105
00106 con_sys.add_low_level_constraints();
00107
00108 dimension_type n_rows = con_sys.num_rows() - 1;
00109 con_sys[0].swap(con_sys[n_rows]);
00110 con_sys.set_sorted(false);
00111
00112 con_sys.set_index_first_pending_row(n_rows);
00113 con_sys.erase_to_end(n_rows);
00114
00115
00116 set_constraints_up_to_date();
00117 assert(OK());
00118 }
00119
00120 template <typename Box>
00121 void
00122 Polyhedron::shrink_bounding_box(Box& box, Complexity_Class complexity) const {
00123 bool reduce_complexity = (complexity != ANY_COMPLEXITY);
00124 if (!reduce_complexity
00125 || (!has_something_pending() && constraints_are_minimized())) {
00126
00127
00128 if (is_universe())
00129 return;
00130 }
00131 if (reduce_complexity) {
00132 if (marked_empty()
00133 || (generators_are_up_to_date() && gen_sys.num_rows() == 0)) {
00134 box.set_empty();
00135 return;
00136 }
00137 else if (constraints_are_up_to_date()) {
00138
00139 for (Constraint_System::const_iterator i = con_sys.begin(),
00140 cs_end = con_sys.end(); i != cs_end; ++i)
00141 if (i->is_inconsistent()) {
00142 box.set_empty();
00143 return;
00144 }
00145
00146
00147 if (complexity == SIMPLEX_COMPLEXITY
00148
00149 && is_necessarily_closed()) {
00150 LP_Problem lp(con_sys);
00151 if (!lp.is_satisfiable()) {
00152 box.set_empty();
00153 return;
00154 }
00155 }
00156 }
00157 }
00158 else
00159
00160
00161 if (is_empty()) {
00162 box.set_empty();
00163 return;
00164 }
00165
00166 if (space_dim == 0)
00167 return;
00168
00169
00170
00171
00172 std::vector<LBoundary>
00173 lower_bound(space_dim,
00174 LBoundary(ERational(PLUS_INFINITY), LBoundary::OPEN));
00175
00176 std::vector<UBoundary>
00177 upper_bound(space_dim,
00178 UBoundary(ERational(MINUS_INFINITY), UBoundary::OPEN));
00179
00180 if (!reduce_complexity && has_something_pending())
00181 process_pending();
00182
00183
00184
00185
00186 if (reduce_complexity &&
00187 (!generators_are_up_to_date() || has_pending_constraints())) {
00188
00189 assert(constraints_are_up_to_date());
00190
00191
00192
00193
00194 Constraint_System cs(con_sys);
00195 if (cs.num_pending_rows() > 0)
00196 cs.unset_pending_rows();
00197 if (has_pending_constraints() || !constraints_are_minimized())
00198 cs.simplify();
00199
00200 const Constraint_System::const_iterator cs_begin = cs.begin();
00201 const Constraint_System::const_iterator cs_end = cs.end();
00202
00203 for (Constraint_System::const_iterator i = cs_begin; i != cs_end; ++i) {
00204 dimension_type varid = space_dim;
00205 const Constraint& c = *i;
00206
00207 if (c.is_inconsistent()) {
00208 box.set_empty();
00209 return;
00210 }
00211 for (dimension_type j = space_dim; j-- > 0; ) {
00212
00213
00214 if (c.coefficient(Variable(j)) != 0)
00215 if (varid != space_dim) {
00216 varid = space_dim;
00217 break;
00218 }
00219 else
00220 varid = j;
00221 }
00222 if (varid != space_dim) {
00223 Coefficient_traits::const_reference d = c.coefficient(Variable(varid));
00224 Coefficient_traits::const_reference n = c.inhomogeneous_term();
00225
00226
00227
00228
00229
00230 mpq_class q;
00231 assign_r(q.get_num(), n, ROUND_NOT_NEEDED);
00232 assign_r(q.get_den(), d, ROUND_NOT_NEEDED);
00233 q.canonicalize();
00234
00235 q = -q;
00236 const ERational r(q, ROUND_NOT_NEEDED);
00237 const Constraint::Type c_type = c.type();
00238 switch (c_type) {
00239 case Constraint::EQUALITY:
00240 lower_bound[varid] = LBoundary(r, LBoundary::CLOSED);
00241 upper_bound[varid] = UBoundary(r, UBoundary::CLOSED);
00242 break;
00243 case Constraint::NONSTRICT_INEQUALITY:
00244 case Constraint::STRICT_INEQUALITY:
00245 if (d > 0)
00246
00247
00248 lower_bound[varid]
00249 = LBoundary(r, (c_type == Constraint::NONSTRICT_INEQUALITY
00250 ? LBoundary::CLOSED
00251 : LBoundary::OPEN));
00252 else {
00253
00254
00255
00256 assert(d < 0);
00257 upper_bound[varid]
00258 = UBoundary(r, (c_type == Constraint::NONSTRICT_INEQUALITY
00259 ? UBoundary::CLOSED
00260 : UBoundary::OPEN));
00261 }
00262 break;
00263 }
00264 }
00265 }
00266 }
00267 else {
00268
00269
00270
00271
00272
00273
00274 const Generator_System& gs = gen_sys;
00275
00276
00277
00278 for (Generator_System::const_iterator i = gs.begin(),
00279 gs_end = gs.end(); i != gs_end; ++i) {
00280
00281 const Generator& g = *i;
00282 Generator::Type g_type = g.type();
00283 switch (g_type) {
00284 case Generator::LINE:
00285
00286
00287 for (dimension_type j = space_dim; j-- > 0; )
00288 if (g.coefficient(Variable(j)) != 0) {
00289 lower_bound[j] = LBoundary(ERational(MINUS_INFINITY),
00290 LBoundary::OPEN);
00291 upper_bound[j] = UBoundary(ERational(PLUS_INFINITY),
00292 UBoundary::OPEN);
00293 }
00294 break;
00295 case Generator::RAY:
00296
00297
00298 for (dimension_type j = space_dim; j-- > 0; ) {
00299 int sign = sgn(g.coefficient(Variable(j)));
00300 if (sign < 0)
00301 lower_bound[j] = LBoundary(ERational(MINUS_INFINITY),
00302 LBoundary::OPEN);
00303 else if (sign > 0)
00304 upper_bound[j] = UBoundary(ERational(PLUS_INFINITY),
00305 UBoundary::OPEN);
00306 }
00307 break;
00308 case Generator::POINT:
00309 case Generator::CLOSURE_POINT:
00310 {
00311 Coefficient_traits::const_reference d = g.divisor();
00312 for (dimension_type j = space_dim; j-- > 0; ) {
00313 Coefficient_traits::const_reference n = g.coefficient(Variable(j));
00314 mpq_class q;
00315 assign_r(q.get_num(), n, ROUND_NOT_NEEDED);
00316 assign_r(q.get_den(), d, ROUND_NOT_NEEDED);
00317 q.canonicalize();
00318 const ERational r(q, ROUND_NOT_NEEDED);
00319 LBoundary lb(r,(g_type == Generator::CLOSURE_POINT
00320 ? LBoundary::OPEN
00321 : LBoundary::CLOSED));
00322 if (lb < lower_bound[j])
00323 lower_bound[j] = lb;
00324 UBoundary ub(r, (g_type == Generator::CLOSURE_POINT
00325 ? UBoundary::OPEN
00326 : UBoundary::CLOSED));
00327 if (ub > upper_bound[j])
00328 upper_bound[j] = ub;
00329 }
00330 }
00331 break;
00332 }
00333 }
00334 }
00335
00336 TEMP_INTEGER(n);
00337 TEMP_INTEGER(d);
00338
00339
00340 for (dimension_type j = space_dim; j-- > 0; ) {
00341
00342 const LBoundary& lb = lower_bound[j];
00343 const ERational& lr = lb.bound();
00344 if (!is_plus_infinity(lr) && !is_minus_infinity(lr)) {
00345 n = raw_value(lr).get_num();
00346 d = raw_value(lr).get_den();
00347 box.raise_lower_bound(j, lb.is_closed(), n, d);
00348 }
00349
00350
00351 const UBoundary& ub = upper_bound[j];
00352 const ERational& ur = ub.bound();
00353 if (!is_plus_infinity(ur) && !is_minus_infinity(ur)) {
00354 n = raw_value(ur).get_num();
00355 d = raw_value(ur).get_den();
00356 box.lower_upper_bound(j, ub.is_closed(), n, d);
00357 }
00358 }
00359 }
00360
00361 template <typename Partial_Function>
00362 void
00363 Polyhedron::map_space_dimensions(const Partial_Function& pfunc) {
00364 if (space_dim == 0)
00365 return;
00366
00367 if (pfunc.has_empty_codomain()) {
00368
00369 if (marked_empty()
00370 || (has_pending_constraints()
00371 && !remove_pending_to_obtain_generators())
00372 || (!generators_are_up_to_date() && !update_generators())) {
00373
00374 space_dim = 0;
00375 con_sys.clear();
00376 }
00377 else
00378
00379 set_zero_dim_univ();
00380
00381 assert(OK());
00382 return;
00383 }
00384
00385 const dimension_type new_space_dimension = pfunc.max_in_codomain() + 1;
00386
00387 if (new_space_dimension == space_dim) {
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 std::vector<dimension_type> cycles;
00404 cycles.reserve(space_dim + space_dim/2);
00405
00406
00407 std::deque<bool> visited(space_dim);
00408
00409 for (dimension_type i = space_dim; i-- > 0; ) {
00410 if (!visited[i]) {
00411 dimension_type j = i;
00412 do {
00413 visited[j] = true;
00414
00415 dimension_type k = 0;
00416 if (!pfunc.maps(j, k))
00417 throw_invalid_argument("map_space_dimensions(pfunc)",
00418 " pfunc is inconsistent");
00419 if (k == j)
00420
00421 goto skip;
00422
00423 cycles.push_back(j+1);
00424
00425 j = k;
00426 } while (!visited[j]);
00427
00428 cycles.push_back(0);
00429 skip:
00430 ;
00431 }
00432 }
00433
00434
00435 if (cycles.empty())
00436 return;
00437
00438
00439
00440
00441 if (constraints_are_up_to_date())
00442 con_sys.permute_columns(cycles);
00443
00444 if (generators_are_up_to_date())
00445 gen_sys.permute_columns(cycles);
00446
00447 assert(OK());
00448 return;
00449 }
00450
00451
00452
00453
00454
00455 const Generator_System& old_gensys = generators();
00456
00457 if (old_gensys.num_rows() == 0) {
00458
00459 Polyhedron new_polyhedron(topology(), new_space_dimension, EMPTY);
00460 std::swap(*this, new_polyhedron);
00461 assert(OK());
00462 return;
00463 }
00464
00465
00466 std::vector<dimension_type> pfunc_maps(space_dim, not_a_dimension());
00467 for (dimension_type j = space_dim; j-- > 0; ) {
00468 dimension_type pfunc_j;
00469 if (pfunc.maps(j, pfunc_j))
00470 pfunc_maps[j] = pfunc_j;
00471 }
00472
00473 Generator_System new_gensys;
00474 for (Generator_System::const_iterator i = old_gensys.begin(),
00475 old_gensys_end = old_gensys.end(); i != old_gensys_end; ++i) {
00476 const Generator& old_g = *i;
00477 Linear_Expression e(0 * Variable(new_space_dimension-1));
00478 bool all_zeroes = true;
00479 for (dimension_type j = space_dim; j-- > 0; ) {
00480 if (old_g.coefficient(Variable(j)) != 0
00481 && pfunc_maps[j] != not_a_dimension()) {
00482 e += Variable(pfunc_maps[j]) * old_g.coefficient(Variable(j));
00483 all_zeroes = false;
00484 }
00485 }
00486 switch (old_g.type()) {
00487 case Generator::LINE:
00488 if (!all_zeroes)
00489 new_gensys.insert(line(e));
00490 break;
00491 case Generator::RAY:
00492 if (!all_zeroes)
00493 new_gensys.insert(ray(e));
00494 break;
00495 case Generator::POINT:
00496
00497 new_gensys.insert(point(e, old_g.divisor()));
00498 break;
00499 case Generator::CLOSURE_POINT:
00500
00501 new_gensys.insert(closure_point(e, old_g.divisor()));
00502 break;
00503 }
00504 }
00505 Polyhedron new_polyhedron(topology(), new_gensys);
00506 std::swap(*this, new_polyhedron);
00507 assert(OK(true));
00508 }
00509
00510 }
00511
00512 #endif // !defined(PPL_Polyhedron_templates_hh)