1 /* Polyhedron class implementation
2 (non-inline private or protected functions).
3 Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
4 Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
5
6 This file is part of the Parma Polyhedra Library (PPL).
7
8 The PPL is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The PPL is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software Foundation,
20 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
21
22 For the most up-to-date information see the Parma Polyhedra Library
23 site: http://bugseng.com/products/ppl/ . */
24
25 #include "ppl-config.h"
26 #include "Polyhedron_defs.hh"
27 #include "Scalar_Products_defs.hh"
28 #include "Scalar_Products_inlines.hh"
29 #include "Linear_Form_defs.hh"
30 #include "C_Integer.hh"
31 #include "assertions.hh"
32 #include <string>
33 #include <iostream>
34 #include <sstream>
35 #include <stdexcept>
36
37 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
38 /*! \ingroup PPL_defines
39 \brief
40 Controls the laziness level of the implementation.
41
42 Temporarily used in a few of the function implementations to
43 switch to an even more lazy algorithm. To be removed as soon as
44 we collect enough information to decide which is the better
45 implementation alternative.
46 */
47 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
48 #define BE_LAZY 1
49
50 namespace PPL = Parma_Polyhedra_Library;
51
Polyhedron(const Topology topol,const dimension_type num_dimensions,const Degenerate_Element kind)52 PPL::Polyhedron::Polyhedron(const Topology topol,
53 const dimension_type num_dimensions,
54 const Degenerate_Element kind)
55 : con_sys(topol, default_con_sys_repr),
56 gen_sys(topol, default_gen_sys_repr),
57 sat_c(),
58 sat_g() {
59 // Protecting against space dimension overflow is up to the caller.
60 PPL_ASSERT(num_dimensions <= max_space_dimension());
61
62 if (kind == EMPTY) {
63 status.set_empty();
64 }
65 else if (num_dimensions > 0) {
66 con_sys.add_low_level_constraints();
67 con_sys.adjust_topology_and_space_dimension(topol, num_dimensions);
68 set_constraints_minimized();
69 }
70 space_dim = num_dimensions;
71 PPL_ASSERT_HEAVY(OK());
72 }
73
Polyhedron(const Polyhedron & y,Complexity_Class)74 PPL::Polyhedron::Polyhedron(const Polyhedron& y, Complexity_Class)
75 : con_sys(y.topology(), default_con_sys_repr),
76 gen_sys(y.topology(), default_gen_sys_repr),
77 status(y.status),
78 space_dim(y.space_dim) {
79 // Being a protected method, we simply assert that topologies do match.
80 PPL_ASSERT(topology() == y.topology());
81 if (y.constraints_are_up_to_date()) {
82 con_sys.assign_with_pending(y.con_sys);
83 }
84 if (y.generators_are_up_to_date()) {
85 gen_sys.assign_with_pending(y.gen_sys);
86 }
87 if (y.sat_c_is_up_to_date()) {
88 sat_c = y.sat_c;
89 }
90 if (y.sat_g_is_up_to_date()) {
91 sat_g = y.sat_g;
92 }
93 }
94
Polyhedron(const Topology topol,const Constraint_System & cs)95 PPL::Polyhedron::Polyhedron(const Topology topol, const Constraint_System& cs)
96 : con_sys(topol, default_con_sys_repr),
97 gen_sys(topol, default_gen_sys_repr),
98 sat_c(),
99 sat_g() {
100 // Protecting against space dimension overflow is up to the caller.
101 PPL_ASSERT(cs.space_dimension() <= max_space_dimension());
102
103 // TODO: this implementation is just an executable specification.
104 Constraint_System cs_copy = cs;
105
106 // Try to adapt `cs_copy' to the required topology.
107 const dimension_type cs_copy_space_dim = cs_copy.space_dimension();
108 if (!cs_copy.adjust_topology_and_space_dimension(topol, cs_copy_space_dim)) {
109 throw_topology_incompatible((topol == NECESSARILY_CLOSED)
110 ? "C_Polyhedron(cs)"
111 : "NNC_Polyhedron(cs)", "cs", cs_copy);
112 }
113 // Set the space dimension.
114 space_dim = cs_copy_space_dim;
115
116 if (space_dim > 0) {
117 // Stealing the rows from `cs_copy'.
118 using std::swap;
119 swap(con_sys, cs_copy);
120 if (con_sys.num_pending_rows() > 0) {
121 // Even though `cs_copy' has pending constraints, since the
122 // generators of the polyhedron are not up-to-date, the
123 // polyhedron cannot have pending constraints. By integrating
124 // the pending part of `con_sys' we may loose sortedness.
125 con_sys.set_sorted(false);
126 con_sys.unset_pending_rows();
127 }
128 con_sys.add_low_level_constraints();
129 set_constraints_up_to_date();
130 }
131 else {
132 // Here `space_dim == 0'.
133 // See if an inconsistent constraint has been passed.
134 for (dimension_type i = cs_copy.num_rows(); i-- > 0; ) {
135 if (cs_copy[i].is_inconsistent()) {
136 // Inconsistent constraint found: the polyhedron is empty.
137 set_empty();
138 break;
139 }
140 }
141 }
142 PPL_ASSERT_HEAVY(OK());
143 }
144
Polyhedron(const Topology topol,Constraint_System & cs,Recycle_Input)145 PPL::Polyhedron::Polyhedron(const Topology topol,
146 Constraint_System& cs,
147 Recycle_Input)
148 : con_sys(topol, default_con_sys_repr),
149 gen_sys(topol, default_gen_sys_repr),
150 sat_c(),
151 sat_g() {
152 // Protecting against space dimension overflow is up to the caller.
153 PPL_ASSERT(cs.space_dimension() <= max_space_dimension());
154
155 // Try to adapt `cs' to the required topology.
156 const dimension_type cs_space_dim = cs.space_dimension();
157 if (!cs.adjust_topology_and_space_dimension(topol, cs_space_dim)) {
158 throw_topology_incompatible((topol == NECESSARILY_CLOSED)
159 ? "C_Polyhedron(cs, recycle)"
160 : "NNC_Polyhedron(cs, recycle)", "cs", cs);
161 }
162
163 // Set the space dimension.
164 space_dim = cs_space_dim;
165
166 if (space_dim > 0) {
167 // Stealing the rows from `cs'.
168 swap(con_sys, cs);
169 if (con_sys.num_pending_rows() > 0) {
170 // Even though `cs' has pending constraints, since the generators
171 // of the polyhedron are not up-to-date, the polyhedron cannot
172 // have pending constraints. By integrating the pending part
173 // of `con_sys' we may loose sortedness.
174 con_sys.unset_pending_rows();
175 con_sys.set_sorted(false);
176 }
177 con_sys.add_low_level_constraints();
178 set_constraints_up_to_date();
179 }
180 else {
181 // Here `space_dim == 0'.
182
183 // See if an inconsistent constraint has been passed.
184 for (dimension_type i = cs.num_rows(); i-- > 0; ) {
185 if (cs[i].is_inconsistent()) {
186 // Inconsistent constraint found: the polyhedron is empty.
187 set_empty();
188 break;
189 }
190 }
191 }
192 PPL_ASSERT_HEAVY(OK());
193 }
194
Polyhedron(const Topology topol,const Generator_System & gs)195 PPL::Polyhedron::Polyhedron(const Topology topol, const Generator_System& gs)
196 : con_sys(topol, default_con_sys_repr),
197 gen_sys(topol, default_gen_sys_repr),
198 sat_c(),
199 sat_g() {
200 // Protecting against space dimension overflow is up to the caller.
201 PPL_ASSERT(gs.space_dimension() <= max_space_dimension());
202
203 // An empty set of generators defines the empty polyhedron.
204 if (gs.has_no_rows()) {
205 space_dim = gs.space_dimension();
206 status.set_empty();
207 PPL_ASSERT_HEAVY(OK());
208 return;
209 }
210
211 // Non-empty valid generator systems have a supporting point, at least.
212 if (!gs.has_points()) {
213 throw_invalid_generators((topol == NECESSARILY_CLOSED)
214 ? "C_Polyhedron(gs)"
215 : "NNC_Polyhedron(gs)", "gs");
216 }
217 // TODO: this implementation is just an executable specification.
218 Generator_System gs_copy = gs;
219
220 const dimension_type gs_copy_space_dim = gs_copy.space_dimension();
221 // Try to adapt `gs_copy' to the required topology.
222 if (!gs_copy.adjust_topology_and_space_dimension(topol, gs_copy_space_dim)) {
223 throw_topology_incompatible((topol == NECESSARILY_CLOSED)
224 ? "C_Polyhedron(gs)"
225 : "NNC_Polyhedron(gs)", "gs", gs_copy);
226 }
227
228 if (gs_copy_space_dim > 0) {
229 // Stealing the rows from `gs_copy'.
230 using std::swap;
231 swap(gen_sys, gs_copy);
232 // In a generator system describing a NNC polyhedron,
233 // for each point we must also have the corresponding closure point.
234 if (topol == NOT_NECESSARILY_CLOSED) {
235 gen_sys.add_corresponding_closure_points();
236 }
237 if (gen_sys.num_pending_rows() > 0) {
238 // Even though `gs_copy' has pending generators, since the
239 // constraints of the polyhedron are not up-to-date, the
240 // polyhedron cannot have pending generators. By integrating the
241 // pending part of `gen_sys' we may lose sortedness.
242 gen_sys.set_sorted(false);
243 gen_sys.unset_pending_rows();
244 }
245 // Generators are now up-to-date.
246 set_generators_up_to_date();
247
248 // Set the space dimension.
249 space_dim = gs_copy_space_dim;
250 PPL_ASSERT_HEAVY(OK());
251 return;
252 }
253
254 // Here `gs_copy.num_rows > 0' and `gs_copy_space_dim == 0':
255 // we already checked for both the topology-compatibility
256 // and the supporting point.
257 space_dim = 0;
258 PPL_ASSERT_HEAVY(OK());
259 }
260
Polyhedron(const Topology topol,Generator_System & gs,Recycle_Input)261 PPL::Polyhedron::Polyhedron(const Topology topol,
262 Generator_System& gs,
263 Recycle_Input)
264 : con_sys(topol, default_con_sys_repr),
265 gen_sys(topol, default_gen_sys_repr),
266 sat_c(),
267 sat_g() {
268 // Protecting against space dimension overflow is up to the caller.
269 PPL_ASSERT(gs.space_dimension() <= max_space_dimension());
270
271 // An empty set of generators defines the empty polyhedron.
272 if (gs.has_no_rows()) {
273 space_dim = gs.space_dimension();
274 status.set_empty();
275 PPL_ASSERT_HEAVY(OK());
276 return;
277 }
278
279 // Non-empty valid generator systems have a supporting point, at least.
280 if (!gs.has_points()) {
281 throw_invalid_generators((topol == NECESSARILY_CLOSED)
282 ? "C_Polyhedron(gs, recycle)"
283 : "NNC_Polyhedron(gs, recycle)", "gs");
284 }
285
286 const dimension_type gs_space_dim = gs.space_dimension();
287 // Try to adapt `gs' to the required topology.
288 if (!gs.adjust_topology_and_space_dimension(topol, gs_space_dim)) {
289 throw_topology_incompatible((topol == NECESSARILY_CLOSED)
290 ? "C_Polyhedron(gs, recycle)"
291 : "NNC_Polyhedron(gs, recycle)", "gs", gs);
292 }
293
294 if (gs_space_dim > 0) {
295 // Stealing the rows from `gs'.
296 swap(gen_sys, gs);
297 // In a generator system describing a NNC polyhedron,
298 // for each point we must also have the corresponding closure point.
299 if (topol == NOT_NECESSARILY_CLOSED) {
300 gen_sys.add_corresponding_closure_points();
301 }
302 if (gen_sys.num_pending_rows() > 0) {
303 // Even though `gs' has pending generators, since the constraints
304 // of the polyhedron are not up-to-date, the polyhedron cannot
305 // have pending generators. By integrating the pending part
306 // of `gen_sys' we may loose sortedness.
307 gen_sys.set_sorted(false);
308 gen_sys.unset_pending_rows();
309 }
310 // Generators are now up-to-date.
311 set_generators_up_to_date();
312
313 // Set the space dimension.
314 space_dim = gs_space_dim;
315 PPL_ASSERT_HEAVY(OK());
316 return;
317 }
318
319 // Here `gs.num_rows > 0' and `gs_space_dim == 0':
320 // we already checked for both the topology-compatibility
321 // and the supporting point.
322 space_dim = 0;
323 PPL_ASSERT_HEAVY(OK());
324 }
325
326 PPL::Polyhedron&
operator =(const Polyhedron & y)327 PPL::Polyhedron::operator=(const Polyhedron& y) {
328 // Being a protected method, we simply assert that topologies do match.
329 PPL_ASSERT(topology() == y.topology());
330 space_dim = y.space_dim;
331 if (y.marked_empty()) {
332 set_empty();
333 }
334 else if (space_dim == 0) {
335 set_zero_dim_univ();
336 }
337 else {
338 status = y.status;
339 if (y.constraints_are_up_to_date()) {
340 con_sys.assign_with_pending(y.con_sys);
341 }
342 if (y.generators_are_up_to_date()) {
343 gen_sys.assign_with_pending(y.gen_sys);
344 }
345 if (y.sat_c_is_up_to_date()) {
346 sat_c = y.sat_c;
347 }
348 if (y.sat_g_is_up_to_date()) {
349 sat_g = y.sat_g;
350 }
351 }
352 return *this;
353 }
354
355 PPL::Polyhedron::Three_Valued_Boolean
quick_equivalence_test(const Polyhedron & y) const356 PPL::Polyhedron::quick_equivalence_test(const Polyhedron& y) const {
357 // Private method: the caller must ensure the following.
358 PPL_ASSERT(topology() == y.topology());
359 PPL_ASSERT(space_dim == y.space_dim);
360 PPL_ASSERT(!marked_empty() && !y.marked_empty() && space_dim > 0);
361
362 const Polyhedron& x = *this;
363
364 if (x.is_necessarily_closed()) {
365 if (!x.has_something_pending() && !y.has_something_pending()) {
366 bool css_normalized = false;
367 if (x.constraints_are_minimized() && y.constraints_are_minimized()) {
368 // Equivalent minimized constraint systems have:
369 // - the same number of constraints; ...
370 if (x.con_sys.num_rows() != y.con_sys.num_rows()) {
371 return Polyhedron::TVB_FALSE;
372 }
373 // - the same number of equalities; ...
374 const dimension_type x_num_equalities = x.con_sys.num_equalities();
375 if (x_num_equalities != y.con_sys.num_equalities()) {
376 return Polyhedron::TVB_FALSE;
377 }
378 // - if there are no equalities, they have the same constraints.
379 // Delay this test: try cheaper tests on generators first.
380 css_normalized = (x_num_equalities == 0);
381 }
382
383 if (x.generators_are_minimized() && y.generators_are_minimized()) {
384 // Equivalent minimized generator systems have:
385 // - the same number of generators; ...
386 if (x.gen_sys.num_rows() != y.gen_sys.num_rows()) {
387 return Polyhedron::TVB_FALSE;
388 }
389 // - the same number of lines; ...
390 const dimension_type x_num_lines = x.gen_sys.num_lines();
391 if (x_num_lines != y.gen_sys.num_lines()) {
392 return Polyhedron::TVB_FALSE;
393 }
394 // - if there are no lines, they have the same generators.
395 if (x_num_lines == 0) {
396 // Sort the two systems and check for syntactic identity.
397 x.obtain_sorted_generators();
398 y.obtain_sorted_generators();
399 if (x.gen_sys == y.gen_sys) {
400 return Polyhedron::TVB_TRUE;
401 }
402 else {
403 return Polyhedron::TVB_FALSE;
404 }
405 }
406 }
407
408 if (css_normalized) {
409 // Sort the two systems and check for identity.
410 x.obtain_sorted_constraints();
411 y.obtain_sorted_constraints();
412 if (x.con_sys == y.con_sys) {
413 return Polyhedron::TVB_TRUE;
414 }
415 else {
416 return Polyhedron::TVB_FALSE;
417 }
418 }
419 }
420 }
421 return Polyhedron::TVB_DONT_KNOW;
422 }
423
424 bool
is_included_in(const Polyhedron & y) const425 PPL::Polyhedron::is_included_in(const Polyhedron& y) const {
426 // Private method: the caller must ensure the following.
427 PPL_ASSERT(topology() == y.topology());
428 PPL_ASSERT(space_dim == y.space_dim);
429 PPL_ASSERT(!marked_empty() && !y.marked_empty() && space_dim > 0);
430
431 const Polyhedron& x = *this;
432
433 // `x' cannot have pending constraints, because we need its generators.
434 if (x.has_pending_constraints() && !x.process_pending_constraints()) {
435 return true;
436 }
437 // `y' cannot have pending generators, because we need its constraints.
438 if (y.has_pending_generators()) {
439 y.process_pending_generators();
440 }
441
442 #if BE_LAZY
443 if (!x.generators_are_up_to_date() && !x.update_generators()) {
444 return true;
445 }
446 if (!y.constraints_are_up_to_date()) {
447 y.update_constraints();
448 }
449 #else
450 if (!x.generators_are_minimized()) {
451 x.minimize();
452 }
453 if (!y.constraints_are_minimized()) {
454 y.minimize();
455 }
456 #endif
457
458 PPL_ASSERT_HEAVY(x.OK());
459 PPL_ASSERT_HEAVY(y.OK());
460
461 const Generator_System& gs = x.gen_sys;
462 const Constraint_System& cs = y.con_sys;
463
464 if (x.is_necessarily_closed()) {
465 // When working with necessarily closed polyhedra,
466 // `x' is contained in `y' if and only if all the generators of `x'
467 // satisfy all the inequalities and saturate all the equalities of `y'.
468 // This comes from the definition of a polyhedron as the set of
469 // vectors satisfying a constraint system and the fact that all
470 // vectors in `x' can be obtained by suitably combining its generators.
471 for (dimension_type i = cs.num_rows(); i-- > 0; ) {
472 const Constraint& c = cs[i];
473 if (c.is_inequality()) {
474 for (dimension_type j = gs.num_rows(); j-- > 0; ) {
475 const Generator& g = gs[j];
476 const int sp_sign = Scalar_Products::sign(c, g);
477 if (g.is_line()) {
478 if (sp_sign != 0) {
479 return false;
480 }
481 }
482 else {
483 // `g' is a ray or a point.
484 if (sp_sign < 0) {
485 return false;
486 }
487 }
488 }
489 }
490 else {
491 // `c' is an equality.
492 for (dimension_type j = gs.num_rows(); j-- > 0; ) {
493 if (Scalar_Products::sign(c, gs[j]) != 0) {
494 return false;
495 }
496 }
497 }
498 }
499 }
500 else {
501 // Here we have an NNC polyhedron: using the reduced scalar product,
502 // which ignores the epsilon coefficient.
503 for (dimension_type i = cs.num_rows(); i-- > 0; ) {
504 const Constraint& c = cs[i];
505 switch (c.type()) {
506 case Constraint::NONSTRICT_INEQUALITY:
507 for (dimension_type j = gs.num_rows(); j-- > 0; ) {
508 const Generator& g = gs[j];
509 const int sp_sign = Scalar_Products::reduced_sign(c, g);
510 if (g.is_line()) {
511 if (sp_sign != 0) {
512 return false;
513 }
514 }
515 else {
516 // `g' is a ray or a point or a closure point.
517 if (sp_sign < 0) {
518 return false;
519 }
520 }
521 }
522 break;
523 case Constraint::EQUALITY:
524 for (dimension_type j = gs.num_rows(); j-- > 0; ) {
525 if (Scalar_Products::reduced_sign(c, gs[j]) != 0) {
526 return false;
527 }
528 }
529 break;
530 case Constraint::STRICT_INEQUALITY:
531 for (dimension_type j = gs.num_rows(); j-- > 0; ) {
532 const Generator& g = gs[j];
533 const int sp_sign = Scalar_Products::reduced_sign(c, g);
534 switch (g.type()) {
535 case Generator::POINT:
536 // If a point violates or saturates a strict inequality
537 // (when ignoring the epsilon coefficients) then it is
538 // not included in the polyhedron.
539 if (sp_sign <= 0) {
540 return false;
541 }
542 break;
543 case Generator::LINE:
544 // Lines have to saturate all constraints.
545 if (sp_sign != 0) {
546 return false;
547 }
548 break;
549 case Generator::RAY:
550 // Intentionally fall through.
551 case Generator::CLOSURE_POINT:
552 // The generator is a ray or closure point: usual test.
553 if (sp_sign < 0) {
554 return false;
555 }
556 break;
557 }
558 }
559 break;
560 }
561 }
562 }
563
564 // Inclusion holds.
565 return true;
566 }
567
568 bool
bounds(const Linear_Expression & expr,const bool from_above) const569 PPL::Polyhedron::bounds(const Linear_Expression& expr,
570 const bool from_above) const {
571 // The dimension of `expr' should not be greater than the dimension
572 // of `*this'.
573 const dimension_type expr_space_dim = expr.space_dimension();
574 if (space_dim < expr_space_dim) {
575 throw_dimension_incompatible((from_above
576 ? "bounds_from_above(e)"
577 : "bounds_from_below(e)"), "e", expr);
578 }
579
580 // A zero-dimensional or empty polyhedron bounds everything.
581 if (space_dim == 0
582 || marked_empty()
583 || (has_pending_constraints() && !process_pending_constraints())
584 || (!generators_are_up_to_date() && !update_generators())) {
585 return true;
586 }
587 // The polyhedron has updated, possibly pending generators.
588 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
589 const Generator& g = gen_sys[i];
590 // Only lines and rays in `*this' can cause `expr' to be unbounded.
591 if (g.is_line_or_ray()) {
592 const int sp_sign = Scalar_Products::homogeneous_sign(expr, g);
593 if (sp_sign != 0
594 && (g.is_line()
595 || (from_above && sp_sign > 0)
596 || (!from_above && sp_sign < 0))) {
597 // `*this' does not bound `expr'.
598 return false;
599 }
600 }
601 }
602 // No sources of unboundedness have been found for `expr'
603 // in the given direction.
604 return true;
605 }
606
607 bool
max_min(const Linear_Expression & expr,const bool maximize,Coefficient & ext_n,Coefficient & ext_d,bool & included,Generator & g) const608 PPL::Polyhedron::max_min(const Linear_Expression& expr,
609 const bool maximize,
610 Coefficient& ext_n, Coefficient& ext_d,
611 bool& included,
612 Generator& g) const {
613 // The dimension of `expr' should not be greater than the dimension
614 // of `*this'.
615 const dimension_type expr_space_dim = expr.space_dimension();
616 if (space_dim < expr_space_dim) {
617 throw_dimension_incompatible((maximize
618 ? "maximize(e, ...)"
619 : "minimize(e, ...)"), "e", expr);
620 }
621
622 // Deal with zero-dim polyhedra first.
623 if (space_dim == 0) {
624 if (marked_empty()) {
625 return false;
626 }
627 else {
628 ext_n = expr.inhomogeneous_term();
629 ext_d = 1;
630 included = true;
631 g = point();
632 return true;
633 }
634 }
635
636 // For an empty polyhedron we simply return false.
637 if (marked_empty()
638 || (has_pending_constraints() && !process_pending_constraints())
639 || (!generators_are_up_to_date() && !update_generators())) {
640 return false;
641 }
642
643 // The polyhedron has updated, possibly pending generators.
644 // The following loop will iterate through the generator
645 // to find the extremum.
646 PPL_DIRTY_TEMP(mpq_class, extremum);
647
648 // True if we have no other candidate extremum to compare with.
649 bool first_candidate = true;
650
651 // To store the position of the current candidate extremum.
652 PPL_UNINITIALIZED(dimension_type, ext_position);
653
654 // Whether the current candidate extremum is included or not.
655 PPL_UNINITIALIZED(bool, ext_included);
656
657 PPL_DIRTY_TEMP_COEFFICIENT(sp);
658 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
659 const Generator& gen_sys_i = gen_sys[i];
660 Scalar_Products::homogeneous_assign(sp, expr, gen_sys_i);
661 // Lines and rays in `*this' can cause `expr' to be unbounded.
662 if (gen_sys_i.is_line_or_ray()) {
663 const int sp_sign = sgn(sp);
664 if (sp_sign != 0
665 && (gen_sys_i.is_line()
666 || (maximize && sp_sign > 0)
667 || (!maximize && sp_sign < 0))) {
668 // `expr' is unbounded in `*this'.
669 return false;
670 }
671 }
672 else {
673 // We have a point or a closure point.
674 PPL_ASSERT(gen_sys_i.is_point() || gen_sys_i.is_closure_point());
675 // Notice that we are ignoring the constant term in `expr' here.
676 // We will add it to the extremum as soon as we find it.
677 PPL_DIRTY_TEMP(mpq_class, candidate);
678 assign_r(candidate.get_num(), sp, ROUND_NOT_NEEDED);
679 assign_r(candidate.get_den(), gen_sys_i.expr.inhomogeneous_term(), ROUND_NOT_NEEDED);
680 candidate.canonicalize();
681 const bool g_is_point = gen_sys_i.is_point();
682 if (first_candidate
683 || (maximize
684 && (candidate > extremum
685 || (g_is_point
686 && !ext_included
687 && candidate == extremum)))
688 || (!maximize
689 && (candidate < extremum
690 || (g_is_point
691 && !ext_included
692 && candidate == extremum)))) {
693 // We have a (new) candidate extremum.
694 first_candidate = false;
695 extremum = candidate;
696 ext_position = i;
697 ext_included = g_is_point;
698 }
699 }
700 }
701
702 // Add in the constant term in `expr'.
703 PPL_DIRTY_TEMP(mpz_class, n);
704 assign_r(n, expr.inhomogeneous_term(), ROUND_NOT_NEEDED);
705 extremum += n;
706
707 // The polyhedron is bounded in the right direction and we have
708 // computed the extremum: write the result into the caller's structures.
709 PPL_ASSERT(!first_candidate);
710 ext_n = extremum.get_num();
711 ext_d = extremum.get_den();
712 included = ext_included;
713 g = gen_sys[ext_position];
714
715 return true;
716 }
717
718 void
set_zero_dim_univ()719 PPL::Polyhedron::set_zero_dim_univ() {
720 status.set_zero_dim_univ();
721 space_dim = 0;
722 con_sys.clear();
723 gen_sys.clear();
724 }
725
726 void
set_empty()727 PPL::Polyhedron::set_empty() {
728 status.set_empty();
729 // The polyhedron is empty: we can thus throw away everything.
730 con_sys.clear();
731 gen_sys.clear();
732 sat_c.clear();
733 sat_g.clear();
734 }
735
736 bool
process_pending_constraints() const737 PPL::Polyhedron::process_pending_constraints() const {
738 PPL_ASSERT(space_dim > 0 && !marked_empty());
739 PPL_ASSERT(has_pending_constraints() && !has_pending_generators());
740
741 Polyhedron& x = const_cast<Polyhedron&>(*this);
742
743 // Integrate the pending part of the system of constraints and minimize.
744 // We need `sat_c' up-to-date and `con_sys' sorted (together with `sat_c').
745 if (!x.sat_c_is_up_to_date()) {
746 x.sat_c.transpose_assign(x.sat_g);
747 }
748 if (!x.con_sys.is_sorted()) {
749 x.obtain_sorted_constraints_with_sat_c();
750 }
751 // We sort in place the pending constraints, erasing those constraints
752 // that also occur in the non-pending part of `con_sys'.
753 x.con_sys.sort_pending_and_remove_duplicates();
754 if (x.con_sys.num_pending_rows() == 0) {
755 // All pending constraints were duplicates.
756 x.clear_pending_constraints();
757 PPL_ASSERT_HEAVY(OK(true));
758 return true;
759 }
760
761 const bool empty = add_and_minimize(true, x.con_sys, x.gen_sys, x.sat_c);
762 PPL_ASSERT(x.con_sys.num_pending_rows() == 0);
763
764 if (empty) {
765 x.set_empty();
766 }
767 else {
768 x.clear_pending_constraints();
769 x.clear_sat_g_up_to_date();
770 x.set_sat_c_up_to_date();
771 }
772 PPL_ASSERT_HEAVY(OK(!empty));
773 return !empty;
774 }
775
776 void
process_pending_generators() const777 PPL::Polyhedron::process_pending_generators() const {
778 PPL_ASSERT(space_dim > 0 && !marked_empty());
779 PPL_ASSERT(has_pending_generators() && !has_pending_constraints());
780
781 Polyhedron& x = const_cast<Polyhedron&>(*this);
782
783 // Integrate the pending part of the system of generators and minimize.
784 // We need `sat_g' up-to-date and `gen_sys' sorted (together with `sat_g').
785 if (!x.sat_g_is_up_to_date()) {
786 x.sat_g.transpose_assign(x.sat_c);
787 }
788 if (!x.gen_sys.is_sorted()) {
789 x.obtain_sorted_generators_with_sat_g();
790 }
791 // We sort in place the pending generators, erasing those generators
792 // that also occur in the non-pending part of `gen_sys'.
793 x.gen_sys.sort_pending_and_remove_duplicates();
794 if (x.gen_sys.num_pending_rows() == 0) {
795 // All pending generators were duplicates.
796 x.clear_pending_generators();
797 PPL_ASSERT_HEAVY(OK(true));
798 return;
799 }
800
801 add_and_minimize(false, x.gen_sys, x.con_sys, x.sat_g);
802 PPL_ASSERT(x.gen_sys.num_pending_rows() == 0);
803
804 x.clear_pending_generators();
805 x.clear_sat_c_up_to_date();
806 x.set_sat_g_up_to_date();
807 }
808
809 void
remove_pending_to_obtain_constraints() const810 PPL::Polyhedron::remove_pending_to_obtain_constraints() const {
811 PPL_ASSERT(has_something_pending());
812
813 Polyhedron& x = const_cast<Polyhedron&>(*this);
814
815 // If the polyhedron has pending constraints, simply unset them.
816 if (x.has_pending_constraints()) {
817 // Integrate the pending constraints, which are possibly not sorted.
818 x.con_sys.set_sorted(false);
819 x.con_sys.unset_pending_rows();
820 x.clear_pending_constraints();
821 x.clear_constraints_minimized();
822 x.clear_generators_up_to_date();
823 }
824 else {
825 PPL_ASSERT(x.has_pending_generators());
826 // We must process the pending generators and obtain the
827 // corresponding system of constraints.
828 x.process_pending_generators();
829 }
830 PPL_ASSERT_HEAVY(OK(true));
831 }
832
833 bool
remove_pending_to_obtain_generators() const834 PPL::Polyhedron::remove_pending_to_obtain_generators() const {
835 PPL_ASSERT(has_something_pending());
836
837 Polyhedron& x = const_cast<Polyhedron&>(*this);
838
839 // If the polyhedron has pending generators, simply unset them.
840 if (x.has_pending_generators()) {
841 // Integrate the pending generators, which are possibly not sorted.
842 x.gen_sys.set_sorted(false);
843 x.gen_sys.unset_pending_rows();
844 x.clear_pending_generators();
845 x.clear_generators_minimized();
846 x.clear_constraints_up_to_date();
847 PPL_ASSERT_HEAVY(OK(true));
848 return true;
849 }
850 else {
851 PPL_ASSERT(x.has_pending_constraints());
852 // We must integrate the pending constraints and obtain the
853 // corresponding system of generators.
854 return x.process_pending_constraints();
855 }
856 }
857
858 void
update_constraints() const859 PPL::Polyhedron::update_constraints() const {
860 PPL_ASSERT(space_dim > 0);
861 PPL_ASSERT(!marked_empty());
862 PPL_ASSERT(generators_are_up_to_date());
863 // We assume the polyhedron has no pending constraints or generators.
864 PPL_ASSERT(!has_something_pending());
865
866 Polyhedron& x = const_cast<Polyhedron&>(*this);
867 minimize(false, x.gen_sys, x.con_sys, x.sat_c);
868 // `sat_c' is the only saturation matrix up-to-date.
869 x.set_sat_c_up_to_date();
870 x.clear_sat_g_up_to_date();
871 // The system of constraints and the system of generators
872 // are minimized.
873 x.set_constraints_minimized();
874 x.set_generators_minimized();
875 }
876
877 bool
update_generators() const878 PPL::Polyhedron::update_generators() const {
879 PPL_ASSERT(space_dim > 0);
880 PPL_ASSERT(!marked_empty());
881 PPL_ASSERT(constraints_are_up_to_date());
882 // We assume the polyhedron has no pending constraints or generators.
883 PPL_ASSERT(!has_something_pending());
884
885 Polyhedron& x = const_cast<Polyhedron&>(*this);
886 // If the system of constraints is not consistent the
887 // polyhedron is empty.
888 const bool empty = minimize(true, x.con_sys, x.gen_sys, x.sat_g);
889 if (empty) {
890 x.set_empty();
891 }
892 else {
893 // `sat_g' is the only saturation matrix up-to-date.
894 x.set_sat_g_up_to_date();
895 x.clear_sat_c_up_to_date();
896 // The system of constraints and the system of generators
897 // are minimized.
898 x.set_constraints_minimized();
899 x.set_generators_minimized();
900 }
901 return !empty;
902 }
903
904 void
update_sat_c() const905 PPL::Polyhedron::update_sat_c() const {
906 PPL_ASSERT(constraints_are_minimized());
907 PPL_ASSERT(generators_are_minimized());
908 PPL_ASSERT(!sat_c_is_up_to_date());
909
910 // We only consider non-pending rows.
911 const dimension_type csr = con_sys.first_pending_row();
912 const dimension_type gsr = gen_sys.first_pending_row();
913 Polyhedron& x = const_cast<Polyhedron&>(*this);
914
915 // The columns of `sat_c' represent the constraints and
916 // its rows represent the generators: resize accordingly.
917 x.sat_c.resize(gsr, csr);
918 for (dimension_type i = gsr; i-- > 0; ) {
919 for (dimension_type j = csr; j-- > 0; ) {
920 const int sp_sign = Scalar_Products::sign(con_sys[j], gen_sys[i]);
921 // The negativity of this scalar product would mean
922 // that the generator `gen_sys[i]' violates the constraint
923 // `con_sys[j]' and it is not possible because both generators
924 // and constraints are up-to-date.
925 PPL_ASSERT(sp_sign >= 0);
926 if (sp_sign > 0) {
927 // `gen_sys[i]' satisfies (without saturate) `con_sys[j]'.
928 x.sat_c[i].set(j);
929 }
930 else {
931 // `gen_sys[i]' saturates `con_sys[j]'.
932 x.sat_c[i].clear(j);
933 }
934 }
935 }
936 x.set_sat_c_up_to_date();
937 }
938
939 void
update_sat_g() const940 PPL::Polyhedron::update_sat_g() const {
941 PPL_ASSERT(constraints_are_minimized());
942 PPL_ASSERT(generators_are_minimized());
943 PPL_ASSERT(!sat_g_is_up_to_date());
944
945 // We only consider non-pending rows.
946 const dimension_type csr = con_sys.first_pending_row();
947 const dimension_type gsr = gen_sys.first_pending_row();
948 Polyhedron& x = const_cast<Polyhedron&>(*this);
949
950 // The columns of `sat_g' represent generators and its
951 // rows represent the constraints: resize accordingly.
952 x.sat_g.resize(csr, gsr);
953 for (dimension_type i = csr; i-- > 0; ) {
954 for (dimension_type j = gsr; j-- > 0; ) {
955 const int sp_sign = Scalar_Products::sign(con_sys[i], gen_sys[j]);
956 // The negativity of this scalar product would mean
957 // that the generator `gen_sys[j]' violates the constraint
958 // `con_sys[i]' and it is not possible because both generators
959 // and constraints are up-to-date.
960 PPL_ASSERT(sp_sign >= 0);
961 if (sp_sign > 0) {
962 // `gen_sys[j]' satisfies (without saturate) `con_sys[i]'.
963 x.sat_g[i].set(j);
964 }
965 else {
966 // `gen_sys[j]' saturates `con_sys[i]'.
967 x.sat_g[i].clear(j);
968 }
969 }
970 }
971 x.set_sat_g_up_to_date();
972 }
973
974 void
obtain_sorted_constraints() const975 PPL::Polyhedron::obtain_sorted_constraints() const {
976 PPL_ASSERT(constraints_are_up_to_date());
977 // `con_sys' will be sorted up to `index_first_pending'.
978 Polyhedron& x = const_cast<Polyhedron&>(*this);
979 if (!x.con_sys.is_sorted()) {
980 if (x.sat_g_is_up_to_date()) {
981 // Sorting constraints keeping `sat_g' consistent.
982 x.con_sys.sort_and_remove_with_sat(x.sat_g);
983 // `sat_c' is not up-to-date anymore.
984 x.clear_sat_c_up_to_date();
985 }
986 else if (x.sat_c_is_up_to_date()) {
987 // Using `sat_c' to obtain `sat_g', then it is like previous case.
988 x.sat_g.transpose_assign(x.sat_c);
989 x.con_sys.sort_and_remove_with_sat(x.sat_g);
990 x.set_sat_g_up_to_date();
991 x.clear_sat_c_up_to_date();
992 }
993 else {
994 // If neither `sat_g' nor `sat_c' are up-to-date,
995 // we just sort the constraints.
996 x.con_sys.sort_rows();
997 }
998 }
999
1000 PPL_ASSERT(con_sys.check_sorted());
1001 }
1002
1003 void
obtain_sorted_generators() const1004 PPL::Polyhedron::obtain_sorted_generators() const {
1005 PPL_ASSERT(generators_are_up_to_date());
1006 // `gen_sys' will be sorted up to `index_first_pending'.
1007 Polyhedron& x = const_cast<Polyhedron&>(*this);
1008 if (!x.gen_sys.is_sorted()) {
1009 if (x.sat_c_is_up_to_date()) {
1010 // Sorting generators keeping 'sat_c' consistent.
1011 x.gen_sys.sort_and_remove_with_sat(x.sat_c);
1012 // `sat_g' is not up-to-date anymore.
1013 x.clear_sat_g_up_to_date();
1014 }
1015 else if (x.sat_g_is_up_to_date()) {
1016 // Obtaining `sat_c' from `sat_g' and proceeding like previous case.
1017 x.sat_c.transpose_assign(x.sat_g);
1018 x.gen_sys.sort_and_remove_with_sat(x.sat_c);
1019 x.set_sat_c_up_to_date();
1020 x.clear_sat_g_up_to_date();
1021 }
1022 else {
1023 // If neither `sat_g' nor `sat_c' are up-to-date, we just sort
1024 // the generators.
1025 x.gen_sys.sort_rows();
1026 }
1027 }
1028
1029 PPL_ASSERT(gen_sys.check_sorted());
1030 }
1031
1032 void
obtain_sorted_constraints_with_sat_c() const1033 PPL::Polyhedron::obtain_sorted_constraints_with_sat_c() const {
1034 PPL_ASSERT(constraints_are_up_to_date());
1035 PPL_ASSERT(constraints_are_minimized());
1036 // `con_sys' will be sorted up to `index_first_pending'.
1037 Polyhedron& x = const_cast<Polyhedron&>(*this);
1038 // At least one of the saturation matrices must be up-to-date.
1039 if (!x.sat_c_is_up_to_date() && !x.sat_g_is_up_to_date()) {
1040 x.update_sat_c();
1041 }
1042 if (x.con_sys.is_sorted()) {
1043 if (x.sat_c_is_up_to_date()) {
1044 // If constraints are already sorted and sat_c is up to
1045 // date there is nothing to do.
1046 return;
1047 }
1048 }
1049 else {
1050 if (!x.sat_g_is_up_to_date()) {
1051 // If constraints are not sorted and sat_g is not up-to-date
1052 // we obtain sat_g from sat_c (that has to be up-to-date)...
1053 x.sat_g.transpose_assign(x.sat_c);
1054 x.set_sat_g_up_to_date();
1055 }
1056 // ... and sort it together with constraints.
1057 x.con_sys.sort_and_remove_with_sat(x.sat_g);
1058 }
1059 // Obtaining sat_c from sat_g.
1060 x.sat_c.transpose_assign(x.sat_g);
1061 x.set_sat_c_up_to_date();
1062 // Constraints are sorted now.
1063 x.con_sys.set_sorted(true);
1064
1065 PPL_ASSERT(con_sys.check_sorted());
1066 }
1067
1068 void
obtain_sorted_generators_with_sat_g() const1069 PPL::Polyhedron::obtain_sorted_generators_with_sat_g() const {
1070 PPL_ASSERT(generators_are_up_to_date());
1071 // `gen_sys' will be sorted up to `index_first_pending'.
1072 Polyhedron& x = const_cast<Polyhedron&>(*this);
1073 // At least one of the saturation matrices must be up-to-date.
1074 if (!x.sat_c_is_up_to_date() && !x.sat_g_is_up_to_date()) {
1075 x.update_sat_g();
1076 }
1077
1078 if (x.gen_sys.is_sorted()) {
1079 if (x.sat_g_is_up_to_date()) {
1080 // If generators are already sorted and sat_g is up to
1081 // date there is nothing to do.
1082 return;
1083 }
1084 }
1085 else {
1086 if (!x.sat_c_is_up_to_date()) {
1087 // If generators are not sorted and sat_c is not up-to-date
1088 // we obtain sat_c from sat_g (that has to be up-to-date)...
1089 x.sat_c.transpose_assign(x.sat_g);
1090 x.set_sat_c_up_to_date();
1091 }
1092 // ... and sort it together with generators.
1093 x.gen_sys.sort_and_remove_with_sat(x.sat_c);
1094 }
1095 // Obtaining sat_g from sat_c.
1096 x.sat_g.transpose_assign(sat_c);
1097 x.set_sat_g_up_to_date();
1098 // Generators are sorted now.
1099 x.gen_sys.set_sorted(true);
1100
1101 PPL_ASSERT(gen_sys.check_sorted());
1102 }
1103
1104 bool
minimize() const1105 PPL::Polyhedron::minimize() const {
1106 // 0-dim space or empty polyhedra are already minimized.
1107 if (marked_empty()) {
1108 return false;
1109 }
1110 if (space_dim == 0) {
1111 return true;
1112 }
1113
1114 // If the polyhedron has something pending, process it.
1115 if (has_something_pending()) {
1116 const bool not_empty = process_pending();
1117 PPL_ASSERT_HEAVY(OK());
1118 return not_empty;
1119 }
1120
1121 // Here there are no pending constraints or generators.
1122 // Is the polyhedron already minimized?
1123 if (constraints_are_minimized() && generators_are_minimized()) {
1124 return true;
1125 }
1126
1127 // If constraints or generators are up-to-date, invoking
1128 // update_generators() or update_constraints(), respectively,
1129 // minimizes both constraints and generators.
1130 // If both are up-to-date it does not matter whether we use
1131 // update_generators() or update_constraints():
1132 // both minimize constraints and generators.
1133 if (constraints_are_up_to_date()) {
1134 // We may discover here that `*this' is empty.
1135 const bool ret = update_generators();
1136 PPL_ASSERT_HEAVY(OK());
1137 return ret;
1138 }
1139 else {
1140 PPL_ASSERT(generators_are_up_to_date());
1141 update_constraints();
1142 PPL_ASSERT_HEAVY(OK());
1143 return true;
1144 }
1145 }
1146
1147 bool
strongly_minimize_constraints() const1148 PPL::Polyhedron::strongly_minimize_constraints() const {
1149 PPL_ASSERT(!is_necessarily_closed());
1150
1151 // From the user perspective, the polyhedron will not change.
1152 Polyhedron& x = const_cast<Polyhedron&>(*this);
1153
1154 // We need `con_sys' (weakly) minimized and `gen_sys' up-to-date.
1155 // `minimize()' will process any pending constraints or generators.
1156 if (!minimize()) {
1157 return false;
1158 }
1159 // If the polyhedron `*this' is zero-dimensional
1160 // at this point it must be a universe polyhedron.
1161 if (x.space_dim == 0) {
1162 return true;
1163 }
1164 // We also need `sat_g' up-to-date.
1165 if (!sat_g_is_up_to_date()) {
1166 PPL_ASSERT(sat_c_is_up_to_date());
1167 x.sat_g.transpose_assign(sat_c);
1168 }
1169
1170 // These Bit_Row's will be later used as masks in order to
1171 // check saturation conditions restricted to particular subsets of
1172 // the generator system.
1173 Bit_Row sat_all_but_rays;
1174 Bit_Row sat_all_but_points;
1175 Bit_Row sat_all_but_closure_points;
1176
1177 const dimension_type gs_rows = gen_sys.num_rows();
1178 const dimension_type n_lines = gen_sys.num_lines();
1179 for (dimension_type i = gs_rows; i-- > n_lines; ) {
1180 switch (gen_sys[i].type()) {
1181 case Generator::RAY:
1182 sat_all_but_rays.set(i);
1183 break;
1184 case Generator::POINT:
1185 sat_all_but_points.set(i);
1186 break;
1187 case Generator::CLOSURE_POINT:
1188 sat_all_but_closure_points.set(i);
1189 break;
1190 case Generator::LINE:
1191 // Found a line with index i >= n_lines !
1192 PPL_UNREACHABLE;
1193 break;
1194 }
1195 }
1196 const Bit_Row
1197 sat_lines_and_rays(sat_all_but_points, sat_all_but_closure_points);
1198 const Bit_Row
1199 sat_lines_and_closure_points(sat_all_but_rays, sat_all_but_points);
1200 const Bit_Row
1201 sat_lines(sat_lines_and_rays, sat_lines_and_closure_points);
1202
1203 // These flags are maintained to later decide if we have to add the
1204 // eps_leq_one constraint and whether or not the constraint system
1205 // was changed.
1206 bool changed = false;
1207 bool found_eps_leq_one = false;
1208
1209 // For all the strict inequalities in `con_sys', check for
1210 // eps-redundancy and eventually move them to the bottom part of the
1211 // system.
1212 Constraint_System& cs = x.con_sys;
1213 Bit_Matrix& sat = x.sat_g;
1214 const Variable eps_var(cs.space_dimension());
1215 // Note that cs.num_rows() is *not* constant because the calls to
1216 // cs.remove_row() decrease it.
1217 for (dimension_type i = 0; i < cs.num_rows(); ) {
1218 if (cs[i].is_strict_inequality()) {
1219 // First, check if it is saturated by no closure points
1220 Bit_Row sat_ci;
1221 sat_ci.union_assign(sat[i], sat_lines_and_closure_points);
1222 if (sat_ci == sat_lines) {
1223 // It is saturated by no closure points.
1224 if (!found_eps_leq_one) {
1225 // Check if it is the eps_leq_one constraint.
1226 const Constraint& c = cs[i];
1227 if (c.expression().all_homogeneous_terms_are_zero()
1228 && (c.expression().inhomogeneous_term() + c.epsilon_coefficient() == 0)) {
1229 // We found the eps_leq_one constraint.
1230 found_eps_leq_one = true;
1231 // Consider next constraint.
1232 ++i;
1233 continue;
1234 }
1235 }
1236 // Here `cs[i]' is not the eps_leq_one constraint,
1237 // so it is eps-redundant.
1238 // Remove it, while keeping `sat_g' consistent.
1239 cs.remove_row(i, false);
1240 swap(sat[i], sat[cs.num_rows()]);
1241 // The constraint system is changed.
1242 changed = true;
1243 // Continue by considering next constraint,
1244 // which is already in place due to the swap.
1245 continue;
1246 }
1247 // Now we check if there exists another strict inequality
1248 // constraint having a superset of its saturators,
1249 // when disregarding points.
1250 sat_ci.union_assign(sat[i], sat_all_but_points);
1251 bool eps_redundant = false;
1252 for (dimension_type j = 0; j < cs.num_rows(); ++j) {
1253 if (i != j && cs[j].is_strict_inequality()
1254 && subset_or_equal(sat[j], sat_ci)) {
1255 // Constraint `cs[i]' is eps-redundant:
1256 // remove it, while keeping `sat_g' consistent.
1257 cs.remove_row(i, false);
1258 swap(sat[i], sat[cs.num_rows()]);
1259 eps_redundant = true;
1260 // The constraint system is changed.
1261 changed = true;
1262 break;
1263 }
1264 }
1265 // Continue with next constraint, which is already in place
1266 // due to the swap if we have found an eps-redundant constraint.
1267 if (!eps_redundant) {
1268 ++i;
1269 }
1270 }
1271 else {
1272 // `cs[i]' is not a strict inequality: consider next constraint.
1273 ++i;
1274 }
1275 }
1276
1277 PPL_ASSERT(cs.num_pending_rows() == 0);
1278
1279 if (changed) {
1280 // If the constraint system has been changed, we have erased
1281 // the epsilon-redundant constraints.
1282
1283 // The generator system is no longer up-to-date.
1284 x.clear_generators_up_to_date();
1285
1286 // If we haven't found an upper bound for the epsilon dimension,
1287 // then we have to check whether such an upper bound is implied
1288 // by the remaining constraints (exploiting the simplex algorithm).
1289 if (!found_eps_leq_one) {
1290 MIP_Problem lp;
1291 // KLUDGE: temporarily mark the constraint system as if it was
1292 // necessarily closed, so that we can interpret the epsilon
1293 // dimension as a standard dimension. Be careful to reset the
1294 // topology of `cs' even on exceptional execution path.
1295 cs.mark_as_necessarily_closed();
1296 try {
1297 lp.add_space_dimensions_and_embed(cs.space_dimension());
1298 lp.add_constraints(cs);
1299 cs.mark_as_not_necessarily_closed();
1300 }
1301 catch (...) {
1302 cs.mark_as_not_necessarily_closed();
1303 throw;
1304 }
1305 // The objective function is `epsilon'.
1306 lp.set_objective_function(Variable(x.space_dim));
1307 lp.set_optimization_mode(MAXIMIZATION);
1308 const MIP_Problem_Status status = lp.solve();
1309 PPL_ASSERT(status != UNFEASIBLE_MIP_PROBLEM);
1310 // If the epsilon dimension is actually unbounded,
1311 // then add the eps_leq_one constraint.
1312 if (status == UNBOUNDED_MIP_PROBLEM) {
1313 cs.insert(Constraint::epsilon_leq_one());
1314 }
1315 }
1316 }
1317
1318 PPL_ASSERT_HEAVY(OK());
1319 return true;
1320 }
1321
1322 bool
strongly_minimize_generators() const1323 PPL::Polyhedron::strongly_minimize_generators() const {
1324 PPL_ASSERT(!is_necessarily_closed());
1325
1326 // From the user perspective, the polyhedron will not change.
1327 Polyhedron& x = const_cast<Polyhedron&>(*this);
1328
1329 // We need `gen_sys' (weakly) minimized and `con_sys' up-to-date.
1330 // `minimize()' will process any pending constraints or generators.
1331 if (!minimize()) {
1332 return false;
1333 }
1334 // If the polyhedron `*this' is zero-dimensional
1335 // at this point it must be a universe polyhedron.
1336 if (x.space_dim == 0) {
1337 return true;
1338 }
1339
1340 // We also need `sat_c' up-to-date.
1341 if (!sat_c_is_up_to_date()) {
1342 PPL_ASSERT(sat_g_is_up_to_date());
1343 x.sat_c.transpose_assign(sat_g);
1344 }
1345
1346 // This Bit_Row will have all and only the indexes
1347 // of strict inequalities set to 1.
1348 Bit_Row sat_all_but_strict_ineq;
1349 const dimension_type cs_rows = con_sys.num_rows();
1350 const dimension_type n_equals = con_sys.num_equalities();
1351 for (dimension_type i = cs_rows; i-- > n_equals; ) {
1352 if (con_sys[i].is_strict_inequality()) {
1353 sat_all_but_strict_ineq.set(i);
1354 }
1355 }
1356
1357 // Will record whether or not we changed the generator system.
1358 bool changed = false;
1359
1360 // For all points in the generator system, check for eps-redundancy
1361 // and eventually move them to the bottom part of the system.
1362 Generator_System& gs = const_cast<Generator_System&>(gen_sys);
1363 Bit_Matrix& sat = const_cast<Bit_Matrix&>(sat_c);
1364 const dimension_type old_gs_rows = gs.num_rows();
1365 dimension_type gs_rows = old_gs_rows;
1366 const dimension_type n_lines = gs.num_lines();
1367 bool gs_sorted = gs.is_sorted();
1368
1369 for (dimension_type i = n_lines; i < gs_rows; ) {
1370 Generator& g = gs.sys.rows[i];
1371 if (g.is_point()) {
1372 // Compute the Bit_Row corresponding to the candidate point
1373 // when strict inequality constraints are ignored.
1374 const Bit_Row sat_gs_i(sat[i], sat_all_but_strict_ineq);
1375 // Check if the candidate point is actually eps-redundant:
1376 // namely, if there exists another point that saturates
1377 // all the non-strict inequalities saturated by the candidate.
1378 bool eps_redundant = false;
1379 for (dimension_type j = n_lines; j < gs_rows; ++j) {
1380 const Generator& g2 = gs.sys.rows[j];
1381 if (i != j && g2.is_point() && subset_or_equal(sat[j], sat_gs_i)) {
1382 // Point `g' is eps-redundant:
1383 // move it to the bottom of the generator system,
1384 // while keeping `sat_c' consistent.
1385 --gs_rows;
1386 swap(g, gs.sys.rows[gs_rows]);
1387 swap(sat[i], sat[gs_rows]);
1388 eps_redundant = true;
1389 changed = true;
1390 break;
1391 }
1392 }
1393 if (!eps_redundant) {
1394 // Let all point encodings have epsilon coordinate 1.
1395 if (g.epsilon_coefficient() != g.expr.inhomogeneous_term()) {
1396 g.set_epsilon_coefficient(g.expr.inhomogeneous_term());
1397 // Enforce normalization.
1398 g.expr.normalize();
1399 PPL_ASSERT(g.OK());
1400 changed = true;
1401 }
1402 // Consider next generator.
1403 ++i;
1404 }
1405 }
1406 else {
1407 // Consider next generator.
1408 ++i;
1409 }
1410 }
1411
1412 // If needed, erase the eps-redundant generators.
1413 if (gs_rows < old_gs_rows) {
1414 gs.sys.rows.resize(gs_rows);
1415 }
1416
1417 if (changed) {
1418 // The generator system is no longer sorted.
1419 gs_sorted = false;
1420 // The constraint system is no longer up-to-date.
1421 x.clear_constraints_up_to_date();
1422 }
1423
1424 gs.sys.index_first_pending = gs.num_rows();
1425 gs.set_sorted(gs_sorted);
1426
1427 PPL_ASSERT(gs.sys.OK());
1428
1429 PPL_ASSERT_HEAVY(OK());
1430 return true;
1431 }
1432
1433 void
refine_no_check(const Constraint & c)1434 PPL::Polyhedron::refine_no_check(const Constraint& c) {
1435 PPL_ASSERT(!marked_empty());
1436 PPL_ASSERT(space_dim >= c.space_dimension());
1437
1438 // Dealing with a zero-dimensional space polyhedron first.
1439 if (space_dim == 0) {
1440 if (c.is_inconsistent()) {
1441 set_empty();
1442 }
1443 return;
1444 }
1445
1446 // The constraints (possibly with pending rows) are required.
1447 if (has_pending_generators()) {
1448 process_pending_generators();
1449 }
1450 else if (!constraints_are_up_to_date()) {
1451 update_constraints();
1452 }
1453
1454 const bool adding_pending = can_have_something_pending();
1455
1456 if (c.is_necessarily_closed() || !is_necessarily_closed()) {
1457 // Since `con_sys' is not empty, the topology and space dimension
1458 // of the inserted constraint are automatically adjusted.
1459 if (adding_pending) {
1460 con_sys.insert_pending(c);
1461 }
1462 else {
1463 con_sys.insert(c);
1464 }
1465 }
1466 else {
1467 // Here we know that the system of constraints has at least a row.
1468 // However, by barely invoking `con_sys.insert(c)' we would
1469 // cause a change in the topology of `con_sys', which is wrong.
1470 // Thus, we insert a "topology corrected" copy of `c'.
1471 const Linear_Expression nc_expr(c.expression());
1472 if (c.is_equality()) {
1473 if (adding_pending) {
1474 con_sys.insert_pending(nc_expr == 0);
1475 }
1476 else {
1477 con_sys.insert(nc_expr == 0);
1478 }
1479 }
1480 else {
1481 if (adding_pending) {
1482 con_sys.insert_pending(nc_expr >= 0);
1483 }
1484 else {
1485 con_sys.insert(nc_expr >= 0);
1486 }
1487 }
1488 }
1489
1490 if (adding_pending) {
1491 set_constraints_pending();
1492 }
1493 else {
1494 // Constraints are not minimized and generators are not up-to-date.
1495 clear_constraints_minimized();
1496 clear_generators_up_to_date();
1497 }
1498
1499 // Note: the constraint system may have become unsatisfiable, thus
1500 // we do not check for satisfiability.
1501 PPL_ASSERT_HEAVY(OK());
1502 }
1503
1504 bool
BHZ09_poly_hull_assign_if_exact(const Polyhedron & y)1505 PPL::Polyhedron::BHZ09_poly_hull_assign_if_exact(const Polyhedron& y) {
1506 Polyhedron& x = *this;
1507
1508 // Private method: the caller must ensure the following.
1509 PPL_ASSERT(x.topology() == y.topology());
1510 PPL_ASSERT(x.space_dim == y.space_dim);
1511
1512 // The zero-dim case is trivial.
1513 if (x.space_dim == 0) {
1514 x.upper_bound_assign(y);
1515 return true;
1516 }
1517
1518 // If `x' or `y' are (known to be) empty, the upper bound is exact.
1519 if (x.marked_empty()) {
1520 x = y;
1521 return true;
1522 }
1523 else if (y.is_empty()) {
1524 return true;
1525 }
1526 else if (x.is_empty()) {
1527 x = y;
1528 return true;
1529 }
1530
1531 if (x.is_necessarily_closed()) {
1532 return x.BHZ09_C_poly_hull_assign_if_exact(y);
1533 }
1534 else {
1535 return x.BHZ09_NNC_poly_hull_assign_if_exact(y);
1536 }
1537 }
1538
1539 bool
BHZ09_C_poly_hull_assign_if_exact(const Polyhedron & y)1540 PPL::Polyhedron::BHZ09_C_poly_hull_assign_if_exact(const Polyhedron& y) {
1541 Polyhedron& x = *this;
1542 // Private method: the caller must ensure the following.
1543 PPL_ASSERT(x.is_necessarily_closed() && y.is_necessarily_closed());
1544 PPL_ASSERT(x.space_dim > 0 && x.space_dim == y.space_dim);
1545 PPL_ASSERT(!x.is_empty() && !y.is_empty());
1546
1547 // Minimization is not really required, but it is probably the best
1548 // way of getting constraints, generators and saturation matrices
1549 // up-to-date. Minimization it also removes redundant
1550 // constraints/generators.
1551 (void) x.minimize();
1552 (void) y.minimize();
1553
1554 // Handle a special case: for topologically closed polyhedra P and Q,
1555 // if the affine dimension of P is greater than that of Q, then
1556 // their upper bound is exact if and only if P includes Q.
1557 const dimension_type x_affine_dim = x.affine_dimension();
1558 const dimension_type y_affine_dim = y.affine_dimension();
1559 if (x_affine_dim > y_affine_dim) {
1560 return y.is_included_in(x);
1561 }
1562 else if (x_affine_dim < y_affine_dim) {
1563 if (x.is_included_in(y)) {
1564 x = y;
1565 return true;
1566 }
1567 else {
1568 return false;
1569 }
1570 }
1571
1572 const Constraint_System& x_cs = x.con_sys;
1573 const Generator_System& x_gs = x.gen_sys;
1574 const Generator_System& y_gs = y.gen_sys;
1575 const dimension_type x_gs_num_rows = x_gs.num_rows();
1576 const dimension_type y_gs_num_rows = y_gs.num_rows();
1577
1578 // Step 1: generators of `x' that are redundant in `y', and vice versa.
1579 Bit_Row x_gs_red_in_y;
1580 dimension_type num_x_gs_red_in_y = 0;
1581 for (dimension_type i = x_gs_num_rows; i-- > 0; ) {
1582 if (y.relation_with(x_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
1583 x_gs_red_in_y.set(i);
1584 ++num_x_gs_red_in_y;
1585 }
1586 }
1587
1588 Bit_Row y_gs_red_in_x;
1589 dimension_type num_y_gs_red_in_x = 0;
1590 for (dimension_type i = y_gs_num_rows; i-- > 0; ) {
1591 if (x.relation_with(y_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
1592 y_gs_red_in_x.set(i);
1593 ++num_y_gs_red_in_x;
1594 }
1595 }
1596 // Step 2: filter away special cases.
1597
1598 // Step 2.1: inclusion tests.
1599 if (num_y_gs_red_in_x == y_gs_num_rows) {
1600 // `y' is included into `x': upper bound `x' is exact.
1601 return true;
1602 }
1603 if (num_x_gs_red_in_y == x_gs_num_rows) {
1604 // `x' is included into `y': upper bound `y' is exact.
1605 x = y;
1606 return true;
1607 }
1608
1609 // Step 2.2: if no generator of `x' is redundant for `y', then
1610 // (as by 2.1 there exists a constraint of `x' non-redundant for `y')
1611 // the upper bound is not exact; the same if exchanging `x' and `y'.
1612 if (num_x_gs_red_in_y == 0 || num_y_gs_red_in_x == 0) {
1613 return false;
1614 }
1615
1616 // Step 3: see if `x' has a non-redundant constraint `c_x' that is not
1617 // satisfied by `y' and a non-redundant generator in `y' (see Step 1)
1618 // saturating `c_x'. If so, the upper bound is not exact.
1619
1620 // Make sure the saturation matrix for `x' is up to date.
1621 // Any sat matrix would do: we choose `sat_g' because it matches
1622 // the two nested loops (constraints on rows and generators on columns).
1623 if (!x.sat_g_is_up_to_date()) {
1624 x.update_sat_g();
1625 }
1626 const Bit_Matrix& x_sat = x.sat_g;
1627
1628 Bit_Row all_ones;
1629 all_ones.set_until(x_gs_num_rows);
1630 Bit_Row row_union;
1631 for (dimension_type i = x_cs.num_rows(); i-- > 0; ) {
1632 const bool included
1633 = y.relation_with(x_cs[i]).implies(Poly_Con_Relation::is_included());
1634 if (!included) {
1635 row_union.union_assign(x_gs_red_in_y, x_sat[i]);
1636 if (row_union != all_ones) {
1637 return false;
1638 }
1639 }
1640 }
1641
1642 // Here we know that the upper bound is exact: compute it.
1643 for (dimension_type j = y_gs_num_rows; j-- > 0; ) {
1644 if (!y_gs_red_in_x[j]) {
1645 add_generator(y_gs[j]);
1646 }
1647 }
1648
1649 PPL_ASSERT_HEAVY(OK());
1650 return true;
1651 }
1652
1653 bool
BHZ09_NNC_poly_hull_assign_if_exact(const Polyhedron & y)1654 PPL::Polyhedron::BHZ09_NNC_poly_hull_assign_if_exact(const Polyhedron& y) {
1655 const Polyhedron& x = *this;
1656 // Private method: the caller must ensure the following.
1657 PPL_ASSERT(!x.is_necessarily_closed() && !y.is_necessarily_closed());
1658 PPL_ASSERT(x.space_dim > 0 && x.space_dim == y.space_dim);
1659 PPL_ASSERT(!x.is_empty() && !y.is_empty());
1660
1661 // Minimization is not really required, but it is probably the best
1662 // way of getting constraints, generators and saturation matrices
1663 // up-to-date. Minimization also removes redundant
1664 // constraints/generators.
1665 (void) x.minimize();
1666 (void) y.minimize();
1667
1668 const Generator_System& x_gs = x.gen_sys;
1669 const Generator_System& y_gs = y.gen_sys;
1670 const dimension_type x_gs_num_rows = x_gs.num_rows();
1671 const dimension_type y_gs_num_rows = y_gs.num_rows();
1672
1673 // Compute generators of `x' that are non-redundant in `y' ...
1674 Bit_Row x_gs_non_redundant_in_y;
1675 Bit_Row x_points_non_redundant_in_y;
1676 Bit_Row x_closure_points;
1677 dimension_type num_x_gs_non_redundant_in_y = 0;
1678 for (dimension_type i = x_gs_num_rows; i-- > 0; ) {
1679 const Generator& x_gs_i = x_gs[i];
1680 if (x_gs_i.is_closure_point()) {
1681 x_closure_points.set(i);
1682 }
1683 if (y.relation_with(x_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
1684 continue;
1685 }
1686 x_gs_non_redundant_in_y.set(i);
1687 ++num_x_gs_non_redundant_in_y;
1688 if (x_gs_i.is_point()) {
1689 x_points_non_redundant_in_y.set(i);
1690 }
1691 }
1692
1693 // If `x' is included into `y', the upper bound `y' is exact.
1694 if (num_x_gs_non_redundant_in_y == 0) {
1695 *this = y;
1696 return true;
1697 }
1698
1699 // ... and vice versa, generators of `y' that are non-redundant in `x'.
1700 Bit_Row y_gs_non_redundant_in_x;
1701 Bit_Row y_points_non_redundant_in_x;
1702 Bit_Row y_closure_points;
1703 dimension_type num_y_gs_non_redundant_in_x = 0;
1704 for (dimension_type i = y_gs_num_rows; i-- > 0; ) {
1705 const Generator& y_gs_i = y_gs[i];
1706 if (y_gs_i.is_closure_point()) {
1707 y_closure_points.set(i);
1708 }
1709 if (x.relation_with(y_gs_i).implies(Poly_Gen_Relation::subsumes())) {
1710 continue;
1711 }
1712 y_gs_non_redundant_in_x.set(i);
1713 ++num_y_gs_non_redundant_in_x;
1714 if (y_gs_i.is_point()) {
1715 y_points_non_redundant_in_x.set(i);
1716 }
1717 }
1718
1719 // If `y' is included into `x', the upper bound `x' is exact.
1720 if (num_y_gs_non_redundant_in_x == 0) {
1721 return true;
1722 }
1723
1724 Bit_Row x_non_points_non_redundant_in_y;
1725 x_non_points_non_redundant_in_y
1726 .difference_assign(x_gs_non_redundant_in_y,
1727 x_points_non_redundant_in_y);
1728
1729 const Constraint_System& x_cs = x.con_sys;
1730 const Constraint_System& y_cs = y.con_sys;
1731 const dimension_type x_cs_num_rows = x_cs.num_rows();
1732 const dimension_type y_cs_num_rows = y_cs.num_rows();
1733
1734 // Filter away the points of `x_gs' that would be redundant
1735 // in the topological closure of `y'.
1736 Bit_Row x_points_non_redundant_in_y_closure;
1737 for (dimension_type i = x_points_non_redundant_in_y.first();
1738 i != C_Integer<unsigned long>::max;
1739 i = x_points_non_redundant_in_y.next(i)) {
1740 const Generator& x_p = x_gs[i];
1741 PPL_ASSERT(x_p.is_point());
1742 // NOTE: we cannot use Constraint_System::relation_with()
1743 // as we need to treat strict inequalities as if they were nonstrict.
1744 for (dimension_type j = y_cs_num_rows; j-- > 0; ) {
1745 const Constraint& y_c = y_cs[j];
1746 const int sp_sign = Scalar_Products::reduced_sign(y_c, x_p);
1747 if (sp_sign < 0 || (y_c.is_equality() && sp_sign > 0)) {
1748 x_points_non_redundant_in_y_closure.set(i);
1749 break;
1750 }
1751 }
1752 }
1753
1754 // Make sure the saturation matrix for `x' is up to date.
1755 // Any sat matrix would do: we choose `sat_g' because it matches
1756 // the two nested loops (constraints on rows and generators on columns).
1757 if (!x.sat_g_is_up_to_date()) {
1758 x.update_sat_g();
1759 }
1760 const Bit_Matrix& x_sat = x.sat_g;
1761
1762 Bit_Row x_cs_condition_3;
1763 Bit_Row x_gs_condition_3;
1764 Bit_Row all_ones;
1765 all_ones.set_until(x_gs_num_rows);
1766 Bit_Row saturators;
1767 Bit_Row tmp_set;
1768 for (dimension_type i = x_cs_num_rows; i-- > 0; ) {
1769 const Constraint& x_c = x_cs[i];
1770 // Skip constraint if it is not violated by `y'.
1771 if (y.relation_with(x_c).implies(Poly_Con_Relation::is_included())) {
1772 continue;
1773 }
1774 saturators.difference_assign(all_ones, x_sat[i]);
1775 // Check condition 1.
1776 tmp_set.intersection_assign(x_non_points_non_redundant_in_y, saturators);
1777 if (!tmp_set.empty()) {
1778 return false;
1779 }
1780 if (x_c.is_strict_inequality()) {
1781 // Postpone check for condition 3.
1782 x_cs_condition_3.set(i);
1783 tmp_set.intersection_assign(x_closure_points, saturators);
1784 x_gs_condition_3.union_assign(x_gs_condition_3, tmp_set);
1785 }
1786 else {
1787 // Check condition 2.
1788 tmp_set.intersection_assign(x_points_non_redundant_in_y_closure,
1789 saturators);
1790 if (!tmp_set.empty()) {
1791 return false;
1792 }
1793 }
1794 }
1795
1796 // Now exchange the roles of `x' and `y'
1797 // (the statement of the NNC theorem in BHZ09 is symmetric).
1798
1799 Bit_Row y_non_points_non_redundant_in_x;
1800 y_non_points_non_redundant_in_x
1801 .difference_assign(y_gs_non_redundant_in_x,
1802 y_points_non_redundant_in_x);
1803
1804 // Filter away the points of `y_gs' that would be redundant
1805 // in the topological closure of `x'.
1806 Bit_Row y_points_non_redundant_in_x_closure;
1807 for (dimension_type i = y_points_non_redundant_in_x.first();
1808 i != C_Integer<unsigned long>::max;
1809 i = y_points_non_redundant_in_x.next(i)) {
1810 const Generator& y_p = y_gs[i];
1811 PPL_ASSERT(y_p.is_point());
1812 // NOTE: we cannot use Constraint_System::relation_with()
1813 // as we need to treat strict inequalities as if they were nonstrict.
1814 for (dimension_type j = x_cs_num_rows; j-- > 0; ) {
1815 const Constraint& x_c = x_cs[j];
1816 const int sp_sign = Scalar_Products::reduced_sign(x_c, y_p);
1817 if (sp_sign < 0 || (x_c.is_equality() && sp_sign > 0)) {
1818 y_points_non_redundant_in_x_closure.set(i);
1819 break;
1820 }
1821 }
1822 }
1823
1824 // Make sure the saturation matrix `sat_g' for `y' is up to date.
1825 if (!y.sat_g_is_up_to_date()) {
1826 y.update_sat_g();
1827 }
1828 const Bit_Matrix& y_sat = y.sat_g;
1829
1830 Bit_Row y_cs_condition_3;
1831 Bit_Row y_gs_condition_3;
1832 all_ones.clear();
1833 all_ones.set_until(y_gs_num_rows);
1834 for (dimension_type i = y_cs_num_rows; i-- > 0; ) {
1835 const Constraint& y_c = y_cs[i];
1836 // Skip constraint if it is not violated by `x'.
1837 if (x.relation_with(y_c).implies(Poly_Con_Relation::is_included())) {
1838 continue;
1839 }
1840 saturators.difference_assign(all_ones, y_sat[i]);
1841 // Check condition 1.
1842 tmp_set.intersection_assign(y_non_points_non_redundant_in_x, saturators);
1843 if (!tmp_set.empty()) {
1844 return false;
1845 }
1846 if (y_c.is_strict_inequality()) {
1847 // Postpone check for condition 3.
1848 y_cs_condition_3.set(i);
1849 tmp_set.intersection_assign(y_closure_points, saturators);
1850 y_gs_condition_3.union_assign(y_gs_condition_3, tmp_set);
1851 }
1852 else {
1853 // Check condition 2.
1854 tmp_set.intersection_assign(y_points_non_redundant_in_x_closure,
1855 saturators);
1856 if (!tmp_set.empty()) {
1857 return false;
1858 }
1859 }
1860 }
1861
1862 // Now considering condition 3.
1863
1864 if (x_cs_condition_3.empty() && y_cs_condition_3.empty()) {
1865 // No test for condition 3 is needed.
1866 // The hull is exact: compute it.
1867 for (dimension_type j = y_gs_num_rows; j-- > 0; ) {
1868 if (y_gs_non_redundant_in_x[j]) {
1869 add_generator(y_gs[j]);
1870 }
1871 }
1872 return true;
1873 }
1874
1875 // We have anyway to compute the upper bound and its constraints too.
1876 Polyhedron ub(x);
1877 for (dimension_type j = y_gs_num_rows; j-- > 0; ) {
1878 if (y_gs_non_redundant_in_x[j]) {
1879 ub.add_generator(y_gs[j]);
1880 }
1881 }
1882 (void) ub.minimize();
1883 PPL_ASSERT(!ub.is_empty());
1884
1885 // NOTE: the following computation of x_gs_condition_3_not_in_y
1886 // (resp., y_gs_condition_3_not_in_x) is not required for correctness.
1887 // It is done so as to later apply a speculative test
1888 // (i.e., a non-conclusive but computationally lighter test).
1889
1890 // Filter away from `x_gs_condition_3' those closure points
1891 // that, when considered as points, would belong to `y',
1892 // i.e., those that violate no strict constraint in `y_cs'.
1893 Bit_Row x_gs_condition_3_not_in_y;
1894 for (dimension_type i = y_cs_num_rows; i-- > 0; ) {
1895 const Constraint& y_c = y_cs[i];
1896 if (y_c.is_strict_inequality()) {
1897 for (dimension_type j = x_gs_condition_3.first();
1898 j != C_Integer<unsigned long>::max; j = x_gs_condition_3.next(j)) {
1899 const Generator& x_cp = x_gs[j];
1900 PPL_ASSERT(x_cp.is_closure_point());
1901 const int sp_sign = Scalar_Products::reduced_sign(y_c, x_cp);
1902 PPL_ASSERT(sp_sign >= 0);
1903 if (sp_sign == 0) {
1904 x_gs_condition_3.clear(j);
1905 x_gs_condition_3_not_in_y.set(j);
1906 }
1907 }
1908 if (x_gs_condition_3.empty()) {
1909 break;
1910 }
1911 }
1912 }
1913 // Symmetrically, filter away from `y_gs_condition_3' those
1914 // closure points that, when considered as points, would belong to `x',
1915 // i.e., those that violate no strict constraint in `x_cs'.
1916 Bit_Row y_gs_condition_3_not_in_x;
1917 for (dimension_type i = x_cs_num_rows; i-- > 0; ) {
1918 if (x_cs[i].is_strict_inequality()) {
1919 const Constraint& x_c = x_cs[i];
1920 for (dimension_type j = y_gs_condition_3.first();
1921 j != C_Integer<unsigned long>::max; j = y_gs_condition_3.next(j)) {
1922 const Generator& y_cp = y_gs[j];
1923 PPL_ASSERT(y_cp.is_closure_point());
1924 const int sp_sign = Scalar_Products::reduced_sign(x_c, y_cp);
1925 PPL_ASSERT(sp_sign >= 0);
1926 if (sp_sign == 0) {
1927 y_gs_condition_3.clear(j);
1928 y_gs_condition_3_not_in_x.set(j);
1929 }
1930 }
1931 if (y_gs_condition_3.empty()) {
1932 break;
1933 }
1934 }
1935 }
1936
1937 // NOTE: here we apply the speculative test.
1938 // Check if there exists a closure point in `x_gs_condition_3_not_in_y'
1939 // or `y_gs_condition_3_not_in_x' that belongs (as point) to the hull.
1940 // If so, the hull is not exact.
1941 const Constraint_System& ub_cs = ub.constraints();
1942 for (dimension_type i = ub_cs.num_rows(); i-- > 0; ) {
1943 if (ub_cs[i].is_strict_inequality()) {
1944 const Constraint& ub_c = ub_cs[i];
1945 for (dimension_type j = x_gs_condition_3_not_in_y.first();
1946 j != C_Integer<unsigned long>::max;
1947 j = x_gs_condition_3_not_in_y.next(j)) {
1948 const Generator& x_cp = x_gs[j];
1949 PPL_ASSERT(x_cp.is_closure_point());
1950 const int sp_sign = Scalar_Products::reduced_sign(ub_c, x_cp);
1951 PPL_ASSERT(sp_sign >= 0);
1952 if (sp_sign == 0) {
1953 x_gs_condition_3_not_in_y.clear(j);
1954 }
1955 }
1956 for (dimension_type j = y_gs_condition_3_not_in_x.first();
1957 j != C_Integer<unsigned long>::max;
1958 j = y_gs_condition_3_not_in_x.next(j)) {
1959 const Generator& y_cp = y_gs[j];
1960 PPL_ASSERT(y_cp.is_closure_point());
1961 const int sp_sign = Scalar_Products::reduced_sign(ub_c, y_cp);
1962 PPL_ASSERT(sp_sign >= 0);
1963 if (sp_sign == 0) {
1964 y_gs_condition_3_not_in_x.clear(j);
1965 }
1966 }
1967 }
1968 }
1969
1970 if (!(x_gs_condition_3_not_in_y.empty()
1971 && y_gs_condition_3_not_in_x.empty())) {
1972 // There exist a closure point satisfying condition 3,
1973 // hence the hull is not exact.
1974 return false;
1975 }
1976 // The speculative test was not successful:
1977 // apply the expensive (but conclusive) test for condition 3.
1978
1979 // Consider strict inequalities in `x' violated by `y'.
1980 for (dimension_type i = x_cs_condition_3.first();
1981 i != C_Integer<unsigned long>::max; i = x_cs_condition_3.next(i)) {
1982 const Constraint& x_cs_i = x_cs[i];
1983 PPL_ASSERT(x_cs_i.is_strict_inequality());
1984 // Build the equality constraint induced by x_cs_i.
1985 const Linear_Expression expr(x_cs_i.expression());
1986 const Constraint eq_i(expr == 0);
1987 PPL_ASSERT(!(ub.relation_with(eq_i)
1988 .implies(Poly_Con_Relation::is_disjoint())));
1989 Polyhedron ub_inters_hyperplane(ub);
1990 ub_inters_hyperplane.add_constraint(eq_i);
1991 Polyhedron y_inters_hyperplane(y);
1992 y_inters_hyperplane.add_constraint(eq_i);
1993 if (!y_inters_hyperplane.contains(ub_inters_hyperplane)) {
1994 // The hull is not exact.
1995 return false;
1996 }
1997 }
1998
1999 // Consider strict inequalities in `y' violated by `x'.
2000 for (dimension_type i = y_cs_condition_3.first();
2001 i != C_Integer<unsigned long>::max; i = y_cs_condition_3.next(i)) {
2002 const Constraint& y_cs_i = y_cs[i];
2003 PPL_ASSERT(y_cs_i.is_strict_inequality());
2004 // Build the equality constraint induced by y_cs_i.
2005 const Constraint eq_i(Linear_Expression(y_cs_i.expression()) == 0);
2006 PPL_ASSERT(!(ub.relation_with(eq_i)
2007 .implies(Poly_Con_Relation::is_disjoint())));
2008 Polyhedron ub_inters_hyperplane(ub);
2009 ub_inters_hyperplane.add_constraint(eq_i);
2010 Polyhedron x_inters_hyperplane(x);
2011 x_inters_hyperplane.add_constraint(eq_i);
2012 if (!x_inters_hyperplane.contains(ub_inters_hyperplane)) {
2013 // The hull is not exact.
2014 return false;
2015 }
2016 }
2017
2018 // The hull is exact.
2019 m_swap(ub);
2020 return true;
2021 }
2022
2023 bool
BFT00_poly_hull_assign_if_exact(const Polyhedron & y)2024 PPL::Polyhedron::BFT00_poly_hull_assign_if_exact(const Polyhedron& y) {
2025 // Declare a const reference to *this (to avoid accidental modifications).
2026 const Polyhedron& x = *this;
2027 // Private method: the caller must ensure the following.
2028 PPL_ASSERT(x.is_necessarily_closed());
2029 PPL_ASSERT(x.topology() == y.topology());
2030 PPL_ASSERT(x.space_dim == y.space_dim);
2031
2032 // The zero-dim case is trivial.
2033 if (x.space_dim == 0) {
2034 upper_bound_assign(y);
2035 return true;
2036 }
2037 // If `x' or `y' is (known to be) empty, the convex union is exact.
2038 if (x.marked_empty()) {
2039 *this = y;
2040 return true;
2041 }
2042 else if (y.is_empty()) {
2043 return true;
2044 }
2045 else if (x.is_empty()) {
2046 *this = y;
2047 return true;
2048 }
2049
2050 // Here both `x' and `y' are known to be non-empty.
2051
2052 // Implementation based on Algorithm 8.1 (page 15) in [BemporadFT00TR],
2053 // generalized so as to also allow for unbounded polyhedra.
2054 // The extension to unbounded polyhedra is obtained by mimicking
2055 // what done in Algorithm 8.2 (page 19) with respect to
2056 // Algorithm 6.2 (page 13).
2057 // We also apply a couple of improvements (see steps 2.1, 3.1, 6.1, 7.1)
2058 // so as to quickly handle special cases and avoid the splitting
2059 // of equalities/lines into pairs of inequalities/rays.
2060
2061 (void) x.minimize();
2062 (void) y.minimize();
2063 const Constraint_System& x_cs = x.con_sys;
2064 const Constraint_System& y_cs = y.con_sys;
2065 const Generator_System& x_gs = x.gen_sys;
2066 const Generator_System& y_gs = y.gen_sys;
2067 const dimension_type x_gs_num_rows = x_gs.num_rows();
2068 const dimension_type y_gs_num_rows = y_gs.num_rows();
2069
2070 // Step 1: generators of `x' that are redundant in `y', and vice versa.
2071 std::vector<bool> x_gs_red_in_y(x_gs_num_rows, false);
2072 dimension_type num_x_gs_red_in_y = 0;
2073 for (dimension_type i = x_gs_num_rows; i-- > 0; ) {
2074 if (y.relation_with(x_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
2075 x_gs_red_in_y[i] = true;
2076 ++num_x_gs_red_in_y;
2077 }
2078 }
2079 std::vector<bool> y_gs_red_in_x(y_gs_num_rows, false);
2080 dimension_type num_y_gs_red_in_x = 0;
2081 for (dimension_type i = y_gs_num_rows; i-- > 0; ) {
2082 if (x.relation_with(y_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
2083 y_gs_red_in_x[i] = true;
2084 ++num_y_gs_red_in_x;
2085 }
2086 }
2087
2088 // Step 2: if no redundant generator has been identified,
2089 // then the union is not convex. CHECKME: why?
2090 if (num_x_gs_red_in_y == 0 && num_y_gs_red_in_x == 0) {
2091 return false;
2092 }
2093
2094 // Step 2.1: while at it, also perform quick inclusion tests.
2095 if (num_y_gs_red_in_x == y_gs_num_rows) {
2096 // `y' is included into `x': union is convex.
2097 return true;
2098 }
2099 if (num_x_gs_red_in_y == x_gs_num_rows) {
2100 // `x' is included into `y': union is convex.
2101 *this = y;
2102 return true;
2103 }
2104
2105 // Here we know that `x' is not included in `y', and vice versa.
2106
2107 // Step 3: constraints of `x' that are satisfied by `y', and vice versa.
2108 const dimension_type x_cs_num_rows = x_cs.num_rows();
2109 std::vector<bool> x_cs_red_in_y(x_cs_num_rows, false);
2110 for (dimension_type i = x_cs_num_rows; i-- > 0; ) {
2111 const Constraint& x_cs_i = x_cs[i];
2112 if (y.relation_with(x_cs_i).implies(Poly_Con_Relation::is_included())) {
2113 x_cs_red_in_y[i] = true;
2114 }
2115 else if (x_cs_i.is_equality()) {
2116 // Step 3.1: `x' has an equality not satisfied by `y':
2117 // union is not convex (recall that `y' does not contain `x').
2118 // NOTE: this would be false for NNC polyhedra.
2119 // Example: x = { A == 0 }, y = { 0 < A <= 1 }.
2120 return false;
2121 }
2122 }
2123 const dimension_type y_cs_num_rows = y_cs.num_rows();
2124 std::vector<bool> y_cs_red_in_x(y_cs_num_rows, false);
2125 for (dimension_type i = y_cs_num_rows; i-- > 0; ) {
2126 const Constraint& y_cs_i = y_cs[i];
2127 if (x.relation_with(y_cs_i).implies(Poly_Con_Relation::is_included())) {
2128 y_cs_red_in_x[i] = true;
2129 }
2130 else if (y_cs_i.is_equality()) {
2131 // Step 3.1: `y' has an equality not satisfied by `x':
2132 // union is not convex (see explanation above).
2133 return false;
2134 }
2135 }
2136
2137 // Loop in steps 5-9: for each pair of non-redundant generators,
2138 // compute their "mid-point" and check if it is both in `x' and `y'.
2139
2140 // Note: reasoning at the polyhedral cone level.
2141 Generator mid_g;
2142
2143 for (dimension_type i = x_gs_num_rows; i-- > 0; ) {
2144 if (x_gs_red_in_y[i]) {
2145 continue;
2146 }
2147 const Generator& x_g = x_gs[i];
2148 const bool x_g_is_line = x_g.is_line_or_equality();
2149 for (dimension_type j = y_gs_num_rows; j-- > 0; ) {
2150 if (y_gs_red_in_x[j]) {
2151 continue;
2152 }
2153 const Generator& y_g = y_gs[j];
2154 const bool y_g_is_line = y_g.is_line_or_equality();
2155
2156 // Step 6: compute mid_g = x_g + y_g.
2157 // NOTE: no need to actually compute the "mid-point",
2158 // since any strictly positive combination would do.
2159 mid_g = x_g;
2160 mid_g.expr += y_g.expr;
2161 // A zero ray is not a well formed generator.
2162 const bool illegal_ray
2163 = (mid_g.expr.inhomogeneous_term() == 0 && mid_g.expr.all_homogeneous_terms_are_zero());
2164 // A zero ray cannot be generated from a line: this holds
2165 // because x_row (resp., y_row) is not subsumed by y (resp., x).
2166 PPL_ASSERT(!(illegal_ray && (x_g_is_line || y_g_is_line)));
2167 if (illegal_ray) {
2168 continue;
2169 }
2170 mid_g.expr.normalize();
2171 if (x_g_is_line) {
2172 if (y_g_is_line) {
2173 // mid_row is a line too: sign normalization is needed.
2174 mid_g.sign_normalize();
2175 }
2176 else {
2177 // mid_row is a ray/point.
2178 mid_g.set_is_ray_or_point_or_inequality();
2179 }
2180 }
2181 PPL_ASSERT(mid_g.OK());
2182
2183 // Step 7: check if mid_g is in the union of x and y.
2184 if (x.relation_with(mid_g) == Poly_Gen_Relation::nothing()
2185 && y.relation_with(mid_g) == Poly_Gen_Relation::nothing()) {
2186 return false;
2187 }
2188 // If either x_g or y_g is a line, we should use its
2189 // negation to produce another generator to be tested too.
2190 // NOTE: exclusive-or is meant.
2191 if (!x_g_is_line && y_g_is_line) {
2192 // Step 6.1: (re-)compute mid_row = x_g - y_g.
2193 mid_g = x_g;
2194 mid_g.expr -= y_g.expr;
2195 mid_g.expr.normalize();
2196 PPL_ASSERT(mid_g.OK());
2197 // Step 7.1: check if mid_g is in the union of x and y.
2198 if (x.relation_with(mid_g) == Poly_Gen_Relation::nothing()
2199 && y.relation_with(mid_g) == Poly_Gen_Relation::nothing()) {
2200 return false;
2201 }
2202 }
2203 else if (x_g_is_line && !y_g_is_line) {
2204 // Step 6.1: (re-)compute mid_row = - x_row + y_row.
2205 mid_g = y_g;
2206 mid_g.expr -= x_g.expr;
2207 mid_g.expr.normalize();
2208 PPL_ASSERT(mid_g.OK());
2209 // Step 7.1: check if mid_g is in the union of x and y.
2210 if (x.relation_with(mid_g) == Poly_Gen_Relation::nothing()
2211 && y.relation_with(mid_g) == Poly_Gen_Relation::nothing()) {
2212 return false;
2213 }
2214 }
2215 }
2216 }
2217
2218 // Here we know that the union of x and y is convex.
2219 // TODO: exploit knowledge on the cardinality of non-redundant
2220 // constraints/generators to improve the convex-hull computation.
2221 // Using generators allows for exploiting incrementality.
2222 for (dimension_type j = 0; j < y_gs_num_rows; ++j) {
2223 if (!y_gs_red_in_x[j]) {
2224 add_generator(y_gs[j]);
2225 }
2226 }
2227 PPL_ASSERT_HEAVY(OK());
2228 return true;
2229 }
2230
2231 void
drop_some_non_integer_points(const Variables_Set * vars_p,Complexity_Class complexity)2232 PPL::Polyhedron::drop_some_non_integer_points(const Variables_Set* vars_p,
2233 Complexity_Class complexity) {
2234 // There is nothing to do for an empty set of variables.
2235 if (vars_p != 0 && vars_p->empty()) {
2236 return;
2237 }
2238
2239 // Any empty polyhedron does not contain integer points.
2240 if (marked_empty()) {
2241 return;
2242 }
2243
2244 // A zero-dimensional, universe polyhedron has, by convention, an
2245 // integer point.
2246 if (space_dim == 0) {
2247 set_empty();
2248 return;
2249 }
2250
2251 // The constraints (possibly with pending rows) are required.
2252 if (has_pending_generators()) {
2253 // Processing of pending generators is exponential in the worst case.
2254 if (complexity != ANY_COMPLEXITY) {
2255 return;
2256 }
2257 else {
2258 process_pending_generators();
2259 }
2260 }
2261 if (!constraints_are_up_to_date()) {
2262 // Constraints update is exponential in the worst case.
2263 if (complexity != ANY_COMPLEXITY) {
2264 return;
2265 }
2266 else {
2267 update_constraints();
2268 }
2269 }
2270 // For NNC polyhedra we need to process any pending constraints.
2271 if (!is_necessarily_closed() && has_pending_constraints()) {
2272 if (complexity != ANY_COMPLEXITY) {
2273 return;
2274 }
2275 else if (!process_pending_constraints()) {
2276 // We just discovered the polyhedron is empty.
2277 return;
2278 }
2279 }
2280
2281 PPL_ASSERT(!has_pending_generators() && constraints_are_up_to_date());
2282 PPL_ASSERT(is_necessarily_closed() || !has_pending_constraints());
2283
2284 bool changed = false;
2285 PPL_DIRTY_TEMP_COEFFICIENT(gcd);
2286
2287 const bool con_sys_was_sorted = con_sys.is_sorted();
2288
2289 Variables_Set other_vars;
2290 if (vars_p != 0) {
2291 // Compute the complement of `*vars_p'.
2292 for (dimension_type i = 0; i < space_dim; ++i) {
2293 if (vars_p->find(i) == vars_p->end()) {
2294 other_vars.insert(Variable(i));
2295 }
2296 }
2297 }
2298
2299 for (dimension_type j = con_sys.sys.rows.size(); j-- > 0; ) {
2300 Constraint& c = con_sys.sys.rows[j];
2301 if (c.is_tautological()) {
2302 continue;
2303 }
2304 if (!other_vars.empty()) {
2305 // Skip constraints having a nonzero coefficient for a variable
2306 // that does not occurr in the input set.
2307 if (!c.expression().all_zeroes(other_vars)) {
2308 continue;
2309 }
2310 }
2311
2312 if (!is_necessarily_closed()) {
2313 // Transform all strict inequalities into non-strict ones,
2314 // with the inhomogeneous term incremented by 1.
2315 if (c.epsilon_coefficient() < 0) {
2316 c.set_epsilon_coefficient(0);
2317 Linear_Expression& e = c.expr;
2318 e.set_inhomogeneous_term(e.inhomogeneous_term() - 1);
2319 // Enforce normalization.
2320 // FIXME: is this really necessary?
2321 e.normalize();
2322 PPL_ASSERT(c.OK());
2323 changed = true;
2324 }
2325 }
2326
2327 // Compute the GCD of all the homogeneous terms.
2328 gcd = c.expression().gcd(1, space_dim + 1);
2329
2330 if (gcd != 0 && gcd != 1) {
2331 PPL_ASSERT(c.expr.inhomogeneous_term() % gcd != 0);
2332
2333 // If we have an equality, the polyhedron becomes empty.
2334 if (c.is_equality()) {
2335 set_empty();
2336 return;
2337 }
2338
2339 // Divide the inhomogeneous coefficients by the GCD.
2340 c.expr.exact_div_assign(gcd, 1, space_dim + 1);
2341
2342 PPL_DIRTY_TEMP_COEFFICIENT(c_0);
2343 c_0 = c.expr.inhomogeneous_term();
2344 const int c_0_sign = sgn(c_0);
2345 c_0 /= gcd;
2346 if (c_0_sign < 0) {
2347 --c_0;
2348 }
2349 c.expr.set_inhomogeneous_term(c_0);
2350 PPL_ASSERT(c.OK());
2351 changed = true;
2352 }
2353 }
2354
2355 con_sys.set_sorted(!changed && con_sys_was_sorted);
2356 PPL_ASSERT(con_sys.sys.OK());
2357
2358 if (changed) {
2359 if (is_necessarily_closed()) {
2360 con_sys.insert(Constraint::zero_dim_positivity());
2361 }
2362 else {
2363 con_sys.insert(Constraint::epsilon_leq_one());
2364 }
2365 // After changing the system of constraints, the generators
2366 // are no longer up-to-date and the constraints are no longer
2367 // minimized.
2368 clear_generators_up_to_date();
2369 clear_constraints_minimized();
2370 }
2371 PPL_ASSERT_HEAVY(OK());
2372 }
2373
2374 void
positive_time_elapse_assign_impl(const Polyhedron & y)2375 PPL::Polyhedron::positive_time_elapse_assign_impl(const Polyhedron& y) {
2376 // Private method: the caller must ensure the following.
2377 PPL_ASSERT(!is_necessarily_closed());
2378
2379 Polyhedron& x = *this;
2380 // Dimension-compatibility checks.
2381 if (x.space_dim != y.space_dim) {
2382 throw_dimension_incompatible("positive_time_elapse_assign(y)", "y", y);
2383 }
2384
2385 // Dealing with the zero-dimensional case.
2386 if (x.space_dim == 0) {
2387 if (y.marked_empty()) {
2388 x.set_empty();
2389 }
2390 return;
2391 }
2392
2393 // If either one of `x' or `y' is empty, the result is empty too.
2394 if (x.marked_empty() || y.marked_empty()
2395 || (x.has_pending_constraints() && !x.process_pending_constraints())
2396 || (!x.generators_are_up_to_date() && !x.update_generators())
2397 || (y.has_pending_constraints() && !y.process_pending_constraints())
2398 || (!y.generators_are_up_to_date() && !y.update_generators())) {
2399 x.set_empty();
2400 return;
2401 }
2402
2403 // At this point both generator systems are up-to-date,
2404 // possibly containing pending generators.
2405
2406 // The new system of generators that will replace the one of x.
2407 Generator_System new_gs(x.gen_sys);
2408 dimension_type num_rows = new_gs.num_rows();
2409
2410 // We are going to do all sorts of funny things with new_gs, so we better
2411 // mark it unsorted.
2412 // Note: `sorted' is an attribute of Linear_System, encapsulated by
2413 // Generator_System; hence, the following is equivalent to
2414 // new_gs.set_sorted(false).
2415 new_gs.sys.set_sorted(false);
2416
2417 // Remove all points from new_gs and put them in 'x_points_gs' for later use.
2418 // Notice that we do not remove the corresponding closure points.
2419 Generator_System x_points_gs;
2420 for (dimension_type i = num_rows; i-- > 0;) {
2421 Generator &g = new_gs.sys.rows[i];
2422 if (g.is_point()) {
2423 x_points_gs.insert(g);
2424 --num_rows;
2425 swap(g, new_gs.sys.rows[num_rows]);
2426 }
2427 }
2428 new_gs.sys.rows.resize(num_rows);
2429
2430 // When there are no pending rows, the pending row index points at
2431 // the smallest non-valid row, i.e., it is equal to the number of rows.
2432 // Hence, each time the system is manually resized, the pending row index
2433 // must be updated.
2434 new_gs.unset_pending_rows();
2435 PPL_ASSERT(new_gs.sys.OK());
2436
2437 // We use a pointer in order to avoid copying the generator
2438 // system when it is not necessary, i.e., when y is an NNC.
2439 const Generator_System *gs = &y.gen_sys;
2440 Generator_System y_gs;
2441
2442 // If y is closed, promote its generator system to not necessarily closed.
2443 if (y.is_necessarily_closed()) {
2444 y_gs = y.gen_sys;
2445 y_gs.convert_into_non_necessarily_closed();
2446 y_gs.add_corresponding_closure_points();
2447 gs = &y_gs;
2448 }
2449
2450 PPL_ASSERT(gs->OK());
2451
2452 const dimension_type gs_num_rows = gs->num_rows();
2453
2454 // For each generator g of y...
2455 for (dimension_type i = gs_num_rows; i-- > 0; ) {
2456 const Generator &g = gs->sys.rows[i];
2457
2458 switch (g.type()) {
2459 case Generator::POINT:
2460 // In principle, points should be added to new_gs as rays.
2461 // However, for each point there is a corresponding closure point in "gs".
2462 // Hence, we leave this operation to closure points.
2463
2464 // Insert into new_gs the sum of g and each point of x.
2465 // For each original point gx of x...
2466 for (dimension_type j = x_points_gs.sys.num_rows(); j-- > 0; ) {
2467 const Generator &gx = x_points_gs.sys.rows[j];
2468 PPL_ASSERT(gx.is_point());
2469 // ...insert the point obtained as the sum of g and gx.
2470 Generator new_g = g; // make a copy
2471 Coefficient new_divisor = g.expr.inhomogeneous_term() * gx.expr.inhomogeneous_term();
2472
2473 new_g.expr.linear_combine(gx.expr, gx.expr.inhomogeneous_term(), g.expr.inhomogeneous_term());
2474 new_g.expr.set_inhomogeneous_term(new_divisor);
2475 if (new_g.is_not_necessarily_closed()) {
2476 new_g.set_epsilon_coefficient(g.epsilon_coefficient());
2477 }
2478 new_g.expr.normalize();
2479 PPL_ASSERT(new_g.OK());
2480 new_gs.insert(new_g);
2481 }
2482 break;
2483 case Generator::CLOSURE_POINT:
2484 // If g is not the origin, insert g into new_gs, as a ray.
2485 if (!g.expr.all_homogeneous_terms_are_zero()) {
2486 // Turn a copy of g into a ray.
2487 Generator g_as_a_ray = g;
2488 g_as_a_ray.expr.set_inhomogeneous_term(0);
2489 g_as_a_ray.expr.normalize();
2490 PPL_ASSERT(g_as_a_ray.OK());
2491 // Insert the ray.
2492 new_gs.insert(g_as_a_ray);
2493 }
2494 break;
2495 case Generator::RAY:
2496 case Generator::LINE:
2497 // Insert g into new_gs.
2498 new_gs.insert(g);
2499 break;
2500 }
2501 }
2502 new_gs.add_corresponding_closure_points();
2503 // Notice: add_corresponding_closure_points adds them as pending.
2504 new_gs.unset_pending_rows();
2505
2506 //Polyhedron new_x(...,new_gs);
2507 //swap(x,new_x);
2508
2509 PPL_ASSERT(new_gs.sys.OK());
2510
2511 // Stealing the rows from `new_gs'.
2512 using std::swap;
2513 swap(gen_sys, new_gs);
2514
2515 gen_sys.set_sorted(false);
2516 clear_generators_minimized();
2517 // Generators are now up-to-date.
2518 set_generators_up_to_date();
2519 // Constraints are not up-to-date.
2520 clear_constraints_up_to_date();
2521
2522 PPL_ASSERT_HEAVY(x.OK(true) && y.OK(true));
2523 }
2524
2525 void
throw_invalid_argument(const char * method,const char * reason) const2526 PPL::Polyhedron::throw_invalid_argument(const char* method,
2527 const char* reason) const {
2528 std::ostringstream s;
2529 s << "PPL::";
2530 if (is_necessarily_closed()) {
2531 s << "C_";
2532 }
2533 else {
2534 s << "NNC_";
2535 }
2536 s << "Polyhedron::" << method << ":" << std::endl
2537 << reason << ".";
2538 throw std::invalid_argument(s.str());
2539 }
2540
2541 void
throw_topology_incompatible(const char * method,const char * ph_name,const Polyhedron & ph) const2542 PPL::Polyhedron::throw_topology_incompatible(const char* method,
2543 const char* ph_name,
2544 const Polyhedron& ph) const {
2545 std::ostringstream s;
2546 s << "PPL::";
2547 if (is_necessarily_closed()) {
2548 s << "C_";
2549 }
2550 else {
2551 s << "NNC_";
2552 }
2553 s << "Polyhedron::" << method << ":" << std::endl
2554 << ph_name << " is a ";
2555 if (ph.is_necessarily_closed()) {
2556 s << "C_";
2557 }
2558 else {
2559 s << "NNC_";
2560 }
2561 s << "Polyhedron." << std::endl;
2562 throw std::invalid_argument(s.str());
2563 }
2564
2565 void
throw_topology_incompatible(const char * method,const char * c_name,const Constraint &) const2566 PPL::Polyhedron::throw_topology_incompatible(const char* method,
2567 const char* c_name,
2568 const Constraint&) const {
2569 PPL_ASSERT(is_necessarily_closed());
2570 std::ostringstream s;
2571 s << "PPL::C_Polyhedron::" << method << ":" << std::endl
2572 << c_name << " is a strict inequality.";
2573 throw std::invalid_argument(s.str());
2574 }
2575
2576 void
throw_topology_incompatible(const char * method,const char * g_name,const Generator &) const2577 PPL::Polyhedron::throw_topology_incompatible(const char* method,
2578 const char* g_name,
2579 const Generator&) const {
2580 PPL_ASSERT(is_necessarily_closed());
2581 std::ostringstream s;
2582 s << "PPL::C_Polyhedron::" << method << ":" << std::endl
2583 << g_name << " is a closure point.";
2584 throw std::invalid_argument(s.str());
2585 }
2586
2587 void
throw_topology_incompatible(const char * method,const char * cs_name,const Constraint_System &) const2588 PPL::Polyhedron::throw_topology_incompatible(const char* method,
2589 const char* cs_name,
2590 const Constraint_System&) const {
2591 PPL_ASSERT(is_necessarily_closed());
2592 std::ostringstream s;
2593 s << "PPL::C_Polyhedron::" << method << ":" << std::endl
2594 << cs_name << " contains strict inequalities.";
2595 throw std::invalid_argument(s.str());
2596 }
2597
2598 void
throw_topology_incompatible(const char * method,const char * gs_name,const Generator_System &) const2599 PPL::Polyhedron::throw_topology_incompatible(const char* method,
2600 const char* gs_name,
2601 const Generator_System&) const {
2602 std::ostringstream s;
2603 s << "PPL::C_Polyhedron::" << method << ":" << std::endl
2604 << gs_name << " contains closure points.";
2605 throw std::invalid_argument(s.str());
2606 }
2607
2608 void
throw_dimension_incompatible(const char * method,const char * other_name,dimension_type other_dim) const2609 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
2610 const char* other_name,
2611 dimension_type other_dim) const {
2612 std::ostringstream s;
2613 s << "PPL::"
2614 << (is_necessarily_closed() ? "C_" : "NNC_")
2615 << "Polyhedron::" << method << ":\n"
2616 << "this->space_dimension() == " << space_dimension() << ", "
2617 << other_name << ".space_dimension() == " << other_dim << ".";
2618 throw std::invalid_argument(s.str());
2619 }
2620
2621 void
throw_dimension_incompatible(const char * method,const char * ph_name,const Polyhedron & ph) const2622 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
2623 const char* ph_name,
2624 const Polyhedron& ph) const {
2625 throw_dimension_incompatible(method, ph_name, ph.space_dimension());
2626 }
2627
2628 void
2629 PPL::Polyhedron
throw_dimension_incompatible(const char * method,const char * le_name,const Linear_Expression & le) const2630 ::throw_dimension_incompatible(const char* method,
2631 const char* le_name,
2632 const Linear_Expression& le) const {
2633 throw_dimension_incompatible(method, le_name, le.space_dimension());
2634 }
2635
2636 void
throw_dimension_incompatible(const char * method,const char * c_name,const Constraint & c) const2637 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
2638 const char* c_name,
2639 const Constraint& c) const {
2640 throw_dimension_incompatible(method, c_name, c.space_dimension());
2641 }
2642
2643 void
throw_dimension_incompatible(const char * method,const char * g_name,const Generator & g) const2644 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
2645 const char* g_name,
2646 const Generator& g) const {
2647 throw_dimension_incompatible(method, g_name, g.space_dimension());
2648 }
2649
2650 void
throw_dimension_incompatible(const char * method,const char * cg_name,const Congruence & cg) const2651 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
2652 const char* cg_name,
2653 const Congruence& cg) const {
2654 throw_dimension_incompatible(method, cg_name, cg.space_dimension());
2655 }
2656
2657 void
throw_dimension_incompatible(const char * method,const char * cs_name,const Constraint_System & cs) const2658 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
2659 const char* cs_name,
2660 const Constraint_System& cs) const {
2661 throw_dimension_incompatible(method, cs_name, cs.space_dimension());
2662 }
2663
2664 void
2665 PPL::Polyhedron
throw_dimension_incompatible(const char * method,const char * gs_name,const Generator_System & gs) const2666 ::throw_dimension_incompatible(const char* method,
2667 const char* gs_name,
2668 const Generator_System& gs) const {
2669 throw_dimension_incompatible(method, gs_name, gs.space_dimension());
2670 }
2671
2672 void
2673 PPL::Polyhedron
throw_dimension_incompatible(const char * method,const char * cgs_name,const Congruence_System & cgs) const2674 ::throw_dimension_incompatible(const char* method,
2675 const char* cgs_name,
2676 const Congruence_System& cgs) const {
2677 throw_dimension_incompatible(method, cgs_name, cgs.space_dimension());
2678 }
2679
2680 void
throw_dimension_incompatible(const char * method,const char * var_name,const Variable var) const2681 PPL::Polyhedron::throw_dimension_incompatible(const char* method,
2682 const char* var_name,
2683 const Variable var) const {
2684 std::ostringstream s;
2685 s << "PPL::";
2686 if (is_necessarily_closed()) {
2687 s << "C_";
2688 }
2689 else {
2690 s << "NNC_";
2691 }
2692 s << "Polyhedron::" << method << ":" << std::endl
2693 << "this->space_dimension() == " << space_dimension() << ", "
2694 << var_name << ".space_dimension() == " << var.space_dimension() << ".";
2695 throw std::invalid_argument(s.str());
2696 }
2697
2698 void
2699 PPL::Polyhedron::
throw_dimension_incompatible(const char * method,dimension_type required_space_dim) const2700 throw_dimension_incompatible(const char* method,
2701 dimension_type required_space_dim) const {
2702 std::ostringstream s;
2703 s << "PPL::";
2704 if (is_necessarily_closed()) {
2705 s << "C_";
2706 }
2707 else {
2708 s << "NNC_";
2709 }
2710 s << "Polyhedron::" << method << ":" << std::endl
2711 << "this->space_dimension() == " << space_dimension()
2712 << ", required space dimension == " << required_space_dim << ".";
2713 throw std::invalid_argument(s.str());
2714 }
2715
2716 PPL::dimension_type
check_space_dimension_overflow(const dimension_type dim,const dimension_type max,const Topology topol,const char * method,const char * reason)2717 PPL::Polyhedron::check_space_dimension_overflow(const dimension_type dim,
2718 const dimension_type max,
2719 const Topology topol,
2720 const char* method,
2721 const char* reason) {
2722 const char* const domain = (topol == NECESSARILY_CLOSED)
2723 ? "PPL::C_Polyhedron::" : "PPL::NNC_Polyhedron::";
2724 return PPL::check_space_dimension_overflow(dim, max, domain, method, reason);
2725 }
2726
2727 PPL::dimension_type
check_space_dimension_overflow(const dimension_type dim,const Topology topol,const char * method,const char * reason)2728 PPL::Polyhedron::check_space_dimension_overflow(const dimension_type dim,
2729 const Topology topol,
2730 const char* method,
2731 const char* reason) {
2732 return check_space_dimension_overflow(dim, Polyhedron::max_space_dimension(),
2733 topol, method, reason);
2734 }
2735
2736 void
throw_invalid_generator(const char * method,const char * g_name) const2737 PPL::Polyhedron::throw_invalid_generator(const char* method,
2738 const char* g_name) const {
2739 std::ostringstream s;
2740 s << "PPL::";
2741 if (is_necessarily_closed()) {
2742 s << "C_";
2743 }
2744 else {
2745 s << "NNC_";
2746 }
2747 s << "Polyhedron::" << method << ":" << std::endl
2748 << "*this is an empty polyhedron and "
2749 << g_name << " is not a point.";
2750 throw std::invalid_argument(s.str());
2751 }
2752
2753 void
throw_invalid_generators(const char * method,const char * gs_name) const2754 PPL::Polyhedron::throw_invalid_generators(const char* method,
2755 const char* gs_name) const {
2756 std::ostringstream s;
2757 s << "PPL::";
2758 if (is_necessarily_closed()) {
2759 s << "C_";
2760 }
2761 else {
2762 s << "NNC_";
2763 }
2764 s << "Polyhedron::" << method << ":" << std::endl
2765 << "*this is an empty polyhedron and" << std::endl
2766 << "the non-empty generator system " << gs_name << " contains no points.";
2767 throw std::invalid_argument(s.str());
2768 }
2769