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 "Generator_System.defs.hh"
00027 #include "Generator_System.inlines.hh"
00028 #include "Congruence_System.inlines.hh"
00029 #include "Congruence_System.defs.hh"
00030 #include "Poly_Con_Relation.defs.hh"
00031 #include "Poly_Gen_Relation.defs.hh"
00032 #include "MIP_Problem.defs.hh"
00033 #include "Variables_Set.defs.hh"
00034 #include "Bit_Row.defs.hh"
00035 #include "Temp.defs.hh"
00036 #include <cassert>
00037 #include <vector>
00038 #include <deque>
00039 #include <iostream>
00040 #include <sstream>
00041 #include <stdexcept>
00042 #include <algorithm>
00043
00044 namespace Parma_Polyhedra_Library {
00045
00046 template <typename T>
00047 BD_Shape<T>::BD_Shape(const Congruence_System& cgs)
00048 : dbm(cgs.space_dimension() + 1),
00049 status(),
00050 redundancy_dbm() {
00051 add_congruences(cgs);
00052 }
00053
00054 template <typename T>
00055 BD_Shape<T>::BD_Shape(const Generator_System& gs)
00056 : dbm(gs.space_dimension() + 1), status(), redundancy_dbm() {
00057 const Generator_System::const_iterator gs_begin = gs.begin();
00058 const Generator_System::const_iterator gs_end = gs.end();
00059 if (gs_begin == gs_end) {
00060
00061 set_empty();
00062 return;
00063 }
00064
00065 const dimension_type space_dim = space_dimension();
00066 DB_Row<N>& dbm_0 = dbm[0];
00067 DIRTY_TEMP(N, tmp);
00068
00069 bool dbm_initialized = false;
00070 bool point_seen = false;
00071
00072 for (Generator_System::const_iterator gs_i = gs_begin;
00073 gs_i != gs_end; ++gs_i) {
00074 const Generator& g = *gs_i;
00075 switch (g.type()) {
00076 case Generator::POINT:
00077 point_seen = true;
00078
00079 case Generator::CLOSURE_POINT:
00080 if (!dbm_initialized) {
00081
00082 dbm_initialized = true;
00083 const Coefficient& d = g.divisor();
00084 for (dimension_type i = space_dim; i > 0; --i) {
00085 const Coefficient& g_i = g.coefficient(Variable(i-1));
00086 DB_Row<N>& dbm_i = dbm[i];
00087 for (dimension_type j = space_dim; j > 0; --j)
00088 if (i != j)
00089 div_round_up(dbm_i[j], g.coefficient(Variable(j-1)) - g_i, d);
00090 div_round_up(dbm_i[0], -g_i, d);
00091 }
00092 for (dimension_type j = space_dim; j > 0; --j)
00093 div_round_up(dbm_0[j], g.coefficient(Variable(j-1)), d);
00094
00095 }
00096 else {
00097
00098
00099 const Coefficient& d = g.divisor();
00100 for (dimension_type i = space_dim; i > 0; --i) {
00101 const Coefficient& g_i = g.coefficient(Variable(i-1));
00102 DB_Row<N>& dbm_i = dbm[i];
00103
00104 for (dimension_type j = space_dim; j > 0; --j) {
00105 div_round_up(tmp, g.coefficient(Variable(j-1)) - g_i, d);
00106 max_assign(dbm_i[j], tmp);
00107 }
00108 div_round_up(tmp, -g_i, d);
00109 max_assign(dbm_i[0], tmp);
00110 }
00111 for (dimension_type j = space_dim; j > 0; --j) {
00112 div_round_up(tmp, g.coefficient(Variable(j-1)), d);
00113 max_assign(dbm_0[j], tmp);
00114 }
00115 }
00116 break;
00117 default:
00118
00119 break;
00120 }
00121 }
00122
00123 if (!point_seen)
00124
00125 throw_generic("BD_Shape(gs)",
00126 "the non-empty generator system gs contains no points.");
00127
00128
00129 for (Generator_System::const_iterator gs_i = gs_begin;
00130 gs_i != gs_end; ++gs_i) {
00131 const Generator& g = *gs_i;
00132 switch (g.type()) {
00133 case Generator::LINE:
00134 for (dimension_type i = space_dim; i > 0; --i) {
00135 const Coefficient& g_i = g.coefficient(Variable(i-1));
00136 DB_Row<N>& dbm_i = dbm[i];
00137
00138 for (dimension_type j = space_dim; j > 0; --j)
00139 if (g_i != g.coefficient(Variable(j-1)))
00140 assign_r(dbm_i[j], PLUS_INFINITY, ROUND_NOT_NEEDED);
00141 if (g_i != 0)
00142 assign_r(dbm_i[0], PLUS_INFINITY, ROUND_NOT_NEEDED);
00143 }
00144 for (dimension_type j = space_dim; j > 0; --j)
00145 if (g.coefficient(Variable(j-1)) != 0)
00146 assign_r(dbm_0[j], PLUS_INFINITY, ROUND_NOT_NEEDED);
00147 break;
00148 case Generator::RAY:
00149 for (dimension_type i = space_dim; i > 0; --i) {
00150 const Coefficient& g_i = g.coefficient(Variable(i-1));
00151 DB_Row<N>& dbm_i = dbm[i];
00152
00153 for (dimension_type j = space_dim; j > 0; --j)
00154 if (g_i < g.coefficient(Variable(j-1)))
00155 assign_r(dbm_i[j], PLUS_INFINITY, ROUND_NOT_NEEDED);
00156 if (g_i < 0)
00157 assign_r(dbm_i[0], PLUS_INFINITY, ROUND_NOT_NEEDED);
00158 }
00159 for (dimension_type j = space_dim; j > 0; --j)
00160 if (g.coefficient(Variable(j-1)) > 0)
00161 assign_r(dbm_0[j], PLUS_INFINITY, ROUND_NOT_NEEDED);
00162 break;
00163 default:
00164
00165 break;
00166 }
00167 }
00168 set_shortest_path_closed();
00169 assert(OK());
00170 }
00171
00172 template <typename T>
00173 BD_Shape<T>::BD_Shape(const Polyhedron& ph, const Complexity_Class complexity)
00174 : dbm(), status(), redundancy_dbm() {
00175 const dimension_type num_dimensions = ph.space_dimension();
00176
00177 if (ph.marked_empty()) {
00178 *this = BD_Shape<T>(num_dimensions, EMPTY);
00179 return;
00180 }
00181
00182 if (num_dimensions == 0) {
00183 *this = BD_Shape<T>(num_dimensions, UNIVERSE);
00184 return;
00185 }
00186
00187
00188
00189 if (complexity == ANY_COMPLEXITY
00190 || (!ph.has_pending_constraints() && ph.generators_are_up_to_date())) {
00191 *this = BD_Shape<T>(ph.generators());
00192 return;
00193 }
00194
00195
00196
00197
00198 assert(ph.constraints_are_up_to_date());
00199
00200 if (!ph.has_something_pending() && ph.constraints_are_minimized()) {
00201
00202
00203 if (ph.is_universe()) {
00204 *this = BD_Shape<T>(num_dimensions, UNIVERSE);
00205 return;
00206 }
00207 }
00208
00209
00210 for (Constraint_System::const_iterator i = ph.con_sys.begin(),
00211 cs_end = ph.con_sys.end(); i != cs_end; ++i)
00212 if (i->is_inconsistent()) {
00213 *this = BD_Shape<T>(num_dimensions, EMPTY);
00214 return;
00215 }
00216
00217
00218
00219 if (complexity == SIMPLEX_COMPLEXITY) {
00220 MIP_Problem lp(num_dimensions);
00221 lp.set_optimization_mode(MAXIMIZATION);
00222
00223 const Constraint_System& ph_cs = ph.constraints();
00224 if (!ph_cs.has_strict_inequalities())
00225 lp.add_constraints(ph_cs);
00226 else
00227
00228 for (Constraint_System::const_iterator i = ph_cs.begin(),
00229 ph_cs_end = ph_cs.end(); i != ph_cs_end; ++i) {
00230 const Constraint& c = *i;
00231 if (c.is_strict_inequality())
00232 lp.add_constraint(Linear_Expression(c) >= 0);
00233 else
00234 lp.add_constraint(c);
00235 }
00236
00237
00238 if (!lp.is_satisfiable()) {
00239 *this = BD_Shape<T>(num_dimensions, EMPTY);
00240 return;
00241 }
00242
00243
00244 *this = BD_Shape<T>(num_dimensions, UNIVERSE);
00245
00246 Generator g(point());
00247 TEMP_INTEGER(num);
00248 TEMP_INTEGER(den);
00249 for (dimension_type i = 1; i <= num_dimensions; ++i) {
00250 Variable x(i-1);
00251
00252 lp.set_objective_function(x);
00253 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00254 g = lp.optimizing_point();
00255 lp.evaluate_objective_function(g, num, den);
00256 div_round_up(dbm[0][i], num, den);
00257 }
00258
00259 for (dimension_type j = 1; j <= num_dimensions; ++j) {
00260 if (i == j)
00261 continue;
00262 Variable y(j-1);
00263 lp.set_objective_function(x - y);
00264 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00265 g = lp.optimizing_point();
00266 lp.evaluate_objective_function(g, num, den);
00267 div_round_up(dbm[j][i], num, den);
00268 }
00269 }
00270
00271 lp.set_objective_function(-x);
00272 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00273 g = lp.optimizing_point();
00274 lp.evaluate_objective_function(g, num, den);
00275 div_round_up(dbm[i][0], num, den);
00276 }
00277 }
00278 set_shortest_path_closed();
00279 assert(OK());
00280 return;
00281 }
00282
00283
00284 assert(complexity == POLYNOMIAL_COMPLEXITY);
00285 *this = BD_Shape<T>(num_dimensions, UNIVERSE);
00286 refine_with_constraints(ph.constraints());
00287 }
00288
00289 template <typename T>
00290 dimension_type
00291 BD_Shape<T>::affine_dimension() const {
00292 const dimension_type space_dim = space_dimension();
00293
00294 if (space_dim == 0)
00295 return 0;
00296
00297
00298
00299 shortest_path_closure_assign();
00300 if (marked_empty())
00301 return 0;
00302
00303
00304
00305
00306 std::vector<dimension_type> predecessor;
00307 compute_predecessors(predecessor);
00308
00309
00310
00311 dimension_type affine_dim = 0;
00312
00313 for (dimension_type i = 1; i <= space_dim; ++i)
00314 if (predecessor[i] == i)
00315 ++affine_dim;
00316
00317 return affine_dim;
00318 }
00319
00320 template <typename T>
00321 Congruence_System
00322 BD_Shape<T>::minimized_congruences() const {
00323
00324
00325 shortest_path_closure_assign();
00326
00327 const dimension_type space_dim = space_dimension();
00328 Congruence_System cgs;
00329 if (space_dim == 0) {
00330 if (marked_empty())
00331 cgs = Congruence_System::zero_dim_empty();
00332 }
00333 else if (marked_empty())
00334 cgs.insert((0*Variable(space_dim-1) %= 1) / 0);
00335 else {
00336
00337
00338 cgs.insert(0*Variable(space_dim-1) == 0);
00339
00340 TEMP_INTEGER(num);
00341 TEMP_INTEGER(den);
00342
00343
00344 std::vector<dimension_type> leaders;
00345 compute_leaders(leaders);
00346
00347
00348 const DB_Row<N>& dbm_0 = dbm[0];
00349 for (dimension_type i = 1; i <= space_dim; ++i) {
00350 const dimension_type leader = leaders[i];
00351 if (i != leader) {
00352
00353 if (leader == 0) {
00354
00355 assert(!is_plus_infinity(dbm_0[i]));
00356 numer_denom(dbm_0[i], num, den);
00357 cgs.insert(den*Variable(i-1) == num);
00358 }
00359 else {
00360
00361 assert(!is_plus_infinity(dbm[i][leader]));
00362 numer_denom(dbm[i][leader], num, den);
00363 cgs.insert(den*Variable(leader-1) - den*Variable(i-1) == num);
00364 }
00365 }
00366 }
00367 }
00368 return cgs;
00369 }
00370
00371 template <typename T>
00372 void
00373 BD_Shape<T>::add_constraint(const Constraint& c) {
00374 const dimension_type c_space_dim = c.space_dimension();
00375
00376 if (c_space_dim > space_dimension())
00377 throw_dimension_incompatible("add_constraint(c)", c);
00378
00379
00380 if (c.is_strict_inequality()) {
00381 if (c.is_inconsistent()) {
00382 set_empty();
00383 return;
00384 }
00385 if (c.is_tautological())
00386 return;
00387
00388 throw_generic("add_constraint(c)", "strict inequalities are not allowed");
00389 }
00390
00391 dimension_type num_vars = 0;
00392 dimension_type i = 0;
00393 dimension_type j = 0;
00394 TEMP_INTEGER(coeff);
00395
00396 if (!extract_bounded_difference(c, c_space_dim, num_vars, i, j, coeff))
00397 throw_generic("add_constraint(c)",
00398 "c is not a bounded difference constraint");
00399
00400 const Coefficient& inhomo = c.inhomogeneous_term();
00401 if (num_vars == 0) {
00402
00403 if (inhomo < 0
00404 || (inhomo != 0 && c.is_equality()))
00405 set_empty();
00406 return;
00407 }
00408
00409
00410
00411 const bool negative = (coeff < 0);
00412 N& x = negative ? dbm[i][j] : dbm[j][i];
00413 N& y = negative ? dbm[j][i] : dbm[i][j];
00414 if (negative)
00415 neg_assign(coeff);
00416
00417 bool changed = false;
00418
00419 DIRTY_TEMP(N, d);
00420 div_round_up(d, inhomo, coeff);
00421 if (x > d) {
00422 x = d;
00423 changed = true;
00424 }
00425
00426 if (c.is_equality()) {
00427
00428 TEMP_INTEGER(minus_c_term);
00429 neg_assign(minus_c_term, inhomo);
00430 div_round_up(d, minus_c_term, coeff);
00431 if (y > d) {
00432 y = d;
00433 changed = true;
00434 }
00435 }
00436
00437
00438
00439 if (changed && marked_shortest_path_closed())
00440 reset_shortest_path_closed();
00441 assert(OK());
00442 }
00443
00444 template <typename T>
00445 void
00446 BD_Shape<T>::add_congruence(const Congruence& cg) {
00447 const dimension_type cg_space_dim = cg.space_dimension();
00448
00449
00450 if (space_dimension() < cg_space_dim)
00451 throw_dimension_incompatible("add_congruence(cg)", cg);
00452
00453
00454 if (cg.is_proper_congruence()) {
00455 if (cg.is_tautological())
00456 return;
00457 if (cg.is_inconsistent()) {
00458 set_empty();
00459 return;
00460 }
00461
00462 throw_generic("add_congruence(cg)",
00463 "cg is a non-trivial, proper congruence");
00464 }
00465
00466 assert(cg.is_equality());
00467 Constraint c(cg);
00468 add_constraint(c);
00469 }
00470
00471 template <typename T>
00472 void
00473 BD_Shape<T>::refine_no_check(const Constraint& c) {
00474 assert(!marked_empty());
00475 const dimension_type c_space_dim = c.space_dimension();
00476 assert(c_space_dim <= space_dimension());
00477
00478 dimension_type num_vars = 0;
00479 dimension_type i = 0;
00480 dimension_type j = 0;
00481 TEMP_INTEGER(coeff);
00482
00483 if (!extract_bounded_difference(c, c_space_dim, num_vars, i, j, coeff))
00484 return;
00485
00486 const Coefficient& inhomo = c.inhomogeneous_term();
00487 if (num_vars == 0) {
00488
00489 if (inhomo < 0
00490 || (c.is_equality() && inhomo != 0)
00491 || (c.is_strict_inequality() && inhomo == 0))
00492 set_empty();
00493 return;
00494 }
00495
00496
00497
00498 const bool negative = (coeff < 0);
00499 N& x = negative ? dbm[i][j] : dbm[j][i];
00500 N& y = negative ? dbm[j][i] : dbm[i][j];
00501 if (negative)
00502 neg_assign(coeff);
00503
00504 bool changed = false;
00505
00506 DIRTY_TEMP(N, d);
00507 div_round_up(d, inhomo, coeff);
00508 if (x > d) {
00509 x = d;
00510 changed = true;
00511 }
00512
00513 if (c.is_equality()) {
00514
00515 TEMP_INTEGER(minus_c_term);
00516 neg_assign(minus_c_term, inhomo);
00517 div_round_up(d, minus_c_term, coeff);
00518 if (y > d) {
00519 y = d;
00520 changed = true;
00521 }
00522 }
00523
00524
00525
00526 if (changed && marked_shortest_path_closed())
00527 reset_shortest_path_closed();
00528 assert(OK());
00529 }
00530
00531 template <typename T>
00532 void
00533 BD_Shape<T>::concatenate_assign(const BD_Shape& y) {
00534 BD_Shape& x = *this;
00535
00536 const dimension_type x_space_dim = x.space_dimension();
00537 const dimension_type y_space_dim = y.space_dimension();
00538
00539
00540
00541 if (y_space_dim == 0 && y.marked_empty()) {
00542 set_empty();
00543 return;
00544 }
00545
00546
00547
00548 if (x_space_dim == 0 && marked_empty()) {
00549 dbm.grow(y_space_dim + 1);
00550 assert(OK());
00551 return;
00552 }
00553
00554
00555
00556
00557
00558
00559
00560 add_space_dimensions_and_embed(y_space_dim);
00561 const dimension_type new_space_dim = x_space_dim + y_space_dim;
00562 for (dimension_type i = x_space_dim + 1; i <= new_space_dim; ++i) {
00563 DB_Row<N>& dbm_i = dbm[i];
00564 dbm_i[0] = y.dbm[i - x_space_dim][0];
00565 dbm[0][i] = y.dbm[0][i - x_space_dim];
00566 for (dimension_type j = x_space_dim + 1; j <= new_space_dim; ++j)
00567 dbm_i[j] = y.dbm[i - x_space_dim][j - x_space_dim];
00568 }
00569
00570 if (marked_shortest_path_closed())
00571 reset_shortest_path_closed();
00572 assert(OK());
00573 }
00574
00575 template <typename T>
00576 bool
00577 BD_Shape<T>::contains(const BD_Shape& y) const {
00578 const BD_Shape<T>& x = *this;
00579 const dimension_type x_space_dim = x.space_dimension();
00580
00581
00582 if (x_space_dim != y.space_dimension())
00583 throw_dimension_incompatible("contains(y)", y);
00584
00585
00586
00587
00588
00589 if (x_space_dim == 0) {
00590 if (!marked_empty())
00591 return true;
00592 else
00593 return y.marked_empty();
00594 }
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 y.shortest_path_closure_assign();
00617
00618
00619 if (y.marked_empty())
00620 return true;
00621
00622
00623
00624 for (dimension_type i = x_space_dim + 1; i-- > 0; ) {
00625 const DB_Row<N>& x_dbm_i = x.dbm[i];
00626 const DB_Row<N>& y_dbm_i = y.dbm[i];
00627 for (dimension_type j = x_space_dim + 1; j-- > 0; )
00628 if (x_dbm_i[j] < y_dbm_i[j])
00629 return false;
00630 }
00631 return true;
00632 }
00633
00634 template <typename T>
00635 bool
00636 BD_Shape<T>::is_disjoint_from(const BD_Shape& y) const {
00637 const dimension_type space_dim = space_dimension();
00638
00639 if (space_dim != y.space_dimension())
00640 throw_dimension_incompatible("is_disjoint_from(y)", y);
00641
00642
00643
00644 shortest_path_closure_assign();
00645 if (marked_empty())
00646 return true;
00647 y.shortest_path_closure_assign();
00648 if (y.marked_empty())
00649 return true;
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663 DIRTY_TEMP(N, tmp);
00664 for (dimension_type i = space_dim+1; i-- > 0; ) {
00665 const DB_Row<N>& x_i = dbm[i];
00666 for (dimension_type j = space_dim+1; j-- > 0; ) {
00667 neg_assign_r(tmp, y.dbm[j][i], ROUND_UP);
00668 if (x_i[j] < tmp)
00669 return true;
00670 }
00671 }
00672
00673 return false;
00674 }
00675
00676 template <typename T>
00677 bool
00678 BD_Shape<T>::is_universe() const {
00679 if (marked_empty())
00680 return false;
00681
00682 const dimension_type space_dim = space_dimension();
00683
00684
00685 if (space_dim == 0)
00686 return true;
00687
00688
00689
00690 for (dimension_type i = space_dim + 1; i-- > 0; ) {
00691 const DB_Row<N>& dbm_i = dbm[i];
00692 for (dimension_type j = space_dim + 1; j-- > 0; )
00693 if (!is_plus_infinity(dbm_i[j]))
00694 return false;
00695 }
00696 return true;
00697 }
00698
00699 template <typename T>
00700 bool
00701 BD_Shape<T>::is_bounded() const {
00702 shortest_path_closure_assign();
00703 const dimension_type space_dim = space_dimension();
00704
00705 if (marked_empty() || space_dim == 0)
00706 return true;
00707
00708
00709
00710 for (dimension_type i = space_dim + 1; i-- > 0; ) {
00711 const DB_Row<N>& dbm_i = dbm[i];
00712 for (dimension_type j = space_dim + 1; j-- > 0; )
00713 if (i != j)
00714 if (is_plus_infinity(dbm_i[j]))
00715 return false;
00716 }
00717
00718 return true;
00719 }
00720
00721 template <typename T>
00722 bool
00723 BD_Shape<T>::contains_integer_point() const {
00724
00725 if (is_empty())
00726 return false;
00727
00728 const dimension_type space_dim = space_dimension();
00729 if (space_dim == 0)
00730 return true;
00731
00732
00733
00734 if (std::numeric_limits<T>::is_integer)
00735 return true;
00736
00737
00738
00739 BD_Shape<mpz_class> bds_z(space_dim);
00740 typedef BD_Shape<mpz_class>::N Z;
00741 bds_z.reset_shortest_path_closed();
00742 DIRTY_TEMP(N, tmp);
00743 bool all_integers = true;
00744 for (dimension_type i = space_dim + 1; i-- > 0; ) {
00745 DB_Row<Z>& z_i = bds_z.dbm[i];
00746 const DB_Row<N>& dbm_i = dbm[i];
00747 for (dimension_type j = space_dim + 1; j-- > 0; ) {
00748 const N& dbm_i_j = dbm_i[j];
00749 if (is_plus_infinity(dbm_i_j))
00750 continue;
00751 if (is_integer(dbm_i_j))
00752 assign_r(z_i[j], dbm_i_j, ROUND_NOT_NEEDED);
00753 else {
00754 all_integers = false;
00755 Z& z_i_j = z_i[j];
00756
00757 neg_assign_r(tmp, dbm_i_j, ROUND_NOT_NEEDED);
00758 assign_r(z_i_j, tmp, ROUND_UP);
00759 neg_assign_r(z_i_j, z_i_j, ROUND_NOT_NEEDED);
00760 }
00761 }
00762 }
00763 return all_integers || !bds_z.is_empty();
00764 }
00765
00766 template <typename T>
00767 bool
00768 BD_Shape<T>::constrains(const Variable var) const {
00769
00770 const dimension_type var_space_dim = var.space_dimension();
00771 if (space_dimension() < var_space_dim)
00772 throw_dimension_incompatible("constrains(v)", "v", var);
00773
00774
00775
00776 if (marked_empty())
00777 return true;
00778
00779
00780 const DB_Row<N>& dbm_v = dbm[var_space_dim];
00781 for (dimension_type i = dbm.num_rows(); i-- > 0; ) {
00782 if (!is_plus_infinity(dbm_v[i])
00783 || !is_plus_infinity(dbm[i][var_space_dim]))
00784 return true;
00785 }
00786
00787
00788
00789 return is_empty();
00790 }
00791
00792 template <typename T>
00793 void
00794 BD_Shape<T>
00795 ::compute_predecessors(std::vector<dimension_type>& predecessor) const {
00796 assert(!marked_empty() && marked_shortest_path_closed());
00797 assert(predecessor.size() == 0);
00798
00799
00800
00801
00802
00803 const dimension_type pred_size = dbm.num_rows();
00804
00805 predecessor.reserve(pred_size);
00806 for (dimension_type i = 0; i < pred_size; ++i)
00807 predecessor.push_back(i);
00808
00809 for (dimension_type i = pred_size; i-- > 1; )
00810 if (i == predecessor[i]) {
00811 const DB_Row<N>& dbm_i = dbm[i];
00812 for (dimension_type j = i; j-- > 0; )
00813 if (j == predecessor[j]
00814 && is_additive_inverse(dbm[j][i], dbm_i[j])) {
00815
00816 predecessor[i] = j;
00817 break;
00818 }
00819 }
00820 }
00821
00822 template <typename T>
00823 void
00824 BD_Shape<T>::compute_leaders(std::vector<dimension_type>& leaders) const {
00825 assert(!marked_empty() && marked_shortest_path_closed());
00826 assert(leaders.size() == 0);
00827
00828 compute_predecessors(leaders);
00829
00830 assert(leaders[0] == 0);
00831 for (dimension_type i = 1, l_size = leaders.size(); i != l_size; ++i) {
00832 const dimension_type l_i = leaders[i];
00833 assert(l_i <= i);
00834 if (l_i != i) {
00835 const dimension_type ll_i = leaders[l_i];
00836 assert(ll_i == leaders[ll_i]);
00837 leaders[i] = ll_i;
00838 }
00839 }
00840 }
00841
00842 template <typename T>
00843 bool
00844 BD_Shape<T>::is_shortest_path_reduced() const {
00845
00846 if (marked_empty())
00847 return true;
00848
00849 const dimension_type space_dim = space_dimension();
00850
00851 if (space_dim == 0)
00852 return true;
00853
00854
00855
00856
00857 if (!marked_shortest_path_reduced())
00858 return false;
00859
00860 const BD_Shape x_copy = *this;
00861 x_copy.shortest_path_closure_assign();
00862
00863 if (x_copy.marked_empty())
00864 return false;
00865
00866
00867 std::vector<dimension_type> leader(space_dim + 1);
00868
00869
00870 for (dimension_type i = space_dim + 1; i-- > 0; )
00871 leader[i] = i;
00872
00873
00874
00875
00876
00877 for (dimension_type i = 0; i < space_dim; ++i) {
00878 const DB_Row<N>& x_copy_dbm_i = x_copy.dbm[i];
00879 for (dimension_type j = i + 1; j <= space_dim; ++j)
00880 if (is_additive_inverse(x_copy.dbm[j][i], x_copy_dbm_i[j]))
00881
00882
00883 leader[j] = leader[i];
00884 }
00885
00886
00887
00888
00889
00890
00891 DIRTY_TEMP(N, c);
00892 for (dimension_type k = 0; k <= space_dim; ++k)
00893 if (leader[k] == k) {
00894 const DB_Row<N>& x_k = x_copy.dbm[k];
00895 for (dimension_type i = 0; i <= space_dim; ++i)
00896 if (leader[i] == i) {
00897 const DB_Row<N>& x_i = x_copy.dbm[i];
00898 const Bit_Row& redundancy_i = redundancy_dbm[i];
00899 const N& x_i_k = x_i[k];
00900 for (dimension_type j = 0; j <= space_dim; ++j)
00901 if (leader[j] == j) {
00902 const N& x_i_j = x_i[j];
00903 if (!is_plus_infinity(x_i_j)) {
00904 add_assign_r(c, x_i_k, x_k[j], ROUND_UP);
00905 if (x_i_j >= c && !redundancy_i[j])
00906 return false;
00907 }
00908 }
00909 }
00910 }
00911
00912
00913
00914
00915
00916 std::vector<dimension_type> var_conn(space_dim + 1);
00917 for (dimension_type i = space_dim + 1; i-- > 0; )
00918 var_conn[i] = space_dim + 1;
00919
00920
00921
00922
00923
00924
00925 for (dimension_type i = 0; i <= space_dim; ++i) {
00926
00927
00928 dimension_type t = 0;
00929 dimension_type ld_i = leader[i];
00930
00931 if (ld_i == i) {
00932 for (dimension_type j = 0; j <= space_dim; ++j) {
00933 dimension_type ld_j = leader[j];
00934
00935
00936 if (j != ld_j)
00937 if (!redundancy_dbm[i][j]) {
00938 if (t == 1)
00939
00940 return false;
00941 else
00942 if (ld_j != i)
00943
00944 return false;
00945 else {
00946 ++t;
00947 var_conn[i] = j;
00948 }
00949 }
00950 }
00951 }
00952
00953 else {
00954 for (dimension_type j = 0; j <= space_dim; ++j) {
00955 if (!redundancy_dbm[i][j]) {
00956 dimension_type ld_j = leader[j];
00957 if (ld_i != ld_j)
00958
00959 return false;
00960 else {
00961 if (t == 1)
00962
00963 return false;
00964 else {
00965 ++t;
00966 var_conn[i] = j;
00967 }
00968 }
00969
00970
00971 if (t == 0)
00972 return false;
00973 }
00974 }
00975 }
00976 }
00977
00978
00979
00980 std::vector<bool> just_checked(space_dim + 1);
00981 for (dimension_type i = space_dim + 1; i-- > 0; )
00982 just_checked[i] = false;
00983
00984
00985
00986 for (dimension_type i = 0; i <= space_dim; ++i) {
00987 bool jc_i = just_checked[i];
00988
00989 if (!jc_i) {
00990 dimension_type v_con = var_conn[i];
00991
00992
00993 if (v_con != space_dim + 1) {
00994
00995
00996 while (v_con != i) {
00997 just_checked[v_con] = true;
00998 v_con = var_conn[v_con];
00999
01000
01001 if (just_checked[v_con])
01002 return false;
01003 }
01004 }
01005 }
01006 just_checked[i] = true;
01007 }
01008
01009
01010 return true;
01011 }
01012
01013 template <typename T>
01014 bool
01015 BD_Shape<T>::bounds(const Linear_Expression& expr,
01016 const bool from_above) const {
01017
01018
01019 const dimension_type expr_space_dim = expr.space_dimension();
01020 const dimension_type space_dim = space_dimension();
01021 if (space_dim < expr_space_dim)
01022 throw_dimension_incompatible((from_above
01023 ? "bounds_from_above(e)"
01024 : "bounds_from_below(e)"), "e", expr);
01025
01026 shortest_path_closure_assign();
01027
01028 if (space_dim == 0 || marked_empty())
01029 return true;
01030
01031
01032
01033 const Constraint& c = from_above ? expr <= 0 : expr >= 0;
01034 const dimension_type c_space_dim = c.space_dimension();
01035 dimension_type num_vars = 0;
01036 dimension_type i = 0;
01037 dimension_type j = 0;
01038 TEMP_INTEGER(coeff);
01039
01040 if (extract_bounded_difference(c, c_space_dim, num_vars, i, j, coeff)) {
01041 if (num_vars == 0)
01042
01043 return true;
01044
01045 const N& x = (coeff < 0) ? dbm[i][j] : dbm[j][i];
01046 return !is_plus_infinity(x);
01047 }
01048 else {
01049
01050 Optimization_Mode mode_bounds
01051 = from_above ? MAXIMIZATION : MINIMIZATION;
01052 MIP_Problem mip(space_dim, constraints(), expr, mode_bounds);
01053
01054 return (mip.solve() == OPTIMIZED_MIP_PROBLEM);
01055 }
01056 }
01057
01058 template <typename T>
01059 bool
01060 BD_Shape<T>::max_min(const Linear_Expression& expr,
01061 const bool maximize,
01062 Coefficient& ext_n, Coefficient& ext_d,
01063 bool& included) const {
01064
01065
01066 const dimension_type space_dim = space_dimension();
01067 const dimension_type expr_space_dim = expr.space_dimension();
01068 if (space_dim < expr_space_dim)
01069 throw_dimension_incompatible((maximize
01070 ? "maximize(e, ...)"
01071 : "minimize(e, ...)"), "e", expr);
01072
01073 if (space_dim == 0) {
01074 if (marked_empty())
01075 return false;
01076 else {
01077 ext_n = expr.inhomogeneous_term();
01078 ext_d = 1;
01079 included = true;
01080 return true;
01081 }
01082 }
01083
01084 shortest_path_closure_assign();
01085
01086 if (marked_empty())
01087 return false;
01088
01089
01090
01091 const Constraint& c = maximize ? expr <= 0 : expr >= 0;
01092 const dimension_type c_space_dim = c.space_dimension();
01093 dimension_type num_vars = 0;
01094 dimension_type i = 0;
01095 dimension_type j = 0;
01096 TEMP_INTEGER(coeff);
01097
01098 if (!extract_bounded_difference(c, c_space_dim, num_vars, i, j, coeff)) {
01099 Optimization_Mode mode_max_min
01100 = maximize ? MAXIMIZATION : MINIMIZATION;
01101 MIP_Problem mip(space_dim, constraints(), expr, mode_max_min);
01102 if (mip.solve() == OPTIMIZED_MIP_PROBLEM) {
01103 mip.optimal_value(ext_n, ext_d);
01104 included = true;
01105 return true;
01106 }
01107 else
01108
01109 return false;
01110 }
01111 else {
01112
01113 if (num_vars == 0) {
01114
01115 ext_n = expr.inhomogeneous_term();
01116 ext_d = 1;
01117 included = true;
01118 return true;
01119 }
01120
01121
01122 const N& x = (coeff < 0) ? dbm[i][j] : dbm[j][i];
01123 if (!is_plus_infinity(x)) {
01124
01125 DIRTY_TEMP(N, d);
01126 const Coefficient& b = expr.inhomogeneous_term();
01127 TEMP_INTEGER(minus_b);
01128 neg_assign(minus_b, b);
01129 const Coefficient& sc_b = maximize ? b : minus_b;
01130 assign_r(d, sc_b, ROUND_UP);
01131
01132
01133 DIRTY_TEMP(N, coeff_expr);
01134 const Coefficient& coeff_i = expr.coefficient(Variable(i-1));
01135 const int sign_i = sgn(coeff_i);
01136 if (sign_i > 0)
01137 assign_r(coeff_expr, coeff_i, ROUND_UP);
01138 else {
01139 TEMP_INTEGER(minus_coeff_i);
01140 neg_assign(minus_coeff_i, coeff_i);
01141 assign_r(coeff_expr, minus_coeff_i, ROUND_UP);
01142 }
01143
01144 add_mul_assign_r(d, coeff_expr, x, ROUND_UP);
01145 numer_denom(d, ext_n, ext_d);
01146 if (!maximize)
01147 neg_assign(ext_n);
01148 included = true;
01149 return true;
01150 }
01151
01152
01153 return false;
01154 }
01155 }
01156
01157 template <typename T>
01158 bool
01159 BD_Shape<T>::max_min(const Linear_Expression& expr,
01160 const bool maximize,
01161 Coefficient& ext_n, Coefficient& ext_d,
01162 bool& included,
01163 Generator& g) const {
01164
01165
01166 const dimension_type space_dim = space_dimension();
01167 const dimension_type expr_space_dim = expr.space_dimension();
01168 if (space_dim < expr_space_dim)
01169 throw_dimension_incompatible((maximize
01170 ? "maximize(e, ...)"
01171 : "minimize(e, ...)"), "e", expr);
01172
01173 if (space_dim == 0) {
01174 if (marked_empty())
01175 return false;
01176 else {
01177 ext_n = expr.inhomogeneous_term();
01178 ext_d = 1;
01179 included = true;
01180 g = point();
01181 return true;
01182 }
01183 }
01184
01185 shortest_path_closure_assign();
01186
01187 if (marked_empty())
01188 return false;
01189
01190 Optimization_Mode mode_max_min
01191 = maximize ? MAXIMIZATION : MINIMIZATION;
01192 MIP_Problem mip(space_dim, constraints(), expr, mode_max_min);
01193 if (mip.solve() == OPTIMIZED_MIP_PROBLEM) {
01194 g = mip.optimizing_point();
01195 mip.evaluate_objective_function(g, ext_n, ext_d);
01196 included = true;
01197 return true;
01198 }
01199
01200 return false;
01201 }
01202
01203 template <typename T>
01204 Poly_Con_Relation
01205 BD_Shape<T>::relation_with(const Congruence& cg) const {
01206 const dimension_type cg_space_dim = cg.space_dimension();
01207 const dimension_type space_dim = space_dimension();
01208
01209
01210 if (cg_space_dim > space_dim)
01211 throw_dimension_incompatible("relation_with(cg)", cg);
01212
01213
01214
01215 if (cg.is_equality()) {
01216 Constraint c(cg);
01217 dimension_type num_vars = 0;
01218 dimension_type i = 0;
01219 dimension_type j = 0;
01220 TEMP_INTEGER(coeff);
01221 if (extract_bounded_difference(c, cg_space_dim, num_vars,
01222 i, j, coeff))
01223 return relation_with(c);
01224 }
01225
01226 shortest_path_closure_assign();
01227
01228 if (marked_empty())
01229 return Poly_Con_Relation::saturates()
01230 && Poly_Con_Relation::is_included()
01231 && Poly_Con_Relation::is_disjoint();
01232
01233 if (space_dim == 0) {
01234 if (cg.is_inconsistent())
01235 return Poly_Con_Relation::is_disjoint();
01236 else if (cg.inhomogeneous_term() % cg.modulus() == 0)
01237 return Poly_Con_Relation::saturates()
01238 && Poly_Con_Relation::is_included();
01239 }
01240
01241 DIRTY_TEMP(Coefficient, min_num);
01242 DIRTY_TEMP(Coefficient, min_den);
01243 bool min_included;
01244 TEMP_INTEGER(mod);
01245 mod = cg.modulus();
01246 Linear_Expression le;
01247 for (dimension_type i = cg_space_dim; i-- > 0; )
01248 le += cg.coefficient(Variable(i)) * Variable(i);
01249 bool bounded_below = minimize(le, min_num, min_den, min_included);
01250
01251 if (!bounded_below)
01252 return Poly_Con_Relation::strictly_intersects();
01253
01254 TEMP_INTEGER(v);
01255 TEMP_INTEGER(lower_num);
01256 TEMP_INTEGER(lower_den);
01257 TEMP_INTEGER(lower);
01258 assign_r(lower_num, min_num, ROUND_NOT_NEEDED);
01259 assign_r(lower_den, min_den, ROUND_NOT_NEEDED);
01260 neg_assign(v, cg.inhomogeneous_term());
01261 lower = lower_num / lower_den;
01262 v += ((lower / mod) * mod);
01263 if (v * lower_den < lower_num)
01264 v += mod;
01265 const Constraint& c(le == v);
01266 return relation_with(c);
01267 }
01268
01269
01270 template <typename T>
01271 Poly_Con_Relation
01272 BD_Shape<T>::relation_with(const Constraint& c) const {
01273 const dimension_type c_space_dim = c.space_dimension();
01274 const dimension_type space_dim = space_dimension();
01275
01276
01277 if (c_space_dim > space_dim)
01278 throw_dimension_incompatible("relation_with(c)", c);
01279
01280 shortest_path_closure_assign();
01281
01282 if (marked_empty())
01283 return Poly_Con_Relation::saturates()
01284 && Poly_Con_Relation::is_included()
01285 && Poly_Con_Relation::is_disjoint();
01286
01287 if (space_dim == 0) {
01288 if ((c.is_equality() && c.inhomogeneous_term() != 0)
01289 || (c.is_inequality() && c.inhomogeneous_term() < 0))
01290 return Poly_Con_Relation::is_disjoint();
01291 else if (c.is_strict_inequality() && c.inhomogeneous_term() == 0)
01292
01293
01294 return Poly_Con_Relation::saturates()
01295 && Poly_Con_Relation::is_disjoint();
01296 else if (c.is_equality() || c.inhomogeneous_term() == 0)
01297 return Poly_Con_Relation::saturates()
01298 && Poly_Con_Relation::is_included();
01299 else
01300
01301
01302
01303 return Poly_Con_Relation::is_included();
01304 }
01305
01306 dimension_type num_vars = 0;
01307 dimension_type i = 0;
01308 dimension_type j = 0;
01309 TEMP_INTEGER(coeff);
01310 if (!extract_bounded_difference(c, c_space_dim, num_vars, i, j, coeff)) {
01311
01312
01313
01314
01315
01316
01317 Linear_Expression le;
01318 for (dimension_type k = c_space_dim; k-- > 0; ) {
01319 Variable vk(k);
01320 le += c.coefficient(vk) * vk;
01321 }
01322 DIRTY_TEMP(Coefficient, max_num);
01323 DIRTY_TEMP(Coefficient, max_den);
01324 bool max_included;
01325 DIRTY_TEMP(Coefficient, min_num);
01326 DIRTY_TEMP(Coefficient, min_den);
01327 bool min_included;
01328 bool bounded_above = maximize(le, max_num, max_den, max_included);
01329 bool bounded_below = minimize(le, min_num, min_den, min_included);
01330 if (!bounded_above) {
01331 if (!bounded_below)
01332 return Poly_Con_Relation::strictly_intersects();
01333 min_num += c.inhomogeneous_term() * min_den;
01334 switch (sgn(min_num)) {
01335 case 1:
01336 if (c.is_equality())
01337 return Poly_Con_Relation::is_disjoint();
01338 return Poly_Con_Relation::is_included();
01339 case 0:
01340 if (c.is_strict_inequality() || c.is_equality())
01341 return Poly_Con_Relation::strictly_intersects();
01342 return Poly_Con_Relation::is_included();
01343 case -1:
01344 return Poly_Con_Relation::strictly_intersects();
01345 }
01346 }
01347 if (!bounded_below) {
01348 max_num += c.inhomogeneous_term() * max_den;
01349 switch (sgn(max_num)) {
01350 case 1:
01351 return Poly_Con_Relation::strictly_intersects();
01352 case 0:
01353 if (c.is_strict_inequality())
01354 return Poly_Con_Relation::is_disjoint();
01355 return Poly_Con_Relation::strictly_intersects();
01356 case -1:
01357 return Poly_Con_Relation::is_disjoint();
01358 }
01359 }
01360 else {
01361 max_num += c.inhomogeneous_term() * max_den;
01362 min_num += c.inhomogeneous_term() * min_den;
01363 switch (sgn(max_num)) {
01364 case 1:
01365 switch (sgn(min_num)) {
01366 case 1:
01367 if (c.is_equality())
01368 return Poly_Con_Relation::is_disjoint();
01369 return Poly_Con_Relation::is_included();
01370 case 0:
01371 if (c.is_equality())
01372 return Poly_Con_Relation::strictly_intersects();
01373 if (c.is_strict_inequality())
01374 return Poly_Con_Relation::strictly_intersects();
01375 return Poly_Con_Relation::is_included();
01376 case -1:
01377 return Poly_Con_Relation::strictly_intersects();
01378 }
01379 case 0:
01380 if (min_num == 0) {
01381 if (c.is_strict_inequality())
01382 return Poly_Con_Relation::is_disjoint()
01383 && Poly_Con_Relation::saturates();
01384 return Poly_Con_Relation::is_included()
01385 && Poly_Con_Relation::saturates();
01386 }
01387 if (c.is_strict_inequality())
01388 return Poly_Con_Relation::is_disjoint();
01389 return Poly_Con_Relation::strictly_intersects();
01390 case -1:
01391 return Poly_Con_Relation::is_disjoint();
01392 }
01393 }
01394 }
01395
01396
01397 if (num_vars == 0) {
01398
01399 switch (sgn(c.inhomogeneous_term())) {
01400 case -1:
01401 return Poly_Con_Relation::is_disjoint();
01402 case 0:
01403 if (c.is_strict_inequality())
01404 return Poly_Con_Relation::saturates()
01405 && Poly_Con_Relation::is_disjoint();
01406 else
01407 return Poly_Con_Relation::saturates()
01408 && Poly_Con_Relation::is_included();
01409 case 1:
01410 if (c.is_equality())
01411 return Poly_Con_Relation::is_disjoint();
01412 else
01413 return Poly_Con_Relation::is_included();
01414 }
01415 }
01416
01417
01418
01419 const bool negative = (coeff < 0);
01420 const N& x = negative ? dbm[i][j] : dbm[j][i];
01421 const N& y = negative ? dbm[j][i] : dbm[i][j];
01422 if (negative)
01423 neg_assign(coeff);
01424
01425
01426
01427
01428
01429
01430
01431
01432 DIRTY_TEMP0(mpq_class, q_x);
01433 DIRTY_TEMP0(mpq_class, q_y);
01434 DIRTY_TEMP0(mpq_class, d);
01435 DIRTY_TEMP0(mpq_class, d1);
01436 DIRTY_TEMP0(mpq_class, c_den);
01437 DIRTY_TEMP0(mpq_class, q_den);
01438 assign_r(c_den, coeff, ROUND_NOT_NEEDED);
01439 assign_r(d, c.inhomogeneous_term(), ROUND_NOT_NEEDED);
01440 neg_assign_r(d1, d, ROUND_NOT_NEEDED);
01441 div_assign_r(d, d, c_den, ROUND_NOT_NEEDED);
01442 div_assign_r(d1, d1, c_den, ROUND_NOT_NEEDED);
01443
01444 if (is_plus_infinity(x)) {
01445 if (!is_plus_infinity(y)) {
01446
01447
01448
01449
01450
01451 TEMP_INTEGER(numer);
01452 TEMP_INTEGER(denom);
01453 numer_denom(y, numer, denom);
01454 assign_r(q_den, denom, ROUND_NOT_NEEDED);
01455 assign_r(q_y, numer, ROUND_NOT_NEEDED);
01456 div_assign_r(q_y, q_y, q_den, ROUND_NOT_NEEDED);
01457 if (q_y < d1)
01458 return Poly_Con_Relation::is_disjoint();
01459 if (q_y == d1 && c.is_strict_inequality())
01460 return Poly_Con_Relation::is_disjoint();
01461 }
01462
01463
01464 return Poly_Con_Relation::strictly_intersects();
01465 }
01466
01467
01468 TEMP_INTEGER(numer);
01469 TEMP_INTEGER(denom);
01470 numer_denom(x, numer, denom);
01471 assign_r(q_den, denom, ROUND_NOT_NEEDED);
01472 assign_r(q_x, numer, ROUND_NOT_NEEDED);
01473 div_assign_r(q_x, q_x, q_den, ROUND_NOT_NEEDED);
01474
01475 if (!is_plus_infinity(y)) {
01476 numer_denom(y, numer, denom);
01477 assign_r(q_den, denom, ROUND_NOT_NEEDED);
01478 assign_r(q_y, numer, ROUND_NOT_NEEDED);
01479 div_assign_r(q_y, q_y, q_den, ROUND_NOT_NEEDED);
01480 if (q_x == d && q_y == d1) {
01481 if (c.is_strict_inequality())
01482 return Poly_Con_Relation::saturates()
01483 && Poly_Con_Relation::is_disjoint();
01484 else
01485 return Poly_Con_Relation::saturates()
01486 && Poly_Con_Relation::is_included();
01487 }
01488
01489
01490
01491 if (q_y < d1)
01492 return Poly_Con_Relation::is_disjoint();
01493 if (q_y == d1 && c.is_strict_inequality())
01494 return Poly_Con_Relation::is_disjoint();
01495 }
01496
01497
01498
01499
01500 if (d > q_x) {
01501 if (c.is_equality())
01502 return Poly_Con_Relation::is_disjoint();
01503 else
01504 return Poly_Con_Relation::is_included();
01505 }
01506
01507 if (d == q_x && c.is_nonstrict_inequality())
01508 return Poly_Con_Relation::is_included();
01509
01510
01511 return Poly_Con_Relation::strictly_intersects();
01512 }
01513
01514 template <typename T>
01515 Poly_Gen_Relation
01516 BD_Shape<T>::relation_with(const Generator& g) const {
01517 const dimension_type space_dim = space_dimension();
01518 const dimension_type g_space_dim = g.space_dimension();
01519
01520
01521 if (space_dim < g_space_dim)
01522 throw_dimension_incompatible("relation_with(g)", g);
01523
01524 shortest_path_closure_assign();
01525
01526 if (marked_empty())
01527 return Poly_Gen_Relation::nothing();
01528
01529
01530
01531 if (space_dim == 0)
01532 return Poly_Gen_Relation::subsumes();
01533
01534 const bool is_line = g.is_line();
01535 const bool is_line_or_ray = g.is_line_or_ray();
01536
01537
01538
01539
01540
01541
01542
01543
01544 TEMP_INTEGER(num);
01545 TEMP_INTEGER(den);
01546 TEMP_INTEGER(product);
01547
01548 for (dimension_type i = 0; i <= space_dim; ++i) {
01549 const Coefficient& g_coeff_y = (i > g_space_dim || i == 0)
01550 ? Coefficient(0) : g.coefficient(Variable(i-1));
01551 const DB_Row<N>& dbm_i = dbm[i];
01552 for (dimension_type j = i + 1; j <= space_dim; ++j) {
01553 const Coefficient& g_coeff_x = (j > g_space_dim)
01554 ? Coefficient(0) : g.coefficient(Variable(j-1));
01555 const N& dbm_ij = dbm_i[j];
01556 const N& dbm_ji = dbm[j][i];
01557 if (is_additive_inverse(dbm_ji, dbm_ij)) {
01558
01559
01560 numer_denom(dbm_ij, num, den);
01561 product = 0;
01562 add_mul_assign(product, den, g_coeff_y);
01563 add_mul_assign(product, -den, g_coeff_x);
01564 if (!is_line_or_ray)
01565 add_mul_assign(product, num, g.divisor());
01566 if (product != 0)
01567 return Poly_Gen_Relation::nothing();
01568 }
01569 else {
01570
01571 if (!is_plus_infinity(dbm_ij)) {
01572
01573
01574 numer_denom(dbm_ij, num, den);
01575 product = 0;
01576 add_mul_assign(product, den, g_coeff_y);
01577 add_mul_assign(product, -den, g_coeff_x);
01578 if (!is_line_or_ray)
01579 add_mul_assign(product, num, g.divisor());
01580 if (is_line) {
01581 if (product != 0)
01582
01583 return Poly_Gen_Relation::nothing();
01584 }
01585 else
01586
01587 if (product < 0)
01588 return Poly_Gen_Relation::nothing();
01589 }
01590
01591 if (!is_plus_infinity(dbm_ji)) {
01592
01593
01594 numer_denom(dbm_ji, num, den);
01595 product = 0;
01596 add_mul_assign(product, den, g_coeff_x);
01597 add_mul_assign(product, -den, g_coeff_y);
01598 if (!is_line_or_ray)
01599 add_mul_assign(product, num, g.divisor());
01600 if (is_line) {
01601 if (product != 0)
01602
01603 return Poly_Gen_Relation::nothing();
01604 }
01605 else
01606
01607 if (product < 0)
01608 return Poly_Gen_Relation::nothing();
01609 }
01610 }
01611 }
01612 }
01613
01614
01615 return Poly_Gen_Relation::subsumes();
01616 }
01617
01618 template <typename T>
01619 void
01620 BD_Shape<T>::shortest_path_closure_assign() const {
01621
01622 if (marked_empty() || marked_shortest_path_closed())
01623 return;
01624 const dimension_type num_dimensions = space_dimension();
01625
01626 if (num_dimensions == 0)
01627 return;
01628
01629
01630
01631 BD_Shape& x = const_cast<BD_Shape<T>&>(*this);
01632
01633
01634 for (dimension_type h = num_dimensions + 1; h-- > 0; ) {
01635 assert(is_plus_infinity(x.dbm[h][h]));
01636 assign_r(x.dbm[h][h], 0, ROUND_NOT_NEEDED);
01637 }
01638
01639 DIRTY_TEMP(N, sum);
01640 for (dimension_type k = num_dimensions + 1; k-- > 0; ) {
01641 const DB_Row<N>& x_dbm_k = x.dbm[k];
01642 for (dimension_type i = num_dimensions + 1; i-- > 0; ) {
01643 DB_Row<N>& x_dbm_i = x.dbm[i];
01644 const N& x_dbm_i_k = x_dbm_i[k];
01645 if (!is_plus_infinity(x_dbm_i_k))
01646 for (dimension_type j = num_dimensions + 1; j-- > 0; ) {
01647 const N& x_dbm_k_j = x_dbm_k[j];
01648 if (!is_plus_infinity(x_dbm_k_j)) {
01649
01650 add_assign_r(sum, x_dbm_i_k, x_dbm_k_j, ROUND_UP);
01651 min_assign(x_dbm_i[j], sum);
01652 }
01653 }
01654 }
01655 }
01656
01657
01658
01659 for (dimension_type h = num_dimensions + 1; h-- > 0; ) {
01660 N& x_dbm_hh = x.dbm[h][h];
01661 if (sgn(x_dbm_hh) < 0) {
01662 x.set_empty();
01663 return;
01664 }
01665 else {
01666 assert(sgn(x_dbm_hh) == 0);
01667
01668 assign_r(x_dbm_hh, PLUS_INFINITY, ROUND_NOT_NEEDED);
01669 }
01670 }
01671
01672
01673 x.set_shortest_path_closed();
01674 }
01675
01676 template <typename T>
01677 void
01678 BD_Shape<T>::shortest_path_reduction_assign() const {
01679
01680 if (marked_shortest_path_reduced())
01681 return;
01682
01683 const dimension_type space_dim = space_dimension();
01684
01685 if (space_dim == 0)
01686 return;
01687
01688
01689 shortest_path_closure_assign();
01690
01691
01692 if (marked_empty())
01693 return;
01694
01695
01696
01697
01698
01699 std::vector<dimension_type> predecessor;
01700 compute_predecessors(predecessor);
01701 std::vector<dimension_type> leaders;
01702 compute_leader_indices(predecessor, leaders);
01703 const dimension_type num_leaders = leaders.size();
01704
01705 Bit_Matrix redundancy(space_dim + 1, space_dim + 1);
01706
01707
01708 Bit_Row& red_0 = redundancy[0];
01709 for (dimension_type j = space_dim + 1; j-- > 0; )
01710 red_0.set(j);
01711 for (dimension_type i = space_dim + 1; i-- > 0; )
01712 redundancy[i] = red_0;
01713
01714
01715
01716 DIRTY_TEMP(N, c);
01717 for (dimension_type l_i = 0; l_i < num_leaders; ++l_i) {
01718 const dimension_type i = leaders[l_i];
01719 const DB_Row<N>& dbm_i = dbm[i];
01720 Bit_Row& redundancy_i = redundancy[i];
01721 for (dimension_type l_j = 0; l_j < num_leaders; ++l_j) {
01722 const dimension_type j = leaders[l_j];
01723 if (redundancy_i[j]) {
01724 const N& dbm_i_j = dbm_i[j];
01725 redundancy_i.clear(j);
01726 for (dimension_type l_k = 0; l_k < num_leaders; ++l_k) {
01727 const dimension_type k = leaders[l_k];
01728 add_assign_r(c, dbm_i[k], dbm[k][j], ROUND_UP);
01729 if (dbm_i_j >= c) {
01730 redundancy_i.set(j);
01731 break;
01732 }
01733 }
01734 }
01735 }
01736 }
01737
01738
01739
01740
01741 std::deque<bool> dealt_with(space_dim + 1, false);
01742 for (dimension_type i = space_dim + 1; i-- > 0; )
01743
01744
01745 if (i != predecessor[i] && !dealt_with[i]) {
01746 dimension_type j = i;
01747 while (true) {
01748 const dimension_type pred_j = predecessor[j];
01749 if (j == pred_j) {
01750
01751 assert(redundancy[i][j]);
01752 redundancy[i].clear(j);
01753
01754
01755 break;
01756 }
01757
01758 assert(redundancy[pred_j][j]);
01759 redundancy[pred_j].clear(j);
01760 dealt_with[pred_j] = true;
01761 j = pred_j;
01762 }
01763 }
01764
01765
01766
01767 BD_Shape<T>& x = const_cast<BD_Shape<T>&>(*this);
01768 std::swap(x.redundancy_dbm, redundancy);
01769 x.set_shortest_path_reduced();
01770
01771 assert(is_shortest_path_reduced());
01772 }
01773
01774 template <typename T>
01775 void
01776 BD_Shape<T>::upper_bound_assign(const BD_Shape& y) {
01777 const dimension_type space_dim = space_dimension();
01778
01779
01780 if (space_dim != y.space_dimension())
01781 throw_dimension_incompatible("upper_bound_assign(y)", y);
01782
01783
01784 y.shortest_path_closure_assign();
01785 if (y.marked_empty())
01786 return;
01787 shortest_path_closure_assign();
01788 if (marked_empty()) {
01789 *this = y;
01790 return;
01791 }
01792
01793
01794
01795 assert(space_dim == 0 || marked_shortest_path_closed());
01796 for (dimension_type i = space_dim + 1; i-- > 0; ) {
01797 DB_Row<N>& dbm_i = dbm[i];
01798 const DB_Row<N>& y_dbm_i = y.dbm[i];
01799 for (dimension_type j = space_dim + 1; j-- > 0; ) {
01800 N& dbm_ij = dbm_i[j];
01801 const N& y_dbm_ij = y_dbm_i[j];
01802 if (dbm_ij < y_dbm_ij)
01803 dbm_ij = y_dbm_ij;
01804 }
01805 }
01806
01807
01808 if (marked_shortest_path_reduced())
01809 reset_shortest_path_reduced();
01810 assert(OK());
01811 }
01812
01813 template <typename T>
01814 void
01815 BD_Shape<T>::difference_assign(const BD_Shape& y) {
01816 const dimension_type space_dim = space_dimension();
01817
01818
01819 if (space_dim != y.space_dimension())
01820 throw_dimension_incompatible("difference_assign(y)", y);
01821
01822 BD_Shape new_bd_shape(space_dim, EMPTY);
01823
01824 BD_Shape& x = *this;
01825
01826 x.shortest_path_closure_assign();
01827
01828
01829 if (x.marked_empty())
01830 return;
01831 y.shortest_path_closure_assign();
01832
01833
01834 if (y.marked_empty())
01835 return;
01836
01837
01838
01839
01840 if (space_dim == 0) {
01841 x.set_empty();
01842 return;
01843 }
01844
01845
01846
01847 if (y.contains(x)) {
01848 x.set_empty();
01849 return;
01850 }
01851
01852
01853
01854
01855 const Constraint_System& y_cs = y.constraints();
01856 for (Constraint_System::const_iterator i = y_cs.begin(),
01857 y_cs_end = y_cs.end(); i != y_cs_end; ++i) {
01858 const Constraint& c = *i;
01859
01860
01861
01862
01863
01864
01865 if (x.relation_with(c).implies(Poly_Con_Relation::is_included()))
01866 continue;
01867 BD_Shape z = x;
01868 const Linear_Expression e = Linear_Expression(c);
01869 z.add_constraint(e <= 0);
01870 if (!z.is_empty())
01871 new_bd_shape.upper_bound_assign(z);
01872 if (c.is_equality()) {
01873 z = x;
01874 z.add_constraint(e >= 0);
01875 if (!z.is_empty())
01876 new_bd_shape.upper_bound_assign(z);
01877 }
01878 }
01879 *this = new_bd_shape;
01880 assert(OK());
01881 }
01882
01883 template <typename T>
01884 bool
01885 BD_Shape<T>::simplify_using_context_assign(const BD_Shape& y) {
01886
01887 used(y);
01888 return true;
01889 }
01890
01891 template <typename T>
01892 void
01893 BD_Shape<T>::add_space_dimensions_and_embed(const dimension_type m) {
01894
01895 if (m == 0)
01896 return;
01897
01898 const dimension_type space_dim = space_dimension();
01899 const dimension_type new_space_dim = space_dim + m;
01900 const bool was_zero_dim_univ = (!marked_empty() && space_dim == 0);
01901
01902
01903
01904
01905 dbm.grow(new_space_dim + 1);
01906
01907
01908
01909 if (marked_shortest_path_reduced())
01910 reset_shortest_path_reduced();
01911
01912
01913
01914 if (was_zero_dim_univ)
01915 set_shortest_path_closed();
01916
01917 assert(OK());
01918 }
01919
01920 template <typename T>
01921 void
01922 BD_Shape<T>::add_space_dimensions_and_project(const dimension_type m) {
01923
01924 if (m == 0)
01925 return;
01926
01927 const dimension_type space_dim = space_dimension();
01928
01929
01930
01931
01932 if (space_dim == 0) {
01933 dbm.grow(m + 1);
01934 if (!marked_empty()) {
01935 for (dimension_type i = m + 1; i-- > 0; ) {
01936 DB_Row<N>& dbm_i = dbm[i];
01937 for (dimension_type j = m + 1; j-- > 0; )
01938 if (i != j)
01939 assign_r(dbm_i[j], 0, ROUND_NOT_NEEDED);
01940 }
01941 set_shortest_path_closed();
01942 }
01943 assert(OK());
01944 return;
01945 }
01946
01947
01948
01949
01950
01951 const dimension_type new_space_dim = space_dim + m;
01952 dbm.grow(new_space_dim + 1);
01953
01954
01955 DB_Row<N>& dbm_0 = dbm[0];
01956 for (dimension_type i = space_dim + 1; i <= new_space_dim; ++i) {
01957 assign_r(dbm[i][0], 0, ROUND_NOT_NEEDED);
01958 assign_r(dbm_0[i], 0, ROUND_NOT_NEEDED);
01959 }
01960
01961 if (marked_shortest_path_closed())
01962 reset_shortest_path_closed();
01963 assert(OK());
01964 }
01965
01966 template <typename T>
01967 void
01968 BD_Shape<T>::remove_space_dimensions(const Variables_Set& to_be_removed) {
01969
01970
01971
01972 if (to_be_removed.empty()) {
01973 assert(OK());
01974 return;
01975 }
01976
01977 const dimension_type old_space_dim = space_dimension();
01978
01979
01980 const dimension_type min_space_dim = to_be_removed.space_dimension();
01981 if (old_space_dim < min_space_dim)
01982 throw_dimension_incompatible("remove_space_dimensions(vs)", min_space_dim);
01983
01984
01985 shortest_path_closure_assign();
01986
01987
01988
01989 const dimension_type new_space_dim = old_space_dim - to_be_removed.size();
01990 if (new_space_dim == 0) {
01991 dbm.resize_no_copy(1);
01992 if (!marked_empty())
01993
01994 set_zero_dim_univ();
01995 assert(OK());
01996 return;
01997 }
01998
01999
02000 if (marked_empty()) {
02001 dbm.resize_no_copy(new_space_dim + 1);
02002 assert(OK());
02003 return;
02004 }
02005
02006
02007
02008 if (marked_shortest_path_reduced())
02009 reset_shortest_path_reduced();
02010
02011
02012
02013
02014 Variables_Set::const_iterator tbr = to_be_removed.begin();
02015 Variables_Set::const_iterator tbr_end = to_be_removed.end();
02016 dimension_type dst = *tbr + 1;
02017 dimension_type src = dst + 1;
02018 for (++tbr; tbr != tbr_end; ++tbr) {
02019 const dimension_type tbr_next = *tbr + 1;
02020
02021
02022 while (src < tbr_next) {
02023 std::swap(dbm[dst], dbm[src]);
02024 for (dimension_type i = old_space_dim + 1; i-- > 0; ) {
02025 DB_Row<N>& dbm_i = dbm[i];
02026 assign_or_swap(dbm_i[dst], dbm_i[src]);
02027 }
02028 ++dst;
02029 ++src;
02030 }
02031 ++src;
02032 }
02033
02034
02035 while (src <= old_space_dim) {
02036 std::swap(dbm[dst], dbm[src]);
02037 for (dimension_type i = old_space_dim + 1; i-- > 0; ) {
02038 DB_Row<N>& dbm_i = dbm[i];
02039 assign_or_swap(dbm_i[dst], dbm_i[src]);
02040 }
02041 ++src;
02042 ++dst;
02043 }
02044
02045
02046 dbm.resize_no_copy(new_space_dim + 1);
02047 assert(OK());
02048 }
02049
02050 template <typename T>
02051 template <typename Partial_Function>
02052 void
02053 BD_Shape<T>::map_space_dimensions(const Partial_Function& pfunc) {
02054 const dimension_type space_dim = space_dimension();
02055
02056 if (space_dim == 0)
02057 return;
02058
02059 if (pfunc.has_empty_codomain()) {
02060
02061 remove_higher_space_dimensions(0);
02062 return;
02063 }
02064
02065 const dimension_type new_space_dim = pfunc.max_in_codomain() + 1;
02066
02067
02068 if (new_space_dim < space_dim)
02069 shortest_path_closure_assign();
02070
02071
02072
02073 if (marked_empty()) {
02074 remove_higher_space_dimensions(new_space_dim);
02075 return;
02076 }
02077
02078
02079
02080 if (marked_shortest_path_reduced())
02081 reset_shortest_path_reduced();
02082
02083
02084 DB_Matrix<N> x(new_space_dim+1);
02085
02086
02087
02088 DB_Row<N>& dbm_0 = dbm[0];
02089 DB_Row<N>& x_0 = x[0];
02090 for (dimension_type j = 1; j <= space_dim; ++j) {
02091 dimension_type new_j;
02092 if (pfunc.maps(j - 1, new_j)) {
02093 assign_or_swap(x_0[new_j + 1], dbm_0[j]);
02094 assign_or_swap(x[new_j + 1][0], dbm[j][0]);
02095 }
02096 }
02097
02098 for (dimension_type i = 1; i <= space_dim; ++i) {
02099 dimension_type new_i;
02100 if (pfunc.maps(i - 1, new_i)) {
02101 DB_Row<N>& dbm_i = dbm[i];
02102 ++new_i;
02103 DB_Row<N>& x_new_i = x[new_i];
02104 for (dimension_type j = i+1; j <= space_dim; ++j) {
02105 dimension_type new_j;
02106 if (pfunc.maps(j - 1, new_j)) {
02107 ++new_j;
02108 assign_or_swap(x_new_i[new_j], dbm_i[j]);
02109 assign_or_swap(x[new_j][new_i], dbm[j][i]);
02110 }
02111 }
02112 }
02113 }
02114
02115 std::swap(dbm, x);
02116 assert(OK());
02117 }
02118
02119 template <typename T>
02120 void
02121 BD_Shape<T>::intersection_assign(const BD_Shape& y) {
02122 const dimension_type space_dim = space_dimension();
02123
02124
02125 if (space_dim != y.space_dimension())
02126 throw_dimension_incompatible("intersection_assign(y)", y);
02127
02128
02129
02130 if (marked_empty())
02131 return;
02132 if (y.marked_empty()) {
02133 set_empty();
02134 return;
02135 }
02136
02137
02138
02139
02140 if (space_dim == 0)
02141 return;
02142
02143
02144
02145 bool changed = false;
02146 for (dimension_type i = space_dim + 1; i-- > 0; ) {
02147 DB_Row<N>& dbm_i = dbm[i];
02148 const DB_Row<N>& y_dbm_i = y.dbm[i];
02149 for (dimension_type j = space_dim + 1; j-- > 0; ) {
02150 N& dbm_ij = dbm_i[j];
02151 const N& y_dbm_ij = y_dbm_i[j];
02152 if (dbm_ij > y_dbm_ij) {
02153 dbm_ij = y_dbm_ij;
02154 changed = true;
02155 }
02156 }
02157 }
02158
02159 if (changed && marked_shortest_path_closed())
02160 reset_shortest_path_closed();
02161 assert(OK());
02162 }
02163
02164 template <typename T>
02165 template <typename Iterator>
02166 void
02167 BD_Shape<T>::CC76_extrapolation_assign(const BD_Shape& y,
02168 Iterator first, Iterator last,
02169 unsigned* tp) {
02170 const dimension_type space_dim = space_dimension();
02171
02172
02173 if (space_dim != y.space_dimension())
02174 throw_dimension_incompatible("CC76_extrapolation_assign(y)", y);
02175
02176 #ifndef NDEBUG
02177 {
02178
02179 const BD_Shape x_copy = *this;
02180 const BD_Shape y_copy = y;
02181 assert(x_copy.contains(y_copy));
02182 }
02183 #endif
02184
02185
02186
02187 if (space_dim == 0)
02188 return;
02189
02190 shortest_path_closure_assign();
02191
02192 if (marked_empty())
02193 return;
02194 y.shortest_path_closure_assign();
02195
02196 if (y.marked_empty())
02197 return;
02198
02199
02200 if (tp != 0 && *tp > 0) {
02201 BD_Shape<T> x_tmp(*this);
02202 x_tmp.CC76_extrapolation_assign(y, first, last, 0);
02203
02204 if (!contains(x_tmp))
02205 --(*tp);
02206 return;
02207 }
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217 for (dimension_type i = space_dim + 1; i-- > 0; ) {
02218 DB_Row<N>& dbm_i = dbm[i];
02219 const DB_Row<N>& y_dbm_i = y.dbm[i];
02220 for (dimension_type j = space_dim + 1; j-- > 0; ) {
02221 N& dbm_ij = dbm_i[j];
02222 const N& y_dbm_ij = y_dbm_i[j];
02223 if (y_dbm_ij < dbm_ij) {
02224 Iterator k = std::lower_bound(first, last, dbm_ij);
02225 if (k != last) {
02226 if (dbm_ij < *k)
02227 assign_r(dbm_ij, *k, ROUND_UP);
02228 }
02229 else
02230 assign_r(dbm_ij, PLUS_INFINITY, ROUND_NOT_NEEDED);
02231 }
02232 }
02233 }
02234 reset_shortest_path_closed();
02235 assert(OK());
02236 }
02237
02238 template <typename T>
02239 void
02240 BD_Shape<T>::get_limiting_shape(const Constraint_System& cs,
02241 BD_Shape& limiting_shape) const {
02242 const dimension_type cs_space_dim = cs.space_dimension();
02243
02244 assert(cs_space_dim <= space_dimension());
02245
02246 shortest_path_closure_assign();
02247 bool changed = false;
02248 TEMP_INTEGER(coeff);
02249 TEMP_INTEGER(minus_c_term);
02250 DIRTY_TEMP(N, d);
02251 DIRTY_TEMP(N, d1);
02252 for (Constraint_System::const_iterator cs_i = cs.begin(),
02253 cs_end = cs.end(); cs_i != cs_end; ++cs_i) {
02254 const Constraint& c = *cs_i;
02255 dimension_type num_vars = 0;
02256 dimension_type i = 0;
02257 dimension_type j = 0;
02258
02259 if (extract_bounded_difference(c, cs_space_dim, num_vars, i, j, coeff)) {
02260
02261
02262 const bool negative = (coeff < 0);
02263 const N& x = negative ? dbm[i][j] : dbm[j][i];
02264 const N& y = negative ? dbm[j][i] : dbm[i][j];
02265 DB_Matrix<N>& ls_dbm = limiting_shape.dbm;
02266 N& ls_x = negative ? ls_dbm[i][j] : ls_dbm[j][i];
02267 N& ls_y = negative ? ls_dbm[j][i] : ls_dbm[i][j];
02268 if (negative)
02269 neg_assign(coeff);
02270
02271 div_round_up(d, c.inhomogeneous_term(), coeff);
02272 if (x <= d) {
02273 if (c.is_inequality()) {
02274 if (ls_x > d) {
02275 ls_x = d;
02276 changed = true;
02277 }
02278 }
02279 else {
02280
02281 neg_assign(minus_c_term, c.inhomogeneous_term());
02282 div_round_up(d1, minus_c_term, coeff);
02283 if (y <= d1)
02284 if((ls_x >= d && ls_y > d1) || (ls_x > d && ls_y >= d1)) {
02285 ls_x = d;
02286 ls_y = d1;
02287 changed = true;
02288 }
02289 }
02290 }
02291 }
02292 }
02293
02294
02295
02296 if (changed && limiting_shape.marked_shortest_path_closed())
02297 limiting_shape.reset_shortest_path_closed();
02298 }
02299
02300 template <typename T>
02301 void
02302 BD_Shape<T>::limited_CC76_extrapolation_assign(const BD_Shape& y,
02303 const Constraint_System& cs,
02304 unsigned* tp) {
02305
02306 const dimension_type space_dim = space_dimension();
02307 if (space_dim != y.space_dimension())
02308 throw_dimension_incompatible("limited_CC76_extrapolation_assign(y, cs)",
02309 y);
02310
02311
02312
02313 const dimension_type cs_space_dim = cs.space_dimension();
02314 if (space_dim < cs_space_dim)
02315 throw_generic("limited_CC76_extrapolation_assign(y, cs)",
02316 "cs is space_dimension incompatible");
02317
02318
02319 if (cs.has_strict_inequalities())
02320 throw_generic("limited_CC76_extrapolation_assign(y, cs)",
02321 "cs has strict inequalities");
02322
02323
02324
02325
02326 if (space_dim == 0)
02327 return;
02328
02329 #ifndef NDEBUG
02330 {
02331
02332 const BD_Shape x_copy = *this;
02333 const BD_Shape y_copy = y;
02334 assert(x_copy.contains(y_copy));
02335 }
02336 #endif
02337
02338
02339 if (marked_empty())
02340 return;
02341
02342 if (y.marked_empty())
02343 return;
02344
02345 BD_Shape<T> limiting_shape(space_dim, UNIVERSE);
02346 get_limiting_shape(cs, limiting_shape);
02347 CC76_extrapolation_assign(y, tp);
02348 intersection_assign(limiting_shape);
02349 }
02350
02351 template <typename T>
02352 void
02353 BD_Shape<T>::BHMZ05_widening_assign(const BD_Shape& y, unsigned* tp) {
02354 const dimension_type space_dim = space_dimension();
02355
02356
02357 if (space_dim != y.space_dimension())
02358 throw_dimension_incompatible("BHMZ05_widening_assign(y)", y);
02359
02360 #ifndef NDEBUG
02361 {
02362
02363 const BD_Shape x_copy = *this;
02364 const BD_Shape y_copy = y;
02365 assert(x_copy.contains(y_copy));
02366 }
02367 #endif
02368
02369
02370 const dimension_type y_affine_dim = y.affine_dimension();
02371
02372
02373
02374 if (y_affine_dim == 0)
02375 return;
02376
02377
02378
02379 const dimension_type x_affine_dim = affine_dimension();
02380 assert(x_affine_dim >= y_affine_dim);
02381 if (x_affine_dim != y_affine_dim)
02382 return;
02383
02384
02385 if (tp != 0 && *tp > 0) {
02386 BD_Shape<T> x_tmp(*this);
02387 x_tmp.BHMZ05_widening_assign(y, 0);
02388
02389 if (!contains(x_tmp))
02390 --(*tp);
02391 return;
02392 }
02393
02394
02395 assert(marked_shortest_path_closed() && y.marked_shortest_path_closed());
02396
02397 y.shortest_path_reduction_assign();
02398
02399
02400 for (dimension_type i = space_dim + 1; i-- > 0; ) {
02401 DB_Row<N>& dbm_i = dbm[i];
02402 const DB_Row<N>& y_dbm_i = y.dbm[i];
02403 const Bit_Row& y_redundancy_i = y.redundancy_dbm[i];
02404 for (dimension_type j = space_dim + 1; j-- > 0; ) {
02405 N& dbm_ij = dbm_i[j];
02406
02407
02408
02409 if (y_redundancy_i[j] || y_dbm_i[j] != dbm_ij)
02410 assign_r(dbm_ij, PLUS_INFINITY, ROUND_NOT_NEEDED);
02411 }
02412 }
02413
02414
02415
02416
02417 reset_shortest_path_closed();
02418 assert(OK());
02419 }
02420
02421 template <typename T>
02422 void
02423 BD_Shape<T>::limited_BHMZ05_extrapolation_assign(const BD_Shape& y,
02424 const Constraint_System& cs,
02425 unsigned* tp) {
02426
02427 const dimension_type space_dim = space_dimension();
02428 if (space_dim != y.space_dimension())
02429 throw_dimension_incompatible("limited_BHMZ05_extrapolation_assign(y, cs)",
02430 y);
02431
02432
02433 const dimension_type cs_space_dim = cs.space_dimension();
02434 if (space_dim < cs_space_dim)
02435 throw_generic("limited_BHMZ05_extrapolation_assign(y, cs)",
02436 "cs is space-dimension incompatible");
02437
02438
02439 if (cs.has_strict_inequalities())
02440 throw_generic("limited_BHMZ05_extrapolation_assign(y, cs)",
02441 "cs has strict inequalities");
02442
02443
02444
02445
02446 if (space_dim == 0)
02447 return;
02448
02449 #ifndef NDEBUG
02450 {
02451
02452 const BD_Shape x_copy = *this;
02453 const BD_Shape y_copy = y;
02454 assert(x_copy.contains(y_copy));
02455 }
02456 #endif
02457
02458
02459 if (marked_empty())
02460 return;
02461
02462 if (y.marked_empty())
02463 return;
02464
02465 BD_Shape<T> limiting_shape(space_dim, UNIVERSE);
02466 get_limiting_shape(cs, limiting_shape);
02467 BHMZ05_widening_assign(y, tp);
02468 intersection_assign(limiting_shape);
02469 }
02470
02471 template <typename T>
02472 void
02473 BD_Shape<T>::CC76_narrowing_assign(const BD_Shape& y) {
02474 const dimension_type space_dim = space_dimension();
02475
02476
02477 if (space_dim != y.space_dimension())
02478 throw_dimension_incompatible("CC76_narrowing_assign(y)", y);
02479
02480 #ifndef NDEBUG
02481 {
02482
02483 const BD_Shape x_copy = *this;
02484 const BD_Shape y_copy = y;
02485 assert(y_copy.contains(x_copy));
02486 }
02487 #endif
02488
02489
02490
02491 if (space_dim == 0)
02492 return;
02493
02494 y.shortest_path_closure_assign();
02495
02496 if (y.marked_empty())
02497 return;
02498 shortest_path_closure_assign();
02499
02500 if (marked_empty())
02501 return;
02502
02503
02504
02505 bool changed = false;
02506 for (dimension_type i = space_dim + 1; i-- > 0; ) {
02507 DB_Row<N>& dbm_i = dbm[i];
02508 const DB_Row<N>& y_dbm_i = y.dbm[i];
02509 for (dimension_type j = space_dim + 1; j-- > 0; ) {
02510 N& dbm_ij = dbm_i[j];
02511 const N& y_dbm_ij = y_dbm_i[j];
02512 if (!is_plus_infinity(dbm_ij)
02513 && !is_plus_infinity(y_dbm_ij)
02514 && dbm_ij != y_dbm_ij) {
02515 dbm_ij = y_dbm_ij;
02516 changed = true;
02517 }
02518 }
02519 }
02520 if (changed && marked_shortest_path_closed())
02521 reset_shortest_path_closed();
02522 assert(OK());
02523 }
02524
02525 template <typename T>
02526 void
02527 BD_Shape<T>
02528 ::deduce_v_minus_u_bounds(const dimension_type v,
02529 const dimension_type last_v,
02530 const Linear_Expression& sc_expr,
02531 Coefficient_traits::const_reference sc_den,
02532 const N& ub_v) {
02533 assert(sc_den > 0);
02534 assert(!is_plus_infinity(ub_v));
02535
02536
02537
02538
02539
02540
02541
02542
02543 DIRTY_TEMP0(mpq_class, mpq_sc_den);
02544 assign_r(mpq_sc_den, sc_den, ROUND_NOT_NEEDED);
02545 const DB_Row<N>& dbm_0 = dbm[0];
02546
02547 DIRTY_TEMP0(mpq_class, minus_lb_u);
02548 DIRTY_TEMP0(mpq_class, q);
02549 DIRTY_TEMP0(mpq_class, ub_u);
02550 DIRTY_TEMP(N, up_approx);
02551
02552 for (dimension_type u = last_v; u > 0; --u)
02553 if (u != v) {
02554 const Coefficient& expr_u = sc_expr.coefficient(Variable(u-1));
02555 if (expr_u > 0) {
02556 if (expr_u >= sc_den)
02557
02558 sub_assign_r(dbm[u][v], ub_v, dbm_0[u], ROUND_UP);
02559 else {
02560 DB_Row<N>& dbm_u = dbm[u];
02561 const N& dbm_u0 = dbm_u[0];
02562 if (!is_plus_infinity(dbm_u0)) {
02563
02564
02565
02566
02567
02568
02569 assign_r(minus_lb_u, dbm_u0, ROUND_NOT_NEEDED);
02570 assign_r(q, expr_u, ROUND_NOT_NEEDED);
02571 div_assign_r(q, q, mpq_sc_den, ROUND_NOT_NEEDED);
02572 assign_r(ub_u, dbm_0[u], ROUND_NOT_NEEDED);
02573
02574 add_assign_r(ub_u, ub_u, minus_lb_u, ROUND_NOT_NEEDED);
02575
02576 sub_mul_assign_r(minus_lb_u, q, ub_u, ROUND_NOT_NEEDED);
02577 assign_r(up_approx, minus_lb_u, ROUND_UP);
02578
02579 add_assign_r(dbm_u[v], ub_v, up_approx, ROUND_UP);
02580 }
02581 }
02582 }
02583 }
02584 }
02585
02586 template <typename T>
02587 void
02588 BD_Shape<T>
02589 ::deduce_u_minus_v_bounds(const dimension_type v,
02590 const dimension_type last_v,
02591 const Linear_Expression& sc_expr,
02592 Coefficient_traits::const_reference sc_den,
02593 const N& minus_lb_v) {
02594 assert(sc_den > 0);
02595 assert(!is_plus_infinity(minus_lb_v));
02596
02597
02598
02599
02600
02601
02602
02603
02604 DIRTY_TEMP0(mpq_class, mpq_sc_den);
02605 assign_r(mpq_sc_den, sc_den, ROUND_NOT_NEEDED);
02606 DB_Row<N>& dbm_0 = dbm[0];
02607 DB_Row<N>& dbm_v = dbm[v];
02608
02609 DIRTY_TEMP0(mpq_class, ub_u);
02610 DIRTY_TEMP0(mpq_class, q);
02611 DIRTY_TEMP0(mpq_class, minus_lb_u);
02612 DIRTY_TEMP(N, up_approx);
02613
02614 for (dimension_type u = last_v; u > 0; --u)
02615 if (u != v) {
02616 const Coefficient& expr_u = sc_expr.coefficient(Variable(u-1));
02617 if (expr_u > 0) {
02618 if (expr_u >= sc_den)
02619
02620
02621 sub_assign_r(dbm_v[u], minus_lb_v, dbm[u][0], ROUND_UP);
02622 else {
02623 const N& dbm_0u = dbm_0[u];
02624 if (!is_plus_infinity(dbm_0u)) {
02625
02626
02627
02628
02629
02630
02631 assign_r(ub_u, dbm_0u, ROUND_NOT_NEEDED);
02632 assign_r(q, expr_u, ROUND_NOT_NEEDED);
02633 div_assign_r(q, q, mpq_sc_den, ROUND_NOT_NEEDED);
02634 assign_r(minus_lb_u, dbm[u][0], ROUND_NOT_NEEDED);
02635
02636 add_assign_r(minus_lb_u, minus_lb_u, ub_u, ROUND_NOT_NEEDED);
02637
02638 sub_mul_assign_r(ub_u, q, minus_lb_u, ROUND_NOT_NEEDED);
02639 assign_r(up_approx, ub_u, ROUND_UP);
02640
02641 add_assign_r(dbm_v[u], up_approx, minus_lb_v, ROUND_UP);
02642 }
02643 }
02644 }
02645 }
02646 }
02647
02648 template <typename T>
02649 void
02650 BD_Shape<T>::forget_all_dbm_constraints(const dimension_type v) {
02651 assert(0 < v && v <= dbm.num_rows());
02652 DB_Row<N>& dbm_v = dbm[v];
02653 for (dimension_type i = dbm.num_rows(); i-- > 0; ) {
02654 assign_r(dbm_v[i], PLUS_INFINITY, ROUND_NOT_NEEDED);
02655 assign_r(dbm[i][v], PLUS_INFINITY, ROUND_NOT_NEEDED);
02656 }
02657 }
02658
02659 template <typename T>
02660 void
02661 BD_Shape<T>::forget_binary_dbm_constraints(const dimension_type v) {
02662 assert(0 < v && v <= dbm.num_rows());
02663 DB_Row<N>& dbm_v = dbm[v];
02664 for (dimension_type i = dbm.num_rows()-1; i > 0; --i) {
02665 assign_r(dbm_v[i], PLUS_INFINITY, ROUND_NOT_NEEDED);
02666 assign_r(dbm[i][v], PLUS_INFINITY, ROUND_NOT_NEEDED);
02667 }
02668 }
02669
02670 template <typename T>
02671 void
02672 BD_Shape<T>::unconstrain(const Variable var) {
02673
02674 const dimension_type dim = var.id();
02675 if (space_dimension() < dim)
02676 throw_dimension_incompatible("unconstrain(var)", dim);
02677
02678
02679
02680 shortest_path_closure_assign();
02681
02682
02683 if (marked_empty())
02684 return;
02685
02686 forget_all_dbm_constraints(dim+1);
02687
02688 reset_shortest_path_reduced();
02689 assert(OK());
02690 }
02691
02692 template <typename T>
02693 void
02694 BD_Shape<T>::unconstrain(const Variables_Set& to_be_unconstrained) {
02695
02696
02697 if (to_be_unconstrained.empty())
02698 return;
02699
02700
02701 const dimension_type min_space_dim = to_be_unconstrained.space_dimension();
02702 if (space_dimension() < min_space_dim)
02703 throw_dimension_incompatible("unconstrain(vs)", min_space_dim);
02704
02705
02706
02707 shortest_path_closure_assign();
02708
02709
02710 if (marked_empty())
02711 return;
02712
02713 for (Variables_Set::const_iterator tbu = to_be_unconstrained.begin(),
02714 tbu_end = to_be_unconstrained.end(); tbu != tbu_end; ++tbu)
02715 forget_all_dbm_constraints(*tbu + 1);
02716
02717 reset_shortest_path_reduced();
02718 assert(OK());
02719 }
02720
02721 template <typename T>
02722 void
02723 BD_Shape<T>::refine(const Variable var,
02724 const Relation_Symbol relsym,
02725 const Linear_Expression& expr,
02726 Coefficient_traits::const_reference denominator) {
02727 assert(denominator != 0);
02728 const dimension_type expr_space_dim = expr.space_dimension();
02729 assert(space_dimension() >= expr_space_dim);
02730 const dimension_type v = var.id() + 1;
02731 assert(v <= space_dimension());
02732 assert(expr.coefficient(var) == 0);
02733 assert(relsym != LESS_THAN && relsym != GREATER_THAN);
02734
02735 const Coefficient& b = expr.inhomogeneous_term();
02736
02737
02738 dimension_type t = 0;
02739
02740 dimension_type w = 0;
02741
02742 for (dimension_type i = expr_space_dim; i-- > 0; )
02743 if (expr.coefficient(Variable(i)) != 0) {
02744 if (t++ == 1)
02745 break;
02746 else
02747 w = i+1;
02748 }
02749
02750
02751
02752
02753
02754
02755 if (t == 1 && expr.coefficient(Variable(w-1)) != denominator)
02756 t = 2;
02757
02758
02759
02760
02761
02762 const DB_Row<N>& dbm_0 = dbm[0];
02763 TEMP_INTEGER(minus_den);
02764 neg_assign(minus_den, denominator);
02765
02766 if (t == 0) {
02767
02768 switch (relsym) {
02769 case EQUAL:
02770
02771 add_dbm_constraint(0, v, b, denominator);
02772 add_dbm_constraint(v, 0, b, minus_den);
02773 break;
02774 case LESS_OR_EQUAL:
02775
02776 add_dbm_constraint(0, v, b, denominator);
02777 break;
02778 case GREATER_OR_EQUAL:
02779
02780
02781 add_dbm_constraint(v, 0, b, minus_den);
02782 break;
02783 default:
02784
02785 throw std::runtime_error("PPL internal error");
02786 }
02787 return;
02788 }
02789
02790 if (t == 1) {
02791
02792 assert(expr.coefficient(Variable(w-1)) == denominator);
02793 DIRTY_TEMP(N, d);
02794 switch (relsym) {
02795 case EQUAL:
02796
02797 div_round_up(d, b, denominator);
02798 add_dbm_constraint(w, v, d);
02799
02800
02801 div_round_up(d, b, minus_den);
02802 add_dbm_constraint(v, w, d);
02803 break;
02804 case LESS_OR_EQUAL:
02805
02806 div_round_up(d, b, denominator);
02807 add_dbm_constraint(w, v, d);
02808 break;
02809 case GREATER_OR_EQUAL:
02810
02811
02812 div_round_up(d, b, minus_den);
02813 add_dbm_constraint(v, w, d);
02814 break;
02815 default:
02816
02817 throw std::runtime_error("PPL internal error");
02818 }
02819 return;
02820 }
02821
02822
02823
02824
02825 const bool is_sc = (denominator > 0);
02826 TEMP_INTEGER(minus_b);
02827 neg_assign(minus_b, b);
02828 const Coefficient& sc_b = is_sc ? b : minus_b;
02829 const Coefficient& minus_sc_b = is_sc ? minus_b : b;
02830 const Coefficient& sc_den = is_sc ? denominator : minus_den;
02831 const Coefficient& minus_sc_den = is_sc ? minus_den : denominator;
02832
02833
02834
02835 Linear_Expression minus_expr;
02836 if (!is_sc)
02837 minus_expr = -expr;
02838 const Linear_Expression& sc_expr = is_sc ? expr : minus_expr;
02839
02840 DIRTY_TEMP(N, sum);
02841
02842 PPL_UNINITIALIZED(dimension_type, pinf_index);
02843
02844 dimension_type pinf_count = 0;
02845
02846
02847
02848 TEMP_INTEGER(minus_sc_i);
02849 DIRTY_TEMP(N, coeff_i);
02850
02851 switch (relsym) {
02852 case EQUAL:
02853 {
02854 DIRTY_TEMP(N, neg_sum);
02855
02856 PPL_UNINITIALIZED(dimension_type, neg_pinf_index);
02857
02858 dimension_type neg_pinf_count = 0;
02859
02860
02861
02862
02863
02864 assign_r(sum, sc_b, ROUND_UP);
02865 assign_r(neg_sum, minus_sc_b, ROUND_UP);
02866
02867
02868
02869
02870 for (dimension_type i = w; i > 0; --i) {
02871 const Coefficient& sc_i = sc_expr.coefficient(Variable(i-1));
02872 const int sign_i = sgn(sc_i);
02873 if (sign_i == 0)
02874 continue;
02875 if (sign_i > 0) {
02876 assign_r(coeff_i, sc_i, ROUND_UP);
02877
02878 if (pinf_count <= 1) {
02879 const N& approx_i = dbm_0[i];
02880 if (!is_plus_infinity(approx_i))
02881 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
02882 else {
02883 ++pinf_count;
02884 pinf_index = i;
02885 }
02886 }
02887
02888 if (neg_pinf_count <= 1) {
02889 const N& approx_minus_i = dbm[i][0];
02890 if (!is_plus_infinity(approx_minus_i))
02891 add_mul_assign_r(neg_sum, coeff_i, approx_minus_i, ROUND_UP);
02892 else {
02893 ++neg_pinf_count;
02894 neg_pinf_index = i;
02895 }
02896 }
02897 }
02898 else if (sign_i < 0) {
02899 neg_assign(minus_sc_i, sc_i);
02900
02901 assign_r(coeff_i, minus_sc_i, ROUND_UP);
02902
02903 if (pinf_count <= 1) {
02904 const N& approx_minus_i = dbm[i][0];
02905 if (!is_plus_infinity(approx_minus_i))
02906 add_mul_assign_r(sum, coeff_i, approx_minus_i, ROUND_UP);
02907 else {
02908 ++pinf_count;
02909 pinf_index = i;
02910 }
02911 }
02912
02913 if (neg_pinf_count <= 1) {
02914 const N& approx_i = dbm_0[i];
02915 if (!is_plus_infinity(approx_i))
02916 add_mul_assign_r(neg_sum, coeff_i, approx_i, ROUND_UP);
02917 else {
02918 ++neg_pinf_count;
02919 neg_pinf_index = i;
02920 }
02921 }
02922 }
02923 }
02924
02925 if (pinf_count > 1 && neg_pinf_count > 1) {
02926 assert(OK());
02927 return;
02928 }
02929
02930
02931 reset_shortest_path_closed();
02932
02933
02934
02935
02936
02937 DIRTY_TEMP(N, down_sc_den);
02938 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
02939 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
02940
02941
02942 if (pinf_count <= 1) {
02943
02944 if (down_sc_den != 1)
02945 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
02946
02947 if (pinf_count == 0) {
02948
02949 dbm[0][v] = sum;
02950
02951 deduce_v_minus_u_bounds(v, w, sc_expr, sc_den, sum);
02952 }
02953 else
02954
02955 if (pinf_index != v
02956 && sc_expr.coefficient(Variable(pinf_index-1)) == sc_den)
02957
02958 dbm[pinf_index][v] = sum;
02959 }
02960
02961
02962 if (neg_pinf_count <= 1) {
02963
02964 if (down_sc_den != 1)
02965 div_assign_r(neg_sum, neg_sum, down_sc_den, ROUND_UP);
02966
02967 if (neg_pinf_count == 0) {
02968
02969 DB_Row<N>& dbm_v = dbm[v];
02970 dbm_v[0] = neg_sum;
02971
02972 deduce_u_minus_v_bounds(v, w, sc_expr, sc_den, neg_sum);
02973 }
02974 else
02975
02976 if (neg_pinf_index != v
02977 && sc_expr.coefficient(Variable(neg_pinf_index-1)) == sc_den)
02978
02979
02980 dbm[v][neg_pinf_index] = neg_sum;
02981 }
02982 }
02983 break;
02984
02985 case LESS_OR_EQUAL:
02986
02987
02988
02989
02990 assign_r(sum, sc_b, ROUND_UP);
02991
02992
02993
02994
02995 for (dimension_type i = w; i > 0; --i) {
02996 const Coefficient& sc_i = sc_expr.coefficient(Variable(i-1));
02997 const int sign_i = sgn(sc_i);
02998 if (sign_i == 0)
02999 continue;
03000
03001 const N& approx_i = (sign_i > 0) ? dbm_0[i] : dbm[i][0];
03002 if (is_plus_infinity(approx_i)) {
03003 if (++pinf_count > 1)
03004 break;
03005 pinf_index = i;
03006 continue;
03007 }
03008 if (sign_i > 0)
03009 assign_r(coeff_i, sc_i, ROUND_UP);
03010 else {
03011 neg_assign(minus_sc_i, sc_i);
03012 assign_r(coeff_i, minus_sc_i, ROUND_UP);
03013 }
03014 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
03015 }
03016
03017
03018 if (sc_den != 1) {
03019
03020
03021
03022
03023 DIRTY_TEMP(N, down_sc_den);
03024 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03025 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03026 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
03027 }
03028
03029 if (pinf_count == 0) {
03030
03031 add_dbm_constraint(0, v, sum);
03032
03033 deduce_v_minus_u_bounds(v, w, sc_expr, sc_den, sum);
03034 }
03035 else if (pinf_count == 1)
03036 if (expr.coefficient(Variable(pinf_index-1)) == denominator)
03037
03038 add_dbm_constraint(pinf_index, v, sum);
03039 break;
03040
03041 case GREATER_OR_EQUAL:
03042
03043
03044
03045
03046
03047 assign_r(sum, minus_sc_b, ROUND_UP);
03048
03049
03050 for (dimension_type i = w; i > 0; --i) {
03051 const Coefficient& sc_i = sc_expr.coefficient(Variable(i-1));
03052 const int sign_i = sgn(sc_i);
03053 if (sign_i == 0)
03054 continue;
03055
03056 const N& approx_i = (sign_i > 0) ? dbm[i][0] : dbm_0[i];
03057 if (is_plus_infinity(approx_i)) {
03058 if (++pinf_count > 1)
03059 break;
03060 pinf_index = i;
03061 continue;
03062 }
03063 if (sign_i > 0)
03064 assign_r(coeff_i, sc_i, ROUND_UP);
03065 else {
03066 neg_assign(minus_sc_i, sc_i);
03067 assign_r(coeff_i, minus_sc_i, ROUND_UP);
03068 }
03069 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
03070 }
03071
03072
03073 if (sc_den != 1) {
03074
03075
03076
03077
03078 DIRTY_TEMP(N, down_sc_den);
03079 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03080 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03081 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
03082 }
03083
03084 if (pinf_count == 0) {
03085
03086 add_dbm_constraint(v, 0, sum);
03087
03088 deduce_u_minus_v_bounds(v, w, sc_expr, sc_den, sum);
03089 }
03090 else if (pinf_count == 1)
03091 if (pinf_index != v
03092 && expr.coefficient(Variable(pinf_index-1)) == denominator)
03093
03094
03095 add_dbm_constraint(v, pinf_index, sum);
03096 break;
03097
03098 default:
03099
03100 throw std::runtime_error("PPL internal error");
03101 }
03102
03103 assert(OK());
03104 }
03105
03106 template <typename T>
03107 void
03108 BD_Shape<T>::affine_image(const Variable var,
03109 const Linear_Expression& expr,
03110 Coefficient_traits::const_reference denominator) {
03111
03112 if (denominator == 0)
03113 throw_generic("affine_image(v, e, d)", "d == 0");
03114
03115
03116
03117
03118 const dimension_type space_dim = space_dimension();
03119 const dimension_type expr_space_dim = expr.space_dimension();
03120 if (space_dim < expr_space_dim)
03121 throw_dimension_incompatible("affine_image(v, e, d)", "e", expr);
03122
03123
03124 const dimension_type v = var.id() + 1;
03125 if (v > space_dim)
03126 throw_dimension_incompatible("affine_image(v, e, d)", var.id());
03127
03128
03129 shortest_path_closure_assign();
03130 if (marked_empty())
03131 return;
03132
03133 const Coefficient& b = expr.inhomogeneous_term();
03134
03135
03136 dimension_type t = 0;
03137
03138 dimension_type w = 0;
03139
03140 for (dimension_type i = expr_space_dim; i-- > 0; )
03141 if (expr.coefficient(Variable(i)) != 0) {
03142 if (t++ == 1)
03143 break;
03144 else
03145 w = i+1;
03146 }
03147
03148
03149
03150
03151
03152
03153
03154
03155 TEMP_INTEGER(minus_den);
03156 neg_assign(minus_den, denominator);
03157
03158 if (t == 0) {
03159
03160
03161 forget_all_dbm_constraints(v);
03162
03163 if (marked_shortest_path_reduced())
03164 reset_shortest_path_reduced();
03165
03166 add_dbm_constraint(0, v, b, denominator);
03167 add_dbm_constraint(v, 0, b, minus_den);
03168 assert(OK());
03169 return;
03170 }
03171
03172 if (t == 1) {
03173
03174 const Coefficient& a = expr.coefficient(Variable(w-1));
03175 if (a == denominator || a == minus_den) {
03176
03177 if (w == v) {
03178
03179 if (a == denominator) {
03180 if (b == 0)
03181
03182 return;
03183 else {
03184
03185
03186 DIRTY_TEMP(N, d);
03187 div_round_up(d, b, denominator);
03188 DIRTY_TEMP(N, c);
03189 div_round_up(c, b, minus_den);
03190 DB_Row<N>& dbm_v = dbm[v];
03191 for (dimension_type i = space_dim + 1; i-- > 0; ) {
03192 N& dbm_vi = dbm_v[i];
03193 add_assign_r(dbm_vi, dbm_vi, c, ROUND_UP);
03194 N& dbm_iv = dbm[i][v];
03195 add_assign_r(dbm_iv, dbm_iv, d, ROUND_UP);
03196 }
03197
03198 }
03199 }
03200 else {
03201
03202
03203 forget_binary_dbm_constraints(v);
03204
03205 std::swap(dbm[v][0], dbm[0][v]);
03206
03207 reset_shortest_path_closed();
03208 if (b != 0) {
03209
03210
03211 DIRTY_TEMP(N, c);
03212 div_round_up(c, b, minus_den);
03213 N& dbm_v0 = dbm[v][0];
03214 add_assign_r(dbm_v0, dbm_v0, c, ROUND_UP);
03215 DIRTY_TEMP(N, d);
03216 div_round_up(d, b, denominator);
03217 N& dbm_0v = dbm[0][v];
03218 add_assign_r(dbm_0v, dbm_0v, d, ROUND_UP);
03219 }
03220 }
03221 }
03222 else {
03223
03224
03225
03226 forget_all_dbm_constraints(v);
03227
03228 if (marked_shortest_path_reduced())
03229 reset_shortest_path_reduced();
03230 if (a == denominator) {
03231
03232 add_dbm_constraint(w, v, b, denominator);
03233 add_dbm_constraint(v, w, b, minus_den);
03234 }
03235 else {
03236
03237
03238
03239 const N& dbm_w0 = dbm[w][0];
03240 if (!is_plus_infinity(dbm_w0)) {
03241
03242 DIRTY_TEMP(N, d);
03243 div_round_up(d, b, denominator);
03244 add_assign_r(dbm[0][v], d, dbm_w0, ROUND_UP);
03245 reset_shortest_path_closed();
03246 }
03247 const N& dbm_0w = dbm[0][w];
03248 if (!is_plus_infinity(dbm_0w)) {
03249
03250 DIRTY_TEMP(N, c);
03251 div_round_up(c, b, minus_den);
03252 add_assign_r(dbm[v][0], dbm_0w, c, ROUND_UP);
03253 reset_shortest_path_closed();
03254 }
03255 }
03256 }
03257 assert(OK());
03258 return;
03259 }
03260 }
03261
03262
03263
03264
03265
03266
03267
03268
03269
03270
03271
03272
03273
03274 const bool is_sc = (denominator > 0);
03275 TEMP_INTEGER(minus_b);
03276 neg_assign(minus_b, b);
03277 const Coefficient& sc_b = is_sc ? b : minus_b;
03278 const Coefficient& minus_sc_b = is_sc ? minus_b : b;
03279 const Coefficient& sc_den = is_sc ? denominator : minus_den;
03280 const Coefficient& minus_sc_den = is_sc ? minus_den : denominator;
03281
03282
03283
03284 Linear_Expression minus_expr;
03285 if (!is_sc)
03286 minus_expr = -expr;
03287 const Linear_Expression& sc_expr = is_sc ? expr : minus_expr;
03288
03289 DIRTY_TEMP(N, pos_sum);
03290 DIRTY_TEMP(N, neg_sum);
03291
03292 PPL_UNINITIALIZED(dimension_type, pos_pinf_index);
03293 PPL_UNINITIALIZED(dimension_type, neg_pinf_index);
03294
03295 dimension_type pos_pinf_count = 0;
03296 dimension_type neg_pinf_count = 0;
03297
03298
03299 assign_r(pos_sum, sc_b, ROUND_UP);
03300 assign_r(neg_sum, minus_sc_b, ROUND_UP);
03301
03302
03303 const DB_Row<N>& dbm_0 = dbm[0];
03304
03305 DIRTY_TEMP(N, coeff_i);
03306 TEMP_INTEGER(minus_sc_i);
03307
03308
03309 for (dimension_type i = w; i > 0; --i) {
03310 const Coefficient& sc_i = sc_expr.coefficient(Variable(i-1));
03311 const int sign_i = sgn(sc_i);
03312 if (sign_i > 0) {
03313 assign_r(coeff_i, sc_i, ROUND_UP);
03314
03315 if (pos_pinf_count <= 1) {
03316 const N& up_approx_i = dbm_0[i];
03317 if (!is_plus_infinity(up_approx_i))
03318 add_mul_assign_r(pos_sum, coeff_i, up_approx_i, ROUND_UP);
03319 else {
03320 ++pos_pinf_count;
03321 pos_pinf_index = i;
03322 }
03323 }
03324
03325 if (neg_pinf_count <= 1) {
03326 const N& up_approx_minus_i = dbm[i][0];
03327 if (!is_plus_infinity(up_approx_minus_i))
03328 add_mul_assign_r(neg_sum, coeff_i, up_approx_minus_i, ROUND_UP);
03329 else {
03330 ++neg_pinf_count;
03331 neg_pinf_index = i;
03332 }
03333 }
03334 }
03335 else if (sign_i < 0) {
03336 neg_assign(minus_sc_i, sc_i);
03337
03338 assign_r(coeff_i, minus_sc_i, ROUND_UP);
03339
03340 if (pos_pinf_count <= 1) {
03341 const N& up_approx_minus_i = dbm[i][0];
03342 if (!is_plus_infinity(up_approx_minus_i))
03343 add_mul_assign_r(pos_sum, coeff_i, up_approx_minus_i, ROUND_UP);
03344 else {
03345 ++pos_pinf_count;
03346 pos_pinf_index = i;
03347 }
03348 }
03349
03350 if (neg_pinf_count <= 1) {
03351 const N& up_approx_i = dbm_0[i];
03352 if (!is_plus_infinity(up_approx_i))
03353 add_mul_assign_r(neg_sum, coeff_i, up_approx_i, ROUND_UP);
03354 else {
03355 ++neg_pinf_count;
03356 neg_pinf_index = i;
03357 }
03358 }
03359 }
03360 }
03361
03362
03363 forget_all_dbm_constraints(v);
03364
03365 if (marked_shortest_path_reduced())
03366 reset_shortest_path_reduced();
03367
03368 if (pos_pinf_count > 1 && neg_pinf_count > 1) {
03369 assert(OK());
03370 return;
03371 }
03372
03373
03374 reset_shortest_path_closed();
03375
03376
03377 if (pos_pinf_count <= 1) {
03378
03379 if (sc_den != 1) {
03380
03381
03382
03383
03384 DIRTY_TEMP(N, down_sc_den);
03385 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03386 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03387 div_assign_r(pos_sum, pos_sum, down_sc_den, ROUND_UP);
03388 }
03389
03390 if (pos_pinf_count == 0) {
03391
03392 dbm[0][v] = pos_sum;
03393
03394 deduce_v_minus_u_bounds(v, w, sc_expr, sc_den, pos_sum);
03395 }
03396 else
03397
03398 if (pos_pinf_index != v
03399 && sc_expr.coefficient(Variable(pos_pinf_index-1)) == sc_den)
03400
03401 dbm[pos_pinf_index][v] = pos_sum;
03402 }
03403
03404
03405 if (neg_pinf_count <= 1) {
03406
03407 if (sc_den != 1) {
03408
03409
03410
03411
03412 DIRTY_TEMP(N, down_sc_den);
03413 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03414 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03415 div_assign_r(neg_sum, neg_sum, down_sc_den, ROUND_UP);
03416 }
03417
03418 if (neg_pinf_count == 0) {
03419
03420 DB_Row<N>& dbm_v = dbm[v];
03421 dbm_v[0] = neg_sum;
03422
03423 deduce_u_minus_v_bounds(v, w, sc_expr, sc_den, neg_sum);
03424 }
03425 else
03426
03427 if (neg_pinf_index != v
03428 && sc_expr.coefficient(Variable(neg_pinf_index-1)) == sc_den)
03429
03430
03431 dbm[v][neg_pinf_index] = neg_sum;
03432 }
03433
03434 assert(OK());
03435 }
03436
03437 template <typename T>
03438 void
03439 BD_Shape<T>::affine_preimage(const Variable var,
03440 const Linear_Expression& expr,
03441 Coefficient_traits::const_reference denominator) {
03442
03443 if (denominator == 0)
03444 throw_generic("affine_preimage(v, e, d)", "d == 0");
03445
03446
03447
03448
03449 const dimension_type space_dim = space_dimension();
03450 const dimension_type expr_space_dim = expr.space_dimension();
03451 if (space_dim < expr_space_dim)
03452 throw_dimension_incompatible("affine_preimage(v, e, d)", "e", expr);
03453
03454
03455
03456 const dimension_type v = var.id() + 1;
03457 if (v > space_dim)
03458 throw_dimension_incompatible("affine_preimage(v, e, d)", var.id());
03459
03460
03461 shortest_path_closure_assign();
03462 if (marked_empty())
03463 return;
03464
03465 const Coefficient& b = expr.inhomogeneous_term();
03466
03467
03468 dimension_type t = 0;
03469
03470 dimension_type j = 0;
03471
03472 for (dimension_type i = expr_space_dim; i-- > 0; )
03473 if (expr.coefficient(Variable(i)) != 0) {
03474 if (t++ == 1)
03475 break;
03476 else
03477 j = i;
03478 }
03479
03480
03481
03482
03483
03484
03485
03486
03487 if (t == 0) {
03488
03489 forget_all_dbm_constraints(v);
03490
03491 if (marked_shortest_path_reduced())
03492 reset_shortest_path_reduced();
03493 assert(OK());
03494 return;
03495 }
03496
03497 if (t == 1) {
03498
03499 const Coefficient& a = expr.coefficient(Variable(j));
03500 if (a == denominator || a == -denominator) {
03501
03502 if (j == var.id())
03503
03504 affine_image(var, denominator*var - b, a);
03505 else {
03506
03507
03508 forget_all_dbm_constraints(v);
03509
03510 if (marked_shortest_path_reduced())
03511 reset_shortest_path_reduced();
03512 assert(OK());
03513 }
03514 return;
03515 }
03516 }
03517
03518
03519
03520
03521
03522 const Coefficient& expr_v = expr.coefficient(var);
03523 if (expr_v != 0) {
03524
03525 Linear_Expression inverse((expr_v + denominator)*var);
03526 inverse -= expr;
03527 affine_image(var, inverse, expr_v);
03528 }
03529 else {
03530
03531 forget_all_dbm_constraints(v);
03532
03533 if (marked_shortest_path_reduced())
03534 reset_shortest_path_reduced();
03535 }
03536 assert(OK());
03537 }
03538
03539 template <typename T>
03540 void
03541 BD_Shape<T>
03542 ::bounded_affine_image(const Variable var,
03543 const Linear_Expression& lb_expr,
03544 const Linear_Expression& ub_expr,
03545 Coefficient_traits::const_reference denominator) {
03546
03547 if (denominator == 0)
03548 throw_generic("bounded_affine_image(v, lb, ub, d)", "d == 0");
03549
03550
03551
03552 const dimension_type bds_space_dim = space_dimension();
03553 const dimension_type v = var.id() + 1;
03554 if (v > bds_space_dim)
03555 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
03556 "v", var);
03557
03558
03559 const dimension_type lb_space_dim = lb_expr.space_dimension();
03560 if (bds_space_dim < lb_space_dim)
03561 throw_dimension_incompatible("bounded_affine_image(v, lb, ub)",
03562 "lb", lb_expr);
03563 const dimension_type ub_space_dim = ub_expr.space_dimension();
03564 if (bds_space_dim < ub_space_dim)
03565 throw_dimension_incompatible("bounded_affine_image(v, lb, ub)",
03566 "ub", ub_expr);
03567
03568
03569 shortest_path_closure_assign();
03570 if (marked_empty())
03571 return;
03572
03573 const Coefficient& b = ub_expr.inhomogeneous_term();
03574
03575
03576 dimension_type t = 0;
03577
03578 dimension_type w = 0;
03579
03580 for (dimension_type i = ub_space_dim; i-- > 0; )
03581 if (ub_expr.coefficient(Variable(i)) != 0) {
03582 if (t++ == 1)
03583 break;
03584 else
03585 w = i+1;
03586 }
03587
03588
03589
03590
03591
03592
03593
03594
03595 TEMP_INTEGER(minus_den);
03596 neg_assign(minus_den, denominator);
03597
03598 if (t == 0) {
03599
03600 generalized_affine_image(var,
03601 GREATER_OR_EQUAL,
03602 lb_expr,
03603 denominator);
03604
03605 add_dbm_constraint(0, v, b, denominator);
03606 assert(OK());
03607 return;
03608 }
03609
03610 if (t == 1) {
03611
03612 const Coefficient& a = ub_expr.coefficient(Variable(w-1));
03613 if (a == denominator || a == minus_den) {
03614
03615 if (w == v) {
03616
03617
03618 const Variable new_var = Variable(bds_space_dim);
03619 add_space_dimensions_and_embed(1);
03620
03621 affine_image(new_var, ub_expr, denominator);
03622
03623 shortest_path_closure_assign();
03624 assert(!marked_empty());
03625
03626 generalized_affine_image(var,
03627 GREATER_OR_EQUAL,
03628 lb_expr,
03629 denominator);
03630
03631 add_constraint(var <= new_var);
03632
03633 remove_higher_space_dimensions(bds_space_dim);
03634 return;
03635 }
03636 else {
03637
03638
03639
03640 generalized_affine_image(var,
03641 GREATER_OR_EQUAL,
03642 lb_expr,
03643 denominator);
03644 if (a == denominator) {
03645
03646 add_dbm_constraint(w, v, b, denominator);
03647 }
03648 else {
03649
03650
03651
03652 const N& dbm_w0 = dbm[w][0];
03653 if (!is_plus_infinity(dbm_w0)) {
03654
03655 DIRTY_TEMP(N, d);
03656 div_round_up(d, b, denominator);
03657 add_assign_r(dbm[0][v], d, dbm_w0, ROUND_UP);
03658 reset_shortest_path_closed();
03659 }
03660 }
03661 assert(OK());
03662 return;
03663 }
03664 }
03665 }
03666
03667
03668
03669
03670
03671
03672
03673
03674
03675
03676 const bool is_sc = (denominator > 0);
03677 TEMP_INTEGER(minus_b);
03678 neg_assign(minus_b, b);
03679 const Coefficient& sc_b = is_sc ? b : minus_b;
03680 const Coefficient& sc_den = is_sc ? denominator : minus_den;
03681 const Coefficient& minus_sc_den = is_sc ? minus_den : denominator;
03682
03683
03684
03685 Linear_Expression minus_expr;
03686 if (!is_sc)
03687 minus_expr = -ub_expr;
03688 const Linear_Expression& sc_expr = is_sc ? ub_expr : minus_expr;
03689
03690 DIRTY_TEMP(N, pos_sum);
03691
03692 PPL_UNINITIALIZED(dimension_type, pos_pinf_index);
03693
03694 dimension_type pos_pinf_count = 0;
03695
03696
03697 assign_r(pos_sum, sc_b, ROUND_UP);
03698
03699
03700 const DB_Row<N>& dbm_0 = dbm[0];
03701
03702 DIRTY_TEMP(N, coeff_i);
03703 TEMP_INTEGER(minus_sc_i);
03704
03705
03706 for (dimension_type i = w; i > 0; --i) {
03707 const Coefficient& sc_i = sc_expr.coefficient(Variable(i-1));
03708 const int sign_i = sgn(sc_i);
03709 if (sign_i > 0) {
03710 assign_r(coeff_i, sc_i, ROUND_UP);
03711
03712 if (pos_pinf_count <= 1) {
03713 const N& up_approx_i = dbm_0[i];
03714 if (!is_plus_infinity(up_approx_i))
03715 add_mul_assign_r(pos_sum, coeff_i, up_approx_i, ROUND_UP);
03716 else {
03717 ++pos_pinf_count;
03718 pos_pinf_index = i;
03719 }
03720 }
03721 }
03722 else if (sign_i < 0) {
03723 neg_assign(minus_sc_i, sc_i);
03724
03725 assign_r(coeff_i, minus_sc_i, ROUND_UP);
03726
03727 if (pos_pinf_count <= 1) {
03728 const N& up_approx_minus_i = dbm[i][0];
03729 if (!is_plus_infinity(up_approx_minus_i))
03730 add_mul_assign_r(pos_sum, coeff_i, up_approx_minus_i, ROUND_UP);
03731 else {
03732 ++pos_pinf_count;
03733 pos_pinf_index = i;
03734 }
03735 }
03736 }
03737 }
03738
03739 generalized_affine_image(var,
03740 GREATER_OR_EQUAL,
03741 lb_expr,
03742 denominator);
03743
03744 if (pos_pinf_count > 1) {
03745 return;
03746 }
03747
03748
03749 reset_shortest_path_closed();
03750
03751
03752 if (pos_pinf_count <= 1) {
03753
03754 if (sc_den != 1) {
03755
03756
03757
03758
03759 DIRTY_TEMP(N, down_sc_den);
03760 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03761 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03762 div_assign_r(pos_sum, pos_sum, down_sc_den, ROUND_UP);
03763 }
03764
03765 if (pos_pinf_count == 0) {
03766
03767 dbm[0][v] = pos_sum;
03768
03769 deduce_v_minus_u_bounds(v, w, sc_expr, sc_den, pos_sum);
03770 }
03771 else
03772
03773 if (pos_pinf_index != v
03774 && sc_expr.coefficient(Variable(pos_pinf_index-1)) == sc_den)
03775
03776 dbm[pos_pinf_index][v] = pos_sum;
03777 }
03778 assert(OK());
03779 }
03780
03781 template <typename T>
03782 void
03783 BD_Shape<T>
03784 ::bounded_affine_preimage(const Variable var,
03785 const Linear_Expression& lb_expr,
03786 const Linear_Expression& ub_expr,
03787 Coefficient_traits::const_reference denominator) {
03788
03789 if (denominator == 0)
03790 throw_generic("bounded_affine_preimage(v, lb, ub, d)", "d == 0");
03791
03792
03793
03794 const dimension_type space_dim = space_dimension();
03795 const dimension_type v = var.id() + 1;
03796 if (v > space_dim)
03797 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
03798 "v", var);
03799
03800
03801 const dimension_type lb_space_dim = lb_expr.space_dimension();
03802 if (space_dim < lb_space_dim)
03803 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub)",
03804 "lb", lb_expr);
03805 const dimension_type ub_space_dim = ub_expr.space_dimension();
03806 if (space_dim < ub_space_dim)
03807 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub)",
03808 "ub", ub_expr);
03809
03810
03811 shortest_path_closure_assign();
03812 if (marked_empty())
03813 return;
03814
03815 if (ub_expr.coefficient(var) == 0) {
03816 refine(var, LESS_OR_EQUAL, ub_expr, denominator);
03817 generalized_affine_preimage(var, GREATER_OR_EQUAL,
03818 lb_expr, denominator);
03819 return;
03820 }
03821 if (lb_expr.coefficient(var) == 0) {
03822 refine(var, GREATER_OR_EQUAL, lb_expr, denominator);
03823 generalized_affine_preimage(var, LESS_OR_EQUAL,
03824 ub_expr, denominator);
03825 return;
03826 }
03827
03828 const Coefficient& lb_expr_v = lb_expr.coefficient(var);
03829
03830
03831 const Variable new_var = Variable(space_dim);
03832 add_space_dimensions_and_embed(1);
03833 const Linear_Expression lb_inverse
03834 = lb_expr - (lb_expr_v + denominator)*var;
03835 TEMP_INTEGER(lb_inverse_den);
03836 neg_assign(lb_inverse_den, lb_expr_v);
03837 affine_image(new_var, lb_inverse, lb_inverse_den);
03838 shortest_path_closure_assign();
03839 assert(!marked_empty());
03840 generalized_affine_preimage(var, LESS_OR_EQUAL,
03841 ub_expr, denominator);
03842 if (sgn(denominator) == sgn(lb_inverse_den))
03843 add_constraint(var >= new_var);
03844 else
03845 add_constraint(var <= new_var);
03846
03847 remove_higher_space_dimensions(space_dim);
03848 }
03849
03850 template <typename T>
03851 void
03852 BD_Shape<T>::generalized_affine_image(const Variable var,
03853 const Relation_Symbol relsym,
03854 const Linear_Expression& expr,
03855 Coefficient_traits::const_reference
03856 denominator) {
03857
03858 if (denominator == 0)
03859 throw_generic("generalized_affine_image(v, r, e, d)", "d == 0");
03860
03861
03862
03863
03864 const dimension_type space_dim = space_dimension();
03865 const dimension_type expr_space_dim = expr.space_dimension();
03866 if (space_dim < expr_space_dim)
03867 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
03868 "e", expr);
03869
03870
03871 const dimension_type v = var.id() + 1;
03872 if (v > space_dim)
03873 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
03874 var.id());
03875
03876
03877 if (relsym == LESS_THAN || relsym == GREATER_THAN)
03878 throw_generic("generalized_affine_image(v, r, e, d)",
03879 "r is a strict relation symbol and "
03880 "*this is a BD_Shape");
03881
03882 if (relsym == NOT_EQUAL)
03883 throw_generic("generalized_affine_image(v, r, e, d)",
03884 "r is the disequality relation symbol and "
03885 "*this is a BD_Shape");
03886
03887 if (relsym == EQUAL) {
03888
03889
03890 affine_image(var, expr, denominator);
03891 return;
03892 }
03893
03894
03895 shortest_path_closure_assign();
03896 if (marked_empty())
03897 return;
03898
03899 const Coefficient& b = expr.inhomogeneous_term();
03900
03901
03902 dimension_type t = 0;
03903
03904 dimension_type w = 0;
03905
03906 for (dimension_type i = expr_space_dim; i-- > 0; )
03907 if (expr.coefficient(Variable(i)) != 0) {
03908 if (t++ == 1)
03909 break;
03910 else
03911 w = i+1;
03912 }
03913
03914
03915
03916
03917
03918
03919
03920
03921 DB_Row<N>& dbm_0 = dbm[0];
03922 DB_Row<N>& dbm_v = dbm[v];
03923 TEMP_INTEGER(minus_den);
03924 neg_assign(minus_den, denominator);
03925
03926 if (t == 0) {
03927
03928
03929 forget_all_dbm_constraints(v);
03930
03931 reset_shortest_path_closed();
03932 switch (relsym) {
03933 case LESS_OR_EQUAL:
03934
03935 add_dbm_constraint(0, v, b, denominator);
03936 break;
03937 case GREATER_OR_EQUAL:
03938
03939
03940 add_dbm_constraint(v, 0, b, minus_den);
03941 break;
03942 default:
03943
03944 throw std::runtime_error("PPL internal error");
03945 }
03946 assert(OK());
03947 return;
03948 }
03949
03950 if (t == 1) {
03951
03952 const Coefficient& a = expr.coefficient(Variable(w-1));
03953 if (a == denominator || a == minus_den) {
03954
03955 DIRTY_TEMP(N, d);
03956 switch (relsym) {
03957 case LESS_OR_EQUAL:
03958 div_round_up(d, b, denominator);
03959 if (w == v) {
03960
03961
03962 reset_shortest_path_closed();
03963 if (a == denominator) {
03964
03965
03966
03967 for (dimension_type i = space_dim + 1; i-- > 0; ) {
03968 N& dbm_iv = dbm[i][v];
03969 add_assign_r(dbm_iv, dbm_iv, d, ROUND_UP);
03970 assign_r(dbm_v[i], PLUS_INFINITY, ROUND_NOT_NEEDED);
03971 }
03972 }
03973 else {
03974
03975
03976
03977 N& dbm_v0 = dbm_v[0];
03978 add_assign_r(dbm_0[v], dbm_v0, d, ROUND_UP);
03979
03980 assign_r(dbm_v0, PLUS_INFINITY, ROUND_NOT_NEEDED);
03981 forget_binary_dbm_constraints(v);
03982 }
03983 }
03984 else {
03985
03986
03987
03988 forget_all_dbm_constraints(v);
03989
03990 if (marked_shortest_path_reduced())
03991 reset_shortest_path_reduced();
03992 if (a == denominator)
03993
03994 add_dbm_constraint(w, v, d);
03995 else {
03996
03997
03998
03999 const N& dbm_w0 = dbm[w][0];
04000 if (!is_plus_infinity(dbm_w0)) {
04001
04002 add_assign_r(dbm_0[v], d, dbm_w0, ROUND_UP);
04003
04004 reset_shortest_path_closed();
04005 }
04006 }
04007 }
04008 break;
04009
04010 case GREATER_OR_EQUAL:
04011 div_round_up(d, b, minus_den);
04012 if (w == v) {
04013
04014
04015 reset_shortest_path_closed();
04016 if (a == denominator) {
04017
04018
04019
04020 for (dimension_type i = space_dim + 1; i-- > 0; ) {
04021 N& dbm_vi = dbm_v[i];
04022 add_assign_r(dbm_vi, dbm_vi, d, ROUND_UP);
04023 assign_r(dbm[i][v], PLUS_INFINITY, ROUND_NOT_NEEDED);
04024 }
04025 }
04026 else {
04027
04028
04029
04030 N& dbm_0v = dbm_0[v];
04031 add_assign_r(dbm_v[0], dbm_0v, d, ROUND_UP);
04032
04033 assign_r(dbm_0v, PLUS_INFINITY, ROUND_NOT_NEEDED);
04034 forget_binary_dbm_constraints(v);
04035 }
04036 }
04037 else {
04038
04039
04040
04041 forget_all_dbm_constraints(v);
04042
04043 if (marked_shortest_path_reduced())
04044 reset_shortest_path_reduced();
04045 if (a == denominator)
04046
04047
04048 add_dbm_constraint(v, w, d);
04049 else {
04050
04051
04052
04053
04054 const N& dbm_0w = dbm_0[w];
04055 if (!is_plus_infinity(dbm_0w)) {
04056
04057 add_assign_r(dbm_v[0], dbm_0w, d, ROUND_UP);
04058
04059 reset_shortest_path_closed();
04060 }
04061 }
04062 }
04063 break;
04064
04065 default:
04066
04067 throw std::runtime_error("PPL internal error");
04068 }
04069 assert(OK());
04070 return;
04071 }
04072 }
04073
04074
04075
04076
04077
04078
04079
04080
04081 const bool is_sc = (denominator > 0);
04082 TEMP_INTEGER(minus_b);
04083 neg_assign(minus_b, b);
04084 const Coefficient& sc_b = is_sc ? b : minus_b;
04085 const Coefficient& minus_sc_b = is_sc ? minus_b : b;
04086 const Coefficient& sc_den = is_sc ? denominator : minus_den;
04087 const Coefficient& minus_sc_den = is_sc ? minus_den : denominator;
04088
04089
04090
04091 Linear_Expression minus_expr;
04092 if (!is_sc)
04093 minus_expr = -expr;
04094 const Linear_Expression& sc_expr = is_sc ? expr : minus_expr;
04095
04096 DIRTY_TEMP(N, sum);
04097
04098 PPL_UNINITIALIZED(dimension_type, pinf_index);
04099
04100 dimension_type pinf_count = 0;
04101
04102
04103 DIRTY_TEMP(N, coeff_i);
04104 TEMP_INTEGER(minus_sc_i);
04105
04106 switch (relsym) {
04107 case LESS_OR_EQUAL:
04108
04109
04110
04111 assign_r(sum, sc_b, ROUND_UP);
04112
04113
04114
04115 for (dimension_type i = w; i > 0; --i) {
04116 const Coefficient& sc_i = sc_expr.coefficient(Variable(i-1));
04117 const int sign_i = sgn(sc_i);
04118 if (sign_i == 0)
04119 continue;
04120
04121 const N& approx_i = (sign_i > 0) ? dbm_0[i] : dbm[i][0];
04122 if (is_plus_infinity(approx_i)) {
04123 if (++pinf_count > 1)
04124 break;
04125 pinf_index = i;
04126 continue;
04127 }
04128 if (sign_i > 0)
04129 assign_r(coeff_i, sc_i, ROUND_UP);
04130 else {
04131 neg_assign(minus_sc_i, sc_i);
04132 assign_r(coeff_i, minus_sc_i, ROUND_UP);
04133 }
04134 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
04135 }
04136
04137
04138 forget_all_dbm_constraints(v);
04139
04140 if (marked_shortest_path_reduced())
04141 reset_shortest_path_reduced();
04142
04143 if (pinf_count > 1) {
04144 assert(OK());
04145 return;
04146 }
04147
04148
04149 if (sc_den != 1) {
04150
04151
04152
04153
04154 DIRTY_TEMP(N, down_sc_den);
04155 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
04156 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
04157 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
04158 }
04159
04160 if (pinf_count == 0) {
04161
04162 add_dbm_constraint(0, v, sum);
04163
04164 deduce_v_minus_u_bounds(v, w, sc_expr, sc_den, sum);
04165 }
04166 else if (pinf_count == 1)
04167 if (pinf_index != v
04168 && expr.coefficient(Variable(pinf_index-1)) == denominator)
04169
04170 add_dbm_constraint(pinf_index, v, sum);
04171 break;
04172
04173 case GREATER_OR_EQUAL:
04174
04175
04176
04177
04178
04179 assign_r(sum, minus_sc_b, ROUND_UP);
04180
04181 for (dimension_type i = expr_space_dim + 1; i > 0; --i) {
04182 const Coefficient& sc_i = sc_expr.coefficient(Variable(i-1));
04183 const int sign_i = sgn(sc_i);
04184 if (sign_i == 0)
04185 continue;
04186
04187 const N& approx_i = (sign_i > 0) ? dbm[i][0] : dbm_0[i];
04188 if (is_plus_infinity(approx_i)) {
04189 if (++pinf_count > 1)
04190 break;
04191 pinf_index = i;
04192 continue;
04193 }
04194 if (sign_i > 0)
04195 assign_r(coeff_i, sc_i, ROUND_UP);
04196 else {
04197 neg_assign(minus_sc_i, sc_i);
04198 assign_r(coeff_i, minus_sc_i, ROUND_UP);
04199 }
04200 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
04201 }
04202
04203
04204 forget_all_dbm_constraints(v);
04205
04206 if (marked_shortest_path_reduced())
04207 reset_shortest_path_reduced();
04208
04209 if (pinf_count > 1) {
04210 assert(OK());
04211 return;
04212 }
04213
04214
04215 if (sc_den != 1) {
04216
04217
04218
04219
04220 DIRTY_TEMP(N, down_sc_den);
04221 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
04222 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
04223 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
04224 }
04225
04226 if (pinf_count == 0) {
04227
04228 add_dbm_constraint(v, 0, sum);
04229
04230 deduce_u_minus_v_bounds(v, w, sc_expr, sc_den, sum);
04231 }
04232 else if (pinf_count == 1)
04233 if (pinf_index != v
04234 && expr.coefficient(Variable(pinf_index-1)) == denominator)
04235
04236
04237 add_dbm_constraint(v, pinf_index, sum);
04238 break;
04239
04240 default:
04241
04242 throw std::runtime_error("PPL internal error");
04243 }
04244 assert(OK());
04245 }
04246
04247 template <typename T>
04248 void
04249 BD_Shape<T>::generalized_affine_image(const Linear_Expression& lhs,
04250 const Relation_Symbol relsym,
04251 const Linear_Expression& rhs) {
04252
04253
04254
04255 const dimension_type space_dim = space_dimension();
04256 const dimension_type lhs_space_dim = lhs.space_dimension();
04257 if (space_dim < lhs_space_dim)
04258 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
04259 "e1", lhs);
04260
04261
04262
04263 const dimension_type rhs_space_dim = rhs.space_dimension();
04264 if (space_dim < rhs_space_dim)
04265 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
04266 "e2", rhs);
04267
04268
04269 if (relsym == LESS_THAN || relsym == GREATER_THAN)
04270 throw_generic("generalized_affine_image(e1, r, e2)",
04271 "r is a strict relation symbol and "
04272 "*this is a BD_Shape");
04273
04274 if (relsym == NOT_EQUAL)
04275 throw_generic("generalized_affine_image(e1, r, e2)",
04276 "r is the disequality relation symbol and "
04277 "*this is a BD_Shape");
04278
04279
04280 shortest_path_closure_assign();
04281 if (marked_empty())
04282 return;
04283
04284
04285
04286 dimension_type t_lhs = 0;
04287
04288 dimension_type j_lhs = 0;
04289
04290 for (dimension_type i = lhs_space_dim; i-- > 0; )
04291 if (lhs.coefficient(Variable(i)) != 0) {
04292 if (t_lhs++ == 1)
04293 break;
04294 else
04295 j_lhs = i;
04296 }
04297
04298 const Coefficient& b_lhs = lhs.inhomogeneous_term();
04299
04300 if (t_lhs == 0) {
04301
04302
04303
04304
04305
04306
04307
04308 switch (relsym) {
04309 case LESS_OR_EQUAL:
04310 refine_no_check(lhs <= rhs);
04311 break;
04312 case EQUAL:
04313 refine_no_check(lhs == rhs);
04314 break;
04315 case GREATER_OR_EQUAL:
04316 refine_no_check(lhs >= rhs);
04317 break;
04318 default:
04319
04320 throw std::runtime_error("PPL internal error");
04321 }
04322 }
04323 else if (t_lhs == 1) {
04324
04325
04326
04327 Variable v(j_lhs);
04328
04329 const Coefficient& den = lhs.coefficient(v);
04330 Relation_Symbol new_relsym = relsym;
04331 if (den < 0) {
04332 if (relsym == LESS_OR_EQUAL)
04333 new_relsym = GREATER_OR_EQUAL;
04334 else if (relsym == GREATER_OR_EQUAL)
04335 new_relsym = LESS_OR_EQUAL;
04336 }
04337 Linear_Expression expr = rhs - b_lhs;
04338 generalized_affine_image(v, new_relsym, expr, den);
04339 }
04340 else {
04341
04342
04343 bool lhs_vars_intersects_rhs_vars = false;
04344 std::vector<Variable> lhs_vars;
04345 for (dimension_type i = lhs_space_dim; i-- > 0; )
04346 if (lhs.coefficient(Variable(i)) != 0) {
04347 lhs_vars.push_back(Variable(i));
04348 if (rhs.coefficient(Variable(i)) != 0)
04349 lhs_vars_intersects_rhs_vars = true;
04350 }
04351
04352 if (!lhs_vars_intersects_rhs_vars) {
04353
04354
04355 for (dimension_type i = lhs_vars.size(); i-- > 0; )
04356 forget_all_dbm_constraints(lhs_vars[i].id() + 1);
04357
04358
04359
04360
04361 switch (relsym) {
04362 case LESS_OR_EQUAL:
04363 refine_no_check(lhs <= rhs);
04364 break;
04365 case EQUAL:
04366 refine_no_check(lhs == rhs);
04367 break;
04368 case GREATER_OR_EQUAL:
04369 refine_no_check(lhs >= rhs);
04370 break;
04371 default:
04372
04373 throw std::runtime_error("PPL internal error");
04374 }
04375 }
04376 else {
04377
04378
04379 #if 1 // Simplified computation (see the TODO note below).
04380
04381 for (dimension_type i = lhs_vars.size(); i-- > 0; )
04382 forget_all_dbm_constraints(lhs_vars[i].id() + 1);
04383
04384 #else // Currently unnecessarily complex computation.
04385
04386
04387
04388
04389
04390 const Variable new_var = Variable(space_dim);
04391 add_space_dimensions_and_embed(1);
04392
04393
04394
04395
04396 affine_image(new_var, rhs);
04397
04398
04399 shortest_path_closure_assign();
04400 assert(!marked_empty());
04401 for (dimension_type i = lhs_vars.size(); i-- > 0; )
04402 forget_all_dbm_constraints(lhs_vars[i].id() + 1);
04403
04404
04405
04406
04407
04408
04409 switch (relsym) {
04410 case LESS_OR_EQUAL:
04411 refine_no_check(lhs <= new_var);
04412 break;
04413 case EQUAL:
04414 refine_no_check(lhs == new_var);
04415 break;
04416 case GREATER_OR_EQUAL:
04417 refine_no_check(lhs >= new_var);
04418 break;
04419 default:
04420
04421 throw std::runtime_error("PPL internal error");
04422 }
04423
04424 remove_higher_space_dimensions(space_dim-1);
04425 #endif // Currently unnecessarily complex computation.
04426 }
04427 }
04428
04429 assert(OK());
04430 }
04431
04432 template <typename T>
04433 void
04434 BD_Shape<T>::generalized_affine_preimage(const Variable var,
04435 const Relation_Symbol relsym,
04436 const Linear_Expression& expr,
04437 Coefficient_traits::const_reference
04438 denominator) {
04439
04440 if (denominator == 0)
04441 throw_generic("generalized_affine_preimage(v, r, e, d)", "d == 0");
04442
04443
04444
04445
04446 const dimension_type space_dim = space_dimension();
04447 const dimension_type expr_space_dim = expr.space_dimension();
04448 if (space_dim < expr_space_dim)
04449 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
04450 "e", expr);
04451
04452
04453 const dimension_type v = var.id() + 1;
04454 if (v > space_dim)
04455 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
04456 var.id());
04457
04458
04459 if (relsym == LESS_THAN || relsym == GREATER_THAN)
04460 throw_generic("generalized_affine_preimage(v, r, e, d)",
04461 "r is a strict relation symbol and "
04462 "*this is a BD_Shape");
04463
04464 if (relsym == NOT_EQUAL)
04465 throw_generic("generalized_affine_preimage(v, r, e, d)",
04466 "r is the disequality relation symbol and "
04467 "*this is a BD_Shape");
04468
04469 if (relsym == EQUAL) {
04470
04471
04472 affine_preimage(var, expr, denominator);
04473 return;
04474 }
04475
04476
04477 shortest_path_closure_assign();
04478 if (marked_empty())
04479 return;
04480
04481
04482
04483 const Coefficient& expr_v = expr.coefficient(var);
04484 if (expr_v != 0) {
04485 const Relation_Symbol reversed_relsym = (relsym == LESS_OR_EQUAL)
04486 ? GREATER_OR_EQUAL : LESS_OR_EQUAL;
04487 const Linear_Expression inverse
04488 = expr - (expr_v + denominator)*var;
04489 TEMP_INTEGER(inverse_den);
04490 neg_assign(inverse_den, expr_v);
04491 const Relation_Symbol inverse_relsym
04492 = (sgn(denominator) == sgn(inverse_den)) ? relsym : reversed_relsym;
04493 generalized_affine_image(var, inverse_relsym, inverse, inverse_den);
04494 return;
04495 }
04496
04497 refine(var, relsym, expr, denominator);
04498
04499 if (is_empty())
04500 return;
04501
04502
04503 forget_all_dbm_constraints(v);
04504
04505 if (marked_shortest_path_reduced())
04506 reset_shortest_path_reduced();
04507 assert(OK());
04508 }
04509
04510 template <typename T>
04511 void
04512 BD_Shape<T>::generalized_affine_preimage(const Linear_Expression& lhs,
04513 const Relation_Symbol relsym,
04514 const Linear_Expression& rhs) {
04515
04516
04517
04518 const dimension_type bds_space_dim = space_dimension();
04519 const dimension_type lhs_space_dim = lhs.space_dimension();
04520 if (bds_space_dim < lhs_space_dim)
04521 throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
04522 "e1", lhs);
04523
04524
04525
04526 const dimension_type rhs_space_dim = rhs.space_dimension();
04527 if (bds_space_dim < rhs_space_dim)
04528 throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
04529 "e2", rhs);
04530
04531
04532 if (relsym == LESS_THAN || relsym == GREATER_THAN)
04533 throw_generic("generalized_affine_preimage(e1, r, e2)",
04534 "r is a strict relation symbol and "
04535 "*this is a BD_Shape");
04536
04537 if (relsym == NOT_EQUAL)
04538 throw_generic("generalized_affine_preimage(e1, r, e2)",
04539 "r is the disequality relation symbol and "
04540 "*this is a BD_Shape");
04541
04542
04543 shortest_path_closure_assign();
04544 if (marked_empty())
04545 return;
04546
04547
04548
04549 dimension_type t_lhs = 0;
04550
04551 dimension_type j_lhs = 0;
04552
04553 for (dimension_type i = lhs_space_dim; i-- > 0; )
04554 if (lhs.coefficient(Variable(i)) != 0) {
04555 if (t_lhs++ == 1)
04556 break;
04557 else
04558 j_lhs = i;
04559 }
04560
04561 const Coefficient& b_lhs = lhs.inhomogeneous_term();
04562
04563 if (t_lhs == 0) {
04564
04565
04566 generalized_affine_image(lhs, relsym, rhs);
04567 return;
04568 }
04569 else if (t_lhs == 1) {
04570
04571
04572
04573 Variable v(j_lhs);
04574
04575 const Coefficient& den = lhs.coefficient(v);
04576 Relation_Symbol new_relsym = relsym;
04577 if (den < 0) {
04578 if (relsym == LESS_OR_EQUAL)
04579 new_relsym = GREATER_OR_EQUAL;
04580 else if (relsym == GREATER_OR_EQUAL)
04581 new_relsym = LESS_OR_EQUAL;
04582 }
04583 Linear_Expression expr = rhs - b_lhs;
04584 generalized_affine_preimage(v, new_relsym, expr, den);
04585 }
04586 else {
04587
04588
04589 bool lhs_vars_intersects_rhs_vars = false;
04590 std::vector<Variable> lhs_vars;
04591 for (dimension_type i = lhs_space_dim; i-- > 0; )
04592 if (lhs.coefficient(Variable(i)) != 0) {
04593 lhs_vars.push_back(Variable(i));
04594 if (rhs.coefficient(Variable(i)) != 0)
04595 lhs_vars_intersects_rhs_vars = true;
04596 }
04597
04598 if (!lhs_vars_intersects_rhs_vars) {
04599
04600
04601
04602
04603
04604
04605 switch (relsym) {
04606 case LESS_OR_EQUAL:
04607 refine_no_check(lhs <= rhs);
04608 break;
04609 case EQUAL:
04610 refine_no_check(lhs == rhs);
04611 break;
04612 case GREATER_OR_EQUAL:
04613 refine_no_check(lhs >= rhs);
04614 break;
04615 default:
04616
04617 throw std::runtime_error("PPL internal error");
04618 }
04619
04620
04621 if (is_empty())
04622 return;
04623
04624 for (dimension_type i = lhs_vars.size(); i-- > 0; )
04625 forget_all_dbm_constraints(lhs_vars[i].id() + 1);
04626 }
04627 else {
04628
04629
04630
04631 const Variable new_var = Variable(bds_space_dim);
04632 add_space_dimensions_and_embed(1);
04633
04634
04635
04636
04637 affine_image(new_var, lhs);
04638
04639
04640 shortest_path_closure_assign();
04641 assert(!marked_empty());
04642 for (dimension_type i = lhs_vars.size(); i-- > 0; )
04643 forget_all_dbm_constraints(lhs_vars[i].id() + 1);
04644
04645
04646
04647
04648
04649
04650
04651 switch (relsym) {
04652 case LESS_OR_EQUAL:
04653 refine_no_check(new_var <= rhs);
04654 break;
04655 case EQUAL:
04656 refine_no_check(new_var == rhs);
04657 break;
04658 case GREATER_OR_EQUAL:
04659 refine_no_check(new_var >= rhs);
04660 break;
04661 default:
04662
04663 throw std::runtime_error("PPL internal error");
04664 }
04665
04666 remove_higher_space_dimensions(bds_space_dim);
04667 }
04668 }
04669
04670 assert(OK());
04671 }
04672
04673 template <typename T>
04674 Constraint_System
04675 BD_Shape<T>::constraints() const {
04676 Constraint_System cs;
04677 const dimension_type space_dim = space_dimension();
04678 if (space_dim == 0) {
04679 if (marked_empty())
04680 cs = Constraint_System::zero_dim_empty();
04681 }
04682 else if (marked_empty())
04683 cs.insert(0*Variable(space_dim-1) <= -1);
04684 else if (marked_shortest_path_reduced())
04685
04686 cs = minimized_constraints();
04687 else {
04688
04689
04690 cs.insert(0*Variable(space_dim-1) <= 0);
04691
04692 TEMP_INTEGER(a);
04693 TEMP_INTEGER(b);
04694
04695 const DB_Row<N>& dbm_0 = dbm[0];
04696 for (dimension_type j = 1; j <= space_dim; ++j) {
04697 const Variable x(j-1);
04698 const N& dbm_0j = dbm_0[j];
04699 const N& dbm_j0 = dbm[j][0];
04700 if (is_additive_inverse(dbm_j0, dbm_0j)) {
04701
04702 numer_denom(dbm_0j, b, a);
04703 cs.insert(a*x == b);
04704 }
04705 else {
04706
04707 if (!is_plus_infinity(dbm_0j)) {
04708 numer_denom(dbm_0j, b, a);
04709 cs.insert(a*x <= b);
04710 }
04711 if (!is_plus_infinity(dbm_j0)) {
04712 numer_denom(dbm_j0, b, a);
04713 cs.insert(-a*x <= b);
04714 }
04715 }
04716 }
04717
04718
04719 for (dimension_type i = 1; i <= space_dim; ++i) {
04720 const Variable y(i-1);
04721 const DB_Row<N>& dbm_i = dbm[i];
04722 for (dimension_type j = i + 1; j <= space_dim; ++j) {
04723 const Variable x(j-1);
04724 const N& dbm_ij = dbm_i[j];
04725 const N& dbm_ji = dbm[j][i];
04726 if (is_additive_inverse(dbm_ji, dbm_ij)) {
04727
04728 numer_denom(dbm_ij, b, a);
04729 cs.insert(a*x - a*y == b);
04730 }
04731 else {
04732
04733 if (!is_plus_infinity(dbm_ij)) {
04734 numer_denom(dbm_ij, b, a);
04735 cs.insert(a*x - a*y <= b);
04736 }
04737 if (!is_plus_infinity(dbm_ji)) {
04738 numer_denom(dbm_ji, b, a);
04739 cs.insert(a*y - a*x <= b);
04740 }
04741 }
04742 }
04743 }
04744 }
04745 return cs;
04746 }
04747
04748 template <typename T>
04749 Constraint_System
04750 BD_Shape<T>::minimized_constraints() const {
04751 shortest_path_reduction_assign();
04752 Constraint_System cs;
04753 const dimension_type space_dim = space_dimension();
04754 if (space_dim == 0) {
04755 if (marked_empty())
04756 cs = Constraint_System::zero_dim_empty();
04757 }
04758 else if (marked_empty())
04759 cs.insert(0*Variable(space_dim-1) <= -1);
04760 else {
04761
04762
04763 cs.insert(0*Variable(space_dim-1) <= 0);
04764
04765 TEMP_INTEGER(num);
04766 TEMP_INTEGER(den);
04767
04768
04769 std::vector<dimension_type> leaders;
04770 compute_leaders(leaders);
04771 std::vector<dimension_type> leader_indices;
04772 compute_leader_indices(leaders, leader_indices);
04773 const dimension_type num_leaders = leader_indices.size();
04774
04775
04776 const DB_Row<N>& dbm_0 = dbm[0];
04777 for (dimension_type i = 1; i <= space_dim; ++i) {
04778 const dimension_type leader = leaders[i];
04779 if (i != leader) {
04780
04781 if (leader == 0) {
04782
04783 assert(!is_plus_infinity(dbm_0[i]));
04784 numer_denom(dbm_0[i], num, den);
04785 cs.insert(den*Variable(i-1) == num);
04786 }
04787 else {
04788
04789 assert(!is_plus_infinity(dbm[i][leader]));
04790 numer_denom(dbm[i][leader], num, den);
04791 cs.insert(den*Variable(leader-1) - den*Variable(i-1) == num);
04792 }
04793 }
04794 }
04795
04796
04797
04798 const Bit_Row& red_0 = redundancy_dbm[0];
04799 for (dimension_type l_i = 1; l_i < num_leaders; ++l_i) {
04800 const dimension_type i = leader_indices[l_i];
04801 if (!red_0[i]) {
04802 numer_denom(dbm_0[i], num, den);
04803 cs.insert(den*Variable(i-1) <= num);
04804 }
04805 if (!redundancy_dbm[i][0]) {
04806 numer_denom(dbm[i][0], num, den);
04807 cs.insert(-den*Variable(i-1) <= num);
04808 }
04809 }
04810
04811 for (dimension_type l_i = 1; l_i < num_leaders; ++l_i) {
04812 const dimension_type i = leader_indices[l_i];
04813 const DB_Row<N>& dbm_i = dbm[i];
04814 const Bit_Row& red_i = redundancy_dbm[i];
04815 for (dimension_type l_j = l_i + 1; l_j < num_leaders; ++l_j) {
04816 const dimension_type j = leader_indices[l_j];
04817 if (!red_i[j]) {
04818 numer_denom(dbm_i[j], num, den);
04819 cs.insert(den*Variable(j-1) - den*Variable(i-1) <= num);
04820 }
04821 if (!redundancy_dbm[j][i]) {
04822 numer_denom(dbm[j][i], num, den);
04823 cs.insert(den*Variable(i-1) - den*Variable(j-1) <= num);
04824 }
04825 }
04826 }
04827 }
04828 return cs;
04829 }
04830
04831 template <typename T>
04832 void
04833 BD_Shape<T>::expand_space_dimension(Variable var, dimension_type m) {
04834 dimension_type old_dim = space_dimension();
04835
04836 if (var.space_dimension() > old_dim)
04837 throw_dimension_incompatible("expand_space_dimension(v, m)", "v", var);
04838
04839
04840
04841 if (m > max_space_dimension() - space_dimension())
04842 throw_generic("expand_dimension(v, m)",
04843 "adding m new space dimensions exceeds "
04844 "the maximum allowed space dimension");
04845
04846
04847 if (m == 0)
04848 return;
04849
04850
04851 add_space_dimensions_and_embed(m);
04852
04853
04854
04855
04856 const dimension_type v_id = var.id() + 1;
04857 const DB_Row<N>& dbm_v = dbm[v_id];
04858 for (dimension_type i = old_dim + 1; i-- > 0; ) {
04859 DB_Row<N>& dbm_i = dbm[i];
04860 const N& dbm_i_v = dbm[i][v_id];
04861 const N& dbm_v_i = dbm_v[i];
04862 for (dimension_type j = old_dim+1; j < old_dim+m+1; ++j) {
04863 dbm_i[j] = dbm_i_v;
04864 dbm[j][i] = dbm_v_i;
04865 }
04866 }
04867
04868
04869 if (marked_shortest_path_closed())
04870 reset_shortest_path_closed();
04871 assert(OK());
04872 }
04873
04874 template <typename T>
04875 void
04876 BD_Shape<T>::fold_space_dimensions(const Variables_Set& to_be_folded,
04877 Variable var) {
04878 const dimension_type space_dim = space_dimension();
04879
04880 if (var.space_dimension() > space_dim)
04881 throw_dimension_incompatible("fold_space_dimensions(tbf, v)",
04882 "v", var);
04883
04884
04885 if (to_be_folded.empty())
04886 return;
04887
04888
04889 if (to_be_folded.space_dimension() > space_dim)
04890 throw_dimension_incompatible("fold_space_dimensions(tbf, ...)",
04891 to_be_folded.space_dimension());
04892
04893
04894 if (to_be_folded.find(var.id()) != to_be_folded.end())
04895 throw_generic("fold_space_dimensions(tbf, v)",
04896 "v should not occur in tbf");
04897
04898 shortest_path_closure_assign();
04899 if (!marked_empty()) {
04900
04901
04902
04903
04904 const dimension_type v_id = var.id() + 1;
04905 DB_Row<N>& dbm_v = dbm[v_id];
04906 for (Variables_Set::const_iterator i = to_be_folded.begin(),
04907 tbf_end = to_be_folded.end(); i != tbf_end; ++i) {
04908 const dimension_type tbf_id = *i + 1;
04909 const DB_Row<N>& dbm_tbf = dbm[tbf_id];
04910 for (dimension_type j = space_dim + 1; j-- > 0; ) {
04911 max_assign(dbm[j][v_id], dbm[j][tbf_id]);
04912 max_assign(dbm_v[j], dbm_tbf[j]);
04913 }
04914 }
04915 }
04916 remove_space_dimensions(to_be_folded);
04917 }
04918
04920 template <typename T>
04921 std::ostream&
04922 IO_Operators::operator<<(std::ostream& s, const BD_Shape<T>& c) {
04923 typedef typename BD_Shape<T>::coefficient_type N;
04924 if (c.is_universe())
04925 s << "true";
04926 else {
04927
04928 dimension_type n = c.space_dimension();
04929 if (c.marked_empty())
04930 s << "false";
04931 else {
04932 DIRTY_TEMP(N, v);
04933 bool first = true;
04934 for (dimension_type i = 0; i <= n; ++i)
04935 for (dimension_type j = i + 1; j <= n; ++j) {
04936 const N& c_i_j = c.dbm[i][j];
04937 const N& c_j_i = c.dbm[j][i];
04938 if (is_additive_inverse(c_j_i, c_i_j)) {
04939
04940 if (first)
04941 first = false;
04942 else
04943 s << ", ";
04944 if (i == 0) {
04945
04946 s << Variable(j - 1);
04947 s << " == " << c_i_j;
04948 }
04949 else {
04950
04951 if (sgn(c_i_j) >= 0) {
04952 s << Variable(j - 1);
04953 s << " - ";
04954 s << Variable(i - 1);
04955 s << " == " << c_i_j;
04956 }
04957 else {
04958 s << Variable(i - 1);
04959 s << " - ";
04960 s << Variable(j - 1);
04961 s << " == " << c_j_i;
04962 }
04963 }
04964 }
04965 else {
04966
04967 if (!is_plus_infinity(c_j_i)) {
04968 if (first)
04969 first = false;
04970 else
04971 s << ", ";
04972 if (i == 0) {
04973
04974 s << Variable(j - 1);
04975 neg_assign_r(v, c_j_i, ROUND_DOWN);
04976 s << " >= " << v;
04977 }
04978 else {
04979
04980 if (sgn(c_j_i) >= 0) {
04981 s << Variable(i - 1);
04982 s << " - ";
04983 s << Variable(j - 1);
04984 s << " <= " << c_j_i;
04985 }
04986 else {
04987 s << Variable(j - 1);
04988 s << " - ";
04989 s << Variable(i - 1);
04990 neg_assign_r(v, c_j_i, ROUND_DOWN);
04991 s << " >= " << v;
04992 }
04993 }
04994 }
04995 if (!is_plus_infinity(c_i_j)) {
04996 if (first)
04997 first = false;
04998 else
04999 s << ", ";
05000 if (i == 0) {
05001
05002 s << Variable(j - 1);
05003 s << " <= " << c_i_j;
05004 }
05005 else {
05006
05007 if (sgn(c_i_j) >= 0) {
05008 s << Variable(j - 1);
05009 s << " - ";
05010 s << Variable(i - 1);
05011 s << " <= " << c_i_j;
05012 }
05013 else {
05014 s << Variable(i - 1);
05015 s << " - ";
05016 s << Variable(j - 1);
05017 neg_assign_r(v, c_i_j, ROUND_DOWN);
05018 s << " >= " << v;
05019 }
05020 }
05021 }
05022 }
05023 }
05024 }
05025 }
05026 return s;
05027 }
05028
05029 template <typename T>
05030 void
05031 BD_Shape<T>::ascii_dump(std::ostream& s) const {
05032 status.ascii_dump(s);
05033 s << "\n";
05034 dbm.ascii_dump(s);
05035 s << "\n";
05036 redundancy_dbm.ascii_dump(s);
05037 }
05038
05039 PPL_OUTPUT_TEMPLATE_DEFINITIONS(T, BD_Shape<T>)
05040
05041 template <typename T>
05042 bool
05043 BD_Shape<T>::ascii_load(std::istream& s) {
05044 if (!status.ascii_load(s))
05045 return false;
05046 if (!dbm.ascii_load(s))
05047 return false;
05048 if (!redundancy_dbm.ascii_load(s))
05049 return false;
05050 return true;
05051 }
05052
05053 template <typename T>
05054 memory_size_type
05055 BD_Shape<T>::external_memory_in_bytes() const {
05056 return dbm.external_memory_in_bytes()
05057 + redundancy_dbm.external_memory_in_bytes();
05058 }
05059
05060 template <typename T>
05061 bool
05062 BD_Shape<T>::OK() const {
05063
05064 if (!dbm.OK())
05065 return false;
05066
05067
05068 if (!status.OK())
05069 return false;
05070
05071
05072 if (marked_empty())
05073 return true;
05074
05075
05076 for (dimension_type i = dbm.num_rows(); i-- > 0; )
05077 for (dimension_type j = dbm.num_rows(); j-- > 0; )
05078 if (is_minus_infinity(dbm[i][j])) {
05079 #ifndef NDEBUG
05080 using namespace Parma_Polyhedra_Library::IO_Operators;
05081 std::cerr << "BD_Shape::dbm[" << i << "][" << j << "] = "
05082 << dbm[i][j] << "!"
05083 << std::endl;
05084 #endif
05085 return false;
05086 }
05087
05088
05089 for (dimension_type i = dbm.num_rows(); i-- > 0; )
05090 if (!is_plus_infinity(dbm[i][i])) {
05091 #ifndef NDEBUG
05092 using namespace Parma_Polyhedra_Library::IO_Operators;
05093 std::cerr << "BD_Shape::dbm[" << i << "][" << i << "] = "
05094 << dbm[i][i] << "! (+inf was expected.)"
05095 << std::endl;
05096 #endif
05097 return false;
05098 }
05099
05100
05101 if (marked_shortest_path_closed()) {
05102 BD_Shape x = *this;
05103 x.reset_shortest_path_closed();
05104 x.shortest_path_closure_assign();
05105 if (x.dbm != dbm) {
05106 #ifndef NDEBUG
05107 std::cerr << "BD_Shape is marked as closed but it is not!"
05108 << std::endl;
05109 #endif
05110 return false;
05111 }
05112 }
05113
05114
05115
05116
05117 if (std::numeric_limits<coefficient_type_base>::is_exact) {
05118
05119
05120 if (marked_shortest_path_reduced()) {
05121
05122 for (dimension_type i = dbm.num_rows(); i-- > 0; )
05123 for (dimension_type j = dbm.num_rows(); j-- > 0; )
05124 if (!redundancy_dbm[i][j] && is_plus_infinity(dbm[i][j])) {
05125 #ifndef NDEBUG
05126 using namespace Parma_Polyhedra_Library::IO_Operators;
05127 std::cerr << "BD_Shape::dbm[" << i << "][" << j << "] = "
05128 << dbm[i][j] << " is marked as non-redundant!"
05129 << std::endl;
05130 #endif
05131 return false;
05132 }
05133
05134 BD_Shape x = *this;
05135 x.reset_shortest_path_reduced();
05136 x.shortest_path_reduction_assign();
05137 if (x.redundancy_dbm != redundancy_dbm) {
05138 #ifndef NDEBUG
05139 std::cerr << "BD_Shape is marked as reduced but it is not!"
05140 << std::endl;
05141 #endif
05142 return false;
05143 }
05144 }
05145 }
05146
05147
05148 return true;
05149 }
05150
05151 template <typename T>
05152 void
05153 BD_Shape<T>::throw_dimension_incompatible(const char* method,
05154 const BD_Shape& y) const {
05155 std::ostringstream s;
05156 s << "PPL::BD_Shape::" << method << ":" << std::endl
05157 << "this->space_dimension() == " << space_dimension()
05158 << ", y->space_dimension() == " << y.space_dimension() << ".";
05159 throw std::invalid_argument(s.str());
05160 }
05161
05162 template <typename T>
05163 void
05164 BD_Shape<T>::throw_dimension_incompatible(const char* method,
05165 dimension_type required_dim) const {
05166 std::ostringstream s;
05167 s << "PPL::BD_Shape::" << method << ":" << std::endl
05168 << "this->space_dimension() == " << space_dimension()
05169 << ", required dimension == " << required_dim << ".";
05170 throw std::invalid_argument(s.str());
05171 }
05172
05173 template <typename T>
05174 void
05175 BD_Shape<T>::throw_dimension_incompatible(const char* method,
05176 const Constraint& c) const {
05177 std::ostringstream s;
05178 s << "PPL::BD_Shape::" << method << ":" << std::endl
05179 << "this->space_dimension() == " << space_dimension()
05180 << ", c->space_dimension == " << c.space_dimension() << ".";
05181 throw std::invalid_argument(s.str());
05182 }
05183
05184 template <typename T>
05185 void
05186 BD_Shape<T>::throw_dimension_incompatible(const char* method,
05187 const Congruence& cg) const {
05188 std::ostringstream s;
05189 s << "PPL::BD_Shape::" << method << ":" << std::endl
05190 << "this->space_dimension() == " << space_dimension()
05191 << ", cg->space_dimension == " << cg.space_dimension() << ".";
05192 throw std::invalid_argument(s.str());
05193 }
05194
05195 template <typename T>
05196 void
05197 BD_Shape<T>::throw_dimension_incompatible(const char* method,
05198 const Generator& g) const {
05199 std::ostringstream s;
05200 s << "PPL::BD_Shape::" << method << ":" << std::endl
05201 << "this->space_dimension() == " << space_dimension()
05202 << ", g->space_dimension == " << g.space_dimension() << ".";
05203 throw std::invalid_argument(s.str());
05204 }
05205
05206 template <typename T>
05207 void
05208 BD_Shape<T>::throw_expression_too_complex(const char* method,
05209 const Linear_Expression& e) {
05210 using namespace IO_Operators;
05211 std::ostringstream s;
05212 s << "PPL::BD_Shape::" << method << ":" << std::endl
05213 << e << " is too complex.";
05214 throw std::invalid_argument(s.str());
05215 }
05216
05217
05218 template <typename T>
05219 void
05220 BD_Shape<T>::throw_dimension_incompatible(const char* method,
05221 const char* name_row,
05222 const Linear_Expression& y) const {
05223 std::ostringstream s;
05224 s << "PPL::BD_Shape::" << method << ":" << std::endl
05225 << "this->space_dimension() == " << space_dimension()
05226 << ", " << name_row << "->space_dimension() == "
05227 << y.space_dimension() << ".";
05228 throw std::invalid_argument(s.str());
05229 }
05230
05231 template <typename T>
05232 void
05233 BD_Shape<T>::throw_generic(const char* method, const char* reason) {
05234 std::ostringstream s;
05235 s << "PPL::BD_Shape::" << method << ":" << std::endl
05236 << reason << ".";
05237 throw std::invalid_argument(s.str());
05238 }
05239
05240 }
05241
05242 #endif // !defined(PPL_BD_Shape_templates_hh)