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
00025 #include "Grid.defs.hh"
00026 #include <cstddef>
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
00037
00038
00039
00040
00041
00042 bool
00043 Grid::lower_triangular(const Congruence_System& sys,
00044 const Dimension_Kinds& dim_kinds) {
00045 dimension_type num_cols = sys.num_columns() - 1;
00046 dimension_type row = sys.num_rows();
00047
00048
00049 if (row > num_cols)
00050 return false;
00051
00052
00053 for (dimension_type dim = 0; dim < num_cols; ++dim) {
00054 if (dim_kinds[dim] == CON_VIRTUAL)
00055 continue;
00056 const Congruence& cg = sys[--row];
00057
00058 if (cg[dim] <= 0)
00059 return false;
00060
00061 dimension_type col = dim;
00062 while (++col < num_cols)
00063 if (cg[col] != 0)
00064 return false;
00065 }
00066
00067
00068 return row == 0;
00069 }
00070
00071
00072
00073
00074
00075
00076
00077 bool
00078 Grid::upper_triangular(const Grid_Generator_System& sys,
00079 const Dimension_Kinds& dim_kinds) {
00080 dimension_type num_cols = sys.space_dimension() + 1;
00081 dimension_type row = sys.num_generators();
00082
00083
00084 if (row > num_cols)
00085 return false;
00086
00087
00088 while (num_cols > 0) {
00089 --num_cols;
00090 if (dim_kinds[num_cols] == GEN_VIRTUAL)
00091 continue;
00092 const Grid_Generator& gen = sys[--row];
00093
00094 if (gen[num_cols] <= 0)
00095 return false;
00096
00097 dimension_type col = num_cols;
00098 while (col-- > 0)
00099 if (gen[col] != 0)
00100 return false;
00101 }
00102
00103
00104 return num_cols == row;
00105 }
00106
00107 inline void
00108 Grid::multiply_grid(const Coefficient& multiplier, Grid_Generator& gen,
00109 Grid_Generator_System& dest, const dimension_type num_rows,
00110 const dimension_type num_dims) {
00111 if (multiplier == 1)
00112 return;
00113
00114 if (gen.is_line())
00115
00116 for (dimension_type column = 0; column < num_dims; ++column)
00117 gen[column] *= multiplier;
00118 else {
00119 assert(gen.is_parameter_or_point());
00120
00121 for (dimension_type index = 0; index < num_rows; ++index) {
00122 Grid_Generator& generator = dest[index];
00123 if (generator.is_parameter_or_point())
00124 for (dimension_type column = 0; column < num_dims; ++column)
00125 generator[column] *= multiplier;
00126 }
00127 }
00128 }
00129
00130 inline void
00131 Grid::multiply_grid(const Coefficient& multiplier, Congruence& cg,
00132 Congruence_System& dest, const dimension_type num_rows,
00133 const dimension_type num_dims) {
00134 if (multiplier == 1)
00135 return;
00136
00137 if (cg.is_proper_congruence())
00138
00139 for (dimension_type index = 0; index < num_rows; ++index) {
00140 Congruence& congruence = dest[index];
00141 if (congruence.is_proper_congruence())
00142 for (dimension_type column = 0; column < num_dims; ++column)
00143 congruence[column] *= multiplier;
00144 }
00145 else {
00146 assert(cg.is_equality());
00147
00148 for (dimension_type column = 0; column < num_dims; ++column)
00149 cg[column] *= multiplier;
00150 }
00151 }
00152
00153
00154
00155
00156
00157 void
00158 Grid::conversion(Grid_Generator_System& source, Congruence_System& dest,
00159 Dimension_Kinds& dim_kinds) {
00160 TRACE(cerr << "============= convert gs to cgs" << endl);
00161 TRACE(cerr << "source:" << endl);
00162 TRACE(source.ascii_dump(cerr));
00163 TRACE(cerr << "dest:" << endl);
00164 TRACE(dest.ascii_dump(cerr));
00165 trace_dim_kinds("gs to cgs ", dim_kinds);
00166
00167
00168
00169
00170 assert(upper_triangular(source, dim_kinds));
00171
00172
00173
00174
00175
00176
00177 dimension_type source_num_rows = 0;
00178
00179 dimension_type dest_num_rows = 0;
00180 TEMP_INTEGER(diagonal_lcm);
00181 diagonal_lcm = 1;
00182 const dimension_type dims = source.space_dimension() + 1;
00183 for (dimension_type dim = 0; dim < dims; ++dim)
00184 if (dim_kinds[dim] == GEN_VIRTUAL)
00185
00186 ++dest_num_rows;
00187 else {
00188 if (dim_kinds[dim] == PARAMETER) {
00189
00190
00191
00192 lcm_assign(diagonal_lcm, diagonal_lcm, source[source_num_rows][dim]);
00193
00194 ++dest_num_rows;
00195 }
00196
00197 ++source_num_rows;
00198 }
00199 TRACE(cerr << "diagonal_lcm: " << diagonal_lcm << endl);
00200 TRACE(cerr << "source_num_rows: " << source_num_rows << endl);
00201 TRACE(cerr << "dest_num_rows: " << dest_num_rows << endl);
00202
00203
00204 if (diagonal_lcm == 0)
00205 throw std::runtime_error("PPL internal error: Grid::conversion:"
00206 " source matrix is singular.");
00207
00208 dest.resize_no_copy(dest_num_rows, dims + 1 );
00209
00210
00211
00212
00213 dimension_type source_index = 0, dest_index = dest_num_rows - 1;
00214 for (dimension_type dim = 0; dim < dims; ++dim) {
00215 TRACE(cerr << "init dim " << dim << endl);
00216 if (dim_kinds[dim] == LINE) {
00217 TRACE(cerr << " line" << endl);
00218 ++source_index;
00219 }
00220 else {
00221 Congruence& cg = dest[dest_index];
00222 for (dimension_type j = 0; j < dim; j++)
00223 cg[j] = 0;
00224 for (dimension_type j = dim + 1; j < dims; j++)
00225 cg[j] = 0;
00226
00227 if (dim_kinds[dim] == GEN_VIRTUAL) {
00228 TRACE(cerr << " gen_virtual" << endl);
00229 cg[dims] = 0;
00230 cg[dim] = 1;
00231 }
00232 else {
00233 assert(dim_kinds[dim] == PARAMETER);
00234 TRACE(cerr << " parameter" << endl);
00235 cg[dims] = 1;
00236 cg[dim] = diagonal_lcm / source[source_index][dim];
00237 ++source_index;
00238 }
00239 --dest_index;
00240 }
00241 }
00242
00243 assert(lower_triangular(dest, dim_kinds));
00244
00245 TRACE(cerr << "dest after init:" << endl);
00246 TRACE(dest.ascii_dump(cerr));
00247
00248
00249
00250
00251
00252
00253
00254
00255 source_index = source_num_rows;
00256 dest_index = 0;
00257
00258 for (dimension_type dim = dims; dim-- > 0; ) {
00259 TRACE(cerr << "dim: " << dim << endl);
00260
00261 if (dim_kinds[dim] != GEN_VIRTUAL) {
00262 --source_index;
00263 TEMP_INTEGER(source_dim);
00264 source_dim = source[source_index][dim];
00265
00266
00267
00268 for (dimension_type row = dest_index; row-- > 0; ) {
00269 TRACE(cerr << " row " << row << endl);
00270 TRACE(dest.ascii_dump(cerr));
00271
00272 Congruence& cg = dest[row];
00273
00274
00275
00276
00277 TEMP_INTEGER(multiplier);
00278 gcd_assign(multiplier, cg[dim], source_dim);
00279 multiplier = source_dim / multiplier;
00280 multiply_grid(multiplier, cg, dest, dest_num_rows, dims);
00281
00282 cg[dim] /= source_dim;
00283 }
00284 TRACE(cerr << "dest after dividing grid:" << endl);
00285 TRACE(dest.ascii_dump(cerr));
00286 }
00287
00288
00289
00290
00291
00292
00293
00294 dimension_type tem_source_index = source_index;
00295 if (dim_kinds[dim] != LINE)
00296 ++dest_index;
00297 for (dimension_type dim_prec = dim; dim_prec-- > 0; ) {
00298 TRACE(cerr << " dim_prec: " << dim_prec);
00299 TRACE(cerr << " dest_index: " << dest_index);
00300 TRACE(cerr << " tem_source_index: " << tem_source_index << endl);
00301 if (dim_kinds[dim_prec] != GEN_VIRTUAL) {
00302 --tem_source_index;
00303 TEMP_INTEGER(source_dim);
00304 source_dim = source[tem_source_index][dim];
00305 TRACE(cerr << " rows:" << endl);
00306
00307
00308
00309
00310
00311
00312
00313
00314 for (dimension_type row = dest_index; row-- > 0; ) {
00315 assert(row < dest_num_rows);
00316 TRACE(cerr << " " << row << endl);
00317 Congruence& cg = dest[row];
00318 cg[dim_prec] -= source_dim * cg[dim];
00319 }
00320 }
00321 }
00322
00323 TRACE(cerr << "dest after processing preceding rows:" << endl);
00324 TRACE(dest.ascii_dump(cerr));
00325 }
00326
00327 Coefficient_traits::const_reference modulus = dest[dest_num_rows - 1][0];
00328 for (dimension_type row = 0; row < dest_num_rows; ++row) {
00329 Congruence& cg = dest[row];
00330 if (cg[dims] > 0)
00331
00332 cg[dims] = modulus;
00333 }
00334 TRACE(cerr << "dest after setting moduli:" << endl);
00335 TRACE(dest.ascii_dump(cerr));
00336
00337 assert(lower_triangular(dest, dim_kinds));
00338
00339 #ifdef STRONG_REDUCTION
00340 for (dimension_type dim = dims, i = 0; dim-- > 0; )
00341 if (dim_kinds[dim] != CON_VIRTUAL)
00342
00343 reduce_reduced<Congruence_System, Congruence>
00344 (dest, dim, i++, 0, dim, dim_kinds, false);
00345 TRACE(cerr << "dest after strong reduction:" << endl);
00346 TRACE(dest.ascii_dump(cerr));
00347 #endif
00348
00349 trace_dim_kinds("gs to cgs end ", dim_kinds);
00350
00351 TRACE(cerr << "------------------- gs to cgs conversion done." << endl);
00352 }
00353
00354 void
00355 Grid::conversion(Congruence_System& source, Grid_Generator_System& dest,
00356 Dimension_Kinds& dim_kinds) {
00357 TRACE(cerr << "============= convert cgs to gs" << endl);
00358 TRACE(cerr << "source:" << endl);
00359 TRACE(source.ascii_dump(cerr));
00360 TRACE(cerr << "dest:" << endl);
00361 TRACE(dest.ascii_dump(cerr));
00362 trace_dim_kinds("cgs to gs ", dim_kinds);
00363
00364
00365
00366
00367 assert(lower_triangular(source, dim_kinds));
00368
00369
00370
00371 dimension_type source_num_rows = 0, dest_num_rows = 0;
00372 TEMP_INTEGER(diagonal_lcm);
00373 diagonal_lcm = 1;
00374 dimension_type dims = source.num_columns() - 1;
00375 for (dimension_type dim = dims; dim-- > 0; )
00376 if (dim_kinds[dim] == CON_VIRTUAL)
00377
00378 ++dest_num_rows;
00379 else {
00380 if (dim_kinds[dim] == PROPER_CONGRUENCE) {
00381
00382
00383
00384 lcm_assign(diagonal_lcm, diagonal_lcm, source[source_num_rows][dim]);
00385
00386 ++dest_num_rows;
00387 }
00388
00389 ++source_num_rows;
00390 }
00391 TRACE(cerr << "diagonal_lcm: " << diagonal_lcm << endl);
00392 TRACE(cerr << "source_num_rows: " << source_num_rows << endl);
00393 TRACE(cerr << "dest_num_rows: " << dest_num_rows << endl);
00394
00395
00396 if (diagonal_lcm == 0)
00397 throw std::runtime_error("PPL internal error: Grid::conversion:"
00398 " source matrix is singular.");
00399
00400 dest.set_index_first_pending_row(dest_num_rows);
00401 dest.resize_no_copy(dest_num_rows, dims + 1 );
00402
00403
00404
00405
00406
00407
00408
00409 dimension_type source_index = 0;
00410
00411 dimension_type dest_index = dest_num_rows - 1;
00412 for (dimension_type dim = dims; dim-- > 0; ) {
00413 TRACE(cerr << "init dim " << dim << endl);
00414 if (dim_kinds[dim] == EQUALITY) {
00415 TRACE(cerr << " equality" << endl);
00416 ++source_index;
00417 }
00418 else {
00419 Grid_Generator& g = dest[dest_index];
00420 for (dimension_type j = 0; j < dim; ++j)
00421 g[j] = 0;
00422 for (dimension_type j = dim + 1; j < dims; ++j)
00423 g[j] = 0;
00424
00425 if (dim_kinds[dim] == CON_VIRTUAL) {
00426 TRACE(cerr << " con_virtual" << endl);
00427 g.set_is_line();
00428 g[dim] = 1;
00429 }
00430 else {
00431 assert(dim_kinds[dim] == PROPER_CONGRUENCE);
00432 TRACE(cerr << " proper_congruence" << endl);
00433 g.set_is_parameter_or_point();
00434 g[dim] = diagonal_lcm / source[source_index][dim];
00435 ++source_index;
00436 }
00437 --dest_index;
00438 }
00439 }
00440
00441 assert(upper_triangular(dest, dim_kinds));
00442
00443 TRACE(cerr << "dest after init:" << endl);
00444 TRACE(dest.ascii_dump(cerr));
00445
00446
00447
00448
00449
00450
00451
00452
00453 source_index = source_num_rows;
00454 dest_index = 0;
00455
00456 for (dimension_type dim = 0; dim < dims; ++dim) {
00457 TRACE(cerr << "dim: " << dim << endl);
00458
00459 if (dim_kinds[dim] != CON_VIRTUAL) {
00460 --source_index;
00461 TEMP_INTEGER(source_dim);
00462 source_dim = source[source_index][dim];
00463
00464
00465
00466 for (dimension_type row = dest_index; row-- > 0; ) {
00467 TRACE(cerr << " row " << row << endl);
00468 TRACE(dest.ascii_dump(cerr));
00469
00470 Grid_Generator& g = dest[row];
00471
00472
00473
00474
00475 TEMP_INTEGER(red_source_dim);
00476 gcd_assign(red_source_dim, g[dim], source_dim);
00477 red_source_dim = source_dim / red_source_dim;
00478 multiply_grid(red_source_dim, g, dest, dest_num_rows,
00479 dims + 1 );
00480
00481 g[dim] /= source_dim;
00482 }
00483 TRACE(cerr << "dest after dividing grid:" << endl);
00484 TRACE(dest.ascii_dump(cerr));
00485 }
00486
00487
00488
00489
00490
00491
00492
00493 dimension_type tem_source_index = source_index;
00494 if (dim_kinds[dim] != EQUALITY)
00495 ++dest_index;
00496 for (dimension_type dim_fol = dim + 1; dim_fol < dims; ++dim_fol) {
00497 TRACE(cerr << " dim_fol: " << dim_fol);
00498 TRACE(cerr << " dest_index: " << dest_index);
00499 TRACE(cerr << " tem_source_index: " << tem_source_index << endl);
00500 if (dim_kinds[dim_fol] != CON_VIRTUAL) {
00501 --tem_source_index;
00502 TEMP_INTEGER(source_dim);
00503 source_dim = source[tem_source_index][dim];
00504 TRACE(cerr << " rows:" << endl);
00505
00506
00507
00508
00509
00510
00511
00512
00513 for (dimension_type row = dest_index; row-- > 0; ) {
00514 assert(row < dest_num_rows);
00515 TRACE(cerr << " " << row << endl);
00516 Grid_Generator& g = dest[row];
00517 g[dim_fol] -= source_dim * g[dim];
00518 }
00519 }
00520 }
00521 TRACE(cerr << "dest after processing preceding rows:" << endl);
00522 TRACE(dest.ascii_dump(cerr));
00523 }
00524
00525 assert(upper_triangular(dest, dim_kinds));
00526
00527 #ifdef STRONG_REDUCTION
00528 for (dimension_type dim = 0, i = 0; dim < dims; ++dim)
00529 if (dim_kinds[dim] != GEN_VIRTUAL)
00530
00531 reduce_reduced<Grid_Generator_System, Grid_Generator>
00532 (dest, dim, i++, dim, dims - 1, dim_kinds);
00533 TRACE(cerr << "dest after strong reduction:" << endl);
00534 TRACE(dest.ascii_dump(cerr));
00535 #endif
00536
00537
00538
00539 Coefficient_traits::const_reference system_divisor = dest[0][0];
00540 for (dimension_type row = 1, dim = 1; dim < dims; ++dim)
00541 switch (dim_kinds[dim]) {
00542 case PARAMETER:
00543 dest[row].divisor() = system_divisor;
00544 case LINE:
00545 ++row;
00546 case GEN_VIRTUAL:
00547 break;
00548 }
00549 TRACE(cerr << "dest after updating param divisors:" << endl);
00550 TRACE(dest.ascii_dump(cerr));
00551
00552 trace_dim_kinds("cgs to gs end ", dim_kinds);
00553
00554 TRACE(cerr << "------------------- cgs to gs conversion done." << endl);
00555 }
00556
00557 #undef TRACE
00558
00559 }