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_Grid_templates_hh
00024 #define PPL_Grid_templates_hh 1
00025
00026 #include "Grid_Generator.defs.hh"
00027 #include "Grid_Generator_System.defs.hh"
00028 #include "Grid_Generator_System.inlines.hh"
00029 #include <algorithm>
00030 #include <deque>
00031
00032 namespace Parma_Polyhedra_Library {
00033
00034 template <typename Interval>
00035 Grid::Grid(const Box<Interval>& box,
00036 Complexity_Class)
00037 : con_sys(),
00038 gen_sys() {
00039 if (box.space_dimension() > max_space_dimension())
00040 throw_space_dimension_overflow("Grid(box, from_bounding_box)",
00041 "the space dimension of box "
00042 "exceeds the maximum allowed "
00043 "space dimension");
00044
00045 space_dim = box.space_dimension();
00046
00047 if (box.is_empty()) {
00048
00049 set_empty();
00050 assert(OK());
00051 return;
00052 }
00053
00054 if (space_dim == 0)
00055 set_zero_dim_univ();
00056 else {
00057
00058 con_sys.increase_space_dimension(space_dim);
00059
00060 TEMP_INTEGER(l_n);
00061 TEMP_INTEGER(l_d);
00062 TEMP_INTEGER(u_n);
00063 TEMP_INTEGER(u_d);
00064 gen_sys.insert(grid_point(0*Variable(space_dim-1)));
00065 Grid_Generator& point = gen_sys[0];
00066 for (dimension_type k = space_dim; k-- > 0; ) {
00067 bool closed = false;
00068
00069 if (box.get_lower_bound(k, closed, l_n, l_d)) {
00070 if (box.get_upper_bound(k, closed, u_n, u_d))
00071 if (l_n * u_d == u_n * l_d) {
00072
00073
00074 con_sys.insert(l_d * Variable(k) == l_n);
00075
00076
00077
00078 const Coefficient& point_divisor = point.divisor();
00079 gcd_assign(u_n, l_d, point_divisor);
00080
00081 exact_div_assign(u_n, point_divisor, u_n);
00082 if (l_d < 0)
00083 neg_assign(u_n);
00084
00085 point.scale_to_divisor(l_d * u_n);
00086
00087 if (l_d < 0)
00088 neg_assign(u_n);
00089
00090 point[k + 1] = l_n * u_n;
00091
00092 continue;
00093 }
00094 }
00095
00096 gen_sys.insert(grid_line(Variable(k)));
00097 }
00098 set_congruences_up_to_date();
00099 set_generators_up_to_date();
00100 gen_sys.unset_pending_rows();
00101 gen_sys.set_sorted(false);
00102 }
00103
00104 assert(OK());
00105 }
00106
00107 template <typename Box>
00108 Grid::Grid(const Box& box, From_Covering_Box)
00109 : con_sys(),
00110 gen_sys() {
00111
00112 if (box.space_dimension() > max_space_dimension())
00113 throw_space_dimension_overflow("Grid(box, from_covering_box)",
00114 "the space dimension of box "
00115 "exceeds the maximum allowed "
00116 "space dimension");
00117
00118 space_dim = box.space_dimension();
00119
00120 TEMP_INTEGER(l_n);
00121 TEMP_INTEGER(l_d);
00122
00123
00124
00125 for (dimension_type k = space_dim; k-- > 0; ) {
00126 bool closed = false;
00127
00128 if (box.get_lower_bound(k, closed, l_n, l_d) && !closed)
00129 throw_invalid_argument("Grid(box, from_covering_box)", "box");
00130 if (box.get_upper_bound(k, closed, l_n, l_d) && !closed)
00131 throw_invalid_argument("Grid(box, from_covering_box)", "box");
00132 }
00133
00134 if (box.is_empty()) {
00135
00136 set_empty();
00137 assert(OK());
00138 return;
00139 }
00140
00141 if (space_dim == 0)
00142 set_zero_dim_univ();
00143 else {
00144
00145 con_sys.increase_space_dimension(space_dim);
00146
00147 TEMP_INTEGER(u_n);
00148 TEMP_INTEGER(u_d);
00149 TEMP_INTEGER(d);
00150 gen_sys.insert(grid_point(0*Variable(space_dim-1)));
00151 Grid_Generator& point = gen_sys[0];
00152 for (dimension_type k = space_dim; k-- > 0; ) {
00153 bool closed = false;
00154
00155 if (box.get_lower_bound(k, closed, l_n, l_d)) {
00156
00157 const Coefficient& point_divisor = point.divisor();
00158 assert(l_d > 0);
00159 assert(point_divisor > 0);
00160
00161 gcd_assign(d, l_d, point_divisor);
00162
00163
00164 exact_div_assign(d, point_divisor, d);
00165
00166 point.scale_to_divisor(l_d * d);
00167
00168
00169 point[k + 1] = l_n * d;
00170
00171 if (box.get_upper_bound(k, closed, u_n, u_d)) {
00172 if (l_n * u_d == u_n * l_d) {
00173
00174
00175 gen_sys.insert(grid_line(Variable(k)));
00176 continue;
00177 }
00178 assert(l_d > 0);
00179 assert(u_d > 0);
00180 gcd_assign(d, l_d, u_d);
00181
00182 exact_div_assign(l_d, l_d, d);
00183 exact_div_assign(d, u_d, d);
00184 l_n *= d;
00185
00186
00187 u_n = (u_n * l_d) - l_n;
00188
00189
00190 l_d *= u_d;
00191
00192 con_sys.insert((l_d * Variable(k) %= l_n) / u_n);
00193 gen_sys.insert(parameter(u_n * Variable(k), l_d));
00194 }
00195 else
00196
00197
00198 con_sys.insert(l_d * Variable(k) == l_n);
00199 }
00200 else
00201 if (box.get_upper_bound(k, closed, u_n, u_d)) {
00202 const Coefficient& point_divisor = point.divisor();
00203 assert(u_d > 0);
00204 assert(point_divisor > 0);
00205
00206 gcd_assign(d, u_d, point_divisor);
00207
00208
00209 exact_div_assign(d, point_divisor, d);
00210
00211 point.scale_to_divisor(u_d * d);
00212
00213
00214 point[k + 1] = u_n * d;
00215
00216
00217 con_sys.insert(u_d * Variable(k) == u_n);
00218 }
00219 else {
00220
00221 set_empty();
00222 assert(OK());
00223 return;
00224 }
00225 }
00226 normalize_divisors(gen_sys);
00227 set_congruences_up_to_date();
00228 set_generators_up_to_date();
00229 gen_sys.set_sorted(false);
00230 gen_sys.unset_pending_rows();
00231 }
00232
00233 assert(OK());
00234 }
00235
00236 template <typename Interval>
00237 void
00238 Grid::get_covering_box(Box<Interval>& box) const {
00239
00240 if (space_dim > box.space_dimension())
00241 throw_dimension_incompatible("get_covering_box(box)", "box",
00242 box.space_dimension());
00243
00244 Box<Interval> new_box(box.space_dimension());
00245
00246 if (marked_empty()) {
00247 box = new_box;
00248 box.set_empty();
00249 return;
00250 }
00251 if (space_dim == 0) {
00252 return;
00253 }
00254 if (!generators_are_up_to_date() && !update_generators()) {
00255
00256 box = new_box;
00257 box.set_empty();
00258 return;
00259 }
00260
00261 assert(!gen_sys.has_no_rows());
00262
00263 dimension_type num_dims = gen_sys.num_columns() - 2 ;
00264 dimension_type num_rows = gen_sys.num_rows();
00265
00266 TEMP_INTEGER(gcd);
00267 TEMP_INTEGER(bound);
00268
00269 if (num_rows > 1) {
00270 Row interval_sizes(num_dims, Row::Flags());
00271 std::vector<bool> interval_emptiness(num_dims, false);
00272
00273
00274
00275
00276
00277 for (dimension_type dim = num_dims; dim-- > 0; )
00278 interval_sizes[dim] = 0;
00279 const Grid_Generator *first_point = NULL;
00280 for (dimension_type row = 0; row < num_rows; ++row) {
00281 Grid_Generator& gen = const_cast<Grid_Generator&>(gen_sys[row]);
00282 if (gen.is_line()) {
00283 for (dimension_type dim = 0; dim < num_dims; ++dim)
00284 if (!interval_emptiness[dim] && gen[dim+1] != 0) {
00285
00286
00287 new_box.add_constraint(Variable(dim) == 0);
00288 interval_emptiness[dim] = true;
00289 }
00290 continue;
00291 }
00292 if (gen.is_point()) {
00293 if (first_point == NULL) {
00294 first_point = &gen_sys[row];
00295 continue;
00296 }
00297 const Grid_Generator& point = *first_point;
00298
00299 dimension_type dim = num_dims;
00300 do {
00301 gen[dim] -= point[dim];
00302 }
00303 while (dim-- > 0);
00304 gen.set_divisor(point.divisor());
00305 }
00306 for (dimension_type dim = num_dims; dim-- > 0; )
00307 if (!interval_emptiness[dim])
00308 gcd_assign(interval_sizes[dim], interval_sizes[dim], gen[dim+1]);
00309 }
00310
00311
00312
00313
00314
00315 const Grid_Generator& point = *first_point;
00316 const Coefficient& divisor = point.divisor();
00317 TEMP_INTEGER(lower_bound);
00318 for (dimension_type dim = num_dims; dim-- > 0; ) {
00319 if (interval_emptiness[dim])
00320 continue;
00321
00322 lower_bound = point[dim+1];
00323
00324
00325
00326 if (interval_sizes[dim] != 0) {
00327
00328
00329 lower_bound %= interval_sizes[dim];
00330
00331
00332
00333 if (lower_bound > 0) {
00334 if (interval_sizes[dim] - lower_bound < lower_bound)
00335 lower_bound -= interval_sizes[dim];
00336 }
00337 else if (lower_bound < 0
00338 && interval_sizes[dim] + lower_bound < - lower_bound)
00339 lower_bound += interval_sizes[dim];
00340
00341
00342 bound = interval_sizes[dim] + lower_bound;
00343 gcd_assign(gcd, bound, divisor);
00344 exact_div_assign(bound, bound, gcd);
00345 exact_div_assign(gcd, divisor, gcd);
00346
00347 new_box.add_constraint(gcd*Variable(dim) <= bound);
00348 }
00349
00350
00351 gcd_assign(gcd, lower_bound, divisor);
00352 exact_div_assign(lower_bound, lower_bound, gcd);
00353 exact_div_assign(gcd, divisor, gcd);
00354
00355 new_box.add_constraint(gcd*Variable(dim) >= lower_bound);
00356 }
00357 }
00358 else {
00359 const Grid_Generator& point = gen_sys[0];
00360 const Coefficient& divisor = point.divisor();
00361
00362 for (dimension_type dim = num_dims; dim-- > 0; ) {
00363
00364 gcd_assign(gcd, point[dim+1], divisor);
00365 exact_div_assign(bound, point[dim+1], gcd);
00366 exact_div_assign(gcd, divisor, gcd);
00367
00368 new_box.add_constraint(gcd*Variable(dim) >= bound);
00369 }
00370 }
00371
00372 box.swap(new_box);
00373 }
00374
00375 template <typename Partial_Function>
00376 void
00377 Grid::map_space_dimensions(const Partial_Function& pfunc) {
00378 if (space_dim == 0)
00379 return;
00380
00381 if (pfunc.has_empty_codomain()) {
00382
00383 if (marked_empty()
00384 || (!generators_are_up_to_date() && !update_generators())) {
00385
00386 space_dim = 0;
00387 set_empty();
00388 }
00389 else
00390
00391 set_zero_dim_univ();
00392
00393 assert(OK());
00394 return;
00395 }
00396
00397 dimension_type new_space_dimension = pfunc.max_in_codomain() + 1;
00398
00399 if (new_space_dimension == space_dim) {
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415 std::vector<dimension_type> cycles;
00416 cycles.reserve(space_dim + space_dim/2);
00417
00418
00419 std::deque<bool> visited(space_dim);
00420
00421 for (dimension_type i = space_dim; i-- > 0; ) {
00422 if (!visited[i]) {
00423 dimension_type j = i;
00424 do {
00425 visited[j] = true;
00426
00427 dimension_type k = 0;
00428 if (!pfunc.maps(j, k))
00429 throw_invalid_argument("map_space_dimensions(pfunc)",
00430 " pfunc is inconsistent");
00431 if (k == j)
00432
00433 goto skip;
00434
00435 cycles.push_back(j+1);
00436
00437 j = k;
00438 } while (!visited[j]);
00439
00440 cycles.push_back(0);
00441 skip:
00442 ;
00443 }
00444 }
00445
00446
00447 if (cycles.empty())
00448 return;
00449
00450
00451 if (congruences_are_up_to_date()) {
00452 con_sys.permute_columns(cycles);
00453 clear_congruences_minimized();
00454 }
00455
00456 if (generators_are_up_to_date()) {
00457 gen_sys.permute_columns(cycles);
00458 clear_generators_minimized();
00459 }
00460
00461 assert(OK());
00462 return;
00463 }
00464
00465
00466
00467
00468 const Grid_Generator_System& old_gensys = grid_generators();
00469
00470 if (old_gensys.has_no_rows()) {
00471
00472 Grid new_grid(new_space_dimension, EMPTY);
00473 std::swap(*this, new_grid);
00474 assert(OK());
00475 return;
00476 }
00477
00478
00479 std::vector<dimension_type> pfunc_maps(space_dim, not_a_dimension());
00480 for (dimension_type j = space_dim; j-- > 0; ) {
00481 dimension_type pfunc_j;
00482 if (pfunc.maps(j, pfunc_j))
00483 pfunc_maps[j] = pfunc_j;
00484 }
00485
00486 Grid_Generator_System new_gensys;
00487
00488 new_gensys.set_sorted(false);
00489
00490 Grid_Generator_System::const_iterator i;
00491 Grid_Generator_System::const_iterator old_gensys_end = old_gensys.end();
00492 for (i = old_gensys.begin(); i != old_gensys_end; ++i)
00493 if (i->is_point())
00494 break;
00495 assert(i != old_gensys_end);
00496 const Coefficient& system_divisor = i->divisor();
00497 for (i = old_gensys.begin(); i != old_gensys_end; ++i) {
00498 const Grid_Generator& old_g = *i;
00499 Linear_Expression e(0 * Variable(new_space_dimension-1));
00500 bool all_zeroes = true;
00501 for (dimension_type j = space_dim; j-- > 0; ) {
00502 if (old_g.coefficient(Variable(j)) != 0
00503 && pfunc_maps[j] != not_a_dimension()) {
00504 e += Variable(pfunc_maps[j]) * old_g.coefficient(Variable(j));
00505 all_zeroes = false;
00506 }
00507 }
00508 switch (old_g.type()) {
00509 case Grid_Generator::LINE:
00510 if (!all_zeroes)
00511 new_gensys.insert(grid_line(e));
00512 break;
00513 case Grid_Generator::PARAMETER:
00514 if (!all_zeroes)
00515 new_gensys.insert(parameter(e, system_divisor));
00516 break;
00517 case Grid_Generator::POINT:
00518 new_gensys.insert(grid_point(e, old_g.divisor()));
00519 break;
00520 case Grid_Generator::CLOSURE_POINT:
00521 default:
00522 assert(0);
00523 }
00524 }
00525
00526 Grid new_grid(new_gensys);
00527 std::swap(*this, new_grid);
00528
00529 assert(OK(true));
00530 }
00531
00532 #ifdef STRONG_REDUCTION
00533 template <typename M, typename R>
00534 void
00535 Grid::reduce_reduced(M& sys,
00536 const dimension_type dim,
00537 const dimension_type pivot_index,
00538 const dimension_type start,
00539 const dimension_type end,
00540 const Dimension_Kinds& dim_kinds,
00541 const bool generators) {
00542 R& pivot = sys[pivot_index];
00543
00544 const Coefficient& pivot_dim = pivot[dim];
00545
00546 if (pivot_dim == 0)
00547 return;
00548
00549 TEMP_INTEGER(pivot_dim_half);
00550 pivot_dim_half = (pivot_dim + 1) / 2;
00551 Dimension_Kind row_kind = dim_kinds[dim];
00552 Dimension_Kind line_or_equality, virtual_kind;
00553 int jump;
00554 if (generators) {
00555 line_or_equality = LINE;
00556 virtual_kind = GEN_VIRTUAL;
00557 jump = -1;
00558 }
00559 else {
00560 line_or_equality = EQUALITY;
00561 virtual_kind = CON_VIRTUAL;
00562 jump = 1;
00563 }
00564
00565 TEMP_INTEGER(num_rows_to_subtract);
00566 TEMP_INTEGER(row_dim_remainder);
00567 for (dimension_type row_index = pivot_index, kinds_index = dim + jump;
00568 row_index-- > 0;
00569 kinds_index += jump) {
00570
00571 while (dim_kinds[kinds_index] == virtual_kind)
00572 kinds_index += jump;
00573
00574
00575 if (row_kind == line_or_equality
00576 || (row_kind == PARAMETER
00577 && dim_kinds[kinds_index] == PARAMETER)) {
00578 R& row = sys[row_index];
00579
00580 const Coefficient& row_dim = row[dim];
00581
00582 num_rows_to_subtract = row_dim / pivot_dim;
00583
00584
00585
00586
00587
00588 row_dim_remainder = row_dim % pivot_dim;
00589 if (row_dim_remainder < 0) {
00590 if (row_dim_remainder <= -pivot_dim_half)
00591 --num_rows_to_subtract;
00592 }
00593 else if (row_dim_remainder > 0 && row_dim_remainder > pivot_dim_half)
00594 ++num_rows_to_subtract;
00595
00596
00597
00598
00599
00600
00601 if (num_rows_to_subtract != 0)
00602 for (dimension_type col = start; col <= end; ++col)
00603 sub_mul_assign(row[col], num_rows_to_subtract, pivot[col]);
00604 }
00605 }
00606 }
00607 #endif // STRONG_REDUCTION
00608
00609 }
00610
00611 #endif // !defined(PPL_Grid_templates_hh)