00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <config.h>
00024 #include <cassert>
00025
00026 #include "Grid.defs.hh"
00027
00028 namespace Parma_Polyhedra_Library {
00029
00030 #define TRACE(x)
00031
00032
00033 TRACE(using std::endl);
00034 TRACE(using std::cerr);
00035
00036 #ifdef STRONG_REDUCTION
00037 template <typename M, typename R>
00038 void
00039 Grid::reduce_reduced(M& sys, dimension_type dim, dimension_type pivot_index,
00040 dimension_type start, dimension_type end,
00041 Dimension_Kinds& dim_kinds, bool generators) {
00042 R& pivot = sys[pivot_index];
00043
00044 TEMP_INTEGER(pivot_dim);
00045 pivot_dim = pivot[dim];
00046
00047 if (pivot_dim == 0)
00048 return;
00049
00050 TEMP_INTEGER(pivot_dim_half);
00051 pivot_dim_half = (pivot_dim + 1) / 2;
00052 Dimension_Kind row_kind = dim_kinds[dim];
00053 Dimension_Kind line_or_equality, virtual_kind;
00054 int jump;
00055 if (generators) {
00056 line_or_equality = LINE;
00057 virtual_kind = GEN_VIRTUAL;
00058 jump = -1;
00059 } else {
00060 line_or_equality = EQUALITY;
00061 virtual_kind = CON_VIRTUAL;
00062 jump = 1;
00063 }
00064
00065 for (dimension_type row_index = pivot_index, kinds_index = dim + jump;
00066 row_index-- > 0;
00067 kinds_index += jump) {
00068
00069 while (dim_kinds[kinds_index] == virtual_kind)
00070 kinds_index += jump;
00071
00072 if (row_kind == line_or_equality
00073 || (row_kind == PARAMETER
00074 && dim_kinds[kinds_index] == PARAMETER)) {
00075 R& row = sys[row_index];
00076
00077 TEMP_INTEGER(row_dim);
00078 row_dim = row[dim];
00079
00080 TEMP_INTEGER(num_rows_to_subtract);
00081 num_rows_to_subtract = row_dim / pivot_dim;
00082
00083
00084
00085
00086
00087 Coefficient& row_dim_rem = row_dim;
00088 row_dim_rem %= pivot_dim;
00089 if (row_dim_rem < 0) {
00090 if (row_dim_rem <= -pivot_dim_half)
00091 --num_rows_to_subtract;
00092 }
00093 else if (row_dim_rem > 0 && row_dim_rem > pivot_dim_half)
00094 num_rows_to_subtract++;
00095
00096
00097
00098
00099
00100
00101 if (num_rows_to_subtract != 0)
00102 for (dimension_type col = start; col <= end; ++col)
00103 row[col] -= num_rows_to_subtract * pivot[col];
00104 }
00105 }
00106 }
00107 #endif // STRONG_REDUCTION
00108
00109 inline void
00110 Grid::reduce_line_with_line(Grid_Generator& row, Grid_Generator& pivot,
00111 dimension_type column) {
00112 TRACE(cerr << "reduce_line_with_line" << endl);
00113
00114 TEMP_INTEGER(gcd);
00115 gcd_assign(gcd, pivot[column], row[column]);
00116
00117 TEMP_INTEGER(red_pivot_col);
00118 TEMP_INTEGER(red_row_col);
00119 red_pivot_col = pivot[column] / gcd;
00120 red_row_col = row[column] / gcd;
00121
00122
00123 row[column] = 0;
00124 for (dimension_type col = pivot.size() - 2 ;
00125 col > column;
00126 --col)
00127 row[col] = (red_pivot_col * row[col]) - (red_row_col * pivot[col]);
00128 }
00129
00130 inline void
00131 Grid::reduce_equality_with_equality(Congruence& row, Congruence& pivot,
00132 dimension_type column) {
00133 TRACE(cerr << "reduce_equality_with_equality" << endl);
00134
00135 assert(row.modulus() == 0 && pivot.modulus() == 0);
00136
00137 TEMP_INTEGER(gcd);
00138 gcd_assign(gcd, pivot[column], row[column]);
00139
00140 TEMP_INTEGER(red_pivot_col);
00141 TEMP_INTEGER(red_row_col);
00142 red_pivot_col = pivot[column] / gcd;
00143 red_row_col = row[column] / gcd;
00144
00145
00146 row[column] = 0;
00147 for (dimension_type col = 0; col < column; ++col)
00148 row[col] = (red_pivot_col * row[col]) - (red_row_col * pivot[col]);
00149 }
00150
00151 template <typename R>
00152 void
00153 Grid::reduce_pc_with_pc(R& row, R& pivot,
00154 dimension_type column,
00155 dimension_type start,
00156 dimension_type end) {
00157 TEMP_INTEGER(gcd);
00158 TEMP_INTEGER(s);
00159 TEMP_INTEGER(t);
00160 gcdext_assign(gcd, pivot[column], row[column], s, t);
00161
00162 TRACE(cerr << " gcd " << gcd << ", s " << s << ", t " << t << endl);
00163
00164
00165 TEMP_INTEGER(red_pivot_col);
00166 TEMP_INTEGER(red_row_col);
00167 red_pivot_col = pivot[column] / gcd;
00168 red_row_col = row[column] / gcd;
00169 TRACE(cerr << " red_pivot_col " << red_pivot_col
00170 << ", red_row_col " << red_row_col << endl);
00171
00172
00173
00174
00175
00176
00177 assert(pivot.size() > 0);
00178 assert(row.size() > 0);
00179 pivot[column] = gcd;
00180 row[column] = 0;
00181 for (dimension_type col = start; col < end; ++col) {
00182 TEMP_INTEGER(pivot_col);
00183 TEMP_INTEGER(row_col);
00184 pivot_col = pivot[col];
00185 row_col = row[col];
00186 pivot[col] = (s * pivot_col) + (t * row_col);
00187 row[col] = (red_pivot_col * row_col) - (red_row_col * pivot_col);
00188 }
00189 }
00190
00191 void
00192 Grid::reduce_parameter_with_line(Grid_Generator& row,
00193 Grid_Generator& pivot,
00194 dimension_type column,
00195 Grid_Generator_System& sys) {
00196
00197
00198 TRACE(cerr << "reduce_parameter_with_line" << endl);
00199
00200 dimension_type num_cols = sys.num_columns() - 1 ;
00201
00202
00203
00204 if (row[column] == pivot[column]) {
00205 for (dimension_type col = 0; col < num_cols; ++col)
00206 row[col] -= pivot[col];
00207 return;
00208 }
00209
00210 TEMP_INTEGER(gcd);
00211 gcd_assign(gcd, pivot[column], row[column]);
00212
00213 TEMP_INTEGER(red_pivot_col);
00214 TEMP_INTEGER(red_row_col);
00215 red_pivot_col = pivot[column] / gcd;
00216 red_row_col = row[column] / gcd;
00217
00218
00219
00220
00221 #ifdef STRONG_REDUCTION
00222
00223
00224
00225 if (red_pivot_col < 0) {
00226 neg_assign(red_pivot_col);
00227 neg_assign(red_row_col);
00228 }
00229 #endif
00230 for (dimension_type index = 0; index < sys.num_generators(); ++index) {
00231 Grid_Generator& row = sys[index];
00232 if (row.is_parameter_or_point())
00233 for (dimension_type col = 0; col < num_cols; ++col)
00234 row[col] *= red_pivot_col;
00235 }
00236
00237
00238 row[column] = 0;
00239 for (dimension_type col = num_cols - 1; col > column; --col)
00240 row[col] -= red_row_col * pivot[col];
00241 }
00242
00243 void
00244 Grid::reduce_congruence_with_equality(Congruence& row,
00245 Congruence& pivot,
00246 dimension_type column,
00247 Congruence_System& sys) {
00248
00249
00250 TRACE(cerr << "reduce_congruence_with_equality" << endl);
00251 assert(row.modulus() > 0 && pivot.modulus() == 0);
00252
00253 dimension_type num_cols = sys.num_columns();
00254
00255
00256
00257 if (row[column] == pivot[column]) {
00258 for (dimension_type col = 0; col < num_cols; ++col)
00259 row[col] -= pivot[col];
00260 return;
00261 }
00262
00263 TEMP_INTEGER(gcd);
00264 gcd_assign(gcd, pivot[column], row[column]);
00265 TEMP_INTEGER(red_pivot_col);
00266 TEMP_INTEGER(red_row_a);
00267 red_pivot_col = pivot[column] / gcd;
00268 red_row_a = row[column] / gcd;
00269
00270
00271
00272 if (red_pivot_col < 0) {
00273 neg_assign(red_pivot_col);
00274 neg_assign(red_row_a);
00275 }
00276
00277
00278
00279 for (dimension_type index = 0; index < sys.num_rows(); ++index) {
00280 Congruence& row = sys[index];
00281 if (row.is_proper_congruence())
00282 for (dimension_type col = 0; col < num_cols; ++col)
00283 row[col] *= red_pivot_col;
00284 }
00285 --num_cols;
00286 row[column] = 0;
00287
00288
00289 for (dimension_type col = 0; col < column; ++col)
00290 row[col] -= red_row_a * pivot[col];
00291 }
00292
00293 #ifndef NDEBUG
00294 template <typename M, typename R>
00295 bool
00296 Grid::rows_are_zero(M& system, dimension_type first,
00297 dimension_type last, dimension_type row_size) {
00298 while (first <= last) {
00299 R& row = system[first++];
00300 for (dimension_type col = 0; col < row_size; ++col)
00301 if (row[col] != 0)
00302 return false;
00303 }
00304 return true;
00305 }
00306 #endif
00307
00308 void
00309 Grid::simplify(Grid_Generator_System& sys, Dimension_Kinds& dim_kinds) {
00310 TRACE(cerr << "==== simplify (reduce) gs:" << endl);
00311 TRACE(cerr << "sys:" << endl);
00312 TRACE(sys.ascii_dump(cerr));
00313 assert(sys.num_generators() > 0);
00314 assert(sys.num_columns() > 0);
00315
00316
00317
00318
00319 dimension_type num_cols = sys.num_columns() - 1 ;
00320
00321 if (dim_kinds.size() != num_cols)
00322 dim_kinds.resize(num_cols);
00323
00324 dimension_type num_rows = sys.num_generators();
00325 TRACE(cerr << " num_rows " << num_rows << endl);
00326
00327
00328
00329
00330
00331 dimension_type pivot_index = 0;
00332 for (dimension_type dim = 0; dim < num_cols; ++dim) {
00333 TRACE(cerr << "dim " << dim << endl);
00334 trace_dim_kinds(" ", dim_kinds);
00335
00336
00337 dimension_type row_index = pivot_index;
00338 TRACE(cerr << " row_index " << row_index << endl);
00339
00340
00341 while (row_index < num_rows && sys[row_index][dim] == 0) {
00342 TRACE(cerr << " .");
00343 ++row_index;
00344 }
00345 TRACE(cerr << endl);
00346
00347 if (row_index == num_rows) {
00348
00349 TRACE(cerr << " Marking virtual row" << endl);
00350 dim_kinds[dim] = GEN_VIRTUAL;
00351 }
00352 else {
00353 if (row_index != pivot_index)
00354 std::swap(sys[row_index], sys[pivot_index]);
00355 Grid_Generator& pivot = sys[pivot_index];
00356 bool pivot_is_line = pivot.is_line();
00357
00358
00359
00360 TRACE(cerr << " Reducing all following rows" << endl);
00361 while (row_index < num_rows - 1) {
00362 ++row_index;
00363 TRACE(cerr << " row_index " << row_index << endl);
00364
00365 Grid_Generator& row = sys[row_index];
00366
00367 if (row[dim] == 0)
00368 continue;
00369
00370 if (row.is_line())
00371 if (pivot_is_line)
00372 reduce_line_with_line(row, pivot, dim);
00373 else {
00374 assert(pivot.is_parameter_or_point());
00375 std::swap(row, pivot);
00376 pivot_is_line = true;
00377 reduce_parameter_with_line(row, pivot, dim, sys);
00378 }
00379 else {
00380 assert(row.is_parameter_or_point());
00381 if (pivot_is_line)
00382 reduce_parameter_with_line(row, pivot, dim, sys);
00383 else {
00384 assert(pivot.is_parameter_or_point());
00385 reduce_pc_with_pc(row, pivot, dim, dim + 1, num_cols);
00386 }
00387 }
00388 }
00389
00390 if (pivot_is_line)
00391 dim_kinds[dim] = LINE;
00392 else {
00393 assert(pivot.is_parameter_or_point());
00394 dim_kinds[dim] = PARAMETER;
00395 }
00396
00397 #ifdef STRONG_REDUCTION
00398
00399 if (pivot[dim] < 0)
00400 pivot.negate(dim, num_cols - 1);
00401 TRACE(cerr << " rr pivot_index " << pivot_index << endl);
00402 TRACE(sys.ascii_dump(cerr));
00403
00404 reduce_reduced<Grid_Generator_System, Grid_Generator>
00405 (sys, dim, pivot_index, dim, num_cols - 1, dim_kinds);
00406 #endif
00407
00408 ++pivot_index;
00409 }
00410 TRACE(sys.ascii_dump(cerr));
00411 }
00412 trace_dim_kinds("gs simpl end ", dim_kinds);
00413
00414
00415 if (num_rows > pivot_index) {
00416 TRACE(cerr << "clipping trailing" << endl);
00417 #ifndef NDEBUG
00418 bool ret = rows_are_zero<Grid_Generator_System,Grid_Generator>
00419 (sys,
00420 pivot_index,
00421 sys.num_generators() - 1,
00422 sys.num_columns() - 1);
00423 assert(ret == true);
00424 #endif
00425 sys.erase_to_end(pivot_index);
00426 }
00427
00428 sys.unset_pending_rows();
00429
00430
00431
00432 TRACE(cerr << "updating param divisors" << endl);
00433 Coefficient_traits::const_reference system_divisor = sys[0][0];
00434 for (dimension_type row = 1, dim = 1, num_cols = sys.num_columns() - 1;
00435 dim < num_cols;
00436 ++dim)
00437 switch (dim_kinds[dim]) {
00438 case PARAMETER:
00439 sys[row].divisor() = system_divisor;
00440 case LINE:
00441 ++row;
00442 case GEN_VIRTUAL:
00443 break;
00444 }
00445
00446 assert(sys.OK());
00447
00448 TRACE(cerr << "---- simplify (reduce) gs done." << endl);
00449 }
00450
00451 bool
00452 Grid::simplify(Congruence_System& sys, Dimension_Kinds& dim_kinds) {
00453 TRACE(cerr << "======== simplify (reduce) cgs:" << endl);
00454 TRACE(cerr << "sys:" << endl);
00455 TRACE(sys.ascii_dump(cerr));
00456 assert(sys.num_columns() > 2);
00457
00458
00459
00460
00461
00462 sys.normalize_moduli();
00463
00464 dimension_type num_cols = sys.num_columns() - 1 ;
00465
00466 if (dim_kinds.size() != num_cols)
00467 dim_kinds.resize(num_cols);
00468
00469 dimension_type num_rows = sys.num_rows();
00470 TRACE(cerr << " num_rows " << num_rows << endl);
00471
00472
00473
00474
00475
00476 dimension_type pivot_index = 0;
00477 for (dimension_type dim = num_cols; dim-- > 0; ) {
00478 TRACE(cerr << "dim " << dim << endl);
00479 trace_dim_kinds(" ", dim_kinds);
00480
00481
00482 dimension_type row_index = pivot_index;
00483 TRACE(cerr << " row_index " << row_index << endl);
00484
00485
00486 while (row_index < num_rows && sys[row_index][dim] == 0) {
00487 TRACE(cerr << " .");
00488 ++row_index;
00489 }
00490 TRACE(cerr << endl);
00491
00492 if (row_index == num_rows) {
00493
00494
00495 TRACE(cerr << " Marking virtual row" << endl);
00496 dim_kinds[dim] = CON_VIRTUAL;
00497 }
00498 else {
00499
00500 if (row_index != pivot_index)
00501 std::swap(sys[row_index], sys[pivot_index]);
00502 Congruence& pivot = sys[pivot_index];
00503 bool pivot_is_equality = pivot.is_equality();
00504
00505
00506
00507 TRACE(cerr << " Reducing all following rows" << endl);
00508 while (row_index < num_rows - 1) {
00509 ++row_index;
00510 TRACE(cerr << " row_index " << row_index << endl);
00511
00512 Congruence& row = sys[row_index];
00513
00514 if (row[dim] == 0)
00515 continue;
00516
00517 if (row.is_equality())
00518 if (pivot_is_equality)
00519 reduce_equality_with_equality(row, pivot, dim);
00520 else {
00521 assert(pivot.is_proper_congruence());
00522 std::swap(row, pivot);
00523 pivot_is_equality = true;
00524 reduce_congruence_with_equality(row, pivot, dim, sys);
00525 }
00526 else {
00527 assert(row.is_proper_congruence());
00528 if (pivot_is_equality)
00529 reduce_congruence_with_equality(row, pivot, dim, sys);
00530 else {
00531 assert(pivot.is_proper_congruence());
00532 reduce_pc_with_pc(row, pivot, dim, 0, dim);
00533 }
00534 }
00535 }
00536
00537 if (pivot_is_equality)
00538 dim_kinds[dim] = EQUALITY;
00539 else {
00540 assert(pivot.is_proper_congruence());
00541 dim_kinds[dim] = PROPER_CONGRUENCE;
00542 }
00543
00544 #ifdef STRONG_REDUCTION
00545
00546 if (pivot[dim] < 0)
00547 pivot.negate(0, dim);
00548
00549 reduce_reduced<Congruence_System, Congruence>
00550 (sys, dim, pivot_index, 0, dim, dim_kinds, false);
00551 #endif
00552 ++pivot_index;
00553 }
00554 TRACE(sys.ascii_dump(cerr));
00555 }
00556
00557 dimension_type& reduced_num_rows = pivot_index;
00558
00559
00560 if (num_rows > 1 && num_rows > reduced_num_rows) {
00561 TRACE(cerr << "clipping trailing" << endl);
00562 #ifndef NDEBUG
00563 bool ret = rows_are_zero<Congruence_System,Congruence>
00564 (sys,
00565 reduced_num_rows,
00566 num_rows - 1,
00567 num_cols);
00568 assert(ret == true);
00569 #endif
00570
00571
00572 if (reduced_num_rows > 0)
00573 sys.erase_to_end(reduced_num_rows);
00574 else
00575 sys.erase_to_end(1);
00576 }
00577
00578 assert(sys.num_rows() == reduced_num_rows
00579 || (sys.num_rows() == 1 && reduced_num_rows == 0));
00580
00581 if (reduced_num_rows > 0) {
00582
00583
00584 Congruence& last_row = sys[reduced_num_rows - 1];
00585 if (dim_kinds[0] == PROPER_CONGRUENCE) {
00586 if (last_row.inhomogeneous_term() % last_row.modulus() != 0) {
00587
00588 last_row.set_is_equality();
00589 dim_kinds[0] = EQUALITY;
00590 goto return_empty;
00591 }
00592 }
00593 else if (dim_kinds[0] == EQUALITY) {
00594
00595
00596
00597 return_empty:
00598 last_row[0] = 1;
00599 dim_kinds.resize(1);
00600 std::swap(sys.rows[0], sys.rows.back());
00601 sys.erase_to_end(1);
00602
00603 trace_dim_kinds("cgs simpl end ", dim_kinds);
00604 assert(sys.OK());
00605 TRACE(cerr << "---- simplify (reduce) cgs done (empty)." << endl);
00606 return true;
00607 }
00608 }
00609 else if (num_rows > 0) {
00610 assert(sys.num_rows() == 1);
00611
00612
00613 dim_kinds[0] = PROPER_CONGRUENCE;
00614 sys[0][num_cols] = 1;
00615 reduced_num_rows = 1;
00616 }
00617
00618
00619 dimension_type mod_index = num_cols;
00620 if (dim_kinds[0] == CON_VIRTUAL) {
00621
00622 dim_kinds[0] = PROPER_CONGRUENCE;
00623 sys.add_zero_rows(1, Linear_Row::Flags(NECESSARILY_CLOSED,
00624 Linear_Row::RAY_OR_POINT_OR_INEQUALITY));
00625 Congruence& new_last_row = sys[reduced_num_rows];
00626 new_last_row[mod_index] = 1;
00627
00628 dimension_type row_index = reduced_num_rows;
00629 while (row_index-- > 0) {
00630 Congruence& row = sys[row_index];
00631 if (row[mod_index] > 0) {
00632 new_last_row[mod_index] = row[mod_index];
00633 break;
00634 }
00635 }
00636 new_last_row[0] = new_last_row[mod_index];
00637 #ifdef STRONG_REDUCTION
00638 ++reduced_num_rows;
00639 #endif
00640 }
00641 else {
00642 Congruence& last_row = sys[reduced_num_rows - 1];
00643 last_row[0] = last_row[mod_index];
00644 }
00645
00646 #ifdef STRONG_REDUCTION
00647
00648 reduce_reduced<Congruence_System, Congruence>
00649 (sys, 0, reduced_num_rows - 1, 0, 0, dim_kinds, false);
00650 #endif
00651
00652 trace_dim_kinds("cgs simpl end ", dim_kinds);
00653 assert(sys.OK());
00654 TRACE(cerr << "---- simplify (reduce) cgs done." << endl);
00655 return false;
00656 }
00657
00658 #undef TRACE
00659
00660 }