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