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_Octagonal_Shape_templates_hh
00024 #define PPL_Octagonal_Shape_templates_hh 1
00025
00026 #include "Generator_System.defs.hh"
00027 #include "Generator_System.inlines.hh"
00028 #include "Congruence_System.defs.hh"
00029 #include "Congruence_System.inlines.hh"
00030 #include <cassert>
00031 #include <vector>
00032 #include <deque>
00033 #include <string>
00034 #include <iostream>
00035 #include <sstream>
00036 #include <stdexcept>
00037 #include <algorithm>
00038
00039 namespace Parma_Polyhedra_Library {
00040
00041 template <typename T>
00042 Octagonal_Shape<T>::Octagonal_Shape(const Polyhedron& ph,
00043 const Complexity_Class complexity)
00044 : matrix(0), space_dim(0), status() {
00045 const dimension_type num_dimensions = ph.space_dimension();
00046
00047 if (ph.marked_empty()) {
00048 *this = Octagonal_Shape(num_dimensions, EMPTY);
00049 return;
00050 }
00051
00052 if (num_dimensions == 0) {
00053 *this = Octagonal_Shape(num_dimensions, UNIVERSE);
00054 return;
00055 }
00056
00057
00058
00059 if (complexity == ANY_COMPLEXITY
00060 || (!ph.has_pending_constraints() && ph.generators_are_up_to_date())) {
00061 *this = Octagonal_Shape(ph.generators());
00062 return;
00063 }
00064
00065
00066
00067
00068 assert(ph.constraints_are_up_to_date());
00069
00070 if (!ph.has_something_pending() && ph.constraints_are_minimized()) {
00071
00072
00073 if (ph.is_universe()) {
00074 *this = Octagonal_Shape(num_dimensions, UNIVERSE);
00075 return;
00076 }
00077 }
00078
00079
00080 for (Constraint_System::const_iterator i = ph.con_sys.begin(),
00081 cs_end = ph.con_sys.end(); i != cs_end; ++i)
00082 if (i->is_inconsistent()) {
00083 *this = Octagonal_Shape(num_dimensions, EMPTY);
00084 return;
00085 }
00086
00087
00088
00089 if (complexity == SIMPLEX_COMPLEXITY) {
00090 MIP_Problem lp(num_dimensions);
00091 lp.set_optimization_mode(MAXIMIZATION);
00092
00093 const Constraint_System& ph_cs = ph.constraints();
00094 if (!ph_cs.has_strict_inequalities())
00095 lp.add_constraints(ph_cs);
00096 else
00097
00098 for (Constraint_System::const_iterator i = ph_cs.begin(),
00099 ph_cs_end = ph_cs.end(); i != ph_cs_end; ++i) {
00100 const Constraint& c = *i;
00101 if (c.is_strict_inequality())
00102 lp.add_constraint(Linear_Expression(c) >= 0);
00103 else
00104 lp.add_constraint(c);
00105 }
00106
00107
00108 if (!lp.is_satisfiable()) {
00109 *this = Octagonal_Shape<T>(num_dimensions, EMPTY);
00110 return;
00111 }
00112
00113
00114 *this = Octagonal_Shape<T>(num_dimensions, UNIVERSE);
00115
00116 Generator g(point());
00117 TEMP_INTEGER(num);
00118 TEMP_INTEGER(den);
00119 for (dimension_type i = 0; i < num_dimensions; ++i) {
00120 Variable x(i);
00121
00122 lp.set_objective_function(x);
00123 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00124 g = lp.optimizing_point();
00125 lp.evaluate_objective_function(g, num, den);
00126 num *= 2;
00127 div_round_up(matrix[2*i+1][2*i], num, den);
00128 }
00129
00130 for (dimension_type j = 0; j < i; ++j) {
00131 Variable y(j);
00132 lp.set_objective_function(x + y);
00133 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00134 g = lp.optimizing_point();
00135 lp.evaluate_objective_function(g, num, den);
00136 div_round_up(matrix[2*i+1][2*j], num, den);
00137 }
00138 }
00139
00140 for (dimension_type j = 0; j < num_dimensions; ++j) {
00141 if (i == j)
00142 continue;
00143 Variable y(j);
00144 lp.set_objective_function(x - y);
00145 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00146 g = lp.optimizing_point();
00147 lp.evaluate_objective_function(g, num, den);
00148 div_round_up((i < j ? matrix[2*j][2*i] : matrix[2*i+1][2*j+1]),
00149 num, den);
00150 }
00151 }
00152
00153 for (dimension_type j = 0; j < num_dimensions; ++j) {
00154 if (i == j)
00155 continue;
00156 Variable y(j);
00157 lp.set_objective_function(x - y);
00158 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00159 g = lp.optimizing_point();
00160 lp.evaluate_objective_function(g, num, den);
00161 div_round_up((i < j ? matrix[2*j][2*i] : matrix[2*i+1][2*j+1]),
00162 num, den);
00163 }
00164 }
00165
00166 for (dimension_type j = 0; j < i; ++j) {
00167 Variable y(j);
00168 lp.set_objective_function(-x - y);
00169 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00170 g = lp.optimizing_point();
00171 lp.evaluate_objective_function(g, num, den);
00172 div_round_up(matrix[2*i][2*j+1], num, den);
00173 }
00174 }
00175
00176 lp.set_objective_function(-x);
00177 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
00178 g = lp.optimizing_point();
00179 lp.evaluate_objective_function(g, num, den);
00180 num *= 2;
00181 div_round_up(matrix[2*i][2*i+1], num, den);
00182 }
00183 }
00184 set_strongly_closed();
00185 assert(OK());
00186 return;
00187 }
00188
00189
00190 assert(complexity == POLYNOMIAL_COMPLEXITY);
00191 *this = Octagonal_Shape(num_dimensions, UNIVERSE);
00192 refine_with_constraints(ph.constraints());
00193 }
00194
00195 template <typename T>
00196 Octagonal_Shape<T>::Octagonal_Shape(const Generator_System& gs)
00197 : matrix(gs.space_dimension()),
00198 space_dim(gs.space_dimension()),
00199 status() {
00200 const Generator_System::const_iterator gs_begin = gs.begin();
00201 const Generator_System::const_iterator gs_end = gs.end();
00202 if (gs_begin == gs_end) {
00203
00204 set_empty();
00205 return;
00206 }
00207
00208 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
00209 typename OR_Matrix<N>::row_iterator mat_begin = matrix.row_begin();
00210
00211 DIRTY_TEMP(N, tmp);
00212 bool mat_initialized = false;
00213 bool point_seen = false;
00214
00215 for (Generator_System::const_iterator k = gs_begin; k != gs_end; ++k) {
00216 const Generator& g = *k;
00217 switch (g.type()) {
00218 case Generator::POINT:
00219 point_seen = true;
00220
00221 case Generator::CLOSURE_POINT:
00222 if (!mat_initialized) {
00223
00224 mat_initialized = true;
00225 const Coefficient& d = g.divisor();
00226 for (dimension_type i = 0; i < space_dim; ++i) {
00227 const Coefficient& g_i = g.coefficient(Variable(i));
00228 const dimension_type di = 2*i;
00229 Row_Reference x_i = *(mat_begin+di);
00230 Row_Reference x_ii = *(mat_begin+di+1);
00231 for (dimension_type j = 0; j < i; ++j) {
00232 const Coefficient& g_j = g.coefficient(Variable(j));
00233 const dimension_type dj = 2*j;
00234
00235
00236
00237
00238 div_round_up(x_i[dj], g_j - g_i, d);
00239 div_round_up(x_ii[dj+1], g_i - g_j, d);
00240
00241 div_round_up(x_i[dj+1], -g_j - g_i, d);
00242 div_round_up(x_ii[dj], g_i + g_j, d);
00243 }
00244
00245 div_round_up(x_i[di+1], -g_i - g_i, d);
00246 div_round_up(x_ii[di], g_i + g_i, d);
00247 }
00248 }
00249 else {
00250
00251
00252 const Coefficient& d = g.divisor();
00253 for (dimension_type i = 0; i < space_dim; ++i) {
00254 const Coefficient& g_i = g.coefficient(Variable(i));
00255 const dimension_type di = 2*i;
00256 Row_Reference x_i = *(mat_begin+di);
00257 Row_Reference x_ii = *(mat_begin+di+1);
00258 for (dimension_type j = 0; j < i; ++j) {
00259 const Coefficient& g_j = g.coefficient(Variable(j));
00260 const dimension_type dj = 2*j;
00261
00262
00263
00264
00265 div_round_up(tmp, g_j - g_i, d);
00266 max_assign(x_i[dj], tmp);
00267 div_round_up(tmp, g_i - g_j, d);
00268 max_assign(x_ii[dj+1], tmp);
00269
00270 div_round_up(tmp, -g_j - g_i, d);
00271 max_assign(x_i[dj+1], tmp);
00272 div_round_up(tmp, g_i + g_j, d);
00273 max_assign(x_ii[dj], tmp);
00274 }
00275
00276 div_round_up(tmp, -g_i - g_i, d);
00277 max_assign(x_i[di+1], tmp);
00278 div_round_up(tmp, g_i + g_i, d);
00279 max_assign(x_ii[di], tmp);
00280 }
00281 }
00282 break;
00283 default:
00284
00285 break;
00286 }
00287 }
00288
00289 if (!point_seen)
00290
00291 throw_generic("Octagonal_Shape(gs)",
00292 "the non-empty generator system gs contains no points.");
00293
00294
00295 for (Generator_System::const_iterator k = gs_begin; k != gs_end; ++k) {
00296 const Generator& g = *k;
00297 switch (g.type()) {
00298 case Generator::LINE:
00299 for (dimension_type i = 0; i < space_dim; ++i) {
00300 const Coefficient& g_i = g.coefficient(Variable(i));
00301 const dimension_type di = 2*i;
00302 Row_Reference x_i = *(mat_begin+di);
00303 Row_Reference x_ii = *(mat_begin+di+1);
00304 for (dimension_type j = 0; j < i; ++j) {
00305 const Coefficient& g_j = g.coefficient(Variable(j));
00306 const dimension_type dj = 2*j;
00307
00308 if (g_i != g_j) {
00309
00310 assign_r(x_i[dj], PLUS_INFINITY, ROUND_NOT_NEEDED);
00311 assign_r(x_ii[dj+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
00312 }
00313 if (g_i != -g_j) {
00314
00315 assign_r(x_i[dj+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
00316 assign_r(x_ii[dj], PLUS_INFINITY, ROUND_NOT_NEEDED);
00317 }
00318 }
00319 if (g_i != 0) {
00320
00321 assign_r(x_i[di+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
00322 assign_r(x_ii[di], PLUS_INFINITY, ROUND_NOT_NEEDED);
00323 }
00324 }
00325 break;
00326 case Generator::RAY:
00327 for (dimension_type i = 0; i < space_dim; ++i) {
00328 const Coefficient& g_i = g.coefficient(Variable(i));
00329 const dimension_type di = 2*i;
00330 Row_Reference x_i = *(mat_begin+di);
00331 Row_Reference x_ii = *(mat_begin+di+1);
00332 for (dimension_type j = 0; j < i; ++j) {
00333 const Coefficient& g_j = g.coefficient(Variable(j));
00334 const dimension_type dj = 2*j;
00335
00336
00337 if (g_i < g_j)
00338
00339 assign_r(x_i[dj], PLUS_INFINITY, ROUND_NOT_NEEDED);
00340 if (g_i > g_j)
00341
00342 assign_r(x_ii[dj+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
00343 if (g_i < -g_j)
00344
00345 assign_r(x_i[dj+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
00346 if (g_i > -g_j)
00347
00348 assign_r(x_ii[dj], PLUS_INFINITY, ROUND_NOT_NEEDED);
00349 }
00350
00351 if (g_i < 0)
00352
00353 assign_r(x_i[di+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
00354 if (g_i > 0)
00355
00356 assign_r(x_ii[di], PLUS_INFINITY, ROUND_NOT_NEEDED);
00357 }
00358 break;
00359 default:
00360
00361 break;
00362 }
00363 }
00364 set_strongly_closed();
00365 assert(OK());
00366 }
00367
00368 template <typename T>
00369 void
00370 Octagonal_Shape<T>::add_constraint(const Constraint& c) {
00371 const dimension_type c_space_dim = c.space_dimension();
00372
00373 if (c_space_dim > space_dim)
00374 throw_dimension_incompatible("add_constraint(c)", c);
00375
00376
00377 if (c.is_strict_inequality()) {
00378 if (c.is_inconsistent()) {
00379 set_empty();
00380 return;
00381 }
00382 if (c.is_tautological())
00383 return;
00384
00385 throw_generic("add_constraint(c)", "strict inequalities are not allowed");
00386 }
00387
00388 dimension_type num_vars = 0;
00389 dimension_type i = 0;
00390 dimension_type j = 0;
00391 TEMP_INTEGER(coeff);
00392 TEMP_INTEGER(term);
00393
00394 if (!extract_octagonal_difference(c, c_space_dim, num_vars,
00395 i, j, coeff, term))
00396 throw_generic("add_constraint(c)",
00397 "c is not an octagonal constraint");
00398
00399 if (num_vars == 0) {
00400
00401 if (c.inhomogeneous_term() < 0
00402 || (c.is_equality() && c.inhomogeneous_term() != 0))
00403 set_empty();
00404 return;
00405 }
00406
00407
00408 typename OR_Matrix<N>::row_iterator i_iter = matrix.row_begin() + i;
00409 typename OR_Matrix<N>::row_reference_type m_i = *i_iter;
00410 N& m_i_j = m_i[j];
00411
00412 if (coeff < 0)
00413 neg_assign(coeff);
00414
00415 bool is_oct_changed = false;
00416
00417 DIRTY_TEMP(N, d);
00418 div_round_up(d, term, coeff);
00419 if (m_i_j > d) {
00420 m_i_j = d;
00421 is_oct_changed = true;
00422 }
00423
00424 if (c.is_equality()) {
00425
00426 if (i%2 == 0)
00427 ++i_iter;
00428 else
00429 --i_iter;
00430
00431 typename OR_Matrix<N>::row_reference_type m_ci = *i_iter;
00432 dimension_type cj = coherent_index(j);
00433 N& m_ci_cj = m_ci[cj];
00434
00435 neg_assign(term);
00436 div_round_up(d, term, coeff);
00437 if (m_ci_cj > d) {
00438 m_ci_cj = d;
00439 is_oct_changed = true;
00440 }
00441 }
00442
00443
00444 if (is_oct_changed && marked_strongly_closed())
00445 reset_strongly_closed();
00446 assert(OK());
00447 }
00448
00449 template <typename T>
00450 void
00451 Octagonal_Shape<T>::add_congruence(const Congruence& cg) {
00452 const dimension_type cg_space_dim = cg.space_dimension();
00453
00454
00455 if (space_dimension() < cg_space_dim)
00456 throw_dimension_incompatible("add_congruence(cg)", cg);
00457
00458
00459 if (cg.is_proper_congruence()) {
00460 if (cg.is_tautological())
00461 return;
00462 if (cg.is_inconsistent()) {
00463 set_empty();
00464 return;
00465 }
00466
00467 throw_generic("add_congruence(cg)",
00468 "cg is a non-trivial, proper congruence");
00469 }
00470
00471 assert(cg.is_equality());
00472 Constraint c(cg);
00473 add_constraint(c);
00474 }
00475
00476 template <typename T>
00477 void
00478 Octagonal_Shape<T>::refine_no_check(const Constraint& c) {
00479 assert(!marked_empty());
00480 const dimension_type c_space_dim = c.space_dimension();
00481 assert(c_space_dim <= space_dim);
00482
00483 dimension_type num_vars = 0;
00484 dimension_type i = 0;
00485 dimension_type j = 0;
00486 TEMP_INTEGER(coeff);
00487 TEMP_INTEGER(term);
00488
00489 if (!extract_octagonal_difference(c, c_space_dim, num_vars,
00490 i, j, coeff, term))
00491 return;
00492
00493 if (num_vars == 0) {
00494 const Coefficient& c_inhomo = c.inhomogeneous_term();
00495
00496 if (c_inhomo < 0
00497 || (c_inhomo != 0 && c.is_equality())
00498 || (c_inhomo == 0 && c.is_strict_inequality()))
00499 set_empty();
00500 return;
00501 }
00502
00503
00504 typename OR_Matrix<N>::row_iterator i_iter = matrix.row_begin() + i;
00505 typename OR_Matrix<N>::row_reference_type m_i = *i_iter;
00506 N& m_i_j = m_i[j];
00507
00508 if (coeff < 0)
00509 neg_assign(coeff);
00510
00511 bool is_oct_changed = false;
00512
00513 DIRTY_TEMP(N, d);
00514 div_round_up(d, term, coeff);
00515 if (m_i_j > d) {
00516 m_i_j = d;
00517 is_oct_changed = true;
00518 }
00519
00520 if (c.is_equality()) {
00521
00522 if (i%2 == 0)
00523 ++i_iter;
00524 else
00525 --i_iter;
00526
00527 typename OR_Matrix<N>::row_reference_type m_ci = *i_iter;
00528 dimension_type cj = coherent_index(j);
00529 N& m_ci_cj = m_ci[cj];
00530
00531 neg_assign(term);
00532 div_round_up(d, term, coeff);
00533 if (m_ci_cj > d) {
00534 m_ci_cj = d;
00535 is_oct_changed = true;
00536 }
00537 }
00538
00539
00540 if (is_oct_changed && marked_strongly_closed())
00541 reset_strongly_closed();
00542 assert(OK());
00543 }
00544
00545 template <typename T>
00546 dimension_type
00547 Octagonal_Shape<T>::affine_dimension() const {
00548 const dimension_type n_rows = matrix.num_rows();
00549
00550 if (n_rows == 0)
00551 return 0;
00552
00553
00554
00555 strong_closure_assign();
00556 if (marked_empty())
00557 return 0;
00558
00559
00560
00561
00562
00563 std::vector<dimension_type> leaders;
00564 compute_leaders(leaders);
00565
00566
00567
00568 dimension_type affine_dim = 0;
00569 for (dimension_type i = 0; i < n_rows; i += 2)
00570
00571 if (leaders[i] == i && leaders[i+1] == i+1)
00572 ++affine_dim;
00573
00574 return affine_dim;
00575 }
00576
00577 template <typename T>
00578 Congruence_System
00579 Octagonal_Shape<T>::minimized_congruences() const {
00580
00581
00582 strong_closure_assign();
00583 const dimension_type space_dim = space_dimension();
00584 Congruence_System cgs;
00585 if (space_dim == 0) {
00586 if (marked_empty())
00587 cgs = Congruence_System::zero_dim_empty();
00588 }
00589 else if (marked_empty())
00590 cgs.insert((0*Variable(space_dim-1) %= 1) / 0);
00591 else {
00592
00593
00594 cgs.insert(0*Variable(space_dim-1) == 0);
00595
00596
00597
00598
00599 std::vector<dimension_type> leaders;
00600 compute_leaders(leaders);
00601
00602 TEMP_INTEGER(num);
00603 TEMP_INTEGER(den);
00604 for (dimension_type i = 0, i_end = 2*space_dim; i != i_end; i += 2) {
00605 const dimension_type lead_i = leaders[i];
00606 if (i == lead_i) {
00607 if (leaders[i+1] == i)
00608
00609 goto singular;
00610 else
00611
00612 continue;
00613 }
00614 else {
00615
00616 if (leaders[i+1] == lead_i)
00617
00618 goto singular;
00619 else
00620
00621 goto non_singular;
00622 }
00623
00624 singular:
00625
00626
00627 {
00628 const Variable x(i/2);
00629 const N& c_ii_i = matrix[i+1][i];
00630 #ifndef NDEBUG
00631 const N& c_i_ii = matrix[i][i+1];
00632 assert(is_additive_inverse(c_i_ii, c_ii_i));
00633 #endif
00634 numer_denom(c_ii_i, num, den);
00635 den *= 2;
00636 cgs.insert(den*x == num);
00637 }
00638 continue;
00639
00640 non_singular:
00641
00642
00643 {
00644 const N& c_i_li = matrix[i][lead_i];
00645 #ifndef NDEBUG
00646 const N& c_ii_lii = matrix[i+1][coherent_index(lead_i)];
00647 assert(is_additive_inverse(c_ii_lii, c_i_li));
00648 #endif
00649 const Variable x(lead_i/2);
00650 const Variable y(i/2);
00651 numer_denom(c_i_li, num, den);
00652 if (lead_i % 2 == 0)
00653 cgs.insert(den*x - den*y == num);
00654 else
00655 cgs.insert(den*x + den*y + num == 0);
00656 }
00657 continue;
00658 }
00659 }
00660 return cgs;
00661 }
00662
00663 template <typename T>
00664 void
00665 Octagonal_Shape<T>::concatenate_assign(const Octagonal_Shape& y) {
00666
00667
00668 if (y.space_dim == 0) {
00669 if (y.marked_empty())
00670 set_empty();
00671 return;
00672 }
00673
00674
00675
00676 if (space_dim == 0 && marked_empty()) {
00677 add_space_dimensions_and_embed(y.space_dim);
00678 return;
00679 }
00680
00681
00682
00683 dimension_type old_num_rows = matrix.num_rows();
00684
00685
00686
00687
00688
00689
00690 add_space_dimensions_and_embed(y.space_dim);
00691 typename OR_Matrix<N>::const_element_iterator
00692 y_it = y.matrix.element_begin();
00693 for(typename OR_Matrix<N>::row_iterator i = matrix.row_begin()+old_num_rows,
00694 matrix_row_end = matrix.row_end(); i != matrix_row_end; ++i) {
00695 typename OR_Matrix<N>::row_reference_type r = *i;
00696 dimension_type rs_i = i.row_size();
00697 for (dimension_type j = old_num_rows; j < rs_i; ++j, ++y_it)
00698 r[j] = *y_it;
00699 }
00700
00701
00702 if (marked_strongly_closed())
00703 reset_strongly_closed();
00704 assert(OK());
00705 }
00706
00707 template <typename T>
00708 bool
00709 Octagonal_Shape<T>::contains(const Octagonal_Shape& y) const {
00710
00711 if (space_dim != y.space_dim)
00712 throw_dimension_incompatible("contains(y)", y);
00713
00714
00715
00716
00717
00718 if (space_dim == 0) {
00719 if (!marked_empty())
00720 return true;
00721 else
00722 return y.marked_empty();
00723 }
00724
00725
00726 y.strong_closure_assign();
00727
00728 if (y.marked_empty())
00729 return true;
00730
00731
00732
00733 for (typename OR_Matrix<N>::const_element_iterator
00734 i = matrix.element_begin(), j = y.matrix.element_begin(),
00735 matrix_element_end = matrix.element_end(); i != matrix_element_end; ++i, ++j)
00736 if (*i < *j)
00737 return false;
00738 return true;
00739 }
00740
00741 template <typename T>
00742 bool
00743 Octagonal_Shape<T>::is_disjoint_from(const Octagonal_Shape& y) const {
00744
00745 if (space_dim != y.space_dim)
00746 throw_dimension_incompatible("is_disjoint_from(y)", y);
00747
00748
00749 strong_closure_assign();
00750 if (marked_empty())
00751 return true;
00752 y.strong_closure_assign();
00753 if (y.marked_empty())
00754 return true;
00755
00756
00757
00758
00759
00760
00761
00762 const dimension_type n_rows = matrix.num_rows();
00763
00764 typedef typename OR_Matrix<N>::const_row_iterator Row_Iterator;
00765 typedef typename OR_Matrix<N>::const_row_reference_type Row_Reference;
00766
00767 const Row_Iterator m_begin = matrix.row_begin();
00768 const Row_Iterator m_end = matrix.row_end();
00769
00770 const Row_Iterator y_begin = y.matrix.row_begin();
00771 const Row_Iterator y_end = y.matrix.row_end();
00772
00773 DIRTY_TEMP(N, neg_y_ci_cj);
00774 for (Row_Iterator i_iter = m_begin; i_iter != m_end; ++i_iter) {
00775 const dimension_type i = i_iter.index();
00776 const dimension_type ci = coherent_index(i);
00777 const dimension_type rs_i = i_iter.row_size();
00778 Row_Reference m_i = *i_iter;
00779 Row_Reference m_ci = (i%2) ? *(i_iter-1) : *(i_iter+1);
00780 for (dimension_type j = 0; j < n_rows; ++j) {
00781 const dimension_type cj = coherent_index(j);
00782 Row_Reference m_cj = *(m_begin + cj);
00783 const N& m_i_j = (j < rs_i) ? m_i[j] : m_cj[ci];
00784 Row_Reference y_ci = *(y_begin + ci);
00785 Row_Reference y_j = *(y_begin + j);
00786 const N& y_ci_cj = (j < rs_i) ? y_ci[cj] : y_j[i];
00787 neg_assign_r(neg_y_ci_cj, y_ci_cj, ROUND_UP);
00788 if (m_i_j < neg_y_ci_cj)
00789 return true;
00790 }
00791 }
00792 return false;
00793 }
00794
00795 template <typename T>
00796 bool
00797 Octagonal_Shape<T>::is_universe() const {
00798
00799 if (marked_empty())
00800 return false;
00801
00802
00803
00804 if (space_dim == 0)
00805 return true;
00806
00807
00808 for (typename OR_Matrix<N>::const_element_iterator
00809 i = matrix.element_begin(), matrix_element_end = matrix.element_end();
00810 i != matrix_element_end;
00811 ++i)
00812 if (!is_plus_infinity(*i))
00813 return false;
00814
00815 return true;
00816 }
00817
00818 template <typename T>
00819 bool
00820 Octagonal_Shape<T>::is_bounded() const {
00821 strong_closure_assign();
00822
00823 if (marked_empty() || space_dim == 0)
00824 return true;
00825
00826
00827 for (typename OR_Matrix<N>::const_row_iterator i = matrix.row_begin(),
00828 matrix_row_end = matrix.row_end(); i != matrix_row_end; ++i) {
00829 typename OR_Matrix<N>::const_row_reference_type x_i = *i;
00830 const dimension_type i_index = i.index();
00831 for (dimension_type j = i.row_size(); j-- > 0; )
00832 if (i_index != j)
00833 if (is_plus_infinity(x_i[j]))
00834 return false;
00835 }
00836
00837 return true;
00838 }
00839
00840 template <typename T>
00841 bool
00842 Octagonal_Shape<T>::contains_integer_point() const {
00843
00844 if (is_empty())
00845 return false;
00846 const dimension_type space_dim = space_dimension();
00847 if (space_dim == 0)
00848 return true;
00849
00850
00851
00852 if (std::numeric_limits<T>::is_integer)
00853 return !tight_coherence_would_make_empty();
00854
00855
00856
00857
00858 Octagonal_Shape<mpz_class> oct_z(space_dim);
00859 oct_z.reset_strongly_closed();
00860
00861 typedef Octagonal_Shape<mpz_class>::N Z;
00862 DIRTY_TEMP(N, tmp);
00863 bool all_integers = true;
00864 typename OR_Matrix<N>::const_element_iterator x_i = matrix.element_begin();
00865 for (typename OR_Matrix<Z>::element_iterator
00866 z_i = oct_z.matrix.element_begin(),
00867 z_end = oct_z.matrix.element_end(); z_i != z_end; ++z_i, ++x_i) {
00868 const N& d = *x_i;
00869 if (is_plus_infinity(d))
00870 continue;
00871 if (is_integer(d))
00872 assign_r(*z_i, d, ROUND_NOT_NEEDED);
00873 else {
00874 all_integers = false;
00875 Z& d_z = *z_i;
00876
00877 neg_assign_r(tmp, d, ROUND_NOT_NEEDED);
00878 assign_r(d_z, tmp, ROUND_UP);
00879 neg_assign_r(d_z, d_z, ROUND_NOT_NEEDED);
00880 }
00881 }
00882
00883 if (all_integers)
00884
00885 oct_z.set_strongly_closed();
00886 else {
00887
00888 oct_z.strong_closure_assign();
00889 if (oct_z.marked_empty())
00890 return false;
00891 }
00892 return !oct_z.tight_coherence_would_make_empty();
00893 }
00894
00895 template <typename T>
00896 bool
00897 Octagonal_Shape<T>::constrains(const Variable var) const {
00898
00899 const dimension_type var_space_dim = var.space_dimension();
00900 if (space_dimension() < var_space_dim)
00901 throw_dimension_incompatible("constrains(v)", "v", var);
00902
00903
00904
00905 if (marked_empty())
00906 return true;
00907
00908
00909 const dimension_type n_v = 2*(var_space_dim - 1);
00910 typename OR_Matrix<N>::const_row_iterator m_iter = matrix.row_begin() + n_v;
00911 typename OR_Matrix<N>::const_row_reference_type r_v = *m_iter;
00912 typename OR_Matrix<N>::const_row_reference_type r_cv = *(++m_iter);
00913 for (dimension_type h = m_iter.row_size(); h-- > 0; ) {
00914 if (!is_plus_infinity(r_v[h]) || !is_plus_infinity(r_cv[h]))
00915 return true;
00916 }
00917 ++m_iter;
00918 for (typename OR_Matrix<N>::const_row_iterator m_end = matrix.row_end();
00919 m_iter != m_end; ++m_iter) {
00920 typename OR_Matrix<N>::const_row_reference_type r = *m_iter;
00921 if (!is_plus_infinity(r[n_v]) || !is_plus_infinity(r[n_v+1]))
00922 return true;
00923 }
00924
00925
00926
00927 return is_empty();
00928 }
00929
00930 template <typename T>
00931 bool
00932 Octagonal_Shape<T>::is_strong_coherent() const {
00933
00934
00935 const dimension_type num_rows = matrix.num_rows();
00936
00937
00938 DIRTY_TEMP(N, semi_sum);
00939
00940
00941
00942
00943
00944 for (dimension_type i = num_rows; i-- > 0; ) {
00945 typename OR_Matrix<N>::const_row_iterator iter = matrix.row_begin() + i;
00946 typename OR_Matrix<N>::const_row_reference_type m_i = *iter;
00947 const N& m_i_ci = m_i[coherent_index(i)];
00948 for (dimension_type j = matrix.row_size(i); j-- > 0; )
00949
00950 if (i != j) {
00951 const N& m_cj_j = matrix[coherent_index(j)][j];
00952 if (!is_plus_infinity(m_i_ci)
00953 && !is_plus_infinity(m_cj_j)) {
00954
00955
00956 add_assign_r(semi_sum, m_i_ci, m_cj_j, ROUND_UP);
00957 div2exp_assign_r(semi_sum, semi_sum, 1, ROUND_UP);
00958 if (m_i[j] > semi_sum)
00959 return false;
00960 }
00961 }
00962 }
00963 return true;
00964 }
00965
00966 template <typename T>
00967 bool
00968 Octagonal_Shape<T>::is_strongly_reduced() const {
00969
00970
00971
00972 if (marked_empty())
00973 return true;
00974
00975 Octagonal_Shape x = *this;
00976
00977
00978 for (typename OR_Matrix<N>::const_row_iterator iter = matrix.row_begin(),
00979 matrix_row_end = matrix.row_end(); iter != matrix_row_end; ++iter) {
00980 typename OR_Matrix<N>::const_row_reference_type m_i = *iter;
00981 const dimension_type i = iter.index();
00982 for (dimension_type j = iter.row_size(); j-- > 0; ) {
00983 if (!is_plus_infinity(m_i[j])) {
00984 Octagonal_Shape x_copy = *this;
00985 assign_r(x_copy.matrix[i][j], PLUS_INFINITY, ROUND_NOT_NEEDED);
00986 if (x == x_copy)
00987 return false;
00988 }
00989 }
00990 }
00991
00992 return true;
00993 }
00994
00995 template <typename T>
00996 bool
00997 Octagonal_Shape<T>::bounds(const Linear_Expression& expr,
00998 const bool from_above) const {
00999
01000
01001 const dimension_type expr_space_dim = expr.space_dimension();
01002 if (space_dim < expr_space_dim)
01003 throw_dimension_incompatible((from_above
01004 ? "bounds_from_above(e)"
01005 : "bounds_from_below(e)"), "e", expr);
01006 strong_closure_assign();
01007
01008
01009 if (space_dim == 0 || marked_empty())
01010 return true;
01011
01012
01013
01014 const Constraint& c = (from_above) ? expr <= 0 : expr >= 0;
01015 dimension_type num_vars = 0;
01016 dimension_type i = 0;
01017 dimension_type j = 0;
01018 TEMP_INTEGER(coeff);
01019 TEMP_INTEGER(term);
01020 if (extract_octagonal_difference(c, c.space_dimension(), num_vars,
01021 i, j, coeff, term)) {
01022 if (num_vars == 0)
01023 return true;
01024
01025 typename OR_Matrix<N>::const_row_iterator i_iter = matrix.row_begin() + i;
01026 typename OR_Matrix<N>::const_row_reference_type m_i = *i_iter;
01027 return !is_plus_infinity(m_i[j]);
01028 }
01029 else {
01030
01031 Optimization_Mode mode_bounds =
01032 from_above ? MAXIMIZATION : MINIMIZATION;
01033 MIP_Problem mip(space_dim, constraints(), expr, mode_bounds);
01034 return (mip.solve() == OPTIMIZED_MIP_PROBLEM);
01035 }
01036 }
01037
01038 template <typename T>
01039 bool
01040 Octagonal_Shape<T>::max_min(const Linear_Expression& expr,
01041 const bool maximize,
01042 Coefficient& ext_n, Coefficient& ext_d,
01043 bool& included) const {
01044
01045
01046 const dimension_type expr_space_dim = expr.space_dimension();
01047 if (space_dim < expr_space_dim)
01048 throw_dimension_incompatible((maximize
01049 ? "maximize(e, ...)"
01050 : "minimize(e, ...)"), "e", expr);
01051
01052 if (space_dim == 0) {
01053 if (marked_empty())
01054 return false;
01055 else {
01056 ext_n = expr.inhomogeneous_term();
01057 ext_d = 1;
01058 included = true;
01059 return true;
01060 }
01061 }
01062
01063 strong_closure_assign();
01064
01065 if (marked_empty())
01066 return false;
01067
01068
01069
01070 const Constraint& c = (maximize) ? expr <= 0 : expr >= 0;
01071 dimension_type num_vars = 0;
01072 dimension_type i = 0;
01073 dimension_type j = 0;
01074 TEMP_INTEGER(coeff);
01075 TEMP_INTEGER(term);
01076 if (!extract_octagonal_difference(c, c.space_dimension(), num_vars,
01077 i, j, coeff, term)) {
01078
01079 Optimization_Mode max_min = (maximize) ? MAXIMIZATION : MINIMIZATION;
01080 MIP_Problem mip(space_dim, constraints(), expr, max_min);
01081 if (mip.solve() == OPTIMIZED_MIP_PROBLEM) {
01082 mip.optimal_value(ext_n, ext_d);
01083 included = true;
01084 return true;
01085 }
01086 else
01087
01088 return false;
01089 }
01090 else {
01091
01092 if (num_vars == 0) {
01093 ext_n = expr.inhomogeneous_term();
01094 ext_d = 1;
01095 included = true;
01096 return true;
01097 }
01098
01099
01100 typename OR_Matrix<N>::const_row_iterator i_iter = matrix.row_begin() + i;
01101 typename OR_Matrix<N>::const_row_reference_type m_i = *i_iter;
01102 DIRTY_TEMP(N, d);
01103 if (!is_plus_infinity(m_i[j])) {
01104 const Coefficient& b = expr.inhomogeneous_term();
01105 TEMP_INTEGER(minus_b);
01106 neg_assign(minus_b, b);
01107 const Coefficient& sc_b = maximize ? b : minus_b;
01108 assign_r(d, sc_b, ROUND_UP);
01109
01110
01111 DIRTY_TEMP(N, coeff_expr);
01112 const Coefficient& coeff_i = expr.coefficient(Variable(i/2));
01113 const int sign_i = sgn(coeff_i);
01114 if (sign_i > 0)
01115 assign_r(coeff_expr, coeff_i, ROUND_UP);
01116 else {
01117 TEMP_INTEGER(minus_coeff_i);
01118 neg_assign(minus_coeff_i, expr.coefficient(Variable(i/2)));
01119 assign_r(coeff_expr, minus_coeff_i, ROUND_UP);
01120 }
01121
01122 if (num_vars == 1) {
01123 DIRTY_TEMP(N, m_i_j);
01124 div2exp_assign_r(m_i_j, m_i[j], 1, ROUND_UP);
01125 add_mul_assign_r(d, coeff_expr, m_i_j, ROUND_UP);
01126 }
01127 else
01128 add_mul_assign_r(d, coeff_expr, m_i[j], ROUND_UP);
01129 numer_denom(d, ext_n, ext_d);
01130 if (!maximize)
01131 neg_assign(ext_n);
01132 included = true;
01133 return true;
01134 }
01135
01136
01137 return false;
01138 }
01139 }
01140
01141 template <typename T>
01142 bool
01143 Octagonal_Shape<T>::max_min(const Linear_Expression& expr,
01144 const bool maximize,
01145 Coefficient& ext_n, Coefficient& ext_d,
01146 bool& included, Generator& g) const {
01147
01148
01149 const dimension_type expr_space_dim = expr.space_dimension();
01150 if (space_dim < expr_space_dim)
01151 throw_dimension_incompatible((maximize
01152 ? "maximize(e, ...)"
01153 : "minimize(e, ...)"), "e", expr);
01154
01155 if (space_dim == 0) {
01156 if (marked_empty())
01157 return false;
01158 else {
01159 ext_n = expr.inhomogeneous_term();
01160 ext_d = 1;
01161 included = true;
01162 g = point();
01163 return true;
01164 }
01165 }
01166
01167 strong_closure_assign();
01168
01169 if (marked_empty())
01170 return false;
01171 if (!is_universe()) {
01172
01173
01174 Optimization_Mode max_min = (maximize) ? MAXIMIZATION : MINIMIZATION;
01175 MIP_Problem mip(space_dim, constraints(), expr, max_min);
01176 if (mip.solve() == OPTIMIZED_MIP_PROBLEM) {
01177 g = mip.optimizing_point();
01178 mip.evaluate_objective_function(g, ext_n, ext_d);
01179 included = true;
01180 return true;
01181 }
01182 }
01183
01184 return false;
01185 }
01186
01187 template <typename T>
01188 Poly_Con_Relation
01189 Octagonal_Shape<T>::relation_with(const Congruence& cg) const {
01190 dimension_type cg_space_dim = cg.space_dimension();
01191
01192
01193 if (cg_space_dim > space_dim) {
01194 throw_dimension_incompatible("relation_with(cg)", cg);
01195 }
01196
01197
01198
01199 if (cg.is_equality()) {
01200 Constraint c(cg);
01201 return relation_with(c);
01202 }
01203
01204 strong_closure_assign();
01205
01206 if (marked_empty())
01207 return Poly_Con_Relation::saturates()
01208 && Poly_Con_Relation::is_included()
01209 && Poly_Con_Relation::is_disjoint();
01210
01211 if (space_dim == 0) {
01212 if (cg.is_inconsistent())
01213 return Poly_Con_Relation::is_disjoint();
01214 else if (cg.inhomogeneous_term() % cg.modulus() == 0)
01215 return Poly_Con_Relation::saturates()
01216 && Poly_Con_Relation::is_included();
01217 }
01218
01219 DIRTY_TEMP(Coefficient, min_num);
01220 DIRTY_TEMP(Coefficient, min_den);
01221 bool min_included;
01222 TEMP_INTEGER(mod);
01223 mod = cg.modulus();
01224 Linear_Expression le;
01225 for (dimension_type i = cg_space_dim; i-- > 0; )
01226 le += cg.coefficient(Variable(i)) * Variable(i);
01227 bool bounded_below = minimize(le, min_num, min_den, min_included);
01228
01229 if (!bounded_below)
01230 return Poly_Con_Relation::strictly_intersects();
01231
01232 TEMP_INTEGER(v);
01233 TEMP_INTEGER(lower_num);
01234 TEMP_INTEGER(lower_den);
01235 TEMP_INTEGER(lower);
01236 assign_r(lower_num, min_num, ROUND_NOT_NEEDED);
01237 assign_r(lower_den, min_den, ROUND_NOT_NEEDED);
01238 neg_assign(v, cg.inhomogeneous_term());
01239 lower = lower_num / lower_den;
01240 v += ((lower / mod) * mod);
01241 if (v * lower_den < lower_num)
01242 v += mod;
01243 const Constraint& c(le == v);
01244 return relation_with(c);
01245 }
01246
01247 template <typename T>
01248 Poly_Con_Relation
01249 Octagonal_Shape<T>::relation_with(const Constraint& c) const {
01250 dimension_type c_space_dim = c.space_dimension();
01251
01252
01253 if (c_space_dim > space_dim)
01254 throw_dimension_incompatible("relation_with(c)", c);
01255
01256
01257 strong_closure_assign();
01258
01259 if (marked_empty())
01260 return Poly_Con_Relation::saturates()
01261 && Poly_Con_Relation::is_included()
01262 && Poly_Con_Relation::is_disjoint();
01263
01264 if (space_dim == 0) {
01265
01266 if ((c.is_equality() && c.inhomogeneous_term() != 0)
01267 || (c.is_inequality() && c.inhomogeneous_term() < 0))
01268 return Poly_Con_Relation::is_disjoint();
01269 else if (c.is_strict_inequality() && c.inhomogeneous_term() == 0)
01270
01271
01272 return Poly_Con_Relation::saturates()
01273 && Poly_Con_Relation::is_disjoint();
01274
01275
01276 else if (c.is_equality() || c.inhomogeneous_term() == 0)
01277 return Poly_Con_Relation::saturates()
01278 && Poly_Con_Relation::is_included();
01279 else
01280
01281
01282
01283 return Poly_Con_Relation::is_included();
01284 }
01285
01286 dimension_type num_vars = 0;
01287 dimension_type i = 0;
01288 dimension_type j = 0;
01289 TEMP_INTEGER(coeff);
01290 TEMP_INTEGER(c_term);
01291 if (!extract_octagonal_difference(c, c_space_dim, num_vars,
01292 i, j, coeff, c_term)) {
01293
01294
01295
01296
01297
01298
01299 Linear_Expression le;
01300 for (dimension_type k = c_space_dim; k-- > 0; ) {
01301 Variable vk(k);
01302 le += c.coefficient(vk) * vk;
01303 }
01304 DIRTY_TEMP(Coefficient, max_num);
01305 DIRTY_TEMP(Coefficient, max_den);
01306 bool max_included;
01307 DIRTY_TEMP(Coefficient, min_num);
01308 DIRTY_TEMP(Coefficient, min_den);
01309 bool min_included;
01310 bool bounded_above = maximize(le, max_num, max_den, max_included);
01311 bool bounded_below = minimize(le, min_num, min_den, min_included);
01312 if (!bounded_above) {
01313 if (!bounded_below)
01314 return Poly_Con_Relation::strictly_intersects();
01315 min_num += c.inhomogeneous_term() * min_den;
01316 switch (sgn(min_num)) {
01317 case 1:
01318 if (c.is_equality())
01319 return Poly_Con_Relation::is_disjoint();
01320 return Poly_Con_Relation::is_included();
01321 case 0:
01322 if (c.is_strict_inequality() || c.is_equality())
01323 return Poly_Con_Relation::strictly_intersects();
01324 return Poly_Con_Relation::is_included();
01325 case -1:
01326 return Poly_Con_Relation::strictly_intersects();
01327 }
01328 }
01329 if (!bounded_below) {
01330 max_num += c.inhomogeneous_term() * max_den;
01331 switch (sgn(max_num)) {
01332 case 1:
01333 return Poly_Con_Relation::strictly_intersects();
01334 case 0:
01335 if (c.is_strict_inequality())
01336 return Poly_Con_Relation::is_disjoint();
01337 return Poly_Con_Relation::strictly_intersects();
01338 case -1:
01339 return Poly_Con_Relation::is_disjoint();
01340 }
01341 }
01342 else {
01343 max_num += c.inhomogeneous_term() * max_den;
01344 min_num += c.inhomogeneous_term() * min_den;
01345 switch (sgn(max_num)) {
01346 case 1:
01347 switch (sgn(min_num)) {
01348 case 1:
01349 if (c.is_equality())
01350 return Poly_Con_Relation::is_disjoint();
01351 return Poly_Con_Relation::is_included();
01352 case 0:
01353 if (c.is_equality())
01354 return Poly_Con_Relation::strictly_intersects();
01355 if (c.is_strict_inequality())
01356 return Poly_Con_Relation::strictly_intersects();
01357 return Poly_Con_Relation::is_included();
01358 case -1:
01359 return Poly_Con_Relation::strictly_intersects();
01360 }
01361 case 0:
01362 if (min_num == 0) {
01363 if (c.is_strict_inequality())
01364 return Poly_Con_Relation::is_disjoint()
01365 && Poly_Con_Relation::saturates();
01366 return Poly_Con_Relation::is_included()
01367 && Poly_Con_Relation::saturates();
01368 }
01369 if (c.is_strict_inequality())
01370 return Poly_Con_Relation::is_disjoint();
01371 return Poly_Con_Relation::strictly_intersects();
01372 case -1:
01373 return Poly_Con_Relation::is_disjoint();
01374 }
01375 }
01376 }
01377
01378 if (num_vars == 0) {
01379
01380 switch (sgn(c.inhomogeneous_term())) {
01381 case -1:
01382 return Poly_Con_Relation::is_disjoint();
01383 case 0:
01384 if (c.is_strict_inequality())
01385 return Poly_Con_Relation::saturates()
01386 && Poly_Con_Relation::is_disjoint();
01387 else
01388 return Poly_Con_Relation::saturates()
01389 && Poly_Con_Relation::is_included();
01390 case 1:
01391 if (c.is_equality())
01392 return Poly_Con_Relation::is_disjoint();
01393 else
01394 return Poly_Con_Relation::is_included();
01395 }
01396 }
01397
01398
01399 typename OR_Matrix<N>::const_row_iterator i_iter = matrix.row_begin() + i;
01400 typename OR_Matrix<N>::const_row_reference_type m_i = *i_iter;
01401 const N& m_i_j = m_i[j];
01402
01403 if (coeff < 0)
01404 neg_assign(coeff);
01405
01406
01407
01408 if (i%2 == 0)
01409 ++i_iter;
01410 else
01411 --i_iter;
01412 typename OR_Matrix<N>::const_row_reference_type m_ci = *i_iter;
01413 const N& m_ci_cj = m_ci[coherent_index(j)];
01414 TEMP_INTEGER(numer);
01415 TEMP_INTEGER(denom);
01416
01417
01418 DIRTY_TEMP0(mpq_class, q_x);
01419 DIRTY_TEMP0(mpq_class, q_y);
01420 DIRTY_TEMP0(mpq_class, d);
01421 DIRTY_TEMP0(mpq_class, d1);
01422 DIRTY_TEMP0(mpq_class, c_den);
01423 DIRTY_TEMP0(mpq_class, q_den);
01424 assign_r(c_den, coeff, ROUND_NOT_NEEDED);
01425 assign_r(d, c_term, ROUND_NOT_NEEDED);
01426 neg_assign_r(d1, d, ROUND_NOT_NEEDED);
01427 div_assign_r(d, d, c_den, ROUND_NOT_NEEDED);
01428 div_assign_r(d1, d1, c_den, ROUND_NOT_NEEDED);
01429
01430 if (is_plus_infinity(m_i_j)) {
01431 if (!is_plus_infinity(m_ci_cj)) {
01432
01433
01434
01435
01436
01437 numer_denom(m_ci_cj, numer, denom);
01438 assign_r(q_den, denom, ROUND_NOT_NEEDED);
01439 assign_r(q_y, numer, ROUND_NOT_NEEDED);
01440 div_assign_r(q_y, q_y, q_den, ROUND_NOT_NEEDED);
01441 if (q_y < d1)
01442 return Poly_Con_Relation::is_disjoint();
01443 if (q_y == d1 && c.is_strict_inequality())
01444 return Poly_Con_Relation::is_disjoint();
01445 }
01446
01447
01448 return Poly_Con_Relation::strictly_intersects();
01449 }
01450
01451
01452 numer_denom(m_i_j, numer, denom);
01453 assign_r(q_den, denom, ROUND_NOT_NEEDED);
01454 assign_r(q_x, numer, ROUND_NOT_NEEDED);
01455 div_assign_r(q_x, q_x, q_den, ROUND_NOT_NEEDED);
01456
01457 if (!is_plus_infinity(m_ci_cj)) {
01458 numer_denom(m_ci_cj, numer, denom);
01459 assign_r(q_den, denom, ROUND_NOT_NEEDED);
01460 assign_r(q_y, numer, ROUND_NOT_NEEDED);
01461 div_assign_r(q_y, q_y, q_den, ROUND_NOT_NEEDED);
01462 if (q_x == d && q_y == d1) {
01463 if (c.is_strict_inequality())
01464 return Poly_Con_Relation::saturates()
01465 && Poly_Con_Relation::is_disjoint();
01466 else
01467 return Poly_Con_Relation::saturates()
01468 && Poly_Con_Relation::is_included();
01469 }
01470
01471
01472 if (q_y < d1)
01473 return Poly_Con_Relation::is_disjoint();
01474 if (q_y == d1 && c.is_strict_inequality())
01475 return Poly_Con_Relation::is_disjoint();
01476 }
01477
01478
01479
01480
01481 if (d > q_x) {
01482 if (c.is_equality())
01483 return Poly_Con_Relation::is_disjoint();
01484 else
01485 return Poly_Con_Relation::is_included();
01486 }
01487
01488 if (d == q_x && c.is_nonstrict_inequality())
01489 return Poly_Con_Relation::is_included();
01490
01491
01492 return Poly_Con_Relation::strictly_intersects();
01493 }
01494
01495 template <typename T>
01496 Poly_Gen_Relation
01497 Octagonal_Shape<T>::relation_with(const Generator& g) const {
01498 const dimension_type g_space_dim = g.space_dimension();
01499
01500
01501 if (space_dim < g_space_dim)
01502 throw_dimension_incompatible("relation_with(g)", g);
01503
01504
01505
01506 strong_closure_assign();
01507
01508
01509 if (marked_empty())
01510 return Poly_Gen_Relation::nothing();
01511
01512
01513
01514 if (space_dim == 0)
01515 return Poly_Gen_Relation::subsumes();
01516
01517 const bool is_line = g.is_line();
01518 const bool is_line_or_ray = g.is_line_or_ray();
01519
01520
01521
01522
01523
01524
01525
01526 typedef typename OR_Matrix<N>::const_row_iterator Row_Iterator;
01527 typedef typename OR_Matrix<N>::const_row_reference_type Row_Reference;
01528
01529 const Row_Iterator m_begin = matrix.row_begin();
01530 const Row_Iterator m_end = matrix.row_end();
01531
01532 TEMP_INTEGER(num);
01533 TEMP_INTEGER(den);
01534 TEMP_INTEGER(product);
01535
01536
01537 for (Row_Iterator i_iter = m_begin; i_iter != m_end; i_iter += 2) {
01538 dimension_type i = i_iter.index();
01539 Row_Reference m_i = *i_iter;
01540 Row_Reference m_ii = *(i_iter+1);
01541 const N& m_i_ii = m_i[i+1];
01542 const N& m_ii_i = m_ii[i];
01543
01544 const Variable x(i/2);
01545 const Coefficient& g_coeff_x = (x.space_dimension() > g_space_dim)
01546 ? Coefficient(0) : g.coefficient(x);
01547 if (is_additive_inverse(m_i_ii, m_ii_i)) {
01548
01549
01550
01551
01552 numer_denom(m_ii_i, num, den);
01553 den *= 2;
01554 product = den * g_coeff_x;
01555
01556
01557 if (!is_line_or_ray) {
01558 neg_assign(num);
01559 add_mul_assign(product, num, g.divisor());
01560 }
01561 if (product != 0)
01562 return Poly_Gen_Relation::nothing();
01563 }
01564
01565 else {
01566 if (!is_plus_infinity(m_i_ii)) {
01567
01568
01569
01570 numer_denom(m_i_ii, num, den);
01571 den *= -2;
01572 product = den * g_coeff_x;
01573
01574
01575 if (!is_line_or_ray) {
01576 neg_assign(num);
01577 add_mul_assign(product, num, g.divisor());
01578 }
01579 if (is_line && product != 0)
01580 return Poly_Gen_Relation::nothing();
01581 else
01582
01583
01584
01585
01586 if (product > 0)
01587 return Poly_Gen_Relation::nothing();
01588 }
01589 if (!is_plus_infinity(m_ii_i)) {
01590
01591 numer_denom(m_ii_i, num, den);
01592 den *= 2;
01593 product = den * g_coeff_x;
01594
01595
01596 if (!is_line_or_ray) {
01597 neg_assign(num);
01598 add_mul_assign(product, num , g.divisor());
01599 }
01600 if (is_line && product != 0)
01601 return Poly_Gen_Relation::nothing();
01602 else
01603
01604
01605
01606
01607 if (product > 0)
01608 return Poly_Gen_Relation::nothing();
01609 }
01610 }
01611 }
01612
01613
01614 for (Row_Iterator i_iter = m_begin ; i_iter != m_end; i_iter += 2) {
01615 dimension_type i = i_iter.index();
01616 Row_Reference m_i = *i_iter;
01617 Row_Reference m_ii = *(i_iter+1);
01618 for (dimension_type j = 0; j < i; j += 2) {
01619 const N& m_i_j = m_i[j];
01620 const N& m_ii_jj = m_ii[j+1];
01621 const N& m_ii_j = m_ii[j];
01622 const N& m_i_jj = m_i[j+1];
01623 const Variable x(j/2);
01624 const Variable y(i/2);
01625 const Coefficient& g_coeff_x = (x.space_dimension() > g_space_dim)
01626 ? Coefficient(0) : g.coefficient(x);
01627 const Coefficient& g_coeff_y = (y.space_dimension() > g_space_dim)
01628 ? Coefficient(0) : g.coefficient(y);
01629
01630
01631 const bool is_binary_equality = is_additive_inverse(m_ii_jj, m_i_j);
01632 const bool is_a_binary_equality = is_additive_inverse(m_i_jj, m_ii_j);
01633 if (is_binary_equality) {
01634
01635
01636
01637
01638
01639 numer_denom(m_i_j, num, den);
01640 product = den * g_coeff_x;
01641 neg_assign(den);
01642 add_mul_assign(product, den, g_coeff_y);
01643
01644
01645 if (!is_line_or_ray) {
01646 neg_assign(num);
01647 add_mul_assign(product, num, g.divisor());
01648 }
01649 if (product != 0)
01650 return Poly_Gen_Relation::nothing();
01651 }
01652 else {
01653 if (!is_plus_infinity(m_i_j)) {
01654
01655
01656
01657
01658
01659 numer_denom(m_i_j, num, den);
01660 product = den * g_coeff_x;
01661 neg_assign(den);
01662 add_mul_assign(product, den, g_coeff_y);
01663
01664
01665 if (!is_line_or_ray) {
01666 neg_assign(num);
01667 add_mul_assign(product, num, g.divisor());
01668 }
01669 if (is_line && product != 0)
01670 return Poly_Gen_Relation::nothing();
01671 else if (product > 0)
01672 return Poly_Gen_Relation::nothing();
01673 }
01674 if (!is_plus_infinity(m_ii_jj)) {
01675
01676
01677
01678
01679
01680 numer_denom(m_ii_jj, num, den);
01681 product = den * g_coeff_y;
01682 neg_assign(den);
01683 add_mul_assign(product, den, g_coeff_x);
01684
01685
01686 if (!is_line_or_ray) {
01687 neg_assign(num);
01688 add_mul_assign(product, num, g.divisor());
01689 }
01690 if (is_line && product != 0)
01691 return Poly_Gen_Relation::nothing();
01692 else if (product > 0)
01693 return Poly_Gen_Relation::nothing();
01694 }
01695 }
01696
01697 if (is_a_binary_equality) {
01698
01699
01700
01701
01702
01703 numer_denom(m_ii_j, num, den);
01704 product = den * g_coeff_x;
01705 add_mul_assign(product, den, g_coeff_y);
01706
01707
01708 if (!is_line_or_ray) {
01709 neg_assign(num);
01710 add_mul_assign(product, num, g.divisor());
01711 }
01712 if (product != 0)
01713 return Poly_Gen_Relation::nothing();
01714 }
01715 else {
01716 if (!is_plus_infinity(m_i_jj)) {
01717
01718
01719
01720
01721
01722 numer_denom(m_i_jj, num, den);
01723 neg_assign(den);
01724 product = den * g_coeff_x;
01725 add_mul_assign(product, den, g_coeff_y);
01726
01727
01728 if (!is_line_or_ray) {
01729 neg_assign(num);
01730 add_mul_assign(product, num, g.divisor());
01731 }
01732 if (is_line && product != 0)
01733 return Poly_Gen_Relation::nothing();
01734 else if (product > 0)
01735 return Poly_Gen_Relation::nothing();
01736 }
01737 if (!is_plus_infinity(m_ii_j)) {
01738
01739
01740
01741
01742
01743 numer_denom(m_ii_j, num, den);
01744 product = den * g_coeff_x;
01745 add_mul_assign(product, den, g_coeff_y);
01746
01747
01748 if (!is_line_or_ray) {
01749 neg_assign(num);
01750 add_mul_assign(product, num, g.divisor());
01751 }
01752 if (is_line && product != 0)
01753 return Poly_Gen_Relation::nothing();
01754 else if (product > 0)
01755 return Poly_Gen_Relation::nothing();
01756 }
01757 }
01758 }
01759 }
01760
01761
01762 return Poly_Gen_Relation::subsumes();
01763 }
01764
01765 template <typename T>
01766 void
01767 Octagonal_Shape<T>::strong_closure_assign() const {
01768
01769 if (marked_empty() || marked_strongly_closed() || space_dim == 0)
01770 return;
01771
01772
01773
01774 Octagonal_Shape& x = const_cast<Octagonal_Shape<T>&>(*this);
01775
01776 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
01777 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
01778
01779 const dimension_type n_rows = x.matrix.num_rows();
01780 const Row_Iterator m_begin = x.matrix.row_begin();
01781 const Row_Iterator m_end = x.matrix.row_end();
01782
01783
01784 for (Row_Iterator i = m_begin; i != m_end; ++i) {
01785 assert(is_plus_infinity((*i)[i.index()]));
01786 assign_r((*i)[i.index()], 0, ROUND_NOT_NEEDED);
01787 }
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800 typename OR_Matrix<N>::element_iterator iter_ij;
01801 std::vector<N> vec_k(n_rows);
01802 std::vector<N> vec_ck(n_rows);
01803 DIRTY_TEMP(N, sum1);
01804 DIRTY_TEMP(N, sum2);
01805 Row_Reference x_k;
01806 Row_Reference x_ck;
01807 Row_Reference x_i;
01808 Row_Reference x_ci;
01809
01810
01811
01812 for (int twice = 0; twice < 2; ++twice) {
01813
01814 Row_Iterator x_k_iter = m_begin;
01815 Row_Iterator x_i_iter = m_begin;
01816 for (dimension_type k = 0; k < n_rows; k += 2) {
01817 const dimension_type ck = k+1;
01818
01819 iter_ij = x.matrix.element_begin();
01820
01821 x_k = *x_k_iter;
01822 ++x_k_iter;
01823 x_ck = *x_k_iter;
01824 ++x_k_iter;
01825
01826 for (dimension_type i = 0; i <= k; i += 2) {
01827 const dimension_type ci = i+1;
01828
01829 vec_k[i] = x_k[i];
01830
01831 vec_k[ci] = x_k[ci];
01832
01833 vec_ck[i] = x_ck[i];
01834
01835 vec_ck[ci] = x_ck[ci];
01836 }
01837 x_i_iter = x_k_iter;
01838 for (dimension_type i = k+2; i < n_rows; i += 2) {
01839 const dimension_type ci = i+1;
01840 x_i = *x_i_iter;
01841 ++x_i_iter;
01842 x_ci = *x_i_iter;
01843 ++x_i_iter;
01844
01845 vec_k[i] = x_ci[ck];
01846
01847 vec_k[ci] = x_i[ck];
01848
01849 vec_ck[i] = x_ci[k];
01850
01851 vec_ck[ci] = x_i[k];
01852 }
01853
01854 for (dimension_type i = 0; i < n_rows; ++i) {
01855 const dimension_type ci = coherent_index(i);
01856 const N& vec_k_ci = vec_k[ci];
01857 const N& vec_ck_ci = vec_ck[ci];
01858
01859
01860 for (dimension_type j = 0; j <= i; ) {
01861
01862
01863
01864 add_assign_r(sum1, vec_ck_ci, vec_k[j], ROUND_UP);
01865 add_assign_r(sum2, vec_k_ci, vec_ck[j], ROUND_UP);
01866 min_assign(sum1, sum2);
01867 min_assign(*iter_ij, sum1);
01868
01869 ++j;
01870 ++iter_ij;
01871
01872 add_assign_r(sum1, vec_ck_ci, vec_k[j], ROUND_UP);
01873 add_assign_r(sum2, vec_k_ci, vec_ck[j], ROUND_UP);
01874 min_assign(sum1, sum2);
01875 min_assign(*iter_ij, sum1);
01876
01877 ++j;
01878 ++iter_ij;
01879 }
01880 }
01881 }
01882 }
01883
01884
01885
01886 for (Row_Iterator i = m_begin; i != m_end; ++i) {
01887 N& x_i_i = (*i)[i.index()];
01888 if (sgn(x_i_i) < 0) {
01889 x.set_empty();
01890 return;
01891 }
01892 else {
01893 assert(sgn(x_i_i) == 0);
01894
01895 assign_r(x_i_i, PLUS_INFINITY, ROUND_NOT_NEEDED);
01896 }
01897 }
01898
01899
01900 x.strong_coherence_assign();
01901
01902 x.set_strongly_closed();
01903 }
01904
01905 template <typename T>
01906 void
01907 Octagonal_Shape<T>::strong_coherence_assign() {
01908
01909
01910
01911
01912
01913 DIRTY_TEMP(N, semi_sum);
01914 for (typename OR_Matrix<N>::row_iterator i_iter = matrix.row_begin(),
01915 i_end = matrix.row_end(); i_iter != i_end; ++i_iter) {
01916 typename OR_Matrix<N>::row_reference_type x_i = *i_iter;
01917 const dimension_type i = i_iter.index();
01918 const N& x_i_ci = x_i[coherent_index(i)];
01919
01920 if (!is_plus_infinity(x_i_ci))
01921 for (dimension_type j = 0, rs_i = i_iter.row_size(); j < rs_i; ++j)
01922 if (i != j) {
01923 const N& x_cj_j = matrix[coherent_index(j)][j];
01924 if (!is_plus_infinity(x_cj_j)) {
01925 add_assign_r(semi_sum, x_i_ci, x_cj_j, ROUND_UP);
01926 div2exp_assign_r(semi_sum, semi_sum, 1, ROUND_UP);
01927 min_assign(x_i[j], semi_sum);
01928 }
01929 }
01930 }
01931 }
01932
01933 template <typename T>
01934 bool
01935 Octagonal_Shape<T>::tight_coherence_would_make_empty() const {
01936 assert(std::numeric_limits<N>::is_integer);
01937 assert(marked_strongly_closed());
01938 const dimension_type space_dim = space_dimension();
01939 for (dimension_type i = 0; i < 2*space_dim; i += 2) {
01940 const dimension_type ci = i+1;
01941 const N& mat_i_ci = matrix[i][ci];
01942 if (!is_plus_infinity(mat_i_ci)
01943
01944 && !is_even(mat_i_ci)
01945
01946 && is_additive_inverse(mat_i_ci, matrix[ci][i]))
01947 return true;
01948 }
01949 return false;
01950 }
01951
01952 template <typename T>
01953 void
01954 Octagonal_Shape<T>
01955 ::incremental_strong_closure_assign(const Variable var) const {
01956
01957 if (var.id() >= space_dim)
01958 throw_dimension_incompatible("incremental_strong_closure_assign(v)",
01959 var.id());
01960
01961
01962 if (marked_empty() || marked_strongly_closed())
01963 return;
01964
01965
01966 if (space_dim == 0)
01967 return;
01968
01969 Octagonal_Shape& x = const_cast<Octagonal_Shape<T>&>(*this);
01970
01971 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
01972 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
01973
01974 const Row_Iterator m_begin = x.matrix.row_begin();
01975 const Row_Iterator m_end = x.matrix.row_end();
01976
01977
01978 for (Row_Iterator i = m_begin; i != m_end; ++i) {
01979 assert(is_plus_infinity((*i)[i.index()]));
01980 assign_r((*i)[i.index()], 0, ROUND_NOT_NEEDED);
01981 }
01982
01983
01984
01985 const dimension_type v = 2*var.id();
01986 const dimension_type cv = v+1;
01987 Row_Iterator v_iter = m_begin + v;
01988 Row_Iterator cv_iter = v_iter + 1;
01989 Row_Reference x_v = *v_iter;
01990 Row_Reference x_cv = *cv_iter;
01991 const dimension_type rs_v = v_iter.row_size();
01992 const dimension_type n_rows = x.matrix.num_rows();
01993 DIRTY_TEMP(N, sum);
01994 for (Row_Iterator k_iter = m_begin; k_iter != m_end; ++k_iter) {
01995 const dimension_type k = k_iter.index();
01996 const dimension_type ck = coherent_index(k);
01997 const dimension_type rs_k = k_iter.row_size();
01998 Row_Reference x_k = *k_iter;
01999 Row_Reference x_ck = (k%2) ? *(k_iter-1) : *(k_iter+1);
02000
02001 for (Row_Iterator i_iter = m_begin; i_iter != m_end; ++i_iter) {
02002 const dimension_type i = i_iter.index();
02003 const dimension_type ci = coherent_index(i);
02004 const dimension_type rs_i = i_iter.row_size();
02005 Row_Reference x_i = *i_iter;
02006 Row_Reference x_ci = (i%2) ? *(i_iter-1) : *(i_iter+1);
02007
02008 const N& x_i_k = (k < rs_i) ? x_i[k] : x_ck[ci];
02009 if (!is_plus_infinity(x_i_k)) {
02010 const N& x_k_v = (v < rs_k) ? x_k[v] : x_cv[ck];
02011 if (!is_plus_infinity(x_k_v)) {
02012 add_assign_r(sum, x_i_k, x_k_v, ROUND_UP);
02013 N& x_i_v = (v < rs_i) ? x_i[v] : x_cv[ci];
02014 min_assign(x_i_v, sum);
02015 }
02016 const N& x_k_cv = (cv < rs_k) ? x_k[cv] : x_v[ck];
02017 if (!is_plus_infinity(x_k_cv)) {
02018 add_assign_r(sum, x_i_k, x_k_cv, ROUND_UP);
02019 N& x_i_cv = (cv < rs_i) ? x_i[cv] : x_v[ci];
02020 min_assign(x_i_cv, sum);
02021 }
02022 }
02023 const N& x_k_i = (i < rs_k) ? x_k[i] : x_ci[ck];
02024 if (!is_plus_infinity(x_k_i)) {
02025 const N& x_v_k = (k < rs_v) ? x_v[k] : x_ck[cv];
02026 if (!is_plus_infinity(x_v_k)) {
02027 N& x_v_i = (i < rs_v) ? x_v[i] : x_ci[cv];
02028 add_assign_r(sum, x_v_k, x_k_i, ROUND_UP);
02029 min_assign(x_v_i, sum);
02030 }
02031 const N& x_cv_k = (k < rs_v) ? x_cv[k] : x_ck[v];
02032 if (!is_plus_infinity(x_cv_k)) {
02033 N& x_cv_i = (i < rs_v) ? x_cv[i] : x_ci[v];
02034 add_assign_r(sum, x_cv_k, x_k_i, ROUND_UP);
02035 min_assign(x_cv_i, sum);
02036 }
02037 }
02038
02039 }
02040 }
02041
02042
02043
02044 for (Row_Iterator i_iter = m_begin; i_iter != m_end; ++i_iter) {
02045 const dimension_type i = i_iter.index();
02046 const dimension_type ci = coherent_index(i);
02047 const dimension_type rs_i = i_iter.row_size();
02048 Row_Reference x_i = *i_iter;
02049 Row_Reference x_ci = (i%2) ? *(i_iter-1) : *(i_iter+1);
02050 const N& x_i_v = (v < rs_i) ? x_i[v] : x_cv[ci];
02051
02052
02053
02054 for (dimension_type j = 0; j < n_rows; ++j) {
02055 const dimension_type cj = coherent_index(j);
02056 Row_Reference x_cj = *(m_begin+cj);
02057 N& x_i_j = (j < rs_i) ? x_i[j] : x_cj[ci];
02058 if (!is_plus_infinity(x_i_v)) {
02059 const N& x_v_j = (j < rs_v) ? x_v[j] : x_cj[cv];
02060 if (!is_plus_infinity(x_v_j)) {
02061 add_assign_r(sum, x_i_v, x_v_j, ROUND_UP);
02062 min_assign(x_i_j, sum);
02063 }
02064 }
02065 const N& x_i_cv = (cv < rs_i) ? x_i[cv] : x_v[ci];
02066 if (!is_plus_infinity(x_i_cv)) {
02067 const N& x_cv_j = (j < rs_v) ? x_cv[j] : x_cj[v];
02068 if (!is_plus_infinity(x_cv_j)) {
02069 add_assign_r(sum, x_i_cv, x_cv_j, ROUND_UP);
02070 min_assign(x_i_j, sum);
02071 }
02072 }
02073 }
02074 }
02075
02076
02077
02078 for (Row_Iterator i = m_begin; i != m_end; ++i) {
02079 N& x_i_i = (*i)[i.index()];
02080 if (sgn(x_i_i) < 0) {
02081 x.set_empty();
02082 return;
02083 }
02084 else {
02085
02086 assert(sgn(x_i_i) == 0);
02087 assign_r(x_i_i, PLUS_INFINITY, ROUND_NOT_NEEDED);
02088 }
02089 }
02090
02091
02092 x.strong_coherence_assign();
02093
02094 x.set_strongly_closed();
02095 }
02096
02097 template <typename T>
02098 void
02099 Octagonal_Shape<T>
02100 ::compute_successors(std::vector<dimension_type>& successor) const {
02101 assert(!marked_empty() && marked_strongly_closed());
02102 assert(successor.size() == 0);
02103
02104
02105
02106 const dimension_type successor_size = matrix.num_rows();
02107
02108 successor.reserve(successor_size);
02109 for (dimension_type i = 0; i < successor_size; ++i)
02110 successor.push_back(i);
02111
02112 for (dimension_type i = successor_size; i-- > 0; ) {
02113 typename OR_Matrix<N>::const_row_iterator i_iter = matrix.row_begin()+i;
02114 typename OR_Matrix<N>::const_row_reference_type m_i = *i_iter;
02115 typename OR_Matrix<N>::const_row_reference_type m_ci = (i%2) ?
02116 *(i_iter-1) : *(i_iter+1);
02117 for (dimension_type j = 0; j < i; ++j) {
02118
02119 dimension_type cj = coherent_index(j);
02120 if (is_additive_inverse(m_ci[cj], m_i[j]))
02121
02122 successor[j] = i;
02123 }
02124 }
02125 }
02126
02127 template <typename T>
02128 void
02129 Octagonal_Shape<T>
02130 ::compute_leaders(std::vector<dimension_type>& leaders) const {
02131 assert(!marked_empty() && marked_strongly_closed());
02132 assert(leaders.size() == 0);
02133
02134
02135
02136 const dimension_type leader_size = matrix.num_rows();
02137
02138 leaders.reserve(leader_size);
02139 for (dimension_type i = 0; i < leader_size; ++i)
02140 leaders.push_back(i);
02141
02142 for (typename OR_Matrix<N>::const_row_iterator i_iter = matrix.row_begin(),
02143 matrix_row_end = matrix.row_end();
02144 i_iter != matrix_row_end; ++i_iter) {
02145 typename OR_Matrix<N>::const_row_reference_type m_i = *i_iter;
02146 dimension_type i = i_iter.index();
02147 typename OR_Matrix<N>::const_row_reference_type m_ci =
02148 (i%2) ? *(i_iter-1) : *(i_iter+1);
02149 for (dimension_type j = 0; j < i; ++j) {
02150 dimension_type cj = coherent_index(j);
02151 if (is_additive_inverse(m_ci[cj], m_i[j]))
02152
02153 leaders[i] = leaders[j];
02154 }
02155 }
02156 }
02157
02158 template <typename T>
02159 void
02160 Octagonal_Shape<T>
02161 ::compute_leaders(std::vector<dimension_type>& successor,
02162 std::vector<dimension_type>& no_sing_leaders,
02163 bool& exist_sing_class,
02164 dimension_type& sing_leader) const {
02165 assert(!marked_empty() && marked_strongly_closed());
02166 assert(no_sing_leaders.size() == 0);
02167 dimension_type successor_size = successor.size();
02168 std::deque<bool> dealt_with(successor_size, false);
02169 for (dimension_type i = 0; i < successor_size; ++i) {
02170 dimension_type next_i = successor[i];
02171 if (!dealt_with[i]) {
02172
02173
02174 if (next_i == coherent_index(i)) {
02175 exist_sing_class = true;
02176 sing_leader = i;
02177 }
02178 else
02179 no_sing_leaders.push_back(i);
02180 }
02181
02182 dealt_with[next_i] = true;
02183 }
02184 }
02185
02186 template <typename T>
02187 void
02188 Octagonal_Shape<T>::strong_reduction_assign() const {
02189
02190 if (space_dim == 0)
02191 return;
02192
02193
02194 strong_closure_assign();
02195
02196
02197 if (marked_empty())
02198 return;
02199
02200
02201
02202
02203
02204 std::vector<dimension_type> no_sing_leaders;
02205 dimension_type sing_leader = 0;
02206 bool exist_sing_class = false;
02207 std::vector<dimension_type> successor;
02208 compute_successors(successor);
02209 compute_leaders(successor, no_sing_leaders, exist_sing_class, sing_leader);
02210 const dimension_type num_no_sing_leaders = no_sing_leaders.size();
02211
02212 Octagonal_Shape aux(space_dim);
02213
02214
02215
02216 for (dimension_type li = 0; li < num_no_sing_leaders; ++li) {
02217 const dimension_type i = no_sing_leaders[li];
02218 const dimension_type ci = coherent_index(i);
02219 typename OR_Matrix<N>::const_row_reference_type m_i =
02220 *(matrix.row_begin()+i);
02221 typename OR_Matrix<N>::row_reference_type aux_i =
02222 *(aux.matrix.row_begin()+i);
02223 if (i%2 == 0) {
02224
02225
02226
02227
02228
02229 if (i != successor[i]) {
02230 dimension_type j = i;
02231 dimension_type next_j = successor[j];
02232 while (j != next_j) {
02233 aux.matrix[next_j][j] = matrix[next_j][j];
02234 j = next_j;
02235 next_j = successor[j];
02236 }
02237 const dimension_type cj = coherent_index(j);
02238 aux.matrix[cj][ci] = matrix[cj][ci];
02239 }
02240 }
02241
02242 dimension_type rs_li = (li%2) ? li :li+1;
02243
02244 DIRTY_TEMP(N, tmp);
02245 for (dimension_type lj = 0 ; lj <= rs_li; ++lj) {
02246 const dimension_type j = no_sing_leaders[lj];
02247 const dimension_type cj = coherent_index(j);
02248 const N& m_i_j = m_i[j];
02249 const N& m_i_ci = m_i[ci];
02250 bool to_add = true;
02251
02252
02253
02254 if (j != ci) {
02255 add_assign_r(tmp, m_i_ci, matrix[cj][j], ROUND_UP);
02256 div2exp_assign_r(tmp, tmp, 1, ROUND_UP);
02257 if (m_i_j >= tmp) {
02258 to_add = false;
02259 continue;
02260 }
02261 }
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271
02272 for (dimension_type lk = 0; lk < num_no_sing_leaders; ++lk) {
02273 const dimension_type k = no_sing_leaders[lk];
02274 if (k != i && k != j) {
02275 dimension_type ck = coherent_index(k);
02276 if (k < j)
02277
02278 add_assign_r(tmp, m_i[k], matrix[cj][ck], ROUND_UP);
02279 else if (k < i)
02280
02281 add_assign_r(tmp, m_i[k], matrix[k][j], ROUND_UP);
02282 else
02283
02284 add_assign_r(tmp, matrix[ck][ci], matrix[k][j], ROUND_UP);
02285
02286
02287 if (m_i_j >= tmp) {
02288 to_add = false;
02289 break;
02290 }
02291 }
02292 }
02293
02294
02295 if (to_add)
02296 aux_i[j] = m_i_j;
02297 }
02298 }
02299
02300
02301
02302
02303
02304 if (exist_sing_class) {
02305 aux.matrix[sing_leader][sing_leader+1]
02306 = matrix[sing_leader][sing_leader+1];
02307 if (successor[sing_leader+1] != sing_leader+1) {
02308 dimension_type j = sing_leader;
02309 dimension_type next_jj = successor[j+1];
02310 while (next_jj != j+1) {
02311 aux.matrix[next_jj][j] = matrix[next_jj][j];
02312 j = next_jj;
02313 next_jj = successor[j+1];
02314 }
02315 aux.matrix[j+1][j] = matrix[j+1][j];
02316 }
02317 else
02318 aux.matrix[sing_leader+1][sing_leader]
02319 = matrix[sing_leader+1][sing_leader];
02320 }
02321
02322 Octagonal_Shape<T>& x = const_cast<Octagonal_Shape<T>&>(*this);
02323 aux.reset_strongly_closed();
02324
02325 #ifndef NDEBUG
02326 {
02327
02328 const Octagonal_Shape x_copy = *this;
02329 const Octagonal_Shape y_copy = aux;
02330 assert(x_copy == y_copy);
02331 }
02332 #endif
02333
02334 std::swap(x, aux);
02335 assert(is_strongly_reduced());
02336 }
02337
02338
02339 template <typename T>
02340 void
02341 Octagonal_Shape<T>::upper_bound_assign(const Octagonal_Shape& y) {
02342
02343 if (space_dim != y.space_dim)
02344 throw_dimension_incompatible("upper_bound_assign(y)", y);
02345
02346
02347 y.strong_closure_assign();
02348 if (y.marked_empty())
02349 return;
02350 strong_closure_assign();
02351 if (marked_empty()) {
02352 *this = y;
02353 return;
02354 }
02355
02356
02357 typename OR_Matrix<N>::const_element_iterator j = y.matrix.element_begin();
02358 for (typename OR_Matrix<N>::element_iterator i = matrix.element_begin(),
02359 matrix_element_end = matrix.element_end();
02360 i != matrix_element_end; ++i, ++j)
02361 max_assign(*i, *j);
02362
02363
02364 assert(OK());
02365 }
02366
02367 template <typename T>
02368 void
02369 Octagonal_Shape<T>::difference_assign(const Octagonal_Shape& y) {
02370
02371 if (space_dim != y.space_dim)
02372 throw_dimension_incompatible("difference_assign(y)", y);
02373
02374 Octagonal_Shape& x = *this;
02375
02376
02377
02378 x.strong_closure_assign();
02379
02380 if (x.marked_empty())
02381 return;
02382
02383 if (y.marked_empty())
02384 return;
02385
02386
02387
02388
02389 if (x.space_dim == 0) {
02390 x.set_empty();
02391 return;
02392 }
02393
02394
02395
02396 if (y.contains(x)) {
02397 x.set_empty();
02398 return;
02399 }
02400
02401 Octagonal_Shape new_oct(space_dim, EMPTY);
02402
02403
02404
02405 const Constraint_System& y_cs = y.constraints();
02406 for (Constraint_System::const_iterator i = y_cs.begin(),
02407 y_cs_end = y_cs.end(); i != y_cs_end; ++i) {
02408 const Constraint& c = *i;
02409
02410
02411
02412
02413 if (x.relation_with(c).implies(Poly_Con_Relation::is_included()))
02414 continue;
02415 Octagonal_Shape z = x;
02416 const Linear_Expression e = Linear_Expression(c);
02417 z.add_constraint(e <= 0);
02418 if (!z.is_empty())
02419 new_oct.upper_bound_assign(z);
02420 if (c.is_equality()) {
02421 z = x;
02422 z.add_constraint(e >= 0);
02423 if (!z.is_empty())
02424 new_oct.upper_bound_assign(z);
02425 }
02426 }
02427 *this = new_oct;
02428 assert(OK());
02429 }
02430
02431 template <typename T>
02432 bool
02433 Octagonal_Shape<T>
02434 ::simplify_using_context_assign(const Octagonal_Shape& y) {
02435
02436 used(y);
02437 return true;
02438 }
02439
02440 template <typename T>
02441 void
02442 Octagonal_Shape<T>::add_space_dimensions_and_embed(dimension_type m) {
02443
02444 if (m == 0)
02445 return;
02446
02447 const dimension_type new_dim = space_dim + m;
02448 const bool was_zero_dim_univ = !marked_empty() && space_dim == 0;
02449
02450
02451
02452 matrix.grow(new_dim);
02453 space_dim = new_dim;
02454
02455
02456 if (was_zero_dim_univ)
02457 set_strongly_closed();
02458
02459 assert(OK());
02460 }
02461
02462 template <typename T>
02463 void
02464 Octagonal_Shape<T>::add_space_dimensions_and_project(dimension_type m) {
02465
02466 if (m == 0)
02467 return;
02468
02469 const dimension_type n = matrix.num_rows();
02470
02471
02472
02473 add_space_dimensions_and_embed(m);
02474
02475
02476 for (typename OR_Matrix<N>::row_iterator i = matrix.row_begin() + n,
02477 matrix_row_end = matrix.row_end(); i != matrix_row_end; i += 2) {
02478 typename OR_Matrix<N>::row_reference_type x_i = *i;
02479 typename OR_Matrix<N>::row_reference_type x_ci = *(i+1);
02480 const dimension_type ind = i.index();
02481 assign_r(x_i[ind+1], 0, ROUND_NOT_NEEDED);
02482 assign_r(x_ci[ind], 0, ROUND_NOT_NEEDED);
02483 }
02484
02485 if (marked_strongly_closed())
02486 reset_strongly_closed();
02487 assert(OK());
02488 }
02489
02490 template <typename T>
02491 void
02492 Octagonal_Shape<T>
02493 ::remove_space_dimensions(const Variables_Set& to_be_removed) {
02494
02495
02496
02497 if (to_be_removed.empty()) {
02498 assert(OK());
02499 return;
02500 }
02501
02502
02503 const dimension_type min_space_dim = to_be_removed.space_dimension();
02504 if (space_dim < min_space_dim)
02505 throw_dimension_incompatible("remove_space_dimensions(vs)", min_space_dim);
02506
02507 const dimension_type new_space_dim = space_dim - to_be_removed.size();
02508
02509 strong_closure_assign();
02510
02511
02512 if (new_space_dim == 0) {
02513 matrix.shrink(0);
02514 if (!marked_empty())
02515
02516 set_zero_dim_univ();
02517 space_dim = 0;
02518 assert(OK());
02519 return;
02520 }
02521
02522
02523
02524
02525 Variables_Set::const_iterator tbr = to_be_removed.begin();
02526 dimension_type ftr = *tbr;
02527 dimension_type ftr_size = 2*ftr*(ftr+1);
02528 typename OR_Matrix<N>::element_iterator
02529 iter = matrix.element_begin()+ftr_size;
02530
02531 dimension_type i = ftr + 1;
02532 while (i < space_dim) {
02533 if (to_be_removed.count(i) != 0)
02534 ++i;
02535 else {
02536 typename OR_Matrix<N>::row_iterator
02537 row_iter = matrix.row_begin()+2*i;
02538 typename OR_Matrix<N>::row_reference_type
02539 row_ref = *row_iter;
02540 typename OR_Matrix<N>::row_reference_type
02541 row_ref1 = *(++row_iter);
02542
02543
02544
02545
02546
02547
02548 for (dimension_type j = 0; j <= i; ++j)
02549 if (to_be_removed.count(j) == 0) {
02550 assign_or_swap(*(iter++), row_ref[2*j]);
02551 assign_or_swap(*(iter++), row_ref[2*j+1]);
02552 }
02553 for (dimension_type j = 0; j <= i; ++j)
02554 if (to_be_removed.count(j) == 0) {
02555 assign_or_swap(*(iter++), row_ref1[2*j]);
02556 assign_or_swap(*(iter++), row_ref1[2*j+1]);
02557 }
02558 ++i;
02559 }
02560 }
02561
02562 matrix.shrink(new_space_dim);
02563 space_dim = new_space_dim;
02564 assert(OK());
02565 }
02566
02567 template <typename T>
02568 template <typename Partial_Function>
02569 void
02570 Octagonal_Shape<T>::map_space_dimensions(const Partial_Function& pfunc) {
02571 if (space_dim == 0)
02572 return;
02573
02574 if (pfunc.has_empty_codomain()) {
02575
02576 remove_higher_space_dimensions(0);
02577 return;
02578 }
02579
02580 const dimension_type new_space_dim = pfunc.max_in_codomain() + 1;
02581
02582
02583 if (new_space_dim < space_dim)
02584 strong_closure_assign();
02585
02586
02587
02588 if (marked_empty()) {
02589 remove_higher_space_dimensions(new_space_dim);
02590 return;
02591 }
02592
02593
02594 OR_Matrix<N> x(new_space_dim);
02595
02596 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
02597 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
02598
02599 Row_Iterator m_begin = x.row_begin();
02600 Row_Iterator m_end = x.row_end();
02601
02602 for (Row_Iterator i_iter = matrix.row_begin(), i_end = matrix.row_end();
02603 i_iter != i_end; i_iter += 2) {
02604 dimension_type new_i;
02605 dimension_type i = i_iter.index()/2;
02606
02607
02608
02609 if (pfunc.maps(i, new_i)) {
02610 Row_Reference r_i = *i_iter;
02611 Row_Reference r_ii = *(i_iter + 1);
02612 dimension_type double_new_i = 2*new_i;
02613 Row_Iterator x_iter = m_begin + double_new_i;
02614 Row_Reference x_i = *x_iter;
02615 Row_Reference x_ii = *(x_iter + 1);
02616 for (dimension_type j = 0; j <= i; ++j) {
02617 dimension_type new_j;
02618
02619 if (pfunc.maps(j, new_j)) {
02620 dimension_type dj = 2*j;
02621 dimension_type double_new_j = 2*new_j;
02622
02623
02624
02625
02626 if (new_i >= new_j) {
02627 assign_or_swap(x_i[double_new_j], r_i[dj]);
02628 assign_or_swap(x_ii[double_new_j], r_ii[dj]);
02629 assign_or_swap(x_ii[double_new_j+1], r_ii[dj + 1]);
02630 assign_or_swap(x_i[double_new_j+1], r_i[dj + 1]);
02631 }
02632 else {
02633 Row_Iterator xj_iter = m_begin + double_new_j;
02634 Row_Reference x_j = *xj_iter;
02635 Row_Reference x_jj = *(xj_iter + 1);
02636 assign_or_swap(x_jj[double_new_i+1], r_i[dj]);
02637 assign_or_swap(x_jj[double_new_i], r_ii[dj]);
02638 assign_or_swap(x_j[double_new_i+1], r_i[dj+1]);
02639 assign_or_swap(x_j[double_new_i], r_ii[dj+1]);
02640 }
02641
02642 }
02643 }
02644 }
02645 }
02646
02647 std::swap(matrix, x);
02648 space_dim = new_space_dim;
02649 assert(OK());
02650 }
02651
02652 template <typename T>
02653 void
02654 Octagonal_Shape<T>::intersection_assign(const Octagonal_Shape& y) {
02655
02656 if (space_dim != y.space_dim)
02657 throw_dimension_incompatible("intersection_assign(y)", y);
02658
02659
02660 if (marked_empty())
02661 return;
02662 if (y.marked_empty()) {
02663 set_empty();
02664 return;
02665 }
02666
02667
02668
02669 if (space_dim == 0)
02670 return;
02671
02672
02673
02674 bool changed = false;
02675
02676 typename OR_Matrix<N>::const_element_iterator j = y.matrix.element_begin();
02677 for (typename OR_Matrix<N>::element_iterator i = matrix.element_begin(),
02678 matrix_element_end = matrix.element_end();
02679 i != matrix_element_end;
02680 ++i, ++j) {
02681 N& elem = *i;
02682 const N& y_elem = *j;
02683 if (y_elem < elem) {
02684 elem = y_elem;
02685 changed = true;
02686 }
02687 }
02688
02689
02690 if (changed && marked_strongly_closed())
02691 reset_strongly_closed();
02692 assert(OK());
02693 }
02694
02695 template <typename T>
02696 template <typename Iterator>
02697 void
02698 Octagonal_Shape<T>::CC76_extrapolation_assign(const Octagonal_Shape& y,
02699 Iterator first, Iterator last,
02700 unsigned* tp) {
02701
02702 if (space_dim != y.space_dim)
02703 throw_dimension_incompatible("CC76_extrapolation_assign(y)", y);
02704
02705 #ifndef NDEBUG
02706 {
02707
02708 const Octagonal_Shape x_copy = *this;
02709 const Octagonal_Shape y_copy = y;
02710 assert(x_copy.contains(y_copy));
02711 }
02712 #endif
02713
02714
02715
02716 if (space_dim == 0)
02717 return;
02718
02719 strong_closure_assign();
02720
02721 if (marked_empty())
02722 return;
02723 y.strong_closure_assign();
02724
02725 if (y.marked_empty())
02726 return;
02727
02728
02729 if (tp != 0 && *tp > 0) {
02730 Octagonal_Shape x_tmp(*this);
02731 x_tmp.CC76_extrapolation_assign(y, first, last, 0);
02732
02733 if (!contains(x_tmp))
02734 --(*tp);
02735 return;
02736 }
02737
02738
02739
02740
02741
02742
02743
02744
02745
02746 typename OR_Matrix<N>::const_element_iterator j = y.matrix.element_begin();
02747 for (typename OR_Matrix<N>::element_iterator i = matrix.element_begin(),
02748 matrix_element_end = matrix.element_end();
02749 i != matrix_element_end;
02750 ++i, ++j) {
02751 const N& y_elem = *j;
02752 N& elem = *i;
02753 if (y_elem < elem) {
02754 Iterator k = std::lower_bound(first, last, elem);
02755 if (k != last) {
02756 if (elem < *k)
02757 assign_r(elem, *k, ROUND_UP);
02758 }
02759 else
02760 assign_r(elem, PLUS_INFINITY, ROUND_NOT_NEEDED);
02761 }
02762 }
02763
02764 reset_strongly_closed();
02765 assert(OK());
02766 }
02767
02768 template <typename T>
02769 void
02770 Octagonal_Shape<T>
02771 ::get_limiting_octagon(const Constraint_System& cs,
02772 Octagonal_Shape& limiting_octagon) const {
02773 const dimension_type cs_space_dim = cs.space_dimension();
02774
02775 assert(cs_space_dim <= space_dim);
02776
02777 strong_closure_assign();
02778 bool is_oct_changed = false;
02779
02780
02781 TEMP_INTEGER(coeff);
02782 TEMP_INTEGER(term);
02783 DIRTY_TEMP(N, d);
02784
02785 for (Constraint_System::const_iterator cs_i = cs.begin(),
02786 cs_end = cs.end(); cs_i != cs_end; ++cs_i) {
02787 const Constraint& c = *cs_i;
02788 dimension_type num_vars = 0;
02789 dimension_type i = 0;
02790 dimension_type j = 0;
02791
02792 if (!extract_octagonal_difference(c, cs_space_dim, num_vars, i, j,
02793 coeff, term))
02794 continue;
02795
02796 typedef typename OR_Matrix<N>::const_row_iterator Row_iterator;
02797 typedef typename OR_Matrix<N>::const_row_reference_type Row_reference;
02798 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
02799 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
02800 Row_iterator m_begin = matrix.row_begin();
02801
02802 Row_iterator i_iter = m_begin + i;
02803 Row_reference m_i = *i_iter;
02804 OR_Matrix<N>& lo_mat = limiting_octagon.matrix;
02805 Row_Iterator lo_iter = lo_mat.row_begin() + i;
02806 Row_Reference lo_m_i = *lo_iter;
02807 N& lo_m_i_j = lo_m_i[j];
02808 if (coeff < 0)
02809 neg_assign(coeff);
02810
02811 div_round_up(d, term, coeff);
02812 if (m_i[j] <= d)
02813 if (c.is_inequality()) {
02814 if (lo_m_i_j > d) {
02815 lo_m_i_j = d;
02816 is_oct_changed = true;
02817 }
02818 else {
02819
02820 if (i%2 == 0) {
02821 ++i_iter;
02822 ++lo_iter;
02823 }
02824 else {
02825 --i_iter;
02826 --lo_iter;
02827 }
02828 Row_reference m_ci = *i_iter;
02829 Row_Reference lo_m_ci = *lo_iter;
02830
02831 dimension_type cj = coherent_index(j);
02832 N& lo_m_ci_cj = lo_m_ci[cj];
02833 neg_assign(term);
02834 div_round_up(d, term, coeff);
02835 if (m_ci[cj] <= d && lo_m_ci_cj > d) {
02836 lo_m_ci_cj = d;
02837 is_oct_changed = true;
02838 }
02839 }
02840 }
02841 }
02842
02843
02844 if (is_oct_changed && limiting_octagon.marked_strongly_closed())
02845 limiting_octagon.reset_strongly_closed();
02846 }
02847
02848 template <typename T>
02849 void
02850 Octagonal_Shape<T>
02851 ::limited_CC76_extrapolation_assign(const Octagonal_Shape& y,
02852 const Constraint_System& cs,
02853 unsigned* tp) {
02854
02855
02856 if (space_dim != y.space_dim)
02857 throw_dimension_incompatible("limited_CC76_extrapolation_assign(y, cs)",
02858 y);
02859
02860 const dimension_type cs_space_dim = cs.space_dimension();
02861 if (space_dim < cs_space_dim)
02862 throw_constraint_incompatible("limited_CC76_extrapolation_assign(y, cs)");
02863
02864
02865 if (cs.has_strict_inequalities())
02866 throw_constraint_incompatible("limited_CC76_extrapolation_assign(y, cs)");
02867
02868
02869
02870
02871 if (space_dim == 0)
02872 return;
02873
02874 #ifndef NDEBUG
02875 {
02876
02877 const Octagonal_Shape x_copy = *this;
02878 const Octagonal_Shape y_copy = y;
02879 assert(x_copy.contains(y_copy));
02880 }
02881 #endif
02882
02883
02884 if (marked_empty())
02885 return;
02886
02887 if (y.marked_empty())
02888 return;
02889
02890 Octagonal_Shape limiting_octagon(space_dim, UNIVERSE);
02891 get_limiting_octagon(cs, limiting_octagon);
02892 CC76_extrapolation_assign(y, tp);
02893 intersection_assign(limiting_octagon);
02894 }
02895
02896 template <typename T>
02897 void
02898 Octagonal_Shape<T>::BHMZ05_widening_assign(const Octagonal_Shape& y,
02899 unsigned* tp) {
02900
02901 if (space_dim != y.space_dim)
02902 throw_dimension_incompatible("BHMZ05_widening_assign(y)", y);
02903
02904 #ifndef NDEBUG
02905 {
02906
02907 const Octagonal_Shape x_copy = *this;
02908 const Octagonal_Shape y_copy = y;
02909 assert(x_copy.contains(y_copy));
02910 }
02911 #endif
02912
02913
02914 const dimension_type y_affine_dim = y.affine_dimension();
02915
02916
02917
02918 if (y_affine_dim == 0)
02919 return;
02920
02921
02922
02923 const dimension_type x_affine_dim = affine_dimension();
02924 assert(x_affine_dim >= y_affine_dim);
02925 if (x_affine_dim != y_affine_dim)
02926 return;
02927
02928
02929 if (tp != 0 && *tp > 0) {
02930 Octagonal_Shape x_tmp(*this);
02931 x_tmp.BHMZ05_widening_assign(y, 0);
02932
02933 if (!contains(x_tmp))
02934 --(*tp);
02935 return;
02936 }
02937
02938
02939 assert(marked_strongly_closed() && y.marked_strongly_closed());
02940
02941 y.strong_reduction_assign();
02942
02943
02944 typename OR_Matrix<N>::const_element_iterator j = y.matrix.element_begin();
02945 for (typename OR_Matrix<N>::element_iterator i = matrix.element_begin(),
02946 matrix_element_end = matrix.element_end();
02947 i != matrix_element_end;
02948 ++i, ++j) {
02949 N& elem = *i;
02950
02951
02952
02953 if (*j != elem)
02954 assign_r(elem, PLUS_INFINITY, ROUND_NOT_NEEDED);
02955 }
02956 reset_strongly_closed();
02957 assert(OK());
02958 }
02959
02960 template <typename T>
02961 void
02962 Octagonal_Shape<T>
02963 ::limited_BHMZ05_extrapolation_assign(const Octagonal_Shape& y,
02964 const Constraint_System& cs,
02965 unsigned* tp) {
02966
02967
02968 if (space_dim != y.space_dim)
02969 throw_dimension_incompatible("limited_BHMZ05_extrapolation_assign(y, cs)",
02970 y);
02971
02972 const dimension_type cs_space_dim = cs.space_dimension();
02973 if (space_dim < cs_space_dim)
02974 throw_constraint_incompatible("limited_CH78_extrapolation_assign(y, cs)");
02975
02976
02977 if (cs.has_strict_inequalities())
02978 throw_constraint_incompatible("limited_CH78_extrapolation_assign(y, cs)");
02979
02980
02981
02982
02983 if (space_dim == 0)
02984 return;
02985
02986 #ifndef NDEBUG
02987 {
02988
02989 const Octagonal_Shape x_copy = *this;
02990 const Octagonal_Shape y_copy = y;
02991 assert(x_copy.contains(y_copy));
02992 }
02993 #endif
02994
02995
02996 if (marked_empty())
02997 return;
02998
02999 if (y.marked_empty())
03000 return;
03001
03002 Octagonal_Shape limiting_octagon(space_dim, UNIVERSE);
03003 get_limiting_octagon(cs, limiting_octagon);
03004 BHMZ05_widening_assign(y, tp);
03005 intersection_assign(limiting_octagon);
03006 }
03007
03008 template <typename T>
03009 void
03010 Octagonal_Shape<T>::CC76_narrowing_assign(const Octagonal_Shape& y) {
03011
03012 if (space_dim != y.space_dim)
03013 throw_dimension_incompatible("CC76_narrowing_assign(y)", y);
03014
03015 #ifndef NDEBUG
03016 {
03017
03018 const Octagonal_Shape x_copy = *this;
03019 const Octagonal_Shape y_copy = y;
03020 assert(y_copy.contains(x_copy));
03021 }
03022 #endif
03023
03024
03025
03026 if (space_dim == 0)
03027 return;
03028
03029 y.strong_closure_assign();
03030
03031 if (y.marked_empty())
03032 return;
03033 strong_closure_assign();
03034
03035 if (marked_empty())
03036 return;
03037
03038
03039
03040 bool is_oct_changed = false;
03041 typename OR_Matrix<N>::const_element_iterator j = y.matrix.element_begin();
03042 for (typename OR_Matrix<N>::element_iterator i = matrix.element_begin(),
03043 matrix_element_end = matrix.element_end();
03044 i != matrix_element_end;
03045 ++i, ++j) {
03046 if (!is_plus_infinity(*i)
03047 && !is_plus_infinity(*j)
03048 && *i != *j) {
03049 *i = *j;
03050 is_oct_changed = true;
03051 }
03052 }
03053
03054 if (is_oct_changed && marked_strongly_closed())
03055 reset_strongly_closed();
03056 assert(OK());
03057 }
03058
03059 template <typename T>
03060 void
03061 Octagonal_Shape<T>
03062 ::deduce_v_pm_u_bounds(const dimension_type v_id,
03063 const dimension_type last_id,
03064 const Linear_Expression& sc_expr,
03065 Coefficient_traits::const_reference sc_den,
03066 const N& ub_v) {
03067
03068 assert(sc_den > 0);
03069 assert(!is_plus_infinity(ub_v));
03070
03071 DIRTY_TEMP0(mpq_class, mpq_sc_den);
03072 assign_r(mpq_sc_den, sc_den, ROUND_NOT_NEEDED);
03073
03074
03075 const dimension_type n_v = 2*v_id;
03076 typename OR_Matrix<N>::row_reference_type m_cv = matrix[n_v+1];
03077
03078
03079 DIRTY_TEMP(N, half);
03080 DIRTY_TEMP0(mpq_class, minus_lb_u);
03081 DIRTY_TEMP0(mpq_class, q);
03082 DIRTY_TEMP0(mpq_class, minus_q);
03083 DIRTY_TEMP0(mpq_class, ub_u);
03084 DIRTY_TEMP0(mpq_class, lb_u);
03085 DIRTY_TEMP(N, up_approx);
03086 TEMP_INTEGER(minus_expr_u);
03087
03088 for (dimension_type u_id = last_id+1; u_id-- > 0; ) {
03089
03090 if (u_id == v_id)
03091 continue;
03092 const Coefficient& expr_u = sc_expr.coefficient(Variable(u_id));
03093
03094 if (expr_u == 0)
03095 continue;
03096
03097 const dimension_type n_u = u_id*2;
03098
03099 if (expr_u > 0) {
03100 if (expr_u >= sc_den) {
03101
03102
03103
03104
03105 div2exp_assign_r(half, matrix[n_u+1][n_u], 1, ROUND_UP);
03106 N& m_v_minus_u = (n_v < n_u) ? matrix[n_u][n_v] : m_cv[n_u+1];
03107 sub_assign_r(m_v_minus_u, ub_v, half, ROUND_UP);
03108 }
03109 else {
03110
03111 typename OR_Matrix<N>::row_reference_type m_u = matrix[n_u];
03112 const N& m_u_cu = m_u[n_u+1];
03113 if (!is_plus_infinity(m_u_cu)) {
03114
03115
03116
03117
03118 assign_r(minus_lb_u, m_u_cu, ROUND_NOT_NEEDED);
03119 div2exp_assign_r(minus_lb_u, minus_lb_u, 1, ROUND_NOT_NEEDED);
03120 assign_r(q, expr_u, ROUND_NOT_NEEDED);
03121 div_assign_r(q, q, mpq_sc_den, ROUND_NOT_NEEDED);
03122 assign_r(ub_u, matrix[n_u+1][n_u], ROUND_NOT_NEEDED);
03123 div2exp_assign_r(ub_u, ub_u, 1, ROUND_NOT_NEEDED);
03124
03125 add_assign_r(ub_u, ub_u, minus_lb_u, ROUND_NOT_NEEDED);
03126
03127 sub_mul_assign_r(minus_lb_u, q, ub_u, ROUND_NOT_NEEDED);
03128 assign_r(up_approx, minus_lb_u, ROUND_UP);
03129
03130 N& m_v_minus_u = (n_v < n_u) ? m_u[n_v] : m_cv[n_u+1];
03131 add_assign_r(m_v_minus_u, ub_v, up_approx, ROUND_UP);
03132 }
03133 }
03134 }
03135 else {
03136 assert(expr_u < 0);
03137
03138 neg_assign(minus_expr_u, expr_u);
03139 if (minus_expr_u >= sc_den) {
03140
03141
03142
03143
03144 div2exp_assign_r(half, matrix[n_u][n_u+1], 1, ROUND_UP);
03145 N& m_v_plus_u = (n_v < n_u) ? matrix[n_u+1][n_v] : m_cv[n_u];
03146 sub_assign_r(m_v_plus_u, ub_v, half, ROUND_UP);
03147 }
03148 else {
03149
03150 typename OR_Matrix<N>::row_reference_type m_cu = matrix[n_u+1];
03151 const N& m_cu_u = m_cu[n_u];
03152 if (!is_plus_infinity(m_cu_u)) {
03153
03154
03155
03156
03157 assign_r(ub_u, m_cu[n_u], ROUND_NOT_NEEDED);
03158 div2exp_assign_r(ub_u, ub_u, 1, ROUND_NOT_NEEDED);
03159 assign_r(minus_q, minus_expr_u, ROUND_NOT_NEEDED);
03160 div_assign_r(minus_q, minus_q, mpq_sc_den, ROUND_NOT_NEEDED);
03161 assign_r(lb_u, matrix[n_u][n_u+1], ROUND_NOT_NEEDED);
03162 div2exp_assign_r(lb_u, lb_u, 1, ROUND_NOT_NEEDED);
03163 neg_assign_r(lb_u, lb_u, ROUND_NOT_NEEDED);
03164
03165 sub_assign_r(lb_u, lb_u, ub_u, ROUND_NOT_NEEDED);
03166
03167 add_mul_assign_r(ub_u, minus_q, lb_u, ROUND_NOT_NEEDED);
03168 assign_r(up_approx, ub_u, ROUND_UP);
03169
03170 N& m_v_plus_u = (n_v < n_u) ? m_cu[n_v] : m_cv[n_u];
03171 add_assign_r(m_v_plus_u, ub_v, up_approx, ROUND_UP);
03172 }
03173 }
03174 }
03175 }
03176 }
03177
03178 template <typename T>
03179 void
03180 Octagonal_Shape<T>
03181 ::deduce_minus_v_pm_u_bounds(const dimension_type v_id,
03182 const dimension_type last_id,
03183 const Linear_Expression& sc_expr,
03184 Coefficient_traits::const_reference sc_den,
03185 const N& minus_lb_v) {
03186
03187 assert(sc_den > 0);
03188 assert(!is_plus_infinity(minus_lb_v));
03189
03190 DIRTY_TEMP0(mpq_class, mpq_sc_den);
03191 assign_r(mpq_sc_den, sc_den, ROUND_NOT_NEEDED);
03192
03193
03194 const dimension_type n_v = 2*v_id;
03195 typename OR_Matrix<N>::row_reference_type m_v = matrix[n_v];
03196
03197
03198 DIRTY_TEMP(N, half);
03199 DIRTY_TEMP0(mpq_class, ub_u);
03200 DIRTY_TEMP0(mpq_class, q);
03201 DIRTY_TEMP0(mpq_class, minus_lb_u);
03202 DIRTY_TEMP(N, up_approx);
03203 TEMP_INTEGER(minus_expr_u);
03204
03205 for (dimension_type u_id = last_id+1; u_id-- > 0; ) {
03206
03207 if (u_id == v_id)
03208 continue;
03209 const Coefficient& expr_u = sc_expr.coefficient(Variable(u_id));
03210
03211 if (expr_u == 0)
03212 continue;
03213
03214 const dimension_type n_u = u_id*2;
03215
03216 if (expr_u > 0) {
03217 if (expr_u >= sc_den) {
03218
03219
03220
03221
03222
03223 div2exp_assign_r(half, matrix[n_u][n_u+1], 1, ROUND_UP);
03224 N& m_u_minus_v = (n_v < n_u) ? matrix[n_u+1][n_v+1] : m_v[n_u];
03225 sub_assign_r(m_u_minus_v, minus_lb_v, half, ROUND_UP);
03226 }
03227 else {
03228
03229 typename OR_Matrix<N>::row_reference_type m_cu = matrix[n_u+1];
03230 const N& m_cu_u = m_cu[n_u];
03231 if (!is_plus_infinity(m_cu_u)) {
03232
03233
03234
03235
03236 assign_r(ub_u, m_cu[n_u], ROUND_NOT_NEEDED);
03237 div2exp_assign_r(ub_u, ub_u, 1, ROUND_NOT_NEEDED);
03238 assign_r(q, expr_u, ROUND_NOT_NEEDED);
03239 div_assign_r(q, q, mpq_sc_den, ROUND_NOT_NEEDED);
03240 assign_r(minus_lb_u, matrix[n_u][n_u+1], ROUND_NOT_NEEDED);
03241 div2exp_assign_r(minus_lb_u, minus_lb_u, 1, ROUND_NOT_NEEDED);
03242
03243 add_assign_r(minus_lb_u, ub_u, minus_lb_u, ROUND_NOT_NEEDED);
03244
03245 sub_mul_assign_r(ub_u, q, minus_lb_u, ROUND_NOT_NEEDED);
03246 assign_r(up_approx, ub_u, ROUND_UP);
03247
03248 N& m_u_minus_v = (n_v < n_u) ? m_cu[n_v+1] : m_v[n_u];
03249 add_assign_r(m_u_minus_v, minus_lb_v, up_approx, ROUND_UP);
03250 }
03251 }
03252 }
03253 else {
03254 assert(expr_u < 0);
03255
03256 neg_assign(minus_expr_u, expr_u);
03257 if (minus_expr_u >= sc_den) {
03258
03259
03260
03261
03262 div2exp_assign_r(half, matrix[n_u+1][n_u], 1, ROUND_UP);
03263 N& m_minus_v_minus_u = (n_v < n_u) ? matrix[n_u][n_v+1] : m_v[n_u+1];
03264 sub_assign_r(m_minus_v_minus_u, minus_lb_v, half, ROUND_UP);
03265 }
03266 else {
03267
03268 typename OR_Matrix<N>::row_reference_type m_u = matrix[n_u];
03269 const N& m_u_cu = m_u[n_u+1];
03270 if (!is_plus_infinity(m_u_cu)) {
03271
03272
03273
03274
03275 assign_r(ub_u, matrix[n_u+1][n_u], ROUND_NOT_NEEDED);
03276 div2exp_assign_r(ub_u, ub_u, 1, ROUND_NOT_NEEDED);
03277 assign_r(q, expr_u, ROUND_NOT_NEEDED);
03278 div_assign_r(q, q, mpq_sc_den, ROUND_NOT_NEEDED);
03279 assign_r(minus_lb_u, m_u[n_u+1], ROUND_NOT_NEEDED);
03280 div2exp_assign_r(minus_lb_u, minus_lb_u, 1, ROUND_NOT_NEEDED);
03281
03282 add_assign_r(ub_u, ub_u, minus_lb_u, ROUND_NOT_NEEDED);
03283
03284 add_mul_assign_r(minus_lb_u, q, ub_u, ROUND_NOT_NEEDED);
03285 assign_r(up_approx, minus_lb_u, ROUND_UP);
03286
03287 N& m_minus_v_minus_u = (n_v < n_u) ? m_u[n_v+1] : m_v[n_u+1];
03288 add_assign_r(m_minus_v_minus_u, minus_lb_v, up_approx, ROUND_UP);
03289 }
03290 }
03291 }
03292 }
03293 }
03294
03295 template <typename T>
03296 void
03297 Octagonal_Shape<T>
03298 ::forget_all_octagonal_constraints(const dimension_type v_id) {
03299 assert(v_id < space_dim);
03300 const dimension_type n_v = 2*v_id;
03301 typename OR_Matrix<N>::row_iterator m_iter = matrix.row_begin() + n_v;
03302 typename OR_Matrix<N>::row_reference_type r_v = *m_iter;
03303 typename OR_Matrix<N>::row_reference_type r_cv = *(++m_iter);
03304 for (dimension_type h = m_iter.row_size(); h-- > 0; ) {
03305 assign_r(r_v[h], PLUS_INFINITY, ROUND_NOT_NEEDED);
03306 assign_r(r_cv[h], PLUS_INFINITY, ROUND_NOT_NEEDED);
03307 }
03308 ++m_iter;
03309 for (typename OR_Matrix<N>::row_iterator m_end = matrix.row_end();
03310 m_iter != m_end; ++m_iter) {
03311 typename OR_Matrix<N>::row_reference_type r = *m_iter;
03312 assign_r(r[n_v], PLUS_INFINITY, ROUND_NOT_NEEDED);
03313 assign_r(r[n_v+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
03314 }
03315 }
03316
03317 template <typename T>
03318 void
03319 Octagonal_Shape<T>
03320 ::forget_binary_octagonal_constraints(const dimension_type v_id) {
03321 assert(v_id < space_dim);
03322 const dimension_type n_v = 2*v_id;
03323 typename OR_Matrix<N>::row_iterator m_iter = matrix.row_begin() + n_v;
03324 typename OR_Matrix<N>::row_reference_type r_v = *m_iter;
03325 typename OR_Matrix<N>::row_reference_type r_cv = *(++m_iter);
03326 for (dimension_type k = n_v; k-- > 0; ) {
03327 assign_r(r_v[k], PLUS_INFINITY, ROUND_NOT_NEEDED);
03328 assign_r(r_cv[k], PLUS_INFINITY, ROUND_NOT_NEEDED);
03329 }
03330 ++m_iter;
03331 for (typename OR_Matrix<N>::row_iterator m_end = matrix.row_end();
03332 m_iter != m_end; ++m_iter) {
03333 typename OR_Matrix<N>::row_reference_type r = *m_iter;
03334 assign_r(r[n_v], PLUS_INFINITY, ROUND_NOT_NEEDED);
03335 assign_r(r[n_v+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
03336 }
03337 }
03338
03339 template <typename T>
03340 void
03341 Octagonal_Shape<T>::unconstrain(const Variable var) {
03342
03343 const dimension_type dim = var.id();
03344 if (space_dimension() < dim)
03345 throw_dimension_incompatible("unconstrain(var)", dim);
03346
03347
03348 strong_closure_assign();
03349
03350
03351 if (marked_empty())
03352 return;
03353
03354 forget_all_octagonal_constraints(dim);
03355
03356 assert(OK());
03357 }
03358
03359 template <typename T>
03360 void
03361 Octagonal_Shape<T>::unconstrain(const Variables_Set& to_be_unconstrained) {
03362
03363
03364 if (to_be_unconstrained.empty())
03365 return;
03366
03367
03368 const dimension_type min_space_dim = to_be_unconstrained.space_dimension();
03369 if (space_dimension() < min_space_dim)
03370 throw_dimension_incompatible("unconstrain(vs)", min_space_dim);
03371
03372
03373 strong_closure_assign();
03374
03375
03376 if (marked_empty())
03377 return;
03378
03379 for (Variables_Set::const_iterator tbu = to_be_unconstrained.begin(),
03380 tbu_end = to_be_unconstrained.end(); tbu != tbu_end; ++tbu)
03381 forget_all_octagonal_constraints(*tbu);
03382
03383 assert(OK());
03384 }
03385
03386 template <typename T>
03387 void
03388 Octagonal_Shape<T>::refine(const Variable var,
03389 const Relation_Symbol relsym,
03390 const Linear_Expression& expr,
03391 Coefficient_traits::const_reference denominator) {
03392 assert(denominator != 0);
03393 const dimension_type expr_space_dim = expr.space_dimension();
03394 assert(space_dim >= expr_space_dim);
03395 const dimension_type var_id = var.id();
03396 assert(var_id <= space_dim);
03397 assert(expr.coefficient(var) == 0);
03398 assert(relsym != LESS_THAN && relsym != GREATER_THAN);
03399
03400 const Coefficient& b = expr.inhomogeneous_term();
03401
03402
03403 dimension_type t = 0;
03404
03405
03406 dimension_type w_id = 0;
03407
03408
03409 for (dimension_type i = expr_space_dim; i-- > 0; )
03410 if (expr.coefficient(Variable(i)) != 0) {
03411 if (t++ == 1)
03412 break;
03413 else
03414 w_id = i;
03415 }
03416
03417
03418
03419
03420
03421 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
03422 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
03423 typedef typename OR_Matrix<N>::const_row_iterator Row_iterator;
03424 typedef typename OR_Matrix<N>::const_row_reference_type Row_reference;
03425
03426 const Row_Iterator m_begin = matrix.row_begin();
03427 const dimension_type n_var = 2*var_id;
03428 TEMP_INTEGER(minus_den);
03429 neg_assign(minus_den, denominator);
03430
03431
03432
03433
03434
03435
03436 if (t == 1 && expr.coefficient(Variable(w_id)) != denominator
03437 && expr.coefficient(Variable(w_id)) != minus_den)
03438 t = 2;
03439
03440 if (t == 0) {
03441
03442 TEMP_INTEGER(two_b);
03443 two_b = 2*b;
03444 switch (relsym) {
03445 case EQUAL:
03446
03447 add_octagonal_constraint(n_var+1, n_var, two_b, denominator);
03448 add_octagonal_constraint(n_var, n_var+1, two_b, minus_den);
03449 break;
03450 case LESS_OR_EQUAL:
03451
03452 add_octagonal_constraint(n_var+1, n_var, two_b, denominator);
03453 break;
03454 case GREATER_OR_EQUAL:
03455
03456
03457 add_octagonal_constraint(n_var, n_var+1, two_b, minus_den);
03458 break;
03459 default:
03460
03461 throw std::runtime_error("PPL internal error");
03462 }
03463 }
03464 else if (t == 1) {
03465
03466 const Coefficient& w_coeff = expr.coefficient(Variable(w_id));
03467 const dimension_type n_w = 2*w_id;
03468 switch (relsym) {
03469 case EQUAL:
03470 if (w_coeff == denominator)
03471
03472 if (var_id < w_id) {
03473 add_octagonal_constraint(n_w, n_var, b, denominator);
03474 add_octagonal_constraint(n_w+1, n_var+1, b, minus_den);
03475 }
03476 else {
03477 add_octagonal_constraint(n_var+1, n_w+1, b, denominator);
03478 add_octagonal_constraint(n_var, n_w, b, minus_den);
03479 }
03480 else
03481
03482 if (var_id < w_id) {
03483 add_octagonal_constraint(n_w+1, n_var, b, denominator);
03484 add_octagonal_constraint(n_w, n_var+1, b, minus_den);
03485 }
03486 else {
03487 add_octagonal_constraint(n_var+1, n_w, b, denominator);
03488 add_octagonal_constraint(n_var, n_w+1, b, minus_den);
03489 }
03490 break;
03491 case LESS_OR_EQUAL:
03492 {
03493 DIRTY_TEMP(N, d);
03494 div_round_up(d, b, denominator);
03495
03496
03497 if (w_coeff == denominator) {
03498
03499 if (var_id < w_id)
03500 add_octagonal_constraint(n_w, n_var, d);
03501 else
03502 add_octagonal_constraint(n_var+1, n_w+1, d);
03503 }
03504 else if (w_coeff == minus_den) {
03505
03506 if (var_id < w_id)
03507 add_octagonal_constraint(n_w+1, n_var, d);
03508 else
03509 add_octagonal_constraint(n_var+1, n_w, d);
03510 }
03511 break;
03512 }
03513
03514 case GREATER_OR_EQUAL:
03515 {
03516 DIRTY_TEMP(N, d);
03517 div_round_up(d, b, minus_den);
03518
03519
03520 if (w_coeff == denominator) {
03521
03522
03523 if (var_id < w_id)
03524 add_octagonal_constraint(n_w+1, n_var+1, d);
03525 else
03526 add_octagonal_constraint(n_var, n_w, d);
03527 }
03528 else if (w_coeff == minus_den) {
03529
03530
03531 if (var_id < w_id)
03532 add_octagonal_constraint(n_w, n_var+1, d);
03533 else
03534 add_octagonal_constraint(n_var, n_w+1, d);
03535 }
03536 break;
03537 }
03538
03539 default:
03540
03541 throw std::runtime_error("PPL internal error");
03542 }
03543 }
03544 else {
03545
03546
03547 const bool is_sc = (denominator > 0);
03548 TEMP_INTEGER(minus_b);
03549 neg_assign(minus_b, b);
03550 const Coefficient& sc_b = is_sc ? b : minus_b;
03551 const Coefficient& minus_sc_b = is_sc ? minus_b : b;
03552 const Coefficient& sc_den = is_sc ? denominator : minus_den;
03553 const Coefficient& minus_sc_den = is_sc ? minus_den : denominator;
03554
03555
03556
03557 Linear_Expression minus_expr;
03558 if (!is_sc)
03559 minus_expr = -expr;
03560 const Linear_Expression& sc_expr = is_sc ? expr : minus_expr;
03561
03562 DIRTY_TEMP(N, sum);
03563
03564 PPL_UNINITIALIZED(dimension_type, pinf_index);
03565
03566 dimension_type pinf_count = 0;
03567
03568 switch (relsym) {
03569 case EQUAL:
03570 {
03571 DIRTY_TEMP(N, neg_sum);
03572
03573 PPL_UNINITIALIZED(dimension_type, neg_pinf_index);
03574
03575 dimension_type neg_pinf_count = 0;
03576
03577
03578 assign_r(sum, sc_b, ROUND_UP);
03579 assign_r(neg_sum, minus_sc_b, ROUND_UP);
03580
03581
03582 DIRTY_TEMP(N, coeff_i);
03583 DIRTY_TEMP(N, half);
03584 TEMP_INTEGER(minus_sc_i);
03585 DIRTY_TEMP(N, minus_coeff_i);
03586
03587
03588 for (Row_iterator m_iter = m_begin, m_iter_end = m_iter + (2*w_id) + 2;
03589 m_iter != m_iter_end; ) {
03590 const dimension_type n_i = m_iter.index();
03591 const dimension_type id = n_i/2;
03592 Row_reference m_i = *m_iter;
03593 ++m_iter;
03594 Row_reference m_ci = *m_iter;
03595 ++m_iter;
03596 const Coefficient& sc_i = sc_expr.coefficient(Variable(id));
03597 const int sign_i = sgn(sc_i);
03598 if (sign_i > 0) {
03599 assign_r(coeff_i, sc_i, ROUND_UP);
03600
03601 if (pinf_count <= 1) {
03602 const N& double_approx_i = m_ci[n_i];
03603 if (!is_plus_infinity(double_approx_i)) {
03604
03605 div2exp_assign_r(half, double_approx_i, 1, ROUND_UP);
03606 add_mul_assign_r(sum, coeff_i, half, ROUND_UP);
03607 }
03608 else {
03609 ++pinf_count;
03610 pinf_index = id;
03611 }
03612 }
03613
03614 if (neg_pinf_count <= 1) {
03615 const N& double_approx_minus_i = m_i[n_i+1];
03616 if (!is_plus_infinity(double_approx_minus_i)) {
03617
03618 div2exp_assign_r(half, double_approx_minus_i, 1, ROUND_UP);
03619 add_mul_assign_r(neg_sum, coeff_i, half, ROUND_UP);
03620 }
03621 else {
03622 ++neg_pinf_count;
03623 neg_pinf_index = id;
03624 }
03625 }
03626 }
03627 else if (sign_i < 0) {
03628 neg_assign_r(minus_sc_i, sc_i, ROUND_NOT_NEEDED);
03629 assign_r(minus_coeff_i, minus_sc_i, ROUND_UP);
03630
03631 if (pinf_count <= 1) {
03632 const N& double_approx_minus_i = m_i[n_i+1];
03633 if (!is_plus_infinity(double_approx_minus_i)) {
03634
03635 div2exp_assign_r(half, double_approx_minus_i, 1, ROUND_UP);
03636 add_mul_assign_r(sum, minus_coeff_i, half, ROUND_UP);
03637 }
03638 else {
03639 ++pinf_count;
03640 pinf_index = id;
03641 }
03642 }
03643
03644 if (neg_pinf_count <= 1) {
03645 const N& double_approx_i = m_ci[n_i];
03646 if (!is_plus_infinity(double_approx_i)) {
03647
03648 div2exp_assign_r(half, double_approx_i, 1, ROUND_UP);
03649 add_mul_assign_r(neg_sum, minus_coeff_i, half, ROUND_UP);
03650 }
03651 else {
03652 ++neg_pinf_count;
03653 neg_pinf_index = id;
03654 }
03655 }
03656 }
03657 }
03658
03659 if (pinf_count > 1 && neg_pinf_count > 1) {
03660 assert(OK());
03661 return;
03662 }
03663
03664
03665 reset_strongly_closed();
03666
03667
03668 if (pinf_count <= 1) {
03669
03670 if (sc_den != 1) {
03671
03672
03673
03674
03675
03676 DIRTY_TEMP(N, down_sc_den);
03677 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03678 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03679 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
03680 }
03681
03682 if (pinf_count == 0) {
03683
03684 DIRTY_TEMP(N, double_sum);
03685 mul2exp_assign_r(double_sum, sum, 1, ROUND_IGNORE);
03686 matrix[n_var+1][n_var] = double_sum;
03687
03688 deduce_v_pm_u_bounds(var_id, w_id, sc_expr, sc_den, sum);
03689 }
03690 else
03691
03692 if (pinf_index != var_id) {
03693 const Coefficient& ppi =
03694 sc_expr.coefficient(Variable(pinf_index));
03695 if (ppi == sc_den)
03696
03697 if (var_id < pinf_index)
03698 matrix[2*pinf_index][n_var] = sum;
03699 else
03700 matrix[n_var+1][2*pinf_index+1] = sum;
03701 else
03702 if (ppi == minus_sc_den) {
03703
03704 if (var_id < pinf_index)
03705 matrix[2*pinf_index+1][n_var] = sum;
03706 else
03707 matrix[n_var+1][2*pinf_index] = sum;
03708 }
03709 }
03710 }
03711
03712
03713 if (neg_pinf_count <= 1) {
03714
03715 if (sc_den != 1) {
03716
03717
03718
03719
03720
03721 DIRTY_TEMP(N, down_sc_den);
03722 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03723 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03724 div_assign_r(neg_sum, neg_sum, down_sc_den, ROUND_UP);
03725 }
03726
03727 if (neg_pinf_count == 0) {
03728
03729 DIRTY_TEMP(N, double_neg_sum);
03730 mul2exp_assign_r(double_neg_sum, neg_sum, 1, ROUND_IGNORE);
03731 matrix[n_var][n_var+1] = double_neg_sum;
03732
03733 deduce_minus_v_pm_u_bounds(var_id, w_id, sc_expr, sc_den, neg_sum);
03734 }
03735 else
03736
03737 if (neg_pinf_index != var_id) {
03738 const Coefficient& npi =
03739 sc_expr.coefficient(Variable(neg_pinf_index));
03740 if (npi == sc_den)
03741
03742
03743 if (neg_pinf_index < var_id)
03744 matrix[n_var][2*neg_pinf_index] = neg_sum;
03745 else
03746 matrix[2*neg_pinf_index+1][n_var+1] = neg_sum;
03747 else
03748 if (npi == minus_sc_den) {
03749
03750
03751 if (neg_pinf_index < var_id)
03752 matrix[n_var][2*neg_pinf_index+1] = neg_sum;
03753 else
03754 matrix[2*neg_pinf_index][n_var+1] = neg_sum;
03755 }
03756 }
03757 }
03758 break;
03759 }
03760
03761 case LESS_OR_EQUAL:
03762 {
03763
03764
03765
03766
03767 assign_r(sum, sc_b, ROUND_UP);
03768
03769
03770 DIRTY_TEMP(N, coeff_i);
03771 DIRTY_TEMP(N, approx_i);
03772 TEMP_INTEGER(minus_sc_i);
03773
03774
03775 for (Row_Iterator m_iter = m_begin, m_end = m_iter + (2*w_id) + 2;
03776 m_iter != m_end; ) {
03777 const dimension_type n_i = m_iter.index();
03778 const dimension_type id = n_i/2;
03779 Row_Reference m_i = *m_iter;
03780 ++m_iter;
03781 Row_Reference m_ci = *m_iter;
03782 ++m_iter;
03783 const Coefficient& sc_i = sc_expr.coefficient(Variable(id));
03784 const int sign_i = sgn(sc_i);
03785 if (sign_i == 0)
03786 continue;
03787
03788 const N& double_approx_i = (sign_i > 0) ? m_ci[n_i] : m_i[n_i+1];
03789 if (is_plus_infinity(double_approx_i)) {
03790 if (++pinf_count > 1)
03791 break;
03792 pinf_index = id;
03793 continue;
03794 }
03795 if (sign_i > 0)
03796 assign_r(coeff_i, sc_i, ROUND_UP);
03797 else {
03798 neg_assign(minus_sc_i, sc_i);
03799 assign_r(coeff_i, minus_sc_i, ROUND_UP);
03800 }
03801 div2exp_assign_r(approx_i, double_approx_i, 1, ROUND_UP);
03802 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
03803 }
03804
03805 if (sc_den != 1) {
03806
03807
03808
03809
03810 DIRTY_TEMP(N, down_sc_den);
03811 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03812 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03813 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
03814 }
03815
03816 if (pinf_count == 0) {
03817
03818 DIRTY_TEMP(N, double_sum);
03819 mul2exp_assign_r(double_sum, sum, 1, ROUND_IGNORE);
03820 add_octagonal_constraint(n_var+1, n_var, double_sum);
03821
03822 deduce_v_pm_u_bounds(var_id, w_id, sc_expr, sc_den, sum);
03823 }
03824 else if (pinf_count == 1) {
03825 dimension_type pinf_ind = 2*pinf_index;
03826 if (expr.coefficient(Variable(pinf_index)) == denominator ) {
03827
03828 if (var_id < pinf_index)
03829 add_octagonal_constraint(pinf_ind, n_var, sum);
03830 else
03831 add_octagonal_constraint(n_var+1, pinf_ind+1, sum);
03832 }
03833 else {
03834 if (expr.coefficient(Variable(pinf_index)) == minus_den) {
03835
03836 if (var_id < pinf_index)
03837 add_octagonal_constraint(pinf_ind+1, n_var, sum);
03838 else
03839 add_octagonal_constraint(n_var+1, pinf_ind, sum);
03840 }
03841 }
03842 }
03843 break;
03844 }
03845
03846 case GREATER_OR_EQUAL:
03847 {
03848
03849
03850
03851
03852
03853 assign_r(sum, minus_sc_b, ROUND_UP);
03854
03855
03856 DIRTY_TEMP(N, coeff_i);
03857 DIRTY_TEMP(N, approx_i);
03858 TEMP_INTEGER(minus_sc_i);
03859 for (Row_Iterator m_iter = m_begin, m_end = m_iter + (2*w_id) + 2;
03860 m_iter != m_end; ) {
03861 const dimension_type n_i = m_iter.index();
03862 const dimension_type id = n_i/2;
03863 Row_Reference m_i = *m_iter;
03864 ++m_iter;
03865 Row_Reference m_ci = *m_iter;
03866 ++m_iter;
03867 const Coefficient& sc_i = sc_expr.coefficient(Variable(id));
03868 const int sign_i = sgn(sc_i);
03869 if (sign_i == 0)
03870 continue;
03871
03872 const N& double_approx_i = (sign_i > 0) ? m_i[n_i+1] : m_ci[n_i];
03873 if (is_plus_infinity(double_approx_i)) {
03874 if (++pinf_count > 1)
03875 break;
03876 pinf_index = id;
03877 continue;
03878 }
03879 if (sign_i > 0)
03880 assign_r(coeff_i, sc_i, ROUND_UP);
03881 else {
03882 neg_assign(minus_sc_i, sc_i);
03883 assign_r(coeff_i, minus_sc_i, ROUND_UP);
03884 }
03885 div2exp_assign_r(approx_i, double_approx_i, 1, ROUND_UP);
03886 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
03887 }
03888
03889
03890 if (sc_den != 1) {
03891
03892
03893
03894
03895 DIRTY_TEMP(N, down_sc_den);
03896 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
03897 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
03898 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
03899 }
03900
03901 if (pinf_count == 0) {
03902
03903 DIRTY_TEMP(N, double_sum);
03904 mul2exp_assign_r(double_sum, sum, 1, ROUND_IGNORE);
03905 add_octagonal_constraint(n_var, n_var+1, double_sum);
03906
03907 deduce_minus_v_pm_u_bounds(var_id, pinf_index, sc_expr, sc_den, sum);
03908 }
03909 else if (pinf_count == 1) {
03910 dimension_type pinf_ind = 2*pinf_index;
03911 if (expr.coefficient(Variable(pinf_index)) == denominator) {
03912
03913
03914 if (pinf_index < var_id)
03915 add_octagonal_constraint(n_var, pinf_ind, sum);
03916 else
03917 add_octagonal_constraint(pinf_ind+1, n_var, sum);
03918 }
03919 else {
03920 if (expr.coefficient(Variable(pinf_index)) == minus_den) {
03921
03922
03923 if (pinf_index < var_id)
03924 add_octagonal_constraint(n_var, pinf_ind+1, sum);
03925 else
03926 add_octagonal_constraint(pinf_ind, n_var+1, sum);
03927 }
03928 }
03929 }
03930 break;
03931 }
03932
03933 default:
03934
03935 throw std::runtime_error("PPL internal error");
03936 }
03937 }
03938 }
03939
03940 template <typename T>
03941 void
03942 Octagonal_Shape<T>::affine_image(const Variable var,
03943 const Linear_Expression& expr,
03944 Coefficient_traits::const_reference
03945 denominator) {
03946
03947 if (denominator == 0)
03948 throw_generic("affine_image(v, e, d)", "d == 0");
03949
03950
03951
03952
03953 const dimension_type expr_space_dim = expr.space_dimension();
03954 if (space_dim < expr_space_dim)
03955 throw_dimension_incompatible("affine_image(v, e, d)", "e", expr);
03956
03957
03958 const dimension_type var_id = var.id();
03959 if (space_dim < var_id + 1)
03960 throw_dimension_incompatible("affine_image(v, e, d)", var.id()+1);
03961
03962 strong_closure_assign();
03963
03964 if (marked_empty())
03965 return;
03966
03967
03968
03969 dimension_type t = 0;
03970
03971 dimension_type w_id = 0;
03972
03973
03974
03975 for (dimension_type i = expr_space_dim; i-- > 0; )
03976 if (expr.coefficient(Variable(i)) != 0) {
03977 if (t++ == 1)
03978 break;
03979 else
03980 w_id = i;
03981 }
03982
03983 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
03984 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
03985 typedef typename OR_Matrix<N>::const_row_iterator Row_iterator;
03986 typedef typename OR_Matrix<N>::const_row_reference_type Row_reference;
03987
03988 const Row_Iterator m_begin = matrix.row_begin();
03989 const Row_Iterator m_end = matrix.row_end();
03990 const dimension_type n_var = 2*var_id;
03991 const Coefficient& b = expr.inhomogeneous_term();
03992 TEMP_INTEGER(minus_den);
03993 neg_assign_r(minus_den, denominator, ROUND_NOT_NEEDED);
03994
03995
03996
03997
03998
03999
04000
04001
04002
04003
04004 if (t == 0) {
04005
04006
04007 forget_all_octagonal_constraints(var_id);
04008 TEMP_INTEGER(two_b);
04009 two_b = 2*b;
04010
04011 add_octagonal_constraint(n_var+1, n_var, two_b, denominator);
04012 add_octagonal_constraint(n_var, n_var+1, two_b, minus_den);
04013 assert(OK());
04014 return;
04015 }
04016
04017 if (t == 1) {
04018
04019 const Coefficient& w_coeff = expr.coefficient(Variable(w_id));
04020 if (w_coeff == denominator || w_coeff == minus_den) {
04021
04022 if (w_id == var_id) {
04023
04024 if (w_coeff == denominator) {
04025 if (b == 0)
04026
04027 return;
04028 else {
04029
04030
04031 DIRTY_TEMP(N, d);
04032 div_round_up(d, b, denominator);
04033 DIRTY_TEMP(N, minus_d);
04034 div_round_up(minus_d, b, minus_den);
04035 Row_Iterator m_iter = m_begin + n_var;
04036 N& m_v_cv = (*m_iter)[n_var+1];
04037 ++m_iter;
04038 N& m_cv_v = (*m_iter)[n_var];
04039 ++m_iter;
04040
04041 for ( ; m_iter != m_end; ++m_iter) {
04042 Row_Reference m_i = *m_iter;
04043 N& m_i_v = m_i[n_var];
04044 add_assign_r(m_i_v, m_i_v, d, ROUND_UP);
04045 N& m_i_cv = m_i[n_var+1];
04046 add_assign_r(m_i_cv, m_i_cv, minus_d, ROUND_UP);
04047 }
04048
04049 mul2exp_assign_r(d, d, 1, ROUND_IGNORE);
04050 add_assign_r(m_cv_v, m_cv_v, d, ROUND_UP);
04051 mul2exp_assign_r(minus_d, minus_d, 1, ROUND_IGNORE);
04052 add_assign_r(m_v_cv, m_v_cv, minus_d, ROUND_UP);
04053 }
04054 reset_strongly_closed();
04055 }
04056
04057 else {
04058
04059
04060 forget_binary_octagonal_constraints(var_id);
04061 Row_Iterator m_iter = m_begin + n_var;
04062 N& m_v_cv = (*m_iter)[n_var+1];
04063 ++m_iter;
04064 N& m_cv_v = (*m_iter)[n_var];
04065
04066 std::swap(m_v_cv, m_cv_v);
04067
04068 reset_strongly_closed();
04069 if (b != 0) {
04070
04071
04072 DIRTY_TEMP(N, d);
04073 div_round_up(d, b, denominator);
04074 DIRTY_TEMP(N, minus_d);
04075 div_round_up(minus_d, b, minus_den);
04076 ++m_iter;
04077 for ( ; m_iter != m_end; ++m_iter) {
04078 Row_Reference m_i = *m_iter;
04079 N& m_i_v = m_i[n_var];
04080 add_assign_r(m_i_v, m_i_v, d, ROUND_UP);
04081 N& m_i_cv = m_i[n_var+1];
04082 add_assign_r(m_i_cv, m_i_cv, minus_d, ROUND_UP);
04083 }
04084 mul2exp_assign_r(d, d, 1, ROUND_IGNORE);
04085 add_assign_r(m_cv_v, m_cv_v, d, ROUND_UP);
04086 mul2exp_assign_r(minus_d, minus_d, 1, ROUND_IGNORE);
04087 add_assign_r(m_v_cv, m_v_cv, minus_d, ROUND_UP);
04088 }
04089 incremental_strong_closure_assign(var);
04090 }
04091 }
04092
04093 else {
04094
04095
04096
04097 forget_all_octagonal_constraints(var_id);
04098 const dimension_type n_w = 2*w_id;
04099
04100 if (w_coeff == denominator)
04101 if (var_id < w_id) {
04102 add_octagonal_constraint(n_w, n_var, b, denominator);
04103 add_octagonal_constraint(n_w+1, n_var+1, b, minus_den);
04104 }
04105 else {
04106 add_octagonal_constraint(n_var+1, n_w+1, b, denominator);
04107 add_octagonal_constraint(n_var, n_w, b, minus_den);
04108 }
04109 else
04110
04111 if (var_id < w_id) {
04112 add_octagonal_constraint(n_w+1, n_var, b, denominator);
04113 add_octagonal_constraint(n_w, n_var+1, b, minus_den);
04114 }
04115 else {
04116 add_octagonal_constraint(n_var+1, n_w, b, denominator);
04117 add_octagonal_constraint(n_var, n_w+1, b, minus_den);
04118 }
04119 incremental_strong_closure_assign(var);
04120 }
04121 assert(OK());
04122 return;
04123 }
04124 }
04125
04126
04127
04128
04129
04130
04131
04132
04133
04134
04135
04136
04137
04138 const bool is_sc = (denominator > 0);
04139 TEMP_INTEGER(minus_b);
04140 neg_assign_r(minus_b, b, ROUND_NOT_NEEDED);
04141
04142 const Coefficient& sc_b = is_sc ? b : minus_b;
04143 const Coefficient& minus_sc_b = is_sc ? minus_b : b;
04144 const Coefficient& sc_den = is_sc ? denominator : minus_den;
04145 const Coefficient& minus_sc_den = is_sc ? minus_den : denominator;
04146
04147
04148
04149 Linear_Expression minus_expr;
04150 if (!is_sc)
04151 minus_expr = -expr;
04152 const Linear_Expression& sc_expr = is_sc ? expr : minus_expr;
04153
04154 DIRTY_TEMP(N, pos_sum);
04155 DIRTY_TEMP(N, neg_sum);
04156
04157 PPL_UNINITIALIZED(dimension_type, pos_pinf_index);
04158 PPL_UNINITIALIZED(dimension_type, neg_pinf_index);
04159
04160 dimension_type pos_pinf_count = 0;
04161 dimension_type neg_pinf_count = 0;
04162
04163
04164 assign_r(pos_sum, sc_b, ROUND_UP);
04165 assign_r(neg_sum, minus_sc_b, ROUND_UP);
04166
04167
04168 DIRTY_TEMP(N, coeff_i);
04169 DIRTY_TEMP(N, minus_coeff_i);
04170 DIRTY_TEMP(N, half);
04171 TEMP_INTEGER(minus_sc_i);
04172
04173
04174 for (Row_iterator m_iter = m_begin, m_iter_end = m_iter + (2*w_id) + 2;
04175 m_iter != m_iter_end; ) {
04176 const dimension_type n_i = m_iter.index();
04177 const dimension_type id = n_i/2;
04178 Row_reference m_i = *m_iter;
04179 ++m_iter;
04180 Row_reference m_ci = *m_iter;
04181 ++m_iter;
04182 const Coefficient& sc_i = sc_expr.coefficient(Variable(id));
04183 const int sign_i = sgn(sc_i);
04184 if (sign_i > 0) {
04185 assign_r(coeff_i, sc_i, ROUND_UP);
04186
04187 if (pos_pinf_count <= 1) {
04188 const N& double_up_approx_i = m_ci[n_i];
04189 if (!is_plus_infinity(double_up_approx_i)) {
04190
04191 div2exp_assign_r(half, double_up_approx_i, 1, ROUND_UP);
04192 add_mul_assign_r(pos_sum, coeff_i, half, ROUND_UP);
04193 }
04194 else {
04195 ++pos_pinf_count;
04196 pos_pinf_index = id;
04197 }
04198 }
04199
04200 if (neg_pinf_count <= 1) {
04201 const N& double_up_approx_minus_i = m_i[n_i+1];
04202 if (!is_plus_infinity(double_up_approx_minus_i)) {
04203
04204 div2exp_assign_r(half, double_up_approx_minus_i, 1, ROUND_UP);
04205 add_mul_assign_r(neg_sum, coeff_i, half, ROUND_UP);
04206 }
04207 else {
04208 ++neg_pinf_count;
04209 neg_pinf_index = id;
04210 }
04211 }
04212 }
04213 else if (sign_i < 0) {
04214 neg_assign_r(minus_sc_i, sc_i, ROUND_NOT_NEEDED);
04215 assign_r(minus_coeff_i, minus_sc_i, ROUND_UP);
04216
04217 if (pos_pinf_count <= 1) {
04218 const N& double_up_approx_minus_i = m_i[n_i+1];
04219 if (!is_plus_infinity(double_up_approx_minus_i)) {
04220
04221 div2exp_assign_r(half, double_up_approx_minus_i, 1, ROUND_UP);
04222 add_mul_assign_r(pos_sum, minus_coeff_i, half, ROUND_UP);
04223 }
04224 else {
04225 ++pos_pinf_count;
04226 pos_pinf_index = id;
04227 }
04228 }
04229
04230 if (neg_pinf_count <= 1) {
04231 const N& double_up_approx_i = m_ci[n_i];
04232 if (!is_plus_infinity(double_up_approx_i)) {
04233
04234 div2exp_assign_r(half, double_up_approx_i, 1, ROUND_UP);
04235 add_mul_assign_r(neg_sum, minus_coeff_i, half, ROUND_UP);
04236 }
04237 else {
04238 ++neg_pinf_count;
04239 neg_pinf_index = id;
04240 }
04241 }
04242 }
04243 }
04244
04245
04246 forget_all_octagonal_constraints(var_id);
04247
04248 if (pos_pinf_count > 1 && neg_pinf_count > 1) {
04249 assert(OK());
04250 return;
04251 }
04252
04253
04254 reset_strongly_closed();
04255
04256
04257 if (pos_pinf_count <= 1) {
04258
04259 if (sc_den != 1) {
04260
04261
04262
04263
04264 DIRTY_TEMP(N, down_sc_den);
04265 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
04266 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
04267 div_assign_r(pos_sum, pos_sum, down_sc_den, ROUND_UP);
04268 }
04269
04270 if (pos_pinf_count == 0) {
04271
04272 DIRTY_TEMP(N, double_pos_sum);
04273 mul2exp_assign_r(double_pos_sum, pos_sum, 1, ROUND_IGNORE);
04274 matrix[n_var+1][n_var] = double_pos_sum;
04275
04276 deduce_v_pm_u_bounds(var_id, w_id, sc_expr, sc_den, pos_sum);
04277 }
04278 else
04279
04280 if (pos_pinf_index != var_id) {
04281 const Coefficient& ppi = sc_expr.coefficient(Variable(pos_pinf_index));
04282 if (ppi == sc_den)
04283
04284 if (var_id < pos_pinf_index)
04285 matrix[2*pos_pinf_index][n_var] = pos_sum;
04286 else
04287 matrix[n_var+1][2*pos_pinf_index+1] = pos_sum;
04288 else
04289 if (ppi == minus_sc_den) {
04290
04291 if (var_id < pos_pinf_index)
04292 matrix[2*pos_pinf_index+1][n_var] = pos_sum;
04293 else
04294 matrix[n_var+1][2*pos_pinf_index] = pos_sum;
04295 }
04296 }
04297 }
04298
04299
04300 if (neg_pinf_count <= 1) {
04301
04302 if (sc_den != 1) {
04303
04304
04305
04306
04307 DIRTY_TEMP(N, down_sc_den);
04308 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
04309 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
04310 div_assign_r(neg_sum, neg_sum, down_sc_den, ROUND_UP);
04311 }
04312
04313 if (neg_pinf_count == 0) {
04314
04315 DIRTY_TEMP(N, double_neg_sum);
04316 mul2exp_assign_r(double_neg_sum, neg_sum, 1, ROUND_IGNORE);
04317 matrix[n_var][n_var+1] = double_neg_sum;
04318
04319 deduce_minus_v_pm_u_bounds(var_id, w_id, sc_expr, sc_den, neg_sum);
04320 }
04321 else
04322
04323 if (neg_pinf_index != var_id) {
04324 const Coefficient& npi = sc_expr.coefficient(Variable(neg_pinf_index));
04325 if (npi == sc_den)
04326
04327
04328 if (neg_pinf_index < var_id)
04329 matrix[n_var][2*neg_pinf_index] = neg_sum;
04330 else
04331 matrix[2*neg_pinf_index+1][n_var+1] = neg_sum;
04332 else
04333 if (npi == minus_sc_den) {
04334
04335
04336 if (neg_pinf_index < var_id)
04337 matrix[n_var][2*neg_pinf_index+1] = neg_sum;
04338 else
04339 matrix[2*neg_pinf_index][n_var+1] = neg_sum;
04340 }
04341 }
04342 }
04343
04344 incremental_strong_closure_assign(var);
04345 assert(OK());
04346 }
04347
04348 template <typename T>
04349 void
04350 Octagonal_Shape<T>::affine_preimage(const Variable var,
04351 const Linear_Expression& expr,
04352 Coefficient_traits::const_reference
04353 denominator) {
04354
04355
04356 if (denominator == 0)
04357 throw_generic("affine_preimage(v, e, d)", "d == 0");
04358
04359
04360
04361
04362 const dimension_type expr_space_dim = expr.space_dimension();
04363 if (space_dim < expr_space_dim)
04364 throw_dimension_incompatible("affine_preimage(v, e, d)", "e", expr);
04365
04366
04367 dimension_type var_id = var.id();
04368 if (space_dim < var_id + 1)
04369 throw_dimension_incompatible("affine_preimage(v, e, d)", var.id()+1);
04370
04371 strong_closure_assign();
04372
04373 if (marked_empty())
04374 return;
04375
04376 const Coefficient& b = expr.inhomogeneous_term();
04377
04378
04379
04380 dimension_type t = 0;
04381
04382
04383 dimension_type w_id = 0;
04384
04385
04386 for (dimension_type i = expr_space_dim; i-- > 0; )
04387 if (expr.coefficient(Variable(i)) != 0) {
04388 if (t++ == 1)
04389 break;
04390 else
04391 w_id = i;
04392 }
04393
04394
04395
04396
04397
04398
04399
04400
04401
04402
04403 if (t == 0) {
04404
04405 forget_all_octagonal_constraints(var_id);
04406 assert(OK());
04407 return;
04408 }
04409
04410 if (t == 1) {
04411
04412 const Coefficient& w_coeff = expr.coefficient(Variable(w_id));
04413 if (w_coeff == denominator || w_coeff == -denominator) {
04414
04415 if (w_id == var_id) {
04416
04417 affine_image(var, denominator*var - b, w_coeff);
04418 }
04419 else {
04420
04421
04422 forget_all_octagonal_constraints(var_id);
04423 assert(OK());
04424 }
04425 return;
04426 }
04427 }
04428
04429
04430
04431
04432 const Coefficient& coeff_v = expr.coefficient(var);
04433 if (coeff_v != 0) {
04434 if (coeff_v > 0) {
04435
04436 Linear_Expression inverse = ((coeff_v + denominator)*var);
04437 inverse -= expr;
04438 affine_image(var, inverse, coeff_v);
04439 }
04440 else {
04441
04442 TEMP_INTEGER(minus_coeff_v);
04443 neg_assign(minus_coeff_v, coeff_v);
04444 Linear_Expression inverse = ((minus_coeff_v - denominator)*var);
04445 inverse += expr;
04446 affine_image(var, inverse, minus_coeff_v);
04447 }
04448 }
04449 else {
04450
04451 forget_all_octagonal_constraints(var_id);
04452 assert(OK());
04453 }
04454 }
04455
04456 template <typename T>
04457 void
04458 Octagonal_Shape<T>
04459 ::generalized_affine_image(const Variable var,
04460 const Relation_Symbol relsym,
04461 const Linear_Expression& expr ,
04462 Coefficient_traits::const_reference denominator) {
04463
04464 if (denominator == 0)
04465 throw_generic("generalized_affine_image(v, r, e, d)", "d == 0");
04466
04467
04468
04469
04470 const dimension_type expr_space_dim = expr.space_dimension();
04471 if (space_dim < expr_space_dim)
04472 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)", "e",
04473 expr);
04474
04475
04476 dimension_type var_id = var.id();
04477 if (space_dim < var_id + 1)
04478 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
04479 var.id()+1);
04480
04481
04482 if (relsym == LESS_THAN || relsym == GREATER_THAN)
04483 throw_generic("generalized_affine_image(v, r, e, d)",
04484 "r is a strict relation symbol and "
04485 "*this is an Octagonal_Shape");
04486
04487 if (relsym == EQUAL) {
04488
04489
04490 affine_image(var, expr, denominator);
04491 return;
04492 }
04493
04494 strong_closure_assign();
04495
04496 if (marked_empty())
04497 return;
04498
04499
04500
04501 dimension_type t = 0;
04502
04503 dimension_type w_id = 0;
04504
04505
04506
04507 for (dimension_type i = expr_space_dim; i-- > 0; )
04508 if (expr.coefficient(Variable(i)) != 0) {
04509 if (t++ == 1)
04510 break;
04511 else
04512 w_id = i;
04513 }
04514
04515 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
04516 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
04517 typedef typename OR_Matrix<N>::const_row_iterator Row_iterator;
04518 typedef typename OR_Matrix<N>::const_row_reference_type Row_reference;
04519
04520 const Row_Iterator m_begin = matrix.row_begin();
04521 const Row_Iterator m_end = matrix.row_end();
04522 const dimension_type n_var = 2*var_id;
04523 const Coefficient& b = expr.inhomogeneous_term();
04524 TEMP_INTEGER(minus_den);
04525 neg_assign_r(minus_den, denominator, ROUND_NOT_NEEDED);
04526
04527
04528
04529
04530
04531
04532
04533
04534
04535
04536 if (t == 0) {
04537
04538 TEMP_INTEGER(two_b);
04539 two_b = 2*b;
04540
04541 forget_all_octagonal_constraints(var_id);
04542
04543 reset_strongly_closed();
04544 switch (relsym) {
04545 case LESS_OR_EQUAL:
04546
04547 add_octagonal_constraint(n_var+1, n_var, two_b, denominator);
04548 break;
04549 case GREATER_OR_EQUAL:
04550
04551
04552 add_octagonal_constraint(n_var, n_var+1, two_b, minus_den);
04553 break;
04554 default:
04555
04556 throw std::runtime_error("PPL internal error");
04557 }
04558 assert(OK());
04559 return;
04560 }
04561
04562 if (t == 1) {
04563
04564 const Coefficient& w_coeff = expr.coefficient(Variable(w_id));
04565 if (w_coeff == denominator || w_coeff == minus_den) {
04566
04567 switch (relsym) {
04568 case LESS_OR_EQUAL:
04569 {
04570 DIRTY_TEMP(N, d);
04571 div_round_up(d, b, denominator);
04572 if (w_id == var_id) {
04573
04574
04575 reset_strongly_closed();
04576 if (w_coeff == denominator) {
04577
04578
04579
04580 Row_Iterator m_iter = m_begin + n_var;
04581 Row_Reference m_v = *m_iter;
04582 N& m_v_cv = m_v[n_var+1];
04583 ++m_iter;
04584 Row_Reference m_cv = *m_iter;
04585 N& m_cv_v = m_cv[n_var];
04586 ++m_iter;
04587
04588 for ( ; m_iter != m_end; ++m_iter) {
04589 Row_Reference m_i = *m_iter;
04590 N& m_i_v = m_i[n_var];
04591 add_assign_r(m_i_v, m_i_v, d, ROUND_UP);
04592 assign_r(m_i[n_var+1], PLUS_INFINITY, ROUND_NOT_NEEDED);
04593 }
04594 for (dimension_type k = n_var; k-- > 0; ) {
04595 assign_r(m_v[k], PLUS_INFINITY, ROUND_NOT_NEEDED);
04596 add_assign_r(m_cv[k], m_cv[k], d, ROUND_UP);
04597 }
04598 mul2exp_assign_r(d, d, 1, ROUND_IGNORE);
04599 add_assign_r(m_cv_v, m_cv_v, d, ROUND_UP);
04600 assign_r(m_v_cv, PLUS_INFINITY, ROUND_NOT_NEEDED);
04601 }
04602 else {
04603
04604
04605 N& m_v_cv = matrix[n_var][n_var+1];
04606 mul2exp_assign_r(d, d, 1, ROUND_IGNORE);
04607 add_assign_r(matrix[n_var+1][n_var], m_v_cv, d, ROUND_UP);
04608 assign_r(m_v_cv, PLUS_INFINITY, ROUND_NOT_NEEDED);
04609 forget_binary_octagonal_constraints(var_id);
04610 }
04611 }
04612 else {
04613
04614
04615
04616 forget_all_octagonal_constraints(var_id);
04617 const dimension_type n_w = 2*w_id;
04618 if (w_coeff == denominator) {
04619
04620 if (var_id < w_id)
04621 add_octagonal_constraint(n_w, n_var, b, denominator);
04622 else
04623 add_octagonal_constraint(n_var+1, n_w+1, b, denominator);
04624 }
04625 else {
04626
04627 if (var_id < w_id)
04628 add_octagonal_constraint(n_w+1, n_var, b, denominator);
04629 else
04630 add_octagonal_constraint(n_var+1, n_w, b, denominator);
04631 }
04632 }
04633 break;
04634 }
04635
04636 case GREATER_OR_EQUAL:
04637 {
04638 DIRTY_TEMP(N, d);
04639 div_round_up(d, b, minus_den);
04640 if (w_id == var_id) {
04641
04642
04643 reset_strongly_closed();
04644 if (w_coeff == denominator) {
04645
04646
04647
04648 Row_Iterator m_iter = m_begin + n_var;
04649 Row_Reference m_v = *m_iter;
04650 N& m_v_cv = m_v[n_var+1];
04651 ++m_iter;
04652 Row_Reference m_cv = *m_iter;
04653 N& m_cv_v = m_cv[n_var];
04654 ++m_iter;
04655
04656 for ( ; m_iter != m_end; ++m_iter) {
04657 Row_Reference m_i = *m_iter;
04658 assign_r(m_i[n_var], PLUS_INFINITY, ROUND_NOT_NEEDED);
04659 add_assign_r(m_i[n_var+1], m_i[n_var+1], d, ROUND_UP);
04660 }
04661 for (dimension_type k = n_var; k-- > 0; ) {
04662 add_assign_r(m_v[k], m_v[k], d, ROUND_UP);
04663 assign_r(m_cv[k], PLUS_INFINITY, ROUND_NOT_NEEDED);
04664 }
04665 mul2exp_assign_r(d, d, 1, ROUND_IGNORE);
04666 add_assign_r(m_v_cv, m_v_cv, d, ROUND_UP);
04667 assign_r(m_cv_v, PLUS_INFINITY, ROUND_NOT_NEEDED);
04668 }
04669 else {
04670
04671
04672 N& m_cv_v = matrix[n_var+1][n_var];
04673 mul2exp_assign_r(d, d, 1, ROUND_IGNORE);
04674 add_assign_r(matrix[n_var][n_var+1], m_cv_v, d, ROUND_UP);
04675 assign_r(m_cv_v, PLUS_INFINITY, ROUND_NOT_NEEDED);
04676 forget_binary_octagonal_constraints(var_id);
04677 }
04678 }
04679 else {
04680
04681
04682
04683 forget_all_octagonal_constraints(var_id);
04684 const dimension_type n_w = 2*w_id;
04685
04686
04687
04688
04689 if (w_coeff == denominator) {
04690
04691
04692 if (var_id < w_id)
04693 add_octagonal_constraint(n_w+1, n_var+1, b, minus_den);
04694 else
04695 add_octagonal_constraint(n_var, n_w, b, minus_den);
04696 }
04697 else {
04698
04699
04700 if (var_id < w_id)
04701 add_octagonal_constraint(n_w, n_var+1, b, minus_den);
04702 else
04703 add_octagonal_constraint(n_var, n_w+1, b, minus_den);
04704 }
04705 }
04706 break;
04707 }
04708
04709 default:
04710
04711 throw std::runtime_error("PPL internal error");
04712 }
04713 assert(OK());
04714 return;
04715 }
04716 }
04717
04718
04719
04720
04721
04722
04723
04724
04725 const bool is_sc = (denominator > 0);
04726 TEMP_INTEGER(minus_b);
04727 neg_assign(minus_b, b);
04728 const Coefficient& sc_b = is_sc ? b : minus_b;
04729 const Coefficient& minus_sc_b = is_sc ? minus_b : b;
04730 const Coefficient& sc_den = is_sc ? denominator : minus_den;
04731 const Coefficient& minus_sc_den = is_sc ? minus_den : denominator;
04732
04733
04734
04735 Linear_Expression minus_expr;
04736 if (!is_sc)
04737 minus_expr = -expr;
04738 const Linear_Expression& sc_expr = is_sc ? expr : minus_expr;
04739
04740 DIRTY_TEMP(N, sum);
04741
04742 PPL_UNINITIALIZED(dimension_type, pinf_index);
04743
04744 dimension_type pinf_count = 0;
04745
04746 switch (relsym) {
04747 case LESS_OR_EQUAL:
04748 {
04749
04750
04751
04752 assign_r(sum, sc_b, ROUND_UP);
04753
04754 DIRTY_TEMP(N, coeff_i);
04755 DIRTY_TEMP(N, approx_i);
04756 TEMP_INTEGER(minus_sc_i);
04757
04758
04759 for (Row_iterator m_iter = m_begin, m_iter_end = m_iter + (2*w_id) + 2;
04760 m_iter != m_iter_end; ) {
04761 const dimension_type n_i = m_iter.index();
04762 const dimension_type id = n_i/2;
04763 Row_reference m_i = *m_iter;
04764 ++m_iter;
04765 Row_reference m_ci = *m_iter;
04766 ++m_iter;
04767 const Coefficient& sc_i = sc_expr.coefficient(Variable(id));
04768 const int sign_i = sgn(sc_i);
04769 if (sign_i == 0)
04770 continue;
04771
04772 const N& double_approx_i = (sign_i > 0) ? m_ci[n_i] : m_i[n_i+1];
04773 if (is_plus_infinity(double_approx_i)) {
04774 if (++pinf_count > 1)
04775 break;
04776 pinf_index = id;
04777 continue;
04778 }
04779 if (sign_i > 0)
04780 assign_r(coeff_i, sc_i, ROUND_UP);
04781 else {
04782 neg_assign(minus_sc_i, sc_i);
04783 assign_r(coeff_i, minus_sc_i, ROUND_UP);
04784 }
04785 div2exp_assign_r(approx_i, double_approx_i, 1, ROUND_UP);
04786 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
04787 }
04788
04789 forget_all_octagonal_constraints(var_id);
04790 reset_strongly_closed();
04791
04792 if (pinf_count > 1) {
04793 assert(OK());
04794 return;
04795 }
04796
04797
04798 if (sc_den != 1) {
04799
04800
04801
04802
04803
04804 DIRTY_TEMP(N, down_sc_den);
04805 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
04806 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
04807 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
04808 }
04809
04810 if (pinf_count == 0) {
04811
04812 DIRTY_TEMP(N, double_sum);
04813 mul2exp_assign_r(double_sum, sum, 1, ROUND_IGNORE);
04814 matrix[n_var+1][n_var] = double_sum;
04815
04816 deduce_v_pm_u_bounds(var_id, w_id, sc_expr, sc_den, sum);
04817 }
04818 else if (pinf_count == 1)
04819 if (pinf_index != var_id) {
04820 const Coefficient& pi = expr.coefficient(Variable(pinf_index));
04821 if (pi == denominator ) {
04822
04823 if (var_id < pinf_index)
04824 matrix[2*pinf_index][n_var] = sum;
04825 else
04826 matrix[n_var+1][2*pinf_index+1] = sum;
04827 }
04828 else {
04829 if (pi == minus_den) {
04830
04831 if (var_id < pinf_index)
04832 matrix[2*pinf_index+1][n_var] = sum;
04833 else
04834 matrix[n_var+1][2*pinf_index] = sum;
04835 }
04836 }
04837 }
04838 break;
04839 }
04840
04841 case GREATER_OR_EQUAL:
04842 {
04843
04844
04845
04846
04847
04848 assign_r(sum, minus_sc_b, ROUND_UP);
04849 DIRTY_TEMP(N, coeff_i);
04850 TEMP_INTEGER(minus_sc_i);
04851 DIRTY_TEMP(N, approx_i);
04852
04853 for (Row_iterator m_iter = m_begin, m_iter_end = m_iter + (2*w_id) + 2;
04854 m_iter != m_iter_end; ) {
04855 const dimension_type n_i = m_iter.index();
04856 const dimension_type id = n_i/2;
04857 Row_reference m_i = *m_iter;
04858 ++m_iter;
04859 Row_reference m_ci = *m_iter;
04860 ++m_iter;
04861 const Coefficient& sc_i = sc_expr.coefficient(Variable(id));
04862 const int sign_i = sgn(sc_i);
04863 if (sign_i == 0)
04864 continue;
04865
04866 const N& double_approx_i = (sign_i > 0) ? m_i[n_i+1] : m_ci[n_i];
04867 if (is_plus_infinity(double_approx_i)) {
04868 if (++pinf_count > 1)
04869 break;
04870 pinf_index = id;
04871 continue;
04872 }
04873 if (sign_i > 0)
04874 assign_r(coeff_i, sc_i, ROUND_UP);
04875 else {
04876 neg_assign(minus_sc_i, sc_i);
04877 assign_r(coeff_i, minus_sc_i, ROUND_UP);
04878 }
04879 div2exp_assign_r(approx_i, double_approx_i, 1, ROUND_UP);
04880 add_mul_assign_r(sum, coeff_i, approx_i, ROUND_UP);
04881 }
04882
04883
04884 forget_all_octagonal_constraints(var_id);
04885 reset_strongly_closed();
04886
04887 if (pinf_count > 1) {
04888 assert(OK());
04889 return;
04890 }
04891
04892
04893 if (sc_den != 1) {
04894
04895
04896
04897
04898
04899 DIRTY_TEMP(N, down_sc_den);
04900 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
04901 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
04902 div_assign_r(sum, sum, down_sc_den, ROUND_UP);
04903 }
04904
04905 if (pinf_count == 0) {
04906
04907 DIRTY_TEMP(N, double_sum);
04908 mul2exp_assign_r(double_sum, sum, 1, ROUND_IGNORE);
04909 matrix[n_var][n_var+1] = double_sum;
04910
04911 deduce_minus_v_pm_u_bounds(var_id, pinf_index, sc_expr, sc_den, sum);
04912 }
04913 else if (pinf_count == 1)
04914 if (pinf_index != var_id) {
04915 const Coefficient& pi = expr.coefficient(Variable(pinf_index));
04916 if (pi == denominator) {
04917
04918
04919 if (pinf_index < var_id)
04920 matrix[n_var][2*pinf_index] = sum;
04921 else
04922 matrix[2*pinf_index+1][n_var+1] = sum;
04923 }
04924 else {
04925 if (pi == minus_den) {
04926
04927
04928 if (pinf_index < var_id)
04929 matrix[n_var][2*pinf_index+1] = sum;
04930 else
04931 matrix[2*pinf_index][n_var+1] = sum;
04932 }
04933 }
04934 }
04935 break;
04936 }
04937
04938 default:
04939
04940 throw std::runtime_error("PPL internal error");
04941 }
04942 incremental_strong_closure_assign(var);
04943 assert(OK());
04944 }
04945
04946 template <typename T>
04947 void
04948 Octagonal_Shape<T>::generalized_affine_image(const Linear_Expression& lhs,
04949 const Relation_Symbol relsym,
04950 const Linear_Expression& rhs) {
04951
04952
04953
04954 dimension_type lhs_space_dim = lhs.space_dimension();
04955 if (space_dim < lhs_space_dim)
04956 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
04957 "e1", lhs);
04958
04959
04960
04961 const dimension_type rhs_space_dim = rhs.space_dimension();
04962 if (space_dim < rhs_space_dim)
04963 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
04964 "e2", rhs);
04965
04966
04967 if (relsym == LESS_THAN || relsym == GREATER_THAN)
04968 throw_generic("generalized_affine_image(e1, r, e2)",
04969 "r is a strict relation symbol and "
04970 "*this is an Octagonal_Shape");
04971
04972 strong_closure_assign();
04973
04974 if (marked_empty())
04975 return;
04976
04977
04978
04979 dimension_type t_lhs = 0;
04980
04981 dimension_type j_lhs = 0;
04982
04983
04984 for (dimension_type i = lhs_space_dim; i-- > 0; )
04985 if (lhs.coefficient(Variable(i)) != 0) {
04986 if (t_lhs++ == 1)
04987 break;
04988 else
04989 j_lhs = i;
04990 }
04991
04992 const Coefficient& b_lhs = lhs.inhomogeneous_term();
04993
04994 if (t_lhs == 0) {
04995
04996
04997
04998
04999
05000
05001
05002
05003 switch (relsym) {
05004 case LESS_OR_EQUAL:
05005 refine_no_check(lhs <= rhs);
05006 break;
05007 case EQUAL:
05008 refine_no_check(lhs == rhs);
05009 break;
05010 case GREATER_OR_EQUAL:
05011 refine_no_check(lhs >= rhs);
05012 break;
05013 default:
05014
05015 throw std::runtime_error("PPL internal error");
05016 }
05017 }
05018
05019 else if (t_lhs == 1) {
05020
05021
05022
05023 Variable v(j_lhs);
05024
05025 const Coefficient& den = lhs.coefficient(v);
05026 Relation_Symbol new_relsym = relsym;
05027 if (den < 0) {
05028 if (relsym == LESS_OR_EQUAL)
05029 new_relsym = GREATER_OR_EQUAL;
05030 else if (relsym == GREATER_OR_EQUAL)
05031 new_relsym = LESS_OR_EQUAL;
05032 }
05033 Linear_Expression expr = rhs - b_lhs;
05034 generalized_affine_image(v, new_relsym, expr, den);
05035 }
05036 else {
05037
05038
05039 bool lhs_vars_intersects_rhs_vars = false;
05040 std::vector<Variable> lhs_vars;
05041 for (dimension_type i = lhs_space_dim; i-- > 0; )
05042 if (lhs.coefficient(Variable(i)) != 0) {
05043 lhs_vars.push_back(Variable(i));
05044 if (rhs.coefficient(Variable(i)) != 0)
05045 lhs_vars_intersects_rhs_vars = true;
05046 }
05047
05048 if (!lhs_vars_intersects_rhs_vars) {
05049
05050
05051 for (dimension_type i = lhs_vars.size(); i-- > 0; ) {
05052 dimension_type lhs_vars_i = lhs_vars[i].id();
05053 forget_all_octagonal_constraints(lhs_vars_i);
05054 }
05055
05056
05057
05058
05059 switch (relsym) {
05060 case LESS_OR_EQUAL:
05061 refine_no_check(lhs <= rhs);
05062 break;
05063 case EQUAL:
05064 refine_no_check(lhs == rhs);
05065 break;
05066 case GREATER_OR_EQUAL:
05067 refine_no_check(lhs >= rhs);
05068 break;
05069 default:
05070
05071 throw std::runtime_error("PPL internal error");
05072 }
05073 }
05074 else {
05075
05076
05077 #if 1 // Simplified computation (see the TODO note below).
05078
05079 for (dimension_type i = lhs_vars.size(); i-- > 0; ) {
05080 dimension_type lhs_vars_i = lhs_vars[i].id();
05081 forget_all_octagonal_constraints(lhs_vars_i);
05082 }
05083
05084 #else // Currently unnecessarily complex computation.
05085
05086
05087
05088
05089
05090 const Variable new_var = Variable(space_dim);
05091 add_space_dimensions_and_embed(1);
05092
05093
05094
05095
05096 affine_image(new_var, rhs);
05097
05098
05099 strong_closure_assign();
05100 assert(!marked_empty());
05101 for (dimension_type i = lhs_vars.size(); i-- > 0; ) {
05102 dimension_type lhs_vars_i = lhs_vars[i].id();
05103 forget_all_octagonal_constraints(lhs_vars_i);
05104 }
05105
05106
05107
05108
05109
05110
05111 switch (relsym) {
05112 case LESS_OR_EQUAL:
05113 refine_no_check(lhs <= new_var);
05114 break;
05115 case EQUAL:
05116 refine_no_check(lhs == new_var);
05117 break;
05118 case GREATER_OR_EQUAL:
05119 refine_no_check(lhs >= new_var);
05120 break;
05121 default:
05122
05123 throw std::runtime_error("PPL internal error");
05124 }
05125
05126 remove_higher_space_dimensions(space_dim-1);
05127 #endif // Currently unnecessarily complex computation.
05128 }
05129 }
05130
05131 assert(OK());
05132 }
05133
05134 template <typename T>
05135 void
05136 Octagonal_Shape<T>::bounded_affine_image(const Variable var,
05137 const Linear_Expression& lb_expr,
05138 const Linear_Expression& ub_expr,
05139 Coefficient_traits::const_reference
05140 denominator) {
05141
05142 if (denominator == 0)
05143 throw_generic("bounded_affine_image(v, lb, ub, d)", "d == 0");
05144
05145
05146 const dimension_type var_id = var.id();
05147 if (space_dim < var_id + 1)
05148 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
05149 var.id()+1);
05150
05151
05152
05153 const dimension_type lb_space_dim = lb_expr.space_dimension();
05154 if (space_dim < lb_space_dim)
05155 throw_dimension_incompatible("bounded_affine_image(v, lb, ub)",
05156 "lb", lb_expr);
05157 const dimension_type ub_space_dim = ub_expr.space_dimension();
05158 if (space_dim < ub_space_dim)
05159 throw_dimension_incompatible("bounded_affine_image(v, lb, ub)",
05160 "ub", ub_expr);
05161
05162 strong_closure_assign();
05163
05164 if (marked_empty())
05165 return;
05166
05167
05168
05169 dimension_type t = 0;
05170
05171 dimension_type w_id = 0;
05172
05173
05174
05175 for (dimension_type i = lb_space_dim; i-- > 0; )
05176 if (lb_expr.coefficient(Variable(i)) != 0) {
05177 if (t++ == 1)
05178 break;
05179 else
05180 w_id = i;
05181 }
05182
05183 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
05184 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
05185 typedef typename OR_Matrix<N>::const_row_iterator Row_iterator;
05186 typedef typename OR_Matrix<N>::const_row_reference_type Row_reference;
05187
05188 const Row_Iterator m_begin = matrix.row_begin();
05189 const Row_Iterator m_end = matrix.row_end();
05190 const dimension_type n_var = 2*var_id;
05191 const Coefficient& b = lb_expr.inhomogeneous_term();
05192 TEMP_INTEGER(minus_den);
05193 neg_assign_r(minus_den, denominator, ROUND_NOT_NEEDED);
05194
05195
05196
05197
05198
05199
05200
05201
05202
05203
05204 if (t == 0) {
05205
05206 generalized_affine_image(var,
05207 LESS_OR_EQUAL,
05208 ub_expr,
05209 denominator);
05210 TEMP_INTEGER(two_b);
05211 two_b = 2*b;
05212
05213 add_octagonal_constraint(n_var, n_var+1, two_b, minus_den);
05214 assert(OK());
05215 return;
05216 }
05217
05218 if (t == 1) {
05219
05220 const Coefficient& w_coeff = lb_expr.coefficient(Variable(w_id));
05221 if (w_coeff == denominator || w_coeff == minus_den) {
05222
05223 if (w_id == var_id) {
05224
05225
05226 const Variable new_var = Variable(space_dim);
05227 add_space_dimensions_and_embed(1);
05228
05229
05230 affine_image(new_var, lb_expr, denominator);
05231
05232 strong_closure_assign();
05233 assert(!marked_empty());
05234
05235 generalized_affine_image(var,
05236 LESS_OR_EQUAL,
05237 ub_expr,
05238 denominator);
05239
05240 refine_no_check(var >= new_var);
05241
05242 remove_higher_space_dimensions(space_dim-1);
05243 return;
05244 }
05245 else {
05246
05247 generalized_affine_image(var,
05248 LESS_OR_EQUAL,
05249 ub_expr,
05250 denominator);
05251
05252
05253 const dimension_type n_w = 2*w_id;
05254
05255 if (w_coeff == denominator)
05256 if (var_id < w_id)
05257 add_octagonal_constraint(n_w+1, n_var+1, b, minus_den);
05258 else
05259 add_octagonal_constraint(n_var, n_w, b, minus_den);
05260 else {
05261
05262 if (var_id < w_id)
05263 add_octagonal_constraint(n_w, n_var+1, b, minus_den);
05264 else
05265 add_octagonal_constraint(n_var, n_w+1, b, minus_den);
05266 }
05267 assert(OK());
05268 return;
05269 }
05270 }
05271 }
05272
05273
05274
05275
05276
05277
05278
05279
05280
05281
05282
05283
05284
05285 const bool is_sc = (denominator > 0);
05286 TEMP_INTEGER(minus_b);
05287 neg_assign_r(minus_b, b, ROUND_NOT_NEEDED);
05288
05289 const Coefficient& minus_sc_b = is_sc ? minus_b : b;
05290 const Coefficient& sc_den = is_sc ? denominator : minus_den;
05291 const Coefficient& minus_sc_den = is_sc ? minus_den : denominator;
05292
05293
05294
05295 Linear_Expression minus_expr;
05296 if (!is_sc)
05297 minus_expr = -lb_expr;
05298 const Linear_Expression& sc_expr = is_sc ? lb_expr : minus_expr;
05299
05300 DIRTY_TEMP(N, neg_sum);
05301
05302 PPL_UNINITIALIZED(dimension_type, neg_pinf_index);
05303
05304 dimension_type neg_pinf_count = 0;
05305
05306
05307 assign_r(neg_sum, minus_sc_b, ROUND_UP);
05308
05309
05310 DIRTY_TEMP(N, coeff_i);
05311 DIRTY_TEMP(N, minus_coeff_i);
05312 DIRTY_TEMP(N, half);
05313 TEMP_INTEGER(minus_sc_i);
05314
05315
05316 for (Row_iterator m_iter = m_begin, m_iter_end = m_iter + (2*w_id) + 2;
05317 m_iter != m_iter_end; ) {
05318 const dimension_type n_i = m_iter.index();
05319 const dimension_type id = n_i/2;
05320 Row_reference m_i = *m_iter;
05321 ++m_iter;
05322 Row_reference m_ci = *m_iter;
05323 ++m_iter;
05324 const Coefficient& sc_i = sc_expr.coefficient(Variable(id));
05325 const int sign_i = sgn(sc_i);
05326 if (sign_i > 0) {
05327 assign_r(coeff_i, sc_i, ROUND_UP);
05328
05329 if (neg_pinf_count <= 1) {
05330 const N& double_up_approx_minus_i = m_i[n_i+1];
05331 if (!is_plus_infinity(double_up_approx_minus_i)) {
05332
05333 div2exp_assign_r(half, double_up_approx_minus_i, 1, ROUND_UP);
05334 add_mul_assign_r(neg_sum, coeff_i, half, ROUND_UP);
05335 }
05336 else {
05337 ++neg_pinf_count;
05338 neg_pinf_index = id;
05339 }
05340 }
05341 }
05342 else if (sign_i < 0) {
05343 neg_assign_r(minus_sc_i, sc_i, ROUND_NOT_NEEDED);
05344 assign_r(minus_coeff_i, minus_sc_i, ROUND_UP);
05345
05346 if (neg_pinf_count <= 1) {
05347 const N& double_up_approx_i = m_ci[n_i];
05348 if (!is_plus_infinity(double_up_approx_i)) {
05349
05350 div2exp_assign_r(half, double_up_approx_i, 1, ROUND_UP);
05351 add_mul_assign_r(neg_sum, minus_coeff_i, half, ROUND_UP);
05352 }
05353 else {
05354 ++neg_pinf_count;
05355 neg_pinf_index = id;
05356 }
05357 }
05358 }
05359 }
05360
05361
05362 generalized_affine_image(var,
05363 LESS_OR_EQUAL,
05364 ub_expr,
05365 denominator);
05366
05367
05368 if (neg_pinf_count > 1) {
05369 return;
05370 }
05371
05372
05373 reset_strongly_closed();
05374
05375
05376 if (neg_pinf_count <= 1) {
05377
05378 if (sc_den != 1) {
05379
05380
05381
05382
05383 DIRTY_TEMP(N, down_sc_den);
05384 assign_r(down_sc_den, minus_sc_den, ROUND_UP);
05385 neg_assign_r(down_sc_den, down_sc_den, ROUND_UP);
05386 div_assign_r(neg_sum, neg_sum, down_sc_den, ROUND_UP);
05387 }
05388
05389 if (neg_pinf_count == 0) {
05390
05391 DIRTY_TEMP(N, double_neg_sum);
05392 mul2exp_assign_r(double_neg_sum, neg_sum, 1, ROUND_IGNORE);
05393 matrix[n_var][n_var+1] = double_neg_sum;
05394
05395 deduce_minus_v_pm_u_bounds(var_id, w_id, sc_expr, sc_den, neg_sum);
05396 }
05397 else
05398
05399 if (neg_pinf_index != var_id) {
05400 const Coefficient& npi = sc_expr.coefficient(Variable(neg_pinf_index));
05401 if (npi == sc_den)
05402
05403
05404 if (neg_pinf_index < var_id)
05405 matrix[n_var][2*neg_pinf_index] = neg_sum;
05406 else
05407 matrix[2*neg_pinf_index+1][n_var+1] = neg_sum;
05408 else
05409 if (npi == minus_sc_den) {
05410
05411
05412 if (neg_pinf_index < var_id)
05413 matrix[n_var][2*neg_pinf_index+1] = neg_sum;
05414 else
05415 matrix[2*neg_pinf_index][n_var+1] = neg_sum;
05416 }
05417 }
05418 }
05419
05420 assert(OK());
05421 }
05422
05423
05424 template <typename T>
05425 void
05426 Octagonal_Shape<T>
05427 ::generalized_affine_preimage(const Variable var,
05428 const Relation_Symbol relsym,
05429 const Linear_Expression& expr,
05430 Coefficient_traits::const_reference
05431 denominator) {
05432
05433 if (denominator == 0)
05434 throw_generic("generalized_affine_preimage(v, r, e, d)", "d == 0");
05435
05436
05437
05438
05439 const dimension_type expr_space_dim = expr.space_dimension();
05440 if (space_dim < expr_space_dim)
05441 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
05442 "e", expr);
05443
05444
05445 const dimension_type var_id = var.id();
05446 if (space_dim < var_id + 1)
05447 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
05448 var.id()+1);
05449
05450
05451 if (relsym == LESS_THAN || relsym == GREATER_THAN)
05452 throw_generic("generalized_affine_preimage(v, r, e, d)",
05453 "r is a strict relation symbol and "
05454 "*this is an Octagonal_Shape");
05455
05456 if (relsym == EQUAL) {
05457
05458
05459 affine_preimage(var, expr, denominator);
05460 return;
05461 }
05462
05463
05464 strong_closure_assign();
05465 if (marked_empty())
05466 return;
05467
05468
05469
05470 const Coefficient& expr_v = expr.coefficient(var);
05471 if (expr_v != 0) {
05472 const Relation_Symbol reversed_relsym = (relsym == LESS_OR_EQUAL)
05473 ? GREATER_OR_EQUAL : LESS_OR_EQUAL;
05474 const Linear_Expression inverse
05475 = expr - (expr_v + denominator)*var;
05476 TEMP_INTEGER(inverse_den);
05477 neg_assign(inverse_den, expr_v);
05478 const Relation_Symbol inverse_relsym
05479 = (sgn(denominator) == sgn(inverse_den)) ? relsym : reversed_relsym;
05480 generalized_affine_image(var, inverse_relsym, inverse, inverse_den);
05481 return;
05482 }
05483
05484
05485
05486
05487
05488 refine(var, relsym, expr, denominator);
05489
05490
05491 if (is_empty())
05492 return;
05493
05494
05495 forget_all_octagonal_constraints(var_id);
05496 assert(OK());
05497 }
05498
05499 template <typename T>
05500 void
05501 Octagonal_Shape<T>
05502 ::generalized_affine_preimage(const Linear_Expression& lhs,
05503 const Relation_Symbol relsym,
05504 const Linear_Expression& rhs) {
05505
05506
05507
05508 dimension_type lhs_space_dim = lhs.space_dimension();
05509 if (space_dim < lhs_space_dim)
05510 throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
05511 "e1", lhs);
05512
05513
05514
05515 const dimension_type rhs_space_dim = rhs.space_dimension();
05516 if (space_dim < rhs_space_dim)
05517 throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
05518 "e2", rhs);
05519
05520
05521 if (relsym == LESS_THAN || relsym == GREATER_THAN)
05522 throw_generic("generalized_affine_preimage(e1, r, e2)",
05523 "r is a strict relation symbol and "
05524 "*this is an Octagonal_Shape");
05525
05526 strong_closure_assign();
05527
05528 if (marked_empty())
05529 return;
05530
05531
05532
05533 dimension_type t_lhs = 0;
05534
05535 dimension_type j_lhs = 0;
05536
05537
05538 for (dimension_type i = lhs_space_dim; i-- > 0; )
05539 if (lhs.coefficient(Variable(i)) != 0) {
05540 if (t_lhs++ == 1)
05541 break;
05542 else
05543 j_lhs = i;
05544 }
05545
05546 const Coefficient& b_lhs = lhs.inhomogeneous_term();
05547
05548
05549
05550 if (t_lhs == 0) {
05551 generalized_affine_image(lhs, relsym, rhs);
05552 return;
05553 }
05554
05555 else if (t_lhs == 1) {
05556
05557
05558
05559 Variable v(j_lhs);
05560
05561 const Coefficient& den = lhs.coefficient(v);
05562 Relation_Symbol new_relsym = relsym;
05563 if (den < 0) {
05564 if (relsym == LESS_OR_EQUAL)
05565 new_relsym = GREATER_OR_EQUAL;
05566 else if (relsym == GREATER_OR_EQUAL)
05567 new_relsym = LESS_OR_EQUAL;
05568 }
05569 Linear_Expression expr = rhs - b_lhs;
05570 generalized_affine_preimage(v, new_relsym, expr, den);
05571 }
05572
05573 else {
05574
05575
05576 bool lhs_vars_intersects_rhs_vars = false;
05577 std::vector<Variable> lhs_vars;
05578 for (dimension_type i = lhs_space_dim; i-- > 0; )
05579 if (lhs.coefficient(Variable(i)) != 0) {
05580 lhs_vars.push_back(Variable(i));
05581 if (rhs.coefficient(Variable(i)) != 0)
05582 lhs_vars_intersects_rhs_vars = true;
05583 }
05584
05585 if (!lhs_vars_intersects_rhs_vars) {
05586
05587
05588
05589
05590
05591 switch (relsym) {
05592 case LESS_OR_EQUAL:
05593 refine_no_check(lhs <= rhs);
05594 break;
05595 case EQUAL:
05596 refine_no_check(lhs == rhs);
05597 break;
05598 case GREATER_OR_EQUAL:
05599 refine_no_check(lhs >= rhs);
05600 break;
05601 default:
05602
05603 throw std::runtime_error("PPL internal error");
05604 }
05605
05606
05607 if (is_empty())
05608 return;
05609
05610 for (dimension_type i = lhs_vars.size(); i-- > 0; ) {
05611 dimension_type lhs_vars_i = lhs_vars[i].id();
05612 forget_all_octagonal_constraints(lhs_vars_i);
05613 }
05614 }
05615 else {
05616
05617
05618
05619
05620
05621
05622 const Variable new_var = Variable(space_dim);
05623 add_space_dimensions_and_embed(1);
05624
05625
05626
05627
05628 affine_image(new_var, lhs);
05629
05630
05631 strong_closure_assign();
05632 assert(!marked_empty());
05633 for (dimension_type i = lhs_vars.size(); i-- > 0; ) {
05634 dimension_type lhs_vars_i = lhs_vars[i].id();
05635 forget_all_octagonal_constraints(lhs_vars_i);
05636 }
05637
05638
05639
05640
05641
05642
05643
05644 switch (relsym) {
05645 case LESS_OR_EQUAL:
05646 refine_no_check(new_var <= rhs);
05647 break;
05648 case EQUAL:
05649 refine_no_check(new_var == rhs);
05650 break;
05651 case GREATER_OR_EQUAL:
05652 refine_no_check(new_var >= rhs);
05653 break;
05654 default:
05655
05656 throw std::runtime_error("PPL internal error");
05657 }
05658
05659 remove_higher_space_dimensions(space_dim-1);
05660 }
05661 }
05662 assert(OK());
05663 }
05664
05665 template <typename T>
05666 void
05667 Octagonal_Shape<T>::bounded_affine_preimage(const Variable var,
05668 const Linear_Expression& lb_expr,
05669 const Linear_Expression& ub_expr,
05670 Coefficient_traits::const_reference
05671 denominator) {
05672
05673 if (denominator == 0)
05674 throw_generic("bounded_affine_preimage(v, lb, ub, d)", "d == 0");
05675
05676
05677 const dimension_type var_id = var.id();
05678 if (space_dim < var_id + 1)
05679 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
05680 var.id()+1);
05681
05682
05683
05684 const dimension_type lb_space_dim = lb_expr.space_dimension();
05685 if (space_dim < lb_space_dim)
05686 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub)",
05687 "lb", lb_expr);
05688 const dimension_type ub_space_dim = ub_expr.space_dimension();
05689 if (space_dim < ub_space_dim)
05690 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub)",
05691 "ub", ub_expr);
05692
05693 strong_closure_assign();
05694
05695 if (marked_empty())
05696 return;
05697
05698 if (ub_expr.coefficient(var) == 0) {
05699 refine(var, LESS_OR_EQUAL, ub_expr, denominator);
05700 generalized_affine_preimage(var, GREATER_OR_EQUAL,
05701 lb_expr, denominator);
05702 return;
05703 }
05704 if (lb_expr.coefficient(var) == 0) {
05705 refine(var, GREATER_OR_EQUAL, lb_expr, denominator);
05706 generalized_affine_preimage(var, LESS_OR_EQUAL,
05707 ub_expr, denominator);
05708 return;
05709 }
05710
05711 const Coefficient& expr_v = lb_expr.coefficient(var);
05712
05713
05714 const Variable new_var = Variable(space_dim);
05715 add_space_dimensions_and_embed(1);
05716 const Linear_Expression lb_inverse
05717 = lb_expr - (expr_v + denominator)*var;
05718 TEMP_INTEGER(inverse_den);
05719 neg_assign(inverse_den, expr_v);
05720 affine_image(new_var, lb_inverse, inverse_den);
05721 strong_closure_assign();
05722 assert(!marked_empty());
05723 generalized_affine_preimage(var, LESS_OR_EQUAL,
05724 ub_expr, denominator);
05725 if (sgn(denominator) == sgn(inverse_den))
05726 refine_no_check(var >= new_var) ;
05727 else
05728 refine_no_check(var <= new_var);
05729
05730 remove_higher_space_dimensions(space_dim-1);
05731 }
05732
05733 template <typename T>
05734 Constraint_System
05735 Octagonal_Shape<T>::constraints() const {
05736 Constraint_System cs;
05737 if (space_dim == 0) {
05738 if (marked_empty())
05739 cs = Constraint_System::zero_dim_empty();
05740 }
05741 else if (marked_empty())
05742 cs.insert(0*Variable(space_dim-1) <= -1);
05743 else {
05744
05745
05746 cs.insert(0*Variable(space_dim-1) <= 0);
05747
05748 typedef typename OR_Matrix<N>::const_row_iterator Row_Iterator;
05749 typedef typename OR_Matrix<N>::const_row_reference_type Row_Reference;
05750
05751 Row_Iterator m_begin = matrix.row_begin();
05752 Row_Iterator m_end = matrix.row_end();
05753
05754 TEMP_INTEGER(a);
05755 TEMP_INTEGER(b);
05756
05757
05758 for (Row_Iterator i_iter = m_begin; i_iter != m_end; ) {
05759 const dimension_type i = i_iter.index();
05760 const Variable x(i/2);
05761 const N& c_i_ii = (*i_iter)[i+1];
05762 ++i_iter;
05763 const N& c_ii_i = (*i_iter)[i];
05764 ++i_iter;
05765
05766 if (is_additive_inverse(c_i_ii, c_ii_i)) {
05767
05768 numer_denom(c_ii_i, b, a);
05769 a *= 2;
05770 cs.insert(a*x == b);
05771 }
05772 else {
05773
05774 if (!is_plus_infinity(c_i_ii)) {
05775 numer_denom(c_i_ii, b, a);
05776 a *= 2;
05777 cs.insert(-a*x <= b);
05778 }
05779 if (!is_plus_infinity(c_ii_i)) {
05780 numer_denom(c_ii_i, b, a);
05781 a *= 2;
05782 cs.insert(a*x <= b);
05783 }
05784 }
05785 }
05786
05787 for (Row_Iterator i_iter = m_begin; i_iter != m_end; ) {
05788 const dimension_type i = i_iter.index();
05789 Row_Reference r_i = *i_iter;
05790 ++i_iter;
05791 Row_Reference r_ii = *i_iter;
05792 ++i_iter;
05793 const Variable y(i/2);
05794 for (dimension_type j = 0; j < i; j += 2) {
05795 const N& c_i_j = r_i[j];
05796 const N& c_ii_jj = r_ii[j+1];
05797 const Variable x(j/2);
05798 if (is_additive_inverse(c_ii_jj, c_i_j)) {
05799
05800 numer_denom(c_i_j, b, a);
05801 cs.insert(a*x - a*y == b);
05802 }
05803 else {
05804
05805 if (!is_plus_infinity(c_i_j)) {
05806 numer_denom(c_i_j, b, a);
05807 cs.insert(a*x - a*y <= b);
05808 }
05809 if (!is_plus_infinity(c_ii_jj)) {
05810 numer_denom(c_ii_jj, b, a);
05811 cs.insert(a*y - a*x <= b);
05812 }
05813 }
05814
05815 const N& c_ii_j = r_ii[j];
05816 const N& c_i_jj = r_i[j+1];
05817 if (is_additive_inverse(c_i_jj, c_ii_j)) {
05818
05819 numer_denom(c_ii_j, b, a);
05820 cs.insert(a*x + a*y == b);
05821 }
05822 else {
05823
05824 if (!is_plus_infinity(c_i_jj)) {
05825 numer_denom(c_i_jj, b, a);
05826 cs.insert(-a*x - a*y <= b);
05827 }
05828 if (!is_plus_infinity(c_ii_j)) {
05829 numer_denom(c_ii_j, b, a);
05830 cs.insert(a*x + a*y <= b);
05831 }
05832 }
05833 }
05834 }
05835 }
05836 return cs;
05837 }
05838
05839 template <typename T>
05840 void
05841 Octagonal_Shape<T>::expand_space_dimension(Variable var, dimension_type m) {
05842
05843 const dimension_type var_id = var.id();
05844 if (var_id+1 > space_dim)
05845 throw_dimension_incompatible("expand_space_dimension(v, m)", var_id+1);
05846
05847
05848
05849 if (m > max_space_dimension() - space_dim)
05850 throw_generic("expand_dimension(v, m)",
05851 "adding m new space dimensions exceeds "
05852 "the maximum allowed space dimension");
05853
05854
05855 if (m == 0)
05856 return;
05857
05858
05859 const dimension_type old_num_rows = matrix.num_rows();
05860
05861
05862 add_space_dimensions_and_embed(m);
05863
05864
05865
05866
05867 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
05868 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
05869 typedef typename OR_Matrix<N>::const_row_iterator Row_iterator;
05870 typedef typename OR_Matrix<N>::const_row_reference_type Row_reference;
05871
05872 const Row_Iterator m_begin = matrix.row_begin();
05873 const Row_Iterator m_end = matrix.row_end();
05874 const dimension_type n_var = 2*var_id;
05875 Row_iterator v_iter = m_begin + n_var;
05876 Row_reference m_v = *v_iter;
05877 Row_reference m_cv = *(v_iter+1);
05878
05879 for (Row_Iterator i_iter = m_begin + old_num_rows; i_iter != m_end;
05880 i_iter += 2) {
05881 Row_Reference m_i = *i_iter;
05882 Row_Reference m_ci = *(i_iter+1);
05883 const dimension_type i = i_iter.index();
05884 const dimension_type ci = i+1;
05885 m_i[ci] = m_v[n_var+1];
05886 m_ci[i] = m_cv[n_var];
05887 for (dimension_type j = 0; j < n_var; ++j) {
05888 m_i[j] = m_v[j];
05889 m_ci[j] = m_cv[j];
05890 }
05891 for (dimension_type j = n_var+2; j < old_num_rows; ++j) {
05892 Row_Iterator j_iter = m_begin + j;
05893 Row_Reference m_j = *j_iter;
05894 Row_Reference m_cj = (j%2) ? *(j_iter-1) : *(j_iter+1);
05895 m_i[j] = m_cj[n_var+1];
05896 m_ci[j] = m_cj[n_var];
05897 }
05898 }
05899
05900
05901 if (marked_strongly_closed())
05902 reset_strongly_closed();
05903 assert(OK());
05904 }
05905
05906 template <typename T>
05907 void
05908 Octagonal_Shape<T>::fold_space_dimensions(const Variables_Set& to_be_folded,
05909 Variable var) {
05910
05911 if (var.space_dimension() > space_dim)
05912 throw_dimension_incompatible("fold_space_dimensions(tbf, v)", "v", var);
05913
05914
05915 if (to_be_folded.empty())
05916 return;
05917
05918
05919 if (to_be_folded.space_dimension() > space_dim)
05920 throw_dimension_incompatible("fold_space_dimensions(tbf, ...)",
05921 to_be_folded.space_dimension());
05922
05923
05924 if (to_be_folded.find(var.id()) != to_be_folded.end())
05925 throw_generic("fold_space_dimensions(tbf, v)",
05926 "v should not occur in tbf");
05927
05928
05929
05930
05931
05932 typedef typename OR_Matrix<N>::row_iterator Row_Iterator;
05933 typedef typename OR_Matrix<N>::row_reference_type Row_Reference;
05934
05935 const Row_Iterator m_begin = matrix.row_begin();
05936 const Row_Iterator m_end = matrix.row_end();
05937
05938 strong_closure_assign();
05939 const dimension_type n_rows = matrix.num_rows();
05940 const dimension_type n_var = 2*var.id();
05941 Row_Iterator v_iter = m_begin + n_var;
05942 Row_Reference m_v = *v_iter;
05943 Row_Reference m_cv = *(v_iter+1);
05944 for (Variables_Set::const_iterator i = to_be_folded.begin(),
05945 tbf_end = to_be_folded.end(); i != tbf_end; ++i) {
05946 const dimension_type tbf_id = *i;
05947 const dimension_type tbf_var = 2*tbf_id;
05948 Row_Iterator tbf_iter = m_begin + tbf_var;
05949 Row_Reference m_tbf = *tbf_iter;
05950 Row_Reference m_ctbf = *(tbf_iter+1);
05951 max_assign(m_v[n_var+1], m_tbf[tbf_var+1]);
05952 max_assign(m_cv[n_var], m_ctbf[tbf_var]);
05953
05954 const dimension_type min_id = std::min(n_var, tbf_var);
05955 const dimension_type max_id = std::max(n_var, tbf_var);
05956
05957 for (dimension_type j = 0; j < min_id; ++j) {
05958 const dimension_type cj = coherent_index(j);
05959 max_assign(m_v[j], m_tbf[j]);
05960 max_assign(m_cv[j], m_ctbf[j]);
05961 max_assign(m_cv[cj], m_ctbf[cj]);
05962 max_assign(m_v[cj], m_tbf[cj]);
05963 }
05964 for (dimension_type j = min_id+2; j < max_id; ++j) {
05965 const dimension_type cj = coherent_index(j);
05966 Row_Iterator j_iter = m_begin + j;
05967 Row_Reference m_j = *j_iter;
05968 Row_Reference m_cj = (j%2) ? *(j_iter-1) : *(j_iter+1);
05969 if (n_var == min_id) {
05970 max_assign(m_cj[n_var+1], m_tbf[j]);
05971 max_assign(m_cj[n_var], m_ctbf[j]);
05972 max_assign(m_j[n_var], m_ctbf[cj]);
05973 max_assign(m_j[n_var+1], m_tbf[cj]);
05974 }
05975 else {
05976 max_assign(m_v[j], m_cj[tbf_var+1]);
05977 max_assign(m_cv[j], m_cj[tbf_var]);
05978 max_assign(m_cv[cj], m_j[tbf_var]);
05979 max_assign(m_v[cj], m_j[tbf_var+1]);
05980 }
05981 }
05982 for (dimension_type j = max_id+2; j < n_rows; ++j) {
05983 Row_Iterator j_iter = m_begin + j;
05984 Row_Reference m_j = *j_iter;
05985 Row_Reference m_cj = (j%2) ? *(j_iter-1) : *(j_iter+1);
05986 max_assign(m_cj[n_var+1], m_cj[tbf_var+1]);
05987 max_assign(m_cj[n_var], m_cj[tbf_var]);
05988 max_assign(m_j[n_var], m_j[tbf_var]);
05989 max_assign(m_j[n_var+1], m_j[tbf_var+1]);
05990 }
05991 }
05992 remove_space_dimensions(to_be_folded);
05993 }
05994
05996 template <typename T>
05997 std::ostream&
05998 IO_Operators::operator<<(std::ostream& s, const Octagonal_Shape<T>& x) {
05999
06000 if (x.marked_empty()) {
06001 s << "false";
06002 return s;
06003 }
06004 if (x.is_universe()) {
06005 s << "true";
06006 return s;
06007 }
06008
06009 typedef typename Octagonal_Shape<T>::coefficient_type N;
06010 typedef typename OR_Matrix<N>::const_row_iterator Row_Iterator;
06011 typedef typename OR_Matrix<N>::const_row_reference_type Row_Reference;
06012
06013
06014 bool first = true;
06015
06016 Row_Iterator m_begin = x.matrix.row_begin();
06017 Row_Iterator m_end = x.matrix.row_end();
06018
06019
06020 DIRTY_TEMP(N, negation);
06021 DIRTY_TEMP(N, half);
06022
06023
06024 for (Row_Iterator i_iter = m_begin; i_iter != m_end; ) {
06025 const dimension_type i = i_iter.index();
06026 const Variable v_i = Variable(i/2);
06027 const N& x_i_ii = (*i_iter)[i+1];
06028 ++i_iter;
06029 const N& x_ii_i = (*i_iter)[i];
06030 ++i_iter;
06031
06032 if (is_additive_inverse(x_i_ii, x_ii_i)) {
06033
06034 assert(!is_plus_infinity(x_i_ii) && !is_plus_infinity(x_ii_i));
06035 if (first)
06036 first = false;
06037 else
06038 s << ", ";
06039
06040
06041 if (div2exp_assign_r(half, x_ii_i, 1, ROUND_UP) == V_EQ)
06042 s << v_i << " == " << half;
06043 else
06044 s << "2*" << v_i << " == " << x_ii_i;
06045 }
06046 else {
06047
06048 if (!is_plus_infinity(x_i_ii)) {
06049 if (first)
06050 first = false;
06051 else
06052 s << ", ";
06053 neg_assign_r(negation, x_i_ii, ROUND_NOT_NEEDED);
06054
06055
06056 if (div2exp_assign_r(half, negation, 1, ROUND_UP) == V_EQ)
06057 s << v_i << " >= " << half;
06058 else
06059 s << "2*" << v_i << " >= " << negation;
06060 }
06061 if (!is_plus_infinity(x_ii_i)) {
06062 if (first)
06063 first = false;
06064 else
06065 s << ", ";
06066
06067
06068 if (div2exp_assign_r(half, x_ii_i, 1, ROUND_UP) == V_EQ)
06069 s << v_i << " <= " << half;
06070 else
06071 s << "2*" << v_i << " <= " << x_ii_i;
06072 }
06073 }
06074 }
06075
06076
06077
06078 for (Row_Iterator i_iter = m_begin; i_iter != m_end; ) {
06079 const dimension_type i = i_iter.index();
06080 const Variable v_i = Variable(i/2);
06081 Row_Reference r_i = *i_iter;
06082 ++i_iter;
06083 Row_Reference r_ii = *i_iter;
06084 ++i_iter;
06085
06086 for (dimension_type j = 0; j < i; j += 2) {
06087 const Variable v_j = Variable(j/2);
06088
06089 const N& x_ii_jj = r_ii[j+1];
06090 const N& x_i_j = r_i[j];
06091
06092 if (is_additive_inverse(x_ii_jj, x_i_j)) {
06093
06094 assert(!is_plus_infinity(x_i_j) && !is_plus_infinity(x_ii_jj));
06095 if (first)
06096 first = false;
06097 else
06098 s << ", ";
06099 if (sgn(x_i_j) >= 0)
06100 s << v_j << " - " << v_i << " == " << x_i_j;
06101 else
06102 s << v_i << " - " << v_j << " == " << x_ii_jj;
06103 }
06104 else {
06105
06106 if (!is_plus_infinity(x_i_j)) {
06107 if (first)
06108 first = false;
06109 else
06110 s << ", ";
06111 if (sgn(x_i_j) >= 0)
06112 s << v_j << " - " << v_i << " <= " << x_i_j;
06113 else {
06114 neg_assign_r(negation, x_i_j, ROUND_DOWN);
06115 s << v_i << " - " << v_j << " >= " << negation;
06116 }
06117 }
06118 if (!is_plus_infinity(x_ii_jj)) {
06119 if (first)
06120 first = false;
06121 else
06122 s << ", ";
06123 if (sgn(x_ii_jj) >= 0)
06124 s << v_i << " - " << v_j << " <= " << x_ii_jj;
06125 else {
06126 neg_assign_r(negation, x_ii_jj, ROUND_DOWN);
06127 s << v_j << " - " << v_i << " >= " << negation;
06128 }
06129 }
06130 }
06131
06132 const N& x_i_jj = r_i[j+1];
06133 const N& x_ii_j = r_ii[j];
06134
06135 if (is_additive_inverse(x_i_jj, x_ii_j)) {
06136
06137 assert(!is_plus_infinity(x_i_jj) && !is_plus_infinity(x_ii_j));
06138 if (first)
06139 first = false;
06140 else
06141 s << ", ";
06142 s << v_j << " + " << v_i << " == " << x_ii_j;
06143 }
06144 else {
06145
06146 if (!is_plus_infinity(x_i_jj)) {
06147 if (first)
06148 first = false;
06149 else
06150 s << ", ";
06151 neg_assign_r(negation, x_i_jj, ROUND_DOWN);
06152 s << v_j << " + " << v_i << " >= " << negation;
06153 }
06154 if (!is_plus_infinity(x_ii_j)) {
06155 if (first)
06156 first = false;
06157 else
06158 s << ", ";
06159 s << v_j << " + " << v_i << " <= " << x_ii_j;
06160 }
06161 }
06162 }
06163 }
06164 return s;
06165 }
06166
06167 template <typename T>
06168 void
06169 Octagonal_Shape<T>::ascii_dump(std::ostream& s) const {
06170 s << "space_dim "
06171 << space_dim
06172 << "\n";
06173 status.ascii_dump(s);
06174 s << "\n";
06175 matrix.ascii_dump(s);
06176 }
06177
06178 PPL_OUTPUT_TEMPLATE_DEFINITIONS(T, Octagonal_Shape<T>)
06179
06180 template <typename T>
06181 bool
06182 Octagonal_Shape<T>::ascii_load(std::istream& s) {
06183 std::string str;
06184
06185 if (!(s >> str) || str != "space_dim")
06186 return false;
06187
06188 if (!(s >> space_dim))
06189 return false;
06190
06191 if (!status.ascii_load(s))
06192 return false;
06193
06194 if (!matrix.ascii_load(s))
06195 return false;
06196
06197 assert(OK());
06198 return true;
06199 }
06200
06201 template <typename T>
06202 memory_size_type
06203 Octagonal_Shape<T>::external_memory_in_bytes() const {
06204 return matrix.external_memory_in_bytes();
06205 }
06206
06207 template <typename T>
06208 bool
06209 Octagonal_Shape<T>::OK() const {
06210
06211 if (!matrix.OK())
06212 return false;
06213
06214
06215 if (!status.OK())
06216 return false;
06217
06218
06219 if (marked_empty())
06220 return true;
06221
06222
06223 if (space_dim == 0)
06224 return true;
06225
06226
06227 for (typename OR_Matrix<N>::const_row_iterator i = matrix.row_begin(),
06228 matrix_row_end = matrix.row_end(); i != matrix_row_end; ++i) {
06229 typename OR_Matrix<N>::const_row_reference_type x_i = *i;
06230 for (dimension_type j = i.row_size(); j-- > 0; )
06231 if (is_minus_infinity(x_i[j])) {
06232 #ifndef NDEBUG
06233 using namespace Parma_Polyhedra_Library::IO_Operators;
06234 std::cerr << "Octagonal_Shape::"
06235 << "matrix[" << i.index() << "][" << j << "] = "
06236 << x_i[j] << "!"
06237 << std::endl;
06238 #endif
06239 return false;
06240 }
06241 }
06242
06243
06244 for (typename OR_Matrix<N>::const_row_iterator i = matrix.row_begin(),
06245 m_end = matrix.row_end(); i != m_end; ++i) {
06246 typename OR_Matrix<N>::const_row_reference_type r = *i;
06247 const N& m_i_i = r[i.index()];
06248 if (!is_plus_infinity(m_i_i)) {
06249 #ifndef NDEBUG
06250 const dimension_type j = i.index();
06251 using namespace Parma_Polyhedra_Library::IO_Operators;
06252 std::cerr << "Octagonal_Shape::matrix[" << j << "][" << j << "] = "
06253 << m_i_i << "! (+inf was expected.)\n";
06254 #endif
06255 return false;
06256 }
06257 }
06258
06259
06260
06261
06262 if (std::numeric_limits<coefficient_type_base>::is_exact) {
06263
06264
06265 if (marked_strongly_closed()) {
06266 Octagonal_Shape x = *this;
06267 x.reset_strongly_closed();
06268 x.strong_closure_assign();
06269 if (x.matrix != matrix) {
06270 #ifndef NDEBUG
06271 std::cerr << "Octagonal_Shape is marked as strongly closed "
06272 << "but it is not!\n";
06273 #endif
06274 return false;
06275 }
06276 }
06277
06278
06279 if (marked_strongly_closed())
06280 if (!is_strong_coherent()) {
06281 #ifndef NDEBUG
06282 std::cerr << "Octagonal_Shape is not strong-coherent!\n";
06283 #endif
06284 return false;
06285 }
06286 }
06287
06288
06289 return true;
06290 }
06291
06292
06293 template <typename T>
06294 void
06295 Octagonal_Shape<T>
06296 ::throw_dimension_incompatible(const char* method,
06297 const Octagonal_Shape& y) const {
06298 std::ostringstream s;
06299 s << "PPL::Octagonal_Shape::" << method << ":\n"
06300 << "this->space_dimension() == " << space_dimension()
06301 << ", y->space_dimension() == " << y.space_dimension() << ".";
06302 throw std::invalid_argument(s.str());
06303 }
06304
06305 template <typename T>
06306 void
06307 Octagonal_Shape<T>
06308 ::throw_dimension_incompatible(const char* method,
06309 dimension_type required_dim) const {
06310 std::ostringstream s;
06311 s << "PPL::Octagonal_Shape::" << method << ":\n"
06312 << "this->space_dimension() == " << space_dimension()
06313 << ", required dimension == " << required_dim << ".";
06314 throw std::invalid_argument(s.str());
06315 }
06316
06317 template <typename T>
06318 void
06319 Octagonal_Shape<T>::throw_dimension_incompatible(const char* method,
06320 const Constraint& c) const {
06321 std::ostringstream s;
06322 s << "PPL::Octagonal_Shape::" << method << ":\n"
06323 << "this->space_dimension() == " << space_dimension()
06324 << ", c->space_dimension == " << c.space_dimension() << ".";
06325 throw std::invalid_argument(s.str());
06326 }
06327
06328 template <typename T>
06329 void
06330 Octagonal_Shape<T>::throw_dimension_incompatible(const char* method,
06331 const Congruence& cg) const {
06332 std::ostringstream s;
06333 s << "PPL::Octagonal_Shape::" << method << ":\n"
06334 << "this->space_dimension() == " << space_dimension()
06335 << ", cg->space_dimension == " << cg.space_dimension() << ".";
06336 throw std::invalid_argument(s.str());
06337 }
06338
06339 template <typename T>
06340 void
06341 Octagonal_Shape<T>::throw_dimension_incompatible(const char* method,
06342 const Generator& g) const {
06343 std::ostringstream s;
06344 s << "PPL::Octagonal_Shape::" << method << ":\n"
06345 << "this->space_dimension() == " << space_dimension()
06346 << ", g->space_dimension == " << g.space_dimension() << ".";
06347 throw std::invalid_argument(s.str());
06348 }
06349
06350 template <typename T>
06351 void
06352 Octagonal_Shape<T>::throw_constraint_incompatible(const char* method) const {
06353 std::ostringstream s;
06354 s << "PPL::Octagonal_Shape::" << method << ":\n"
06355 << "the constraint is incompatible.";
06356 throw std::invalid_argument(s.str());
06357 }
06358
06359 template <typename T>
06360 void
06361 Octagonal_Shape<T>
06362 ::throw_expression_too_complex(const char* method,
06363 const Linear_Expression& e) const {
06364 using namespace IO_Operators;
06365 std::ostringstream s;
06366 s << "PPL::Octagonal_Shape::" << method << ":\n"
06367 << e << " is too complex.";
06368 throw std::invalid_argument(s.str());
06369 }
06370
06371 template <typename T>
06372 void
06373 Octagonal_Shape<T>
06374 ::throw_dimension_incompatible(const char* method,
06375 const char* name_row,
06376 const Linear_Expression& y) const {
06377 std::ostringstream s;
06378 s << "PPL::Octagonal_Shape::" << method << ":\n"
06379 << "this->space_dimension() == " << space_dimension()
06380 << ", " << name_row << "->space_dimension() == "
06381 << y.space_dimension() << ".";
06382 throw std::invalid_argument(s.str());
06383 }
06384
06385
06386 template <typename T>
06387 void
06388 Octagonal_Shape<T>::throw_generic(const char* method,
06389 const char* reason) const {
06390 std::ostringstream s;
06391 s << "PPL::Octagonal_Shape::" << method << ":\n"
06392 << reason << ".";
06393 throw std::invalid_argument(s.str());
06394 }
06395
06396 }
06397
06398 #endif // !defined(PPL_Octagonal_Shape_templates_hh)