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