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