1 /* Grid class implementation: conversion().
2 Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
3 Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
4
5 This file is part of the Parma Polyhedra Library (PPL).
6
7 The PPL is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 The PPL is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software Foundation,
19 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
20
21 For the most up-to-date information see the Parma Polyhedra Library
22 site: http://bugseng.com/products/ppl/ . */
23
24 #include "ppl-config.h"
25 #include "Grid_defs.hh"
26 #include <cstddef>
27
28 namespace Parma_Polyhedra_Library {
29
30 // X 0 0 0 upside down, so x x x X
31 // x X 0 0 x x X 0
32 // x x X 0 x X 0 0
33 // x x x X X 0 0 0
34 //
35 // Where X is greater than zero and x is an integer.
36 bool
lower_triangular(const Congruence_System & sys,const Dimension_Kinds & dim_kinds)37 Grid::lower_triangular(const Congruence_System& sys,
38 const Dimension_Kinds& dim_kinds) {
39 const dimension_type num_columns = sys.space_dimension() + 1;
40
41 // Check for easy square failure case.
42 if (sys.num_rows() > num_columns) {
43 return false;
44 }
45
46 // Check triangularity.
47
48 dimension_type row = 0;
49 for (dimension_type dim = num_columns; dim-- > 0; ) {
50 if (dim_kinds[dim] == CON_VIRTUAL) {
51 continue;
52 }
53 const Congruence& cg = sys[row];
54 ++row;
55 // Check diagonal.
56 if (cg.expr.get(dim) <= 0) {
57 return false;
58 }
59 if (!cg.expr.all_zeroes(dim + 1, num_columns)) {
60 return false;
61 }
62 }
63
64 // Check squareness.
65 return row == sys.num_rows();
66 }
67
68 // X x x x
69 // 0 X x x
70 // 0 0 X x
71 // 0 0 0 X
72 //
73 // Where X is greater than zero and x is an integer.
74 bool
upper_triangular(const Grid_Generator_System & sys,const Dimension_Kinds & dim_kinds)75 Grid::upper_triangular(const Grid_Generator_System& sys,
76 const Dimension_Kinds& dim_kinds) {
77 dimension_type num_columns = sys.space_dimension() + 1;
78 dimension_type row = sys.num_rows();
79
80 // Check for easy square fail case.
81 if (row > num_columns) {
82 return false;
83 }
84
85 // Check triangularity.
86 while (num_columns > 0) {
87 --num_columns;
88 if (dim_kinds[num_columns] == GEN_VIRTUAL) {
89 continue;
90 }
91 const Grid_Generator& gen = sys[--row];
92 // Check diagonal.
93 if (gen.expr.get(num_columns) <= 0) {
94 return false;
95 }
96 // Check elements preceding diagonal.
97 if (!gen.expr.all_zeroes(0, num_columns)) {
98 return false;
99 }
100 }
101
102 // Check for squareness.
103 return num_columns == row;
104 }
105
106 void
multiply_grid(const Coefficient & multiplier,Grid_Generator & gen,Swapping_Vector<Grid_Generator> & dest_rows,const dimension_type num_rows)107 Grid::multiply_grid(const Coefficient& multiplier, Grid_Generator& gen,
108 Swapping_Vector<Grid_Generator>& dest_rows,
109 const dimension_type num_rows) {
110 if (multiplier == 1) {
111 return;
112 }
113
114 if (gen.is_line()) {
115 // Multiply every element of the line.
116 gen.expr *= multiplier;
117 }
118 else {
119 PPL_ASSERT(gen.is_parameter_or_point());
120 // Multiply every element of every parameter.
121
122 for (dimension_type index = num_rows; index-- > 0; ) {
123 Grid_Generator& generator = dest_rows[index];
124 if (generator.is_parameter_or_point()) {
125 generator.expr *= multiplier;
126 }
127 }
128 }
129 }
130
131 void
multiply_grid(const Coefficient & multiplier,Congruence & cg,Swapping_Vector<Congruence> & dest,const dimension_type num_rows)132 Grid::multiply_grid(const Coefficient& multiplier, Congruence& cg,
133 Swapping_Vector<Congruence>& dest,
134 const dimension_type num_rows) {
135 if (multiplier == 1) {
136 return;
137 }
138
139 if (cg.is_proper_congruence()) {
140 // Multiply every element of every congruence.
141 for (dimension_type index = num_rows; index-- > 0; ) {
142 Congruence& congruence = dest[index];
143 if (congruence.is_proper_congruence()) {
144 congruence.scale(multiplier);
145 }
146 }
147 }
148 else {
149 PPL_ASSERT(cg.is_equality());
150 // Multiply every element of the equality.
151 cg.scale(multiplier);
152 }
153 }
154
155 // TODO: Rename the next two methods to convert and this file to
156 // Grid_convert.cc (and equivalently for Polyhedron) to use
157 // verbs consistently as function and method names.
158
159 void
conversion(Grid_Generator_System & source,Congruence_System & dest,Dimension_Kinds & dim_kinds)160 Grid::conversion(Grid_Generator_System& source, Congruence_System& dest,
161 Dimension_Kinds& dim_kinds) {
162 // Quite similar to the congruence to generator version below.
163 // Changes here may be needed there too.
164
165 PPL_ASSERT(upper_triangular(source, dim_kinds));
166
167 // Initialize matrix row number counters and compute the LCM of the
168 // diagonal entries of the parameters in `source'.
169 //
170 // The top-down order of the generator system rows corresponds to
171 // the left-right order of the dimensions, while the congruence
172 // system rows have a bottom-up ordering.
173 dimension_type dest_num_rows = 0;
174 PPL_DIRTY_TEMP_COEFFICIENT(diagonal_lcm);
175 diagonal_lcm = 1;
176 const dimension_type dims = source.space_dimension() + 1;
177 dimension_type source_index = source.num_rows();
178 for (dimension_type dim = dims; dim-- > 0; ) {
179 if (dim_kinds[dim] == GEN_VIRTUAL) {
180 // Virtual generators map to equalities.
181 ++dest_num_rows;
182 }
183 else {
184 --source_index;
185 if (dim_kinds[dim] == PARAMETER) {
186 // Dimension `dim' has a parameter row at `source_index' in
187 // `source', so include in `diagonal_lcm' the `dim'th element
188 // of that row.
189 lcm_assign(diagonal_lcm, diagonal_lcm, source[source_index].expr.get(dim));
190 // Parameters map to proper congruences.
191 ++dest_num_rows;
192 }
193 // Lines map to virtual congruences.
194 }
195 }
196 PPL_ASSERT(source_index == 0);
197
198 // `source' must be regular.
199 PPL_ASSERT(diagonal_lcm != 0);
200
201 dest.clear();
202 dest.set_space_dimension(dims - 1);
203
204 // In `dest' initialize row types and elements, including setting
205 // the diagonal elements to the inverse ratio of the `source'
206 // diagonal elements.
207 source_index = source.num_rows();
208 for (dimension_type dim = dims; dim-- > 0; ) {
209 if (dim_kinds[dim] == LINE) {
210 --source_index;
211 }
212 else {
213 Linear_Expression le;
214 le.set_space_dimension(dest.space_dimension());
215
216 if (dim_kinds[dim] == GEN_VIRTUAL) {
217 le.set(dim, Coefficient_one());
218 Congruence cg(le, Coefficient_zero(), Recycle_Input());
219 dest.insert_verbatim(cg, Recycle_Input());
220 }
221 else {
222 PPL_ASSERT(dim_kinds[dim] == PARAMETER);
223 --source_index;
224 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
225 exact_div_assign(tmp, diagonal_lcm,
226 source[source_index].expr.get(dim));
227 le.set(dim, tmp);
228 Congruence cg(le, Coefficient_one(), Recycle_Input());
229 dest.insert_verbatim(cg, Recycle_Input());
230 }
231 }
232 }
233
234 PPL_ASSERT(source_index == 0);
235 PPL_ASSERT(lower_triangular(dest, dim_kinds));
236
237 // Convert.
238 //
239 // `source_index' and `dest_index' hold the positions of pivot rows
240 // in `source' and `dest'. The order of the rows in `dest' is the
241 // reverse of the order in `source', so the rows are iterated from
242 // last to first (index 0) in `source' and from first to last in
243 // `dest'.
244 source_index = source.num_rows();
245 dimension_type dest_index = 0;
246 PPL_DIRTY_TEMP_COEFFICIENT(multiplier);
247
248 for (dimension_type dim = dims; dim-- > 0; ) {
249 if (dim_kinds[dim] != GEN_VIRTUAL) {
250 --source_index;
251 const Coefficient& source_dim = source[source_index].expr.get(dim);
252
253 // In the rows in `dest' above `dest_index' divide each element
254 // at column `dim' by `source_dim'.
255 for (dimension_type row = dest_index; row-- > 0; ) {
256 Congruence& cg = dest.rows[row];
257
258 // Multiply the representation of `dest' such that entry `dim'
259 // of `g' is a multiple of `source_dim'. This ensures that
260 // the result of the division that follows is a whole number.
261 gcd_assign(multiplier, cg.expression().get(dim), source_dim);
262 exact_div_assign(multiplier, source_dim, multiplier);
263 multiply_grid(multiplier, cg, dest.rows, dest_num_rows);
264
265 cg.expr.exact_div_assign(source_dim, dim, dim + 1);
266 }
267 }
268
269 // Invert and transpose the source row at `source_index' into the
270 // destination row at `dest_index'.
271 //
272 // Consider each dimension `dim_prec' that precedes `dim', as the
273 // rows in `dest' that follow `dim_index' have zeroes at index
274 // `dim'.
275 dimension_type tmp_source_index = source_index;
276 if (dim_kinds[dim] != LINE) {
277 ++dest_index;
278 }
279 for (dimension_type dim_prec = dim; dim_prec-- > 0; ) {
280 if (dim_kinds[dim_prec] != GEN_VIRTUAL) {
281 --tmp_source_index;
282 const Coefficient& source_dim = source[tmp_source_index].expr.get(dim);
283
284 // In order to compute the transpose of the inverse of
285 // `source', subtract source[tmp_source_index][dim] times the
286 // column vector in `dest' at `dim' from the column vector in
287 // `dest' at `dim_prec'.
288 //
289 // I.e., for each row `dest_index' in `dest' that is above the
290 // row `dest_index', subtract dest[tmp_source_index][dim]
291 // times the entry `dim' from the entry at `dim_prec'.
292 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
293 for (dimension_type row = dest_index; row-- > 0; ) {
294 PPL_ASSERT(row < dest_num_rows);
295 Congruence& cg = dest.rows[row];
296 tmp = cg.expr.get(dim_prec);
297 sub_mul_assign(tmp, source_dim, cg.expression().get(dim));
298 cg.expr.set(dim_prec, tmp);
299 }
300 }
301 }
302 }
303 // Set the modulus in every congruence.
304 Coefficient_traits::const_reference modulus
305 = dest.rows[dest_num_rows - 1].inhomogeneous_term();
306 for (dimension_type row = dest_num_rows; row-- > 0; ) {
307 Congruence& cg = dest.rows[row];
308 if (cg.is_proper_congruence()) {
309 cg.set_modulus(modulus);
310 }
311 }
312
313 PPL_ASSERT(lower_triangular(dest, dim_kinds));
314
315 // Since we are reducing the system to "strong minimal form",
316 // reduce the coefficients in the congruence system
317 // using "diagonal" values.
318 for (dimension_type dim = dims, i = 0; dim-- > 0; ) {
319 if (dim_kinds[dim] != CON_VIRTUAL) {
320 // Factor the "diagonal" congruence out of the preceding rows.
321 reduce_reduced<Congruence_System>
322 (dest.rows, dim, i++, 0, dim, dim_kinds, false);
323 }
324 }
325
326 #ifndef NDEBUG
327 // Make sure that all the rows are now OK.
328 for (dimension_type i = dest.num_rows(); i-- > 0; ) {
329 PPL_ASSERT(dest[i].OK());
330 }
331 #endif
332
333 PPL_ASSERT(dest.OK());
334 }
335
336 void
conversion(Congruence_System & source,Grid_Generator_System & dest,Dimension_Kinds & dim_kinds)337 Grid::conversion(Congruence_System& source, Grid_Generator_System& dest,
338 Dimension_Kinds& dim_kinds) {
339 // Quite similar to the generator to congruence version above.
340 // Changes here may be needed there too.
341
342 PPL_ASSERT(lower_triangular(source, dim_kinds));
343
344 // Initialize matrix row number counters and compute the LCM of the
345 // diagonal entries of the proper congruences in `source'.
346 dimension_type source_num_rows = 0;
347 dimension_type dest_num_rows = 0;
348 PPL_DIRTY_TEMP_COEFFICIENT(diagonal_lcm);
349 diagonal_lcm = 1;
350 const dimension_type dims = source.space_dimension() + 1;
351 for (dimension_type dim = dims; dim-- > 0; ) {
352 if (dim_kinds[dim] == CON_VIRTUAL) {
353 // Virtual congruences map to lines.
354 ++dest_num_rows;
355 }
356 else {
357 if (dim_kinds[dim] == PROPER_CONGRUENCE) {
358 // Dimension `dim' has a proper congruence row at
359 // `source_num_rows' in `source', so include in `diagonal_lcm'
360 // the `dim'th element of that row.
361 lcm_assign(diagonal_lcm, diagonal_lcm, source[source_num_rows].expr.get(dim));
362 // Proper congruences map to parameters.
363 ++dest_num_rows;
364 }
365 // Equalities map to virtual generators.
366 ++source_num_rows;
367 }
368
369 }
370 // `source' must be regular.
371 PPL_ASSERT(diagonal_lcm != 0);
372
373 dest.clear();
374 PPL_ASSERT(dims > 0);
375 dest.set_space_dimension(dims - 1);
376
377 // In `dest' initialize row types and elements, including setting
378 // the diagonal elements to the inverse ratio of the `source'
379 // diagonal elements.
380 //
381 // The top-down order of the congruence system rows corresponds to
382 // the right-left order of the dimensions.
383 dimension_type source_index = source_num_rows;
384 // The generator system has a bottom-up ordering.
385 for (dimension_type dim = 0; dim < dims; ++dim) {
386 if (dim_kinds[dim] == EQUALITY) {
387 --source_index;
388 }
389 else {
390 Grid_Generator g(dest.representation());
391 g.set_topology(dest.topology());
392 g.expr.set_space_dimension(dims);
393
394 if (dim_kinds[dim] == CON_VIRTUAL) {
395 g.set_is_line();
396 g.expr.set(0, Coefficient_zero());
397 g.expr.set(dim, Coefficient_one());
398 }
399 else {
400 PPL_ASSERT(dim_kinds[dim] == PROPER_CONGRUENCE);
401 g.set_is_parameter_or_point();
402 --source_index;
403 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
404 exact_div_assign(tmp, diagonal_lcm,
405 source[source_index].expr.get(dim));
406 g.expr.set(0, Coefficient_zero());
407 g.expr.set(dim, tmp);
408 }
409 // Don't assert g.OK() here, because it may fail.
410 // All the rows in `dest' are checked at the end of this function.
411 dest.insert_verbatim(g);
412 }
413 }
414
415 PPL_ASSERT(dest.num_rows() == dest_num_rows);
416 PPL_ASSERT(dest.first_pending_row() == dest_num_rows);
417 PPL_ASSERT(upper_triangular(dest, dim_kinds));
418
419 // Convert.
420 //
421 // `source_index' and `dest_index' hold the positions of pivot rows
422 // in `source' and `dest'. The order of the rows in `dest' is the
423 // reverse of the order in `source', so the rows are iterated from
424 // last to first (index 0) in `source' and from first to last in
425 // `dest'.
426 source_index = source_num_rows;
427 dimension_type dest_index = 0;
428 PPL_DIRTY_TEMP_COEFFICIENT(reduced_source_dim);
429
430 for (dimension_type dim = 0; dim < dims; ++dim) {
431 if (dim_kinds[dim] != CON_VIRTUAL) {
432 --source_index;
433 Coefficient_traits::const_reference source_dim
434 = source[source_index].expr.get(dim);
435
436 // In the rows in `dest' above `dest_index' divide each element
437 // at column `dim' by `source_dim'.
438
439 for (dimension_type i = dest_index; i-- > 0; ) {
440 Grid_Generator& g = dest.sys.rows[i];
441
442 // Multiply the representation of `dest' such that entry `dim'
443 // of `g' is a multiple of `source_dim'. This ensures that
444 // the result of the division that follows is a whole number.
445 gcd_assign(reduced_source_dim, g.expr.get(dim), source_dim);
446 exact_div_assign(reduced_source_dim, source_dim, reduced_source_dim);
447 multiply_grid(reduced_source_dim, g, dest.sys.rows, dest_num_rows);
448
449 g.expr.exact_div_assign(source_dim, dim, dim + 1);
450 // Don't assert g.OK() here, because it may fail.
451 // All the rows in `dest' are checked at the end of this function.
452 }
453 }
454
455 // Invert and transpose the source row at `source_index' into the
456 // destination row at `dest_index'.
457 //
458 // Consider each dimension `dim_fol' that follows `dim', as the
459 // rows in `dest' that follow row `dest_index' are zero at index
460 // `dim'.
461 dimension_type tmp_source_index = source_index;
462 if (dim_kinds[dim] != EQUALITY) {
463 ++dest_index;
464 }
465 for (dimension_type dim_fol = dim + 1; dim_fol < dims; ++dim_fol) {
466 if (dim_kinds[dim_fol] != CON_VIRTUAL) {
467 --tmp_source_index;
468 Coefficient_traits::const_reference source_dim
469 = source[tmp_source_index].expr.get(dim);
470 // In order to compute the transpose of the inverse of
471 // `source', subtract source[tmp_source_index][dim] times the
472 // column vector in `dest' at `dim' from the column vector in
473 // `dest' at `dim_fol'.
474 //
475 // I.e., for each row `dest_index' in `dest' that is above the
476 // row `dest_index', subtract dest[tmp_source_index][dim]
477 // times the entry `dim' from the entry at `dim_fol'.
478
479 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
480 for (dimension_type i = dest_index; i-- > 0; ) {
481 PPL_ASSERT(i < dest_num_rows);
482 Grid_Generator& row = dest.sys.rows[i];
483 tmp = row.expr.get(dim_fol);
484 sub_mul_assign(tmp, source_dim,
485 row.expr.get(dim));
486 row.expr.set(dim_fol, tmp);
487 // Don't assert row.OK() here, because it may fail.
488 // All the rows in `dest' are checked at the end of this function.
489 }
490 }
491 }
492 }
493
494 PPL_ASSERT(upper_triangular(dest, dim_kinds));
495
496 // Since we are reducing the system to "strong minimal form",
497 // reduce the coordinates in the grid_generator system
498 // using "diagonal" values.
499 for (dimension_type dim = 0, i = 0; dim < dims; ++dim) {
500 if (dim_kinds[dim] != GEN_VIRTUAL) {
501 // Factor the "diagonal" generator out of the preceding rows.
502 reduce_reduced<Grid_Generator_System>
503 (dest.sys.rows, dim, i++, dim, dims - 1, dim_kinds);
504 }
505 }
506
507 // Ensure that the parameter divisors are the same as the divisor of
508 // the point.
509 const Coefficient& system_divisor
510 = dest.sys.rows[0].expr.inhomogeneous_term();
511
512 for (dimension_type i = dest.sys.rows.size() - 1, dim = dims; dim-- > 1; ) {
513 switch (dim_kinds[dim]) {
514 case PARAMETER:
515 dest.sys.rows[i].set_divisor(system_divisor);
516 // Intentionally fall through.
517 case LINE:
518 --i;
519 break;
520 case GEN_VIRTUAL:
521 break;
522 }
523 }
524
525 #ifndef NDEBUG
526 // The previous code can modify the rows' fields, exploiting the friendness.
527 // Check that all rows are OK now.
528 for (dimension_type i = dest.sys.rows.size(); i-- > 0; ) {
529 PPL_ASSERT(dest.sys.rows[i].OK());
530 }
531 #endif
532
533 PPL_ASSERT(dest.sys.OK());
534 }
535
536 } // namespace Parma_Polyhedra_Library
537