1 /* Polyhedron class implementation (non-inline public functions).
2 Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
3 Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
4
5 This file is part of the Parma Polyhedra Library (PPL).
6
7 The PPL is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 The PPL is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software Foundation,
19 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
20
21 For the most up-to-date information see the Parma Polyhedra Library
22 site: http://bugseng.com/products/ppl/ . */
23
24 #include "ppl-config.h"
25 #include "Polyhedron_defs.hh"
26 #include "C_Polyhedron_defs.hh"
27 #include "NNC_Polyhedron_defs.hh"
28 #include "Scalar_Products_defs.hh"
29 #include "Scalar_Products_inlines.hh"
30 #include "MIP_Problem_defs.hh"
31 #include "wrap_assign.hh"
32 #include "assertions.hh"
33 #include <cstdlib>
34 #include <iostream>
35
36 #ifndef ENSURE_SORTEDNESS
37 #define ENSURE_SORTEDNESS 0
38 #endif
39
40 namespace PPL = Parma_Polyhedra_Library;
41
42 PPL::dimension_type* PPL::Polyhedron::simplify_num_saturators_p = 0;
43
44 size_t PPL::Polyhedron::simplify_num_saturators_size = 0;
45
46 void
initialize()47 PPL::Polyhedron::initialize() {
48 PPL_ASSERT(simplify_num_saturators_p == 0);
49 PPL_ASSERT(simplify_num_saturators_size == 0);
50 simplify_num_saturators_p = new dimension_type[simplify_num_saturators_size];
51 }
52
53 void
finalize()54 PPL::Polyhedron::finalize() {
55 delete [] simplify_num_saturators_p;
56 simplify_num_saturators_p = 0;
57 simplify_num_saturators_size = 0;
58 }
59
60 PPL::dimension_type
affine_dimension() const61 PPL::Polyhedron::affine_dimension() const {
62 if (is_empty()) {
63 return 0;
64 }
65
66 const Constraint_System& cs = minimized_constraints();
67 dimension_type d = space_dim;
68 for (Constraint_System::const_iterator i = cs.begin(),
69 cs_end = cs.end(); i != cs_end; ++i) {
70 if (i->is_equality()) {
71 --d;
72 }
73 }
74 return d;
75 }
76
77 const PPL::Constraint_System&
constraints() const78 PPL::Polyhedron::constraints() const {
79 if (marked_empty()) {
80 // We want `con_sys' to only contain the unsatisfiable constraint
81 // of the appropriate dimension.
82 if (con_sys.has_no_rows()) {
83 // The 0-dim unsatisfiable constraint is extended to
84 // the appropriate dimension and then stored in `con_sys'.
85 Constraint_System unsat_cs = Constraint_System::zero_dim_empty();
86 unsat_cs.adjust_topology_and_space_dimension(topology(), space_dim);
87 swap(const_cast<Constraint_System&>(con_sys), unsat_cs);
88 }
89 else {
90 // Checking that `con_sys' contains the right thing.
91 PPL_ASSERT(con_sys.space_dimension() == space_dim);
92 PPL_ASSERT(con_sys.num_rows() == 1);
93 PPL_ASSERT(con_sys[0].is_inconsistent());
94 }
95 return con_sys;
96 }
97
98 if (space_dim == 0) {
99 // Zero-dimensional universe.
100 PPL_ASSERT(con_sys.num_rows() == 0 && con_sys.space_dimension() == 0);
101 return con_sys;
102 }
103
104 // If the polyhedron has pending generators, we process them to obtain
105 // the constraints. No processing is needed if the polyhedron has
106 // pending constraints.
107 if (has_pending_generators()) {
108 process_pending_generators();
109 }
110 else if (!constraints_are_up_to_date()) {
111 update_constraints();
112 }
113
114 // TODO: reconsider whether to really sort constraints at this stage.
115 #if ENSURE_SORTEDNESS
116 // We insist in returning a sorted system of constraints,
117 // but sorting is useless if there are pending constraints.
118 if (!has_pending_constraints())
119 obtain_sorted_constraints();
120 #endif
121 return con_sys;
122 }
123
124 const PPL::Constraint_System&
minimized_constraints() const125 PPL::Polyhedron::minimized_constraints() const {
126 // `minimize()' or `strongly_minimize_constraints()'
127 // will process any pending constraints or generators.
128 if (is_necessarily_closed()) {
129 minimize();
130 }
131 else {
132 strongly_minimize_constraints();
133 }
134 return constraints();
135 }
136
137 const PPL::Generator_System&
generators() const138 PPL::Polyhedron::generators() const {
139 if (marked_empty()) {
140 PPL_ASSERT(gen_sys.has_no_rows());
141 // We want `gen_sys' to have the appropriate space dimension,
142 // even though it is an empty generator system.
143 if (gen_sys.space_dimension() != space_dim) {
144 Generator_System gs;
145 gs.adjust_topology_and_space_dimension(topology(), space_dim);
146 swap(const_cast<Generator_System&>(gen_sys), gs);
147 }
148 return gen_sys;
149 }
150
151 if (space_dim == 0) {
152 PPL_ASSERT(gen_sys.num_rows() == 0 && gen_sys.space_dimension() == 0);
153 return Generator_System::zero_dim_univ();
154 }
155
156 // If the polyhedron has pending constraints, we process them to obtain
157 // the generators (we may discover that the polyhedron is empty).
158 // No processing is needed if the polyhedron has pending generators.
159 if ((has_pending_constraints() && !process_pending_constraints())
160 || (!generators_are_up_to_date() && !update_generators())) {
161 // We have just discovered that `*this' is empty.
162 PPL_ASSERT(gen_sys.has_no_rows());
163 // We want `gen_sys' to have the appropriate space dimension,
164 // even though it is an empty generator system.
165 if (gen_sys.space_dimension() != space_dim) {
166 Generator_System gs;
167 gs.adjust_topology_and_space_dimension(topology(), space_dim);
168 swap(const_cast<Generator_System&>(gen_sys), gs);
169 }
170 return gen_sys;
171 }
172
173 // TODO: reconsider whether to really sort generators at this stage.
174 #if ENSURE_SORTEDNESS
175 // We insist in returning a sorted system of generators,
176 // but sorting is useless if there are pending generators.
177 if (!has_pending_generators())
178 obtain_sorted_generators();
179 #else
180 // In the case of an NNC polyhedron, if the generator system is fully
181 // minimized (i.e., minimized and with no pending generator), then
182 // return a sorted system of generators: this is needed so that the
183 // const_iterator could correctly filter out the matched closure points.
184 if (!is_necessarily_closed()
185 && generators_are_minimized() && !has_pending_generators()) {
186 obtain_sorted_generators();
187 }
188 #endif
189 return gen_sys;
190 }
191
192 const PPL::Generator_System&
minimized_generators() const193 PPL::Polyhedron::minimized_generators() const {
194 // `minimize()' or `strongly_minimize_generators()'
195 // will process any pending constraints or generators.
196 if (is_necessarily_closed()) {
197 minimize();
198 }
199 else {
200 strongly_minimize_generators();
201 }
202 // Note: calling generators() on a strongly minimized NNC generator
203 // system will also ensure sortedness, which is required to correctly
204 // filter away the matched closure points.
205 return generators();
206 }
207
208 PPL::Poly_Con_Relation
relation_with(const Constraint & c) const209 PPL::Polyhedron::relation_with(const Constraint& c) const {
210 // Dimension-compatibility check.
211 if (space_dim < c.space_dimension()) {
212 throw_dimension_incompatible("relation_with(c)", "c", c);
213 }
214
215 if (marked_empty()) {
216 return Poly_Con_Relation::saturates()
217 && Poly_Con_Relation::is_included()
218 && Poly_Con_Relation::is_disjoint();
219 }
220 if (space_dim == 0) {
221 if (c.is_inconsistent()) {
222 if (c.is_strict_inequality() && c.inhomogeneous_term() == 0) {
223 // The constraint 0 > 0 implicitly defines the hyperplane 0 = 0;
224 // thus, the zero-dimensional point also saturates it.
225 return Poly_Con_Relation::saturates()
226 && Poly_Con_Relation::is_disjoint();
227 }
228 else {
229 return Poly_Con_Relation::is_disjoint();
230 }
231 }
232 else if (c.is_equality() || c.inhomogeneous_term() == 0) {
233 return Poly_Con_Relation::saturates()
234 && Poly_Con_Relation::is_included();
235 }
236 else {
237 // The zero-dimensional point saturates
238 // neither the positivity constraint 1 >= 0,
239 // nor the strict positivity constraint 1 > 0.
240 return Poly_Con_Relation::is_included();
241 }
242 }
243
244 if ((has_pending_constraints() && !process_pending_constraints())
245 || (!generators_are_up_to_date() && !update_generators())) {
246 // The polyhedron is empty.
247 return Poly_Con_Relation::saturates()
248 && Poly_Con_Relation::is_included()
249 && Poly_Con_Relation::is_disjoint();
250 }
251 return gen_sys.relation_with(c);
252 }
253
254 PPL::Poly_Gen_Relation
relation_with(const Generator & g) const255 PPL::Polyhedron::relation_with(const Generator& g) const {
256 // Dimension-compatibility check.
257 if (space_dim < g.space_dimension()) {
258 throw_dimension_incompatible("relation_with(g)", "g", g);
259 }
260
261 // The empty polyhedron cannot subsume a generator.
262 if (marked_empty()) {
263 return Poly_Gen_Relation::nothing();
264 }
265
266 // A universe polyhedron in a zero-dimensional space subsumes
267 // all the generators of a zero-dimensional space.
268 if (space_dim == 0) {
269 return Poly_Gen_Relation::subsumes();
270 }
271
272 if (has_pending_generators()) {
273 process_pending_generators();
274 }
275 else if (!constraints_are_up_to_date()) {
276 update_constraints();
277 }
278 return
279 con_sys.satisfies_all_constraints(g)
280 ? Poly_Gen_Relation::subsumes()
281 : Poly_Gen_Relation::nothing();
282 }
283
284 PPL::Poly_Con_Relation
relation_with(const Congruence & cg) const285 PPL::Polyhedron::relation_with(const Congruence& cg) const {
286 const dimension_type cg_space_dim = cg.space_dimension();
287 // Dimension-compatibility check.
288 if (space_dim < cg_space_dim) {
289 throw_dimension_incompatible("relation_with(cg)", "cg", cg);
290 }
291
292 if (cg.is_equality()) {
293 const Constraint c(cg);
294 return relation_with(c);
295 }
296
297 if (marked_empty()) {
298 return Poly_Con_Relation::saturates()
299 && Poly_Con_Relation::is_included()
300 && Poly_Con_Relation::is_disjoint();
301 }
302
303 if (space_dim == 0) {
304 if (cg.is_inconsistent()) {
305 return Poly_Con_Relation::is_disjoint();
306 }
307 else {
308 return Poly_Con_Relation::saturates()
309 && Poly_Con_Relation::is_included();
310 }
311 }
312
313 if ((has_pending_constraints() && !process_pending_constraints())
314 || (!generators_are_up_to_date() && !update_generators())) {
315 // The polyhedron is empty.
316 return Poly_Con_Relation::saturates()
317 && Poly_Con_Relation::is_included()
318 && Poly_Con_Relation::is_disjoint();
319 }
320 // Build the equality corresponding to the congruence (ignoring the modulus).
321 Linear_Expression expr(cg.expression());
322 const Constraint c(expr == 0);
323
324 // The polyhedron is non-empty so that there exists a point.
325 // For an arbitrary generator point, compute the scalar product with
326 // the equality.
327 PPL_DIRTY_TEMP_COEFFICIENT(sp_point);
328 for (Generator_System::const_iterator gs_i = gen_sys.begin(),
329 gs_end = gen_sys.end(); gs_i != gs_end; ++gs_i) {
330 if (gs_i->is_point()) {
331 Scalar_Products::assign(sp_point, c, *gs_i);
332 expr -= sp_point;
333 break;
334 }
335 }
336
337 // Find two hyperplanes that satisfy the congruence and are near to
338 // the generating point (so that the point lies on or between these
339 // two hyperplanes).
340 // Then use the relations between the polyhedron and the halfspaces
341 // corresponding to the hyperplanes to determine the result.
342
343 // Compute the distance from the point to an hyperplane.
344 const Coefficient& modulus = cg.modulus();
345 PPL_DIRTY_TEMP_COEFFICIENT(signed_distance);
346 signed_distance = sp_point % modulus;
347 if (signed_distance == 0) {
348 // The point is lying on the hyperplane.
349 return relation_with(expr == 0);
350 }
351 else {
352 // The point is not lying on the hyperplane.
353 expr += signed_distance;
354 }
355 // Build first halfspace constraint.
356 const bool positive = (signed_distance > 0);
357 const Constraint first_halfspace = positive ? (expr >= 0) : (expr <= 0);
358
359 const Poly_Con_Relation first_rels = relation_with(first_halfspace);
360 PPL_ASSERT(!first_rels.implies(Poly_Con_Relation::saturates())
361 && !first_rels.implies(Poly_Con_Relation::is_disjoint()));
362 if (first_rels.implies(Poly_Con_Relation::strictly_intersects())) {
363 return Poly_Con_Relation::strictly_intersects();
364 }
365
366 // Build second halfspace.
367 if (positive) {
368 expr -= modulus;
369 }
370 else {
371 expr += modulus;
372 }
373 const Constraint second_halfspace = positive ? (expr <= 0) : (expr >= 0);
374
375 PPL_ASSERT(first_rels == Poly_Con_Relation::is_included());
376 const Poly_Con_Relation second_rels = relation_with(second_halfspace);
377 PPL_ASSERT(!second_rels.implies(Poly_Con_Relation::saturates())
378 && !second_rels.implies(Poly_Con_Relation::is_disjoint()));
379 if (second_rels.implies(Poly_Con_Relation::strictly_intersects())) {
380 return Poly_Con_Relation::strictly_intersects();
381 }
382
383 PPL_ASSERT(second_rels == Poly_Con_Relation::is_included());
384 return Poly_Con_Relation::is_disjoint();
385 }
386
387 bool
is_universe() const388 PPL::Polyhedron::is_universe() const {
389 if (marked_empty()) {
390 return false;
391 }
392
393 if (space_dim == 0) {
394 return true;
395 }
396 if (!has_pending_generators() && constraints_are_up_to_date()) {
397 // Search for a constraint that is not a tautology.
398 for (dimension_type i = con_sys.num_rows(); i-- > 0; ) {
399 if (!con_sys[i].is_tautological()) {
400 return false;
401 }
402 }
403 // All the constraints are tautologies.
404 return true;
405 }
406
407 PPL_ASSERT(!has_pending_constraints() && generators_are_up_to_date());
408
409 // Try a fast-fail test.
410 dimension_type num_lines = 0;
411 dimension_type num_rays = 0;
412 const dimension_type first_pending = gen_sys.first_pending_row();
413 for (dimension_type i = first_pending; i-- > 0; ) {
414 switch (gen_sys[i].type()) {
415 case Generator::RAY:
416 ++num_rays;
417 break;
418 case Generator::LINE:
419 ++num_lines;
420 break;
421 default:
422 break;
423 }
424 }
425
426 if (has_pending_generators()) {
427 // The non-pending part of `gen_sys' was minimized:
428 // a success-first test is possible in this case.
429 PPL_ASSERT(generators_are_minimized());
430 if (num_lines == space_dim) {
431 PPL_ASSERT(num_rays == 0);
432 return true;
433 }
434 PPL_ASSERT(num_lines < space_dim);
435 // Now scan the pending generators.
436 dimension_type num_pending_lines = 0;
437 dimension_type num_pending_rays = 0;
438 const dimension_type gs_num_rows = gen_sys.num_rows();
439 for (dimension_type i = first_pending; i < gs_num_rows; ++i) {
440 switch (gen_sys[i].type()) {
441 case Generator::RAY:
442 ++num_pending_rays;
443 break;
444 case Generator::LINE:
445 ++num_pending_lines;
446 break;
447 default:
448 break;
449 }
450 }
451 // If no pending rays and lines were found,
452 // then it is not the universe polyhedron.
453 if (num_pending_rays == 0 && num_pending_lines == 0) {
454 return false;
455 }
456 // Factor away the lines already seen (to be on the safe side,
457 // we assume they are all linearly independent).
458 if (num_lines + num_pending_lines < space_dim) {
459 const dimension_type num_dims_missing
460 = space_dim - (num_lines + num_pending_lines);
461 // In order to span an n dimensional space (where n = num_dims_missing),
462 // at least n+1 rays are needed.
463 if (num_rays + num_pending_rays <= num_dims_missing) {
464 return false;
465 }
466 }
467 }
468 else {
469 // There is nothing pending.
470 if (generators_are_minimized()) {
471 // The exact test is possible.
472 PPL_ASSERT(num_rays == 0 || num_lines < space_dim);
473 return num_lines == space_dim;
474 }
475 else {
476 // Only the fast-fail test can be computed: in order to span
477 // an n dimensional space (where n = space_dim - num_lines),
478 // at least n+1 rays are needed.
479 if (num_lines < space_dim && num_lines + num_rays <= space_dim) {
480 return false;
481 }
482 }
483 }
484
485 // We need the polyhedron in minimal form.
486 if (has_pending_generators()) {
487 process_pending_generators();
488 }
489 else if (!constraints_are_minimized()) {
490 minimize();
491 }
492 if (is_necessarily_closed()) {
493 return con_sys.num_rows() == 1
494 && con_sys[0].is_inequality()
495 && con_sys[0].is_tautological();
496 }
497 else {
498 // NNC polyhedron.
499 if (con_sys.num_rows() != 2
500 || con_sys[0].is_equality()
501 || con_sys[1].is_equality()) {
502 return false;
503 }
504 else {
505 // If the system of constraints contains two rows that
506 // are not equalities, we are sure that they are
507 // epsilon constraints: in this case we know that
508 // the polyhedron is universe.
509 #ifndef NDEBUG
510 obtain_sorted_constraints();
511 const Constraint& eps_leq_one = con_sys[0];
512 const Constraint& eps_geq_zero = con_sys[1];
513 PPL_ASSERT(eps_leq_one.inhomogeneous_term() > 0
514 && eps_leq_one.epsilon_coefficient() < 0
515 && eps_geq_zero.inhomogeneous_term() == 0
516 && eps_geq_zero.epsilon_coefficient() > 0);
517 PPL_ASSERT(eps_leq_one.expression().all_homogeneous_terms_are_zero());
518 PPL_ASSERT(eps_geq_zero.expression().all_homogeneous_terms_are_zero());
519 #endif
520 return true;
521 }
522 }
523 }
524
525 bool
is_bounded() const526 PPL::Polyhedron::is_bounded() const {
527 // A zero-dimensional or empty polyhedron is bounded.
528 if (space_dim == 0
529 || marked_empty()
530 || (has_pending_constraints() && !process_pending_constraints())
531 || (!generators_are_up_to_date() && !update_generators())) {
532 return true;
533 }
534
535 // If the system of generators contains any line or a ray,
536 // then the polyhedron is unbounded.
537 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
538 if (gen_sys[i].is_line_or_ray()) {
539 return false;
540 }
541 }
542
543 // The system of generators is composed only by
544 // points and closure points: the polyhedron is bounded.
545 return true;
546 }
547
548 bool
is_topologically_closed() const549 PPL::Polyhedron::is_topologically_closed() const {
550 // Necessarily closed polyhedra are trivially closed.
551 if (is_necessarily_closed()) {
552 return true;
553 }
554 // Any empty or zero-dimensional polyhedron is closed.
555 if (marked_empty()
556 || space_dim == 0
557 || (has_something_pending() && !process_pending())) {
558 return true;
559 }
560
561 // At this point there are no pending constraints or generators.
562 PPL_ASSERT(!has_something_pending());
563
564 if (generators_are_minimized()) {
565 // A polyhedron is closed if and only if all of its (non-redundant)
566 // closure points are matched by a corresponding point.
567 const dimension_type n_rows = gen_sys.num_rows();
568 const dimension_type n_lines = gen_sys.num_lines();
569 for (dimension_type i = n_rows; i-- > n_lines; ) {
570 const Generator& gen_sys_i = gen_sys[i];
571 if (gen_sys_i.is_closure_point()) {
572 bool gen_sys_i_has_no_matching_point = true;
573 for (dimension_type j = n_rows; j-- > n_lines; ) {
574 const Generator& gen_sys_j = gen_sys[j];
575 if (i != j
576 && gen_sys_j.is_point()
577 && gen_sys_i.is_matching_closure_point(gen_sys_j)) {
578 gen_sys_i_has_no_matching_point = false;
579 break;
580 }
581 }
582 if (gen_sys_i_has_no_matching_point) {
583 return false;
584 }
585 }
586 }
587 // All closure points are matched.
588 return true;
589 }
590
591 // A polyhedron is closed if, after strong minimization
592 // of its constraint system, it has no strict inequalities.
593 strongly_minimize_constraints();
594 return marked_empty() || !con_sys.has_strict_inequalities();
595 }
596
597 bool
contains_integer_point() const598 PPL::Polyhedron::contains_integer_point() const {
599 // Any empty polyhedron does not contain integer points.
600 if (marked_empty()) {
601 return false;
602 }
603 // A zero-dimensional, universe polyhedron has, by convention, an
604 // integer point.
605 if (space_dim == 0) {
606 return true;
607 }
608
609 // CHECKME: do we really want to call conversion to check for emptiness?
610 if (has_pending_constraints() && !process_pending()) {
611 // Empty again.
612 return false;
613 }
614
615 // FIXME: do also exploit info regarding rays and lines, if possible.
616 // Is any integer point already available?
617 PPL_ASSERT(!has_pending_constraints());
618 if (generators_are_up_to_date()) {
619 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
620 if (gen_sys[i].is_point() && gen_sys[i].divisor() == 1) {
621 return true;
622 }
623 }
624 }
625 const Constraint_System& cs = constraints();
626 #if 0 // TEMPORARILY DISABLED.
627 MIP_Problem mip(space_dim,
628 cs.begin(), cs.end(),
629 Variables_Set(Variable(0), Variable(space_dim-1)));
630 #else
631 // FIXME: temporary workaround, to be removed as soon as the MIP
632 // problem class will correctly and precisely handle
633 // ((strict) in-) equality constraints having all integer variables.
634 MIP_Problem mip(space_dim);
635 mip.add_to_integer_space_dimensions(Variables_Set(Variable(0),
636 Variable(space_dim-1)));
637 PPL_DIRTY_TEMP_COEFFICIENT(homogeneous_gcd);
638 PPL_DIRTY_TEMP_COEFFICIENT(gcd);
639 PPL_DIRTY_TEMP(mpq_class, rational_inhomogeneous);
640 PPL_DIRTY_TEMP_COEFFICIENT(tightened_inhomogeneous);
641 for (Constraint_System::const_iterator cs_i = cs.begin(),
642 cs_end = cs.end(); cs_i != cs_end; ++cs_i) {
643 const Constraint& c = *cs_i;
644 const Constraint::Type c_type = c.type();
645 const Coefficient& inhomogeneous = c.inhomogeneous_term();
646 if (c_type == Constraint::STRICT_INEQUALITY) {
647 // CHECKME: should we change the behavior of Linear_Expression(c) ?
648 // Compute the GCD of the coefficients of c
649 // (disregarding the inhomogeneous term and the epsilon dimension).
650 homogeneous_gcd = c.expression().gcd(1, space_dim + 1);
651 if (homogeneous_gcd == 0) {
652 // NOTE: since tautological constraints are already filtered away
653 // by iterators, here we must have an inconsistent constraint.
654 PPL_ASSERT(c.is_inconsistent());
655 return false;
656 }
657 Linear_Expression le(c.expression());
658 if (homogeneous_gcd != 1) {
659 le /= homogeneous_gcd;
660 }
661 // Further tighten the constraint if the inhomogeneous term
662 // was integer, i.e., if `homogeneous_gcd' divides `inhomogeneous'.
663 gcd_assign(gcd, homogeneous_gcd, inhomogeneous);
664 if (gcd == homogeneous_gcd) {
665 le -= 1;
666 }
667 mip.add_constraint(le >= 0);
668 }
669 else {
670 // Equality or non-strict inequality.
671 // If possible, avoid useless gcd computations.
672 if (inhomogeneous == 0) {
673 // The inhomogeneous term cannot be tightened.
674 mip.add_constraint(c);
675 }
676 else {
677 // Compute the GCD of the coefficients of c
678 // (disregarding the inhomogeneous term)
679 // to see whether or not the inhomogeneous term can be tightened.
680 homogeneous_gcd = c.expression().gcd(1, space_dim + 1);
681 if (homogeneous_gcd == 0) {
682 // NOTE: since tautological constraints are already filtered away
683 // by iterators, here we must have an inconsistent constraint.
684 PPL_ASSERT(c.is_inconsistent());
685 return false;
686 }
687 else if (homogeneous_gcd == 1) {
688 // The normalized inhomogeneous term is integer:
689 // add the constraint as-is.
690 mip.add_constraint(c);
691 }
692 else {
693 PPL_ASSERT(homogeneous_gcd > 1);
694 // Here the normalized inhomogeneous term is rational:
695 // the constraint has to be tightened.
696 #ifndef NDEBUG
697 // `homogeneous_gcd' does not divide `inhomogeneous'.
698 // FIXME: add a divisibility test for Coefficient.
699 gcd_assign(gcd, homogeneous_gcd, inhomogeneous);
700 PPL_ASSERT(gcd == 1);
701 #endif
702 if (c.type() == Constraint::EQUALITY) {
703 return false;
704 }
705 // Extract the homogeneous part of the constraint.
706 Linear_Expression le(c.expression());
707 le -= inhomogeneous;
708 // Tighten the inhomogeneous term.
709 assign_r(rational_inhomogeneous.get_num(),
710 inhomogeneous, ROUND_NOT_NEEDED);
711 assign_r(rational_inhomogeneous.get_den(),
712 homogeneous_gcd, ROUND_NOT_NEEDED);
713 // Note: canonicalization is not needed (as gcd == 1).
714 PPL_ASSERT(is_canonical(rational_inhomogeneous));
715 assign_r(tightened_inhomogeneous,
716 rational_inhomogeneous, ROUND_DOWN);
717 tightened_inhomogeneous *= homogeneous_gcd;
718 le += tightened_inhomogeneous;
719 mip.add_constraint(le >= 0);
720 }
721 }
722 }
723 }
724 #endif // TEMPORARY WORKAROUND.
725 return mip.is_satisfiable();
726 }
727
728 bool
constrains(const Variable var) const729 PPL::Polyhedron::constrains(const Variable var) const {
730 // `var' should be one of the dimensions of the polyhedron.
731 const dimension_type var_space_dim = var.space_dimension();
732 if (space_dim < var_space_dim) {
733 throw_dimension_incompatible("constrains(v)", "v", var);
734 }
735
736 // An empty polyhedron constrains all variables.
737 if (marked_empty()) {
738 return true;
739 }
740
741 if (generators_are_up_to_date() && !has_pending_constraints()) {
742 // Since generators are up-to-date and there are no pending
743 // constraints, the generator system (since it is well formed)
744 // contains a point. Hence the polyhedron is not empty.
745 if (constraints_are_up_to_date() && !has_pending_generators()) {
746 // Here a variable is constrained if and only if it is
747 // syntactically constrained.
748 goto syntactic_check;
749 }
750 if (generators_are_minimized()) {
751 // Try a quick, incomplete check for the universe polyhedron:
752 // a universe polyhedron constrains no variable.
753 // Count the number of non-pending
754 // (hence, linearly independent) lines.
755 dimension_type num_lines = 0;
756 const dimension_type first_pending = gen_sys.first_pending_row();
757 for (dimension_type i = first_pending; i-- > 0; ) {
758 if (gen_sys[i].is_line()) {
759 ++num_lines;
760 }
761 }
762
763 if (num_lines == space_dim) {
764 return false;
765 }
766 }
767
768 // Scan generators: perhaps we will find a generator equivalent to
769 // line(var) or a pair of generators equivalent to ray(-var) and
770 // ray(var).
771 bool have_positive_ray = false;
772 bool have_negative_ray = false;
773 const dimension_type var_id = var.id();
774 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
775 const Generator& gen_sys_i = gen_sys[i];
776 if (gen_sys_i.is_line_or_ray()) {
777 const int sign = sgn(gen_sys_i.coefficient(var));
778 if (sign != 0) {
779 if (gen_sys_i.expression().all_zeroes(1, var_id)
780 && gen_sys_i.expression().all_zeroes(var_id + 1, space_dim + 1)) {
781
782 if (gen_sys_i.is_line()) {
783 return true;
784 }
785 if (sign > 0) {
786 if (have_negative_ray) {
787 return true;
788 }
789 else {
790 have_positive_ray = true;
791 }
792 }
793 else if (have_positive_ray) {
794 return true;
795 }
796 else {
797 have_negative_ray = true;
798 }
799 }
800 }
801 }
802 }
803
804 // We are still here: at least we know that, since generators are
805 // up-to-date and there are no pending constraints, then the
806 // generator system (since it is well formed) contains a point.
807 // Hence the polyhedron is not empty.
808 if (has_pending_generators()) {
809 process_pending_generators();
810 }
811 else if (!constraints_are_up_to_date()) {
812 update_constraints();
813 }
814 goto syntactic_check;
815 }
816
817 // We must minimize to detect emptiness and obtain constraints.
818 if (!minimize()) {
819 return true;
820 }
821
822 syntactic_check:
823 for (dimension_type i = con_sys.num_rows(); i-- > 0; ) {
824 if (con_sys[i].coefficient(var) != 0) {
825 return true;
826 }
827 }
828 return false;
829 }
830
831 bool
OK(bool check_not_empty) const832 PPL::Polyhedron::OK(bool check_not_empty) const {
833 #ifndef NDEBUG
834 using std::endl;
835 using std::cerr;
836 #endif
837
838 // Check whether the topologies of `con_sys' and `gen_sys' agree.
839 if (con_sys.topology() != gen_sys.topology()) {
840 #ifndef NDEBUG
841 cerr << "Constraints and generators have different topologies!"
842 << endl;
843 #endif
844 goto bomb;
845 }
846
847 // Check whether the status information is legal.
848 if (!status.OK()) {
849 goto bomb;
850 }
851
852 if (marked_empty()) {
853 if (check_not_empty) {
854 // The caller does not want the polyhedron to be empty.
855 #ifndef NDEBUG
856 cerr << "Empty polyhedron!" << endl;
857 #endif
858 goto bomb;
859 }
860
861 // An empty polyhedron is allowed if the system of constraints
862 // either has no rows or only contains an unsatisfiable constraint
863 // and if it has no pending constraints or generators.
864 if (has_something_pending()) {
865 #ifndef NDEBUG
866 cerr << "The polyhedron is empty, "
867 << "but it has something pending" << endl;
868 #endif
869 goto bomb;
870 }
871 if (con_sys.has_no_rows()) {
872 return true;
873 }
874 else {
875 if (con_sys.space_dimension() != space_dim) {
876 #ifndef NDEBUG
877 cerr << "The polyhedron is in a space of dimension "
878 << space_dim
879 << " while the system of constraints is in a space of dimension "
880 << con_sys.space_dimension()
881 << endl;
882 #endif
883 goto bomb;
884 }
885 if (con_sys.num_rows() != 1) {
886 #ifndef NDEBUG
887 cerr << "The system of constraints for an empty polyhedron "
888 << "has more then one row"
889 << endl;
890 #endif
891 goto bomb;
892 }
893 if (!con_sys[0].is_inconsistent()) {
894 #ifndef NDEBUG
895 cerr << "Empty polyhedron with a satisfiable system of constraints"
896 << endl;
897 #endif
898 goto bomb;
899 }
900 // Here we have only one, inconsistent constraint.
901 return true;
902 }
903 }
904
905 // A zero-dimensional, non-empty polyhedron is legal only if the
906 // system of constraint `con_sys' and the system of generators
907 // `gen_sys' have no rows.
908 if (space_dim == 0) {
909 if (has_something_pending()) {
910 #ifndef NDEBUG
911 cerr << "Zero-dimensional polyhedron with something pending"
912 << endl;
913 #endif
914 goto bomb;
915 }
916 if (!con_sys.has_no_rows() || !gen_sys.has_no_rows()) {
917 #ifndef NDEBUG
918 cerr << "Zero-dimensional polyhedron with a non-empty"
919 << endl
920 << "system of constraints or generators."
921 << endl;
922 #endif
923 goto bomb;
924 }
925 return true;
926 }
927
928 // A polyhedron is defined by a system of constraints
929 // or a system of generators: at least one of them must be up to date.
930 if (!constraints_are_up_to_date() && !generators_are_up_to_date()) {
931 #ifndef NDEBUG
932 cerr << "Polyhedron not empty, not zero-dimensional"
933 << endl
934 << "and with neither constraints nor generators up-to-date!"
935 << endl;
936 #endif
937 goto bomb;
938 }
939
940 // Here we check if the size of the matrices is consistent.
941 // Let us suppose that all the matrices are up-to-date; this means:
942 // `con_sys' : number of constraints x poly_num_columns
943 // `gen_sys' : number of generators x poly_num_columns
944 // `sat_c' : number of generators x number of constraints
945 // `sat_g' : number of constraints x number of generators.
946 if (constraints_are_up_to_date()) {
947 if (con_sys.space_dimension() != space_dim) {
948 #ifndef NDEBUG
949 cerr << "Incompatible size! (con_sys and space_dim)"
950 << endl;
951 #endif
952 goto bomb;
953 }
954 if (sat_c_is_up_to_date()) {
955 if (con_sys.first_pending_row() != sat_c.num_columns()) {
956 #ifndef NDEBUG
957 cerr << "Incompatible size! (con_sys and sat_c)"
958 << endl;
959 #endif
960 goto bomb;
961 }
962 }
963 if (sat_g_is_up_to_date()) {
964 if (con_sys.first_pending_row() != sat_g.num_rows()) {
965 #ifndef NDEBUG
966 cerr << "Incompatible size! (con_sys and sat_g)"
967 << endl;
968 #endif
969 goto bomb;
970 }
971 }
972 if (generators_are_up_to_date()) {
973 if (con_sys.space_dimension() != gen_sys.space_dimension()) {
974 #ifndef NDEBUG
975 cerr << "Incompatible size! (con_sys and gen_sys)"
976 << endl;
977 #endif
978 goto bomb;
979 }
980 }
981 }
982
983 if (generators_are_up_to_date()) {
984 if (gen_sys.space_dimension() != space_dim) {
985 #ifndef NDEBUG
986 cerr << "Incompatible size! (gen_sys and space_dim)"
987 << endl;
988 #endif
989 goto bomb;
990 }
991 if (sat_c_is_up_to_date()) {
992 if (gen_sys.first_pending_row() != sat_c.num_rows()) {
993 #ifndef NDEBUG
994 cerr << "Incompatible size! (gen_sys and sat_c)"
995 << endl;
996 #endif
997 goto bomb;
998 }
999 }
1000 if (sat_g_is_up_to_date()) {
1001 if (gen_sys.first_pending_row() != sat_g.num_columns()) {
1002 #ifndef NDEBUG
1003 cerr << "Incompatible size! (gen_sys and sat_g)"
1004 << endl;
1005 #endif
1006 goto bomb;
1007 }
1008 }
1009 if (gen_sys.first_pending_row() == 0) {
1010 #ifndef NDEBUG
1011 cerr << "Up-to-date generator system with all rows pending!"
1012 << endl;
1013 #endif
1014 goto bomb;
1015 }
1016
1017 // A non-empty system of generators describing a polyhedron
1018 // is valid if and only if it contains a point.
1019 if (!gen_sys.has_no_rows() && !gen_sys.has_points()) {
1020 #ifndef NDEBUG
1021 cerr << "Non-empty generator system declared up-to-date "
1022 << "has no points!"
1023 << endl;
1024 #endif
1025 goto bomb;
1026 }
1027
1028 #if 0 // To be activated when Status keeps strong minimization flags.
1029 //=================================================
1030 // TODO: this test is wrong in the general case.
1031 // However, such an invariant does hold for a
1032 // strongly-minimized Generator_System.
1033 // We will activate this test as soon as the Status
1034 // flags will be able to remember if a system is
1035 // strongly minimized.
1036
1037 // Checking that the number of closure points is always
1038 // greater than the number of points.
1039 if (!is_necessarily_closed()) {
1040 dimension_type num_points = 0;
1041 dimension_type num_closure_points = 0;
1042 dimension_type eps_index = gen_sys.space_dimension() + 1;
1043 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
1044 if (!gen_sys[i].is_line_or_ray()) {
1045 if (gen_sys[i][eps_index] > 0) {
1046 ++num_points;
1047 }
1048 else {
1049 ++num_closure_points;
1050 }
1051 }
1052 }
1053 if (num_points > num_closure_points) {
1054 #ifndef NDEBUG
1055 cerr << "# POINTS > # CLOSURE_POINTS" << endl;
1056 #endif
1057 goto bomb;
1058 }
1059 }
1060 //=================================================
1061 #endif
1062
1063 if (generators_are_minimized()) {
1064 // If the system of generators is minimized, the number of
1065 // lines, rays and points of the polyhedron must be the same as
1066 // of a temporary, minimized one. If this does not happen then
1067 // the polyhedron is not OK.
1068 Constraint_System new_con_sys(topology(), default_con_sys_repr);
1069 Generator_System gs_without_pending = gen_sys;
1070 gs_without_pending.remove_trailing_rows(gs_without_pending.num_rows()
1071 - gen_sys.first_pending_row());
1072 gs_without_pending.unset_pending_rows();
1073 Generator_System copy_of_gen_sys = gs_without_pending;
1074 Bit_Matrix new_sat_c;
1075 minimize(false, copy_of_gen_sys, new_con_sys, new_sat_c);
1076 const dimension_type copy_num_lines = copy_of_gen_sys.num_lines();
1077 if (gs_without_pending.num_rows() != copy_of_gen_sys.num_rows()
1078 || gs_without_pending.num_lines() != copy_num_lines
1079 || gs_without_pending.num_rays() != copy_of_gen_sys.num_rays()) {
1080 #ifndef NDEBUG
1081 cerr << "Generators are declared minimized, but they are not!\n"
1082 << "Here is the minimized form of the generators:\n";
1083 copy_of_gen_sys.ascii_dump(cerr);
1084 cerr << endl;
1085 #endif
1086 goto bomb;
1087 }
1088
1089 // CHECKME : the following observation is not formally true
1090 // for a NNC_Polyhedron. But it may be true for its
1091 // representation ...
1092
1093 // If the corresponding polyhedral cone is _pointed_, then
1094 // a minimal system of generators is unique up to positive scaling.
1095 // We thus verify if the cone is pointed (i.e., there are no lines)
1096 // and, after normalizing and sorting a copy of the system `gen_sys'
1097 // of the polyhedron (we use a copy not to modify the polyhedron's
1098 // system) and the system `copy_of_gen_sys' that has been just
1099 // minimized, we check if the two matrices are identical. If
1100 // they are different it means that the generators of the
1101 // polyhedron are declared minimized, but they are not.
1102 if (copy_num_lines == 0) {
1103 copy_of_gen_sys.strong_normalize();
1104 copy_of_gen_sys.sort_rows();
1105 gs_without_pending.strong_normalize();
1106 gs_without_pending.sort_rows();
1107 if (copy_of_gen_sys != gs_without_pending) {
1108 #ifndef NDEBUG
1109 cerr << "Generators are declared minimized, but they are not!\n"
1110 << "(we are in the case:\n"
1111 << "dimension of lineality space equal to 0)\n"
1112 << "Here is the minimized form of the generators:\n";
1113 copy_of_gen_sys.ascii_dump(cerr);
1114 cerr << endl;
1115 #endif
1116 goto bomb;
1117 }
1118 }
1119 }
1120 }
1121
1122 if (constraints_are_up_to_date()) {
1123 if (con_sys.first_pending_row() == 0) {
1124 #ifndef NDEBUG
1125 cerr << "Up-to-date constraint system with all rows pending!"
1126 << endl;
1127 #endif
1128 goto bomb;
1129 }
1130
1131 // A non-empty system of constraints describing a polyhedron
1132 // must contain a constraint with a non-zero inhomogeneous term;
1133 // such a constraint corresponds to (a combination of other
1134 // constraints with):
1135 // -* the positivity constraint, for necessarily closed polyhedra;
1136 // -* the epsilon <= 1 constraint, for NNC polyhedra.
1137 bool no_positivity_constraint = true;
1138 for (dimension_type i = con_sys.num_rows(); i-- > 0; ) {
1139 if (con_sys[i].inhomogeneous_term() != 0) {
1140 no_positivity_constraint = false;
1141 break;
1142 }
1143 }
1144 if (no_positivity_constraint) {
1145 #ifndef NDEBUG
1146 cerr << "Non-empty constraint system has no positivity constraint"
1147 << endl;
1148 #endif
1149 goto bomb;
1150 }
1151
1152 if (!is_necessarily_closed()) {
1153 // A non-empty system of constraints describing a NNC polyhedron
1154 // must also contain a (combination of) the constraint epsilon >= 0,
1155 // i.e., a constraint with a positive epsilon coefficient.
1156 bool no_epsilon_geq_zero = true;
1157 for (dimension_type i = con_sys.num_rows(); i-- > 0; ) {
1158 if (con_sys[i].epsilon_coefficient() > 0) {
1159 no_epsilon_geq_zero = false;
1160 break;
1161 }
1162 }
1163 if (no_epsilon_geq_zero) {
1164 #ifndef NDEBUG
1165 cerr << "Non-empty constraint system for NNC polyhedron "
1166 << "has no epsilon >= 0 constraint"
1167 << endl;
1168 #endif
1169 goto bomb;
1170 }
1171 }
1172
1173 Constraint_System cs_without_pending = con_sys;
1174 cs_without_pending.remove_trailing_rows(cs_without_pending.num_rows()
1175 - con_sys.first_pending_row());
1176 cs_without_pending.unset_pending_rows();
1177 Constraint_System copy_of_con_sys = cs_without_pending;
1178 bool empty = false;
1179 if (check_not_empty || constraints_are_minimized()) {
1180 Generator_System new_gen_sys(topology(), default_gen_sys_repr);
1181 Bit_Matrix new_sat_g;
1182 empty = minimize(true, copy_of_con_sys, new_gen_sys, new_sat_g);
1183 }
1184
1185 if (empty && check_not_empty) {
1186 #ifndef NDEBUG
1187 cerr << "Unsatisfiable system of constraints!"
1188 << endl;
1189 #endif
1190 goto bomb;
1191 }
1192
1193 if (constraints_are_minimized()) {
1194 // If the constraints are minimized, the number of equalities
1195 // and of inequalities of the system of the polyhedron must be
1196 // the same of the temporary minimized one.
1197 // If it does not happen, the polyhedron is not OK.
1198 if (cs_without_pending.num_rows() != copy_of_con_sys.num_rows()
1199 || cs_without_pending.num_equalities()
1200 != copy_of_con_sys.num_equalities()) {
1201 #ifndef NDEBUG
1202 cerr << "Constraints are declared minimized, but they are not!\n"
1203 << "Here is the minimized form of the constraints:\n";
1204 copy_of_con_sys.ascii_dump(cerr);
1205 cerr << endl;
1206 #endif
1207 goto bomb;
1208 }
1209 // The system `copy_of_con_sys' has the form that is obtained
1210 // after applying methods gauss() and back_substitute().
1211 // A system of constraints can be minimal even if it does not
1212 // have this form. So, to verify if the polyhedron is correct,
1213 // we copy the system `con_sys' in a temporary one and then
1214 // modify it using method simplify() (which calls both gauss()
1215 // and back_substitute()).
1216 // If the temporary system and `copy_of_con_sys' are different,
1217 // the polyhedron is not OK.
1218 copy_of_con_sys.strong_normalize();
1219 copy_of_con_sys.sort_rows();
1220 cs_without_pending.simplify();
1221 cs_without_pending.strong_normalize();
1222 cs_without_pending.sort_rows();
1223 if (cs_without_pending != copy_of_con_sys) {
1224 #ifndef NDEBUG
1225 cerr << "Constraints are declared minimized, but they are not!\n"
1226 << "Here is the minimized form of the constraints:\n";
1227 copy_of_con_sys.ascii_dump(cerr);
1228 cerr << endl;
1229 #endif
1230 goto bomb;
1231 }
1232 }
1233 }
1234
1235 if (sat_c_is_up_to_date()) {
1236 for (dimension_type i = sat_c.num_rows(); i-- > 0; ) {
1237 const Generator tmp_gen = gen_sys[i];
1238 const Bit_Row tmp_sat = sat_c[i];
1239 for (dimension_type j = sat_c.num_columns(); j-- > 0; ) {
1240 const bool sat_j = (Scalar_Products::sign(con_sys[j], tmp_gen) == 0);
1241 if (sat_j == tmp_sat[j]) {
1242 #ifndef NDEBUG
1243 cerr << "sat_c is declared up-to-date, but it is not!"
1244 << endl;
1245 #endif
1246 goto bomb;
1247 }
1248 }
1249 }
1250 }
1251 if (sat_g_is_up_to_date()) {
1252 for (dimension_type i = sat_g.num_rows(); i-- > 0; ) {
1253 const Constraint tmp_con = con_sys[i];
1254 const Bit_Row tmp_sat = sat_g[i];
1255 for (dimension_type j = sat_g.num_columns(); j-- > 0; ) {
1256 const bool sat_j = (Scalar_Products::sign(tmp_con, gen_sys[j]) == 0);
1257 if (sat_j == tmp_sat[j]) {
1258 #ifndef NDEBUG
1259 cerr << "sat_g is declared up-to-date, but it is not!"
1260 << endl;
1261 #endif
1262 goto bomb;
1263 }
1264 }
1265 }
1266 }
1267 if (has_pending_constraints()) {
1268 if (con_sys.num_pending_rows() == 0) {
1269 #ifndef NDEBUG
1270 cerr << "The polyhedron is declared to have pending constraints, "
1271 << "but con_sys has no pending rows!"
1272 << endl;
1273 #endif
1274 goto bomb;
1275 }
1276 }
1277
1278 if (has_pending_generators()) {
1279 if (gen_sys.num_pending_rows() == 0) {
1280 #ifndef NDEBUG
1281 cerr << "The polyhedron is declared to have pending generators, "
1282 << "but gen_sys has no pending rows!"
1283 << endl;
1284 #endif
1285 goto bomb;
1286 }
1287 }
1288
1289 return true;
1290
1291 bomb:
1292 #ifndef NDEBUG
1293 cerr << "Here is the guilty polyhedron:"
1294 << endl;
1295 ascii_dump(cerr);
1296 #endif
1297 return false;
1298 }
1299
1300 void
add_constraint(const Constraint & c)1301 PPL::Polyhedron::add_constraint(const Constraint& c) {
1302 // Topology-compatibility check.
1303 if (c.is_strict_inequality() && is_necessarily_closed()) {
1304 // Trivially true/false strict inequalities are legal.
1305 if (c.is_tautological()) {
1306 return;
1307 }
1308 if (c.is_inconsistent()) {
1309 set_empty();
1310 return;
1311 }
1312 // Here c is a non-trivial strict inequality.
1313 throw_topology_incompatible("add_constraint(c)", "c", c);
1314 }
1315
1316 // Dimension-compatibility check:
1317 // the dimension of `c' can not be greater than space_dim.
1318 if (space_dim < c.space_dimension()) {
1319 throw_dimension_incompatible("add_constraint(c)", "c", c);
1320 }
1321
1322 if (!marked_empty()) {
1323 refine_no_check(c);
1324 }
1325
1326 }
1327
1328 void
add_congruence(const Congruence & cg)1329 PPL::Polyhedron::add_congruence(const Congruence& cg) {
1330 // Dimension-compatibility check:
1331 // the dimension of `cg' can not be greater than space_dim.
1332 if (space_dim < cg.space_dimension()) {
1333 throw_dimension_incompatible("add_congruence(cg)", "cg", cg);
1334 }
1335
1336 // Handle the case of proper congruences first.
1337 if (cg.is_proper_congruence()) {
1338 if (cg.is_tautological()) {
1339 return;
1340 }
1341 if (cg.is_inconsistent()) {
1342 set_empty();
1343 return;
1344 }
1345 // Non-trivial and proper congruences are not allowed.
1346 throw_invalid_argument("add_congruence(cg)",
1347 "cg is a non-trivial, proper congruence");
1348 }
1349
1350 PPL_ASSERT(cg.is_equality());
1351 // Handle empty and 0-dim cases first.
1352 if (marked_empty()) {
1353 return;
1354 }
1355 if (space_dim == 0) {
1356 if (cg.is_inconsistent()) {
1357 set_empty();
1358 }
1359 return;
1360 }
1361
1362 // Add the equality.
1363 Linear_Expression le(cg.expression());
1364 const Constraint c(le, Constraint::EQUALITY, NECESSARILY_CLOSED);
1365 refine_no_check(c);
1366 }
1367
1368 void
add_generator(const Generator & g)1369 PPL::Polyhedron::add_generator(const Generator& g) {
1370 // Topology-compatibility check.
1371 if (g.is_closure_point() && is_necessarily_closed()) {
1372 throw_topology_incompatible("add_generator(g)", "g", g);
1373 }
1374 // Dimension-compatibility check:
1375 // the dimension of `g' can not be greater than space_dim.
1376 const dimension_type g_space_dim = g.space_dimension();
1377 if (space_dim < g_space_dim) {
1378 throw_dimension_incompatible("add_generator(g)", "g", g);
1379 }
1380 // Dealing with a zero-dimensional space polyhedron first.
1381 if (space_dim == 0) {
1382 // It is not possible to create 0-dim rays or lines.
1383 PPL_ASSERT(g.is_point() || g.is_closure_point());
1384 // Closure points can only be inserted in non-empty polyhedra.
1385 if (marked_empty()) {
1386 if (g.type() != Generator::POINT) {
1387 throw_invalid_generator("add_generator(g)", "g");
1388 }
1389 else {
1390 set_zero_dim_univ();
1391 }
1392 }
1393 PPL_ASSERT_HEAVY(OK());
1394 return;
1395 }
1396
1397 if (marked_empty()
1398 || (has_pending_constraints() && !process_pending_constraints())
1399 || (!generators_are_up_to_date() && !update_generators())) {
1400 // Here the polyhedron is empty:
1401 // the specification says we can only insert a point.
1402 if (!g.is_point()) {
1403 throw_invalid_generator("add_generator(g)", "g");
1404 }
1405 if (g.is_necessarily_closed() || !is_necessarily_closed()) {
1406 gen_sys.insert(g);
1407 // Since `gen_sys' was empty, after inserting `g' we have to resize
1408 // the system of generators to have the right dimension.
1409 gen_sys.adjust_topology_and_space_dimension(topology(), space_dim);
1410 if (!is_necessarily_closed()) {
1411 // In the NNC topology, each point has to be matched by
1412 // a corresponding closure point:
1413 // turn the just inserted point into the corresponding
1414 // (normalized) closure point.
1415 gen_sys.sys.rows.back().set_epsilon_coefficient(0);
1416 gen_sys.sys.rows.back().expr.normalize();
1417 PPL_ASSERT(gen_sys.sys.rows.back().OK());
1418 PPL_ASSERT(gen_sys.sys.OK());
1419 // Re-insert the point (which is already normalized).
1420 gen_sys.insert(g);
1421 }
1422 }
1423 else {
1424 // Note: here we have a _legal_ topology mismatch,
1425 // because `g' is NOT a closure point (it is a point!)
1426 // However, by barely invoking `gen_sys.insert(g)' we would
1427 // cause a change in the topology of `gen_sys', which is wrong.
1428 // Thus, we insert a "topology corrected" copy of `g'.
1429 const Linear_Expression nc_expr(g.expression());
1430 gen_sys.insert(Generator::point(nc_expr, g.divisor()));
1431 // Since `gen_sys' was empty, after inserting `g' we have to resize
1432 // the system of generators to have the right dimension.
1433 gen_sys.adjust_topology_and_space_dimension(topology(), space_dim);
1434 }
1435 // No longer empty, generators up-to-date and minimized.
1436 clear_empty();
1437 set_generators_minimized();
1438 }
1439 else {
1440 PPL_ASSERT(generators_are_up_to_date());
1441 const bool has_pending = can_have_something_pending();
1442 if (g.is_necessarily_closed() || !is_necessarily_closed()) {
1443 // Since `gen_sys' is not empty, the topology and space dimension
1444 // of the inserted generator are automatically adjusted.
1445 if (has_pending) {
1446 gen_sys.insert_pending(g);
1447 }
1448 else {
1449 gen_sys.insert(g);
1450 }
1451 if (!is_necessarily_closed() && g.is_point()) {
1452 // In the NNC topology, each point has to be matched by
1453 // a corresponding closure point:
1454 // turn the just inserted point into the corresponding
1455 // (normalized) closure point.
1456 gen_sys.sys.rows.back().set_epsilon_coefficient(0);
1457 gen_sys.sys.rows.back().expr.normalize();
1458 PPL_ASSERT(gen_sys.sys.rows.back().OK());
1459 PPL_ASSERT(gen_sys.sys.OK());
1460 // Re-insert the point (which is already normalized).
1461 if (has_pending) {
1462 gen_sys.insert_pending(g);
1463 }
1464 else {
1465 gen_sys.insert(g);
1466 }
1467 }
1468 }
1469 else {
1470 PPL_ASSERT(!g.is_closure_point());
1471 // Note: here we have a _legal_ topology mismatch, because
1472 // `g' is NOT a closure point.
1473 // However, by barely invoking `gen_sys.insert(g)' we would
1474 // cause a change in the topology of `gen_sys', which is wrong.
1475 // Thus, we insert a "topology corrected" copy of `g'.
1476 const Linear_Expression nc_expr(g.expression());
1477 switch (g.type()) {
1478 case Generator::LINE:
1479 if (has_pending) {
1480 gen_sys.insert_pending(Generator::line(nc_expr));
1481 }
1482 else {
1483 gen_sys.insert(Generator::line(nc_expr));
1484 }
1485 break;
1486 case Generator::RAY:
1487 if (has_pending) {
1488 gen_sys.insert_pending(Generator::ray(nc_expr));
1489 }
1490 else {
1491 gen_sys.insert(Generator::ray(nc_expr));
1492 }
1493 break;
1494 case Generator::POINT:
1495 if (has_pending) {
1496 gen_sys.insert_pending(Generator::point(nc_expr, g.divisor()));
1497 }
1498 else {
1499 gen_sys.insert(Generator::point(nc_expr, g.divisor()));
1500 }
1501 break;
1502 case Generator::CLOSURE_POINT:
1503 PPL_UNREACHABLE;
1504 break;
1505 }
1506 }
1507
1508 if (has_pending) {
1509 set_generators_pending();
1510 }
1511 else {
1512 // After adding the new generator,
1513 // constraints are no longer up-to-date.
1514 clear_generators_minimized();
1515 clear_constraints_up_to_date();
1516 }
1517 }
1518 PPL_ASSERT_HEAVY(OK());
1519 }
1520
1521 void
add_recycled_constraints(Constraint_System & cs)1522 PPL::Polyhedron::add_recycled_constraints(Constraint_System& cs) {
1523 // Topology compatibility check.
1524 if (is_necessarily_closed() && cs.has_strict_inequalities()) {
1525 // We check if _all_ strict inequalities in cs are trivially false.
1526 // (The iterators already filter away trivially true constraints.)
1527 for (Constraint_System::const_iterator i = cs.begin(),
1528 i_end = cs.end(); i != i_end; ++i) {
1529 if (i->is_strict_inequality()
1530 && !i->is_inconsistent()) {
1531 throw_topology_incompatible("add_recycled_constraints(cs)",
1532 "cs", cs);
1533 }
1534 }
1535 // If we reach this point, all strict inequalities were inconsistent.
1536 set_empty();
1537 return;
1538 }
1539
1540 // Dimension-compatibility check:
1541 // the dimension of `cs' can not be greater than space_dim.
1542 const dimension_type cs_space_dim = cs.space_dimension();
1543 if (space_dim < cs_space_dim) {
1544 throw_dimension_incompatible("add_recycled_constraints(cs)", "cs", cs);
1545 }
1546 // Adding no constraints is a no-op.
1547 if (cs.has_no_rows()) {
1548 return;
1549 }
1550
1551 if (space_dim == 0) {
1552 // In a 0-dimensional space the constraints are
1553 // tautologies (e.g., 0 == 0 or 1 >= 0 or 1 > 0) or
1554 // inconsistent (e.g., 1 == 0 or -1 >= 0 or 0 > 0).
1555 // In a system of constraints `begin()' and `end()' are equal
1556 // if and only if the system only contains tautologies.
1557 if (cs.begin() != cs.end()) {
1558 // There is a constraint, it must be inconsistent,
1559 // the polyhedron is empty.
1560 status.set_empty();
1561 }
1562 return;
1563 }
1564
1565 if (marked_empty()) {
1566 return;
1567 }
1568
1569 // The constraints (possibly with pending rows) are required.
1570 if (has_pending_generators()) {
1571 process_pending_generators();
1572 }
1573 else if (!constraints_are_up_to_date()) {
1574 update_constraints();
1575 }
1576
1577 // Adjust `cs' to the right topology and space dimension.
1578 // NOTE: we already checked for topology compatibility.
1579 cs.adjust_topology_and_space_dimension(topology(), space_dim);
1580
1581 const bool adding_pending = can_have_something_pending();
1582
1583 // Here we do not require `con_sys' to be sorted.
1584 // also, we _recycle_ (instead of copying) the coefficients of `cs'.
1585 if (adding_pending) {
1586 con_sys.insert_pending(cs, Recycle_Input());
1587 set_constraints_pending();
1588 }
1589 else {
1590 con_sys.insert(cs, Recycle_Input());
1591 // Constraints are not minimized and generators are not up-to-date.
1592 clear_constraints_minimized();
1593 clear_generators_up_to_date();
1594 }
1595
1596 // Note: the constraint system may have become unsatisfiable, thus
1597 // we do not check for satisfiability.
1598 PPL_ASSERT_HEAVY(OK());
1599 }
1600
1601 void
add_constraints(const Constraint_System & cs)1602 PPL::Polyhedron::add_constraints(const Constraint_System& cs) {
1603 // TODO: this is just an executable specification.
1604 Constraint_System cs_copy = cs;
1605 add_recycled_constraints(cs_copy);
1606 }
1607
1608 void
add_recycled_generators(Generator_System & gs)1609 PPL::Polyhedron::add_recycled_generators(Generator_System& gs) {
1610 // Topology compatibility check.
1611 if (is_necessarily_closed() && gs.has_closure_points()) {
1612 throw_topology_incompatible("add_recycled_generators(gs)", "gs", gs);
1613 }
1614 // Dimension-compatibility check:
1615 // the dimension of `gs' can not be greater than space_dim.
1616 const dimension_type gs_space_dim = gs.space_dimension();
1617 if (space_dim < gs_space_dim) {
1618 throw_dimension_incompatible("add_recycled_generators(gs)", "gs", gs);
1619 }
1620
1621 // Adding no generators is a no-op.
1622 if (gs.has_no_rows()) {
1623 return;
1624 }
1625
1626 // Adding valid generators to a zero-dimensional polyhedron
1627 // transform it in the zero-dimensional universe polyhedron.
1628 if (space_dim == 0) {
1629 if (marked_empty() && !gs.has_points()) {
1630 throw_invalid_generators("add_recycled_generators(gs)", "gs");
1631 }
1632 set_zero_dim_univ();
1633 PPL_ASSERT_HEAVY(OK(true));
1634 return;
1635 }
1636
1637 // Adjust `gs' to the right topology and dimensions.
1638 // NOTE: we already checked for topology compatibility.
1639 gs.adjust_topology_and_space_dimension(topology(), space_dim);
1640 // For NNC polyhedra, each point must be matched by
1641 // the corresponding closure point.
1642 if (!is_necessarily_closed()) {
1643 gs.add_corresponding_closure_points();
1644 }
1645
1646 // The generators (possibly with pending rows) are required.
1647 if ((has_pending_constraints() && !process_pending_constraints())
1648 || (!generators_are_up_to_date() && !minimize())) {
1649 // We have just discovered that `*this' is empty.
1650 // So `gs' must contain at least one point.
1651 if (!gs.has_points()) {
1652 throw_invalid_generators("add_recycled_generators(gs)", "gs");
1653 }
1654 // The polyhedron is no longer empty and generators are up-to-date.
1655 swap(gen_sys, gs);
1656 if (gen_sys.num_pending_rows() > 0) {
1657 // Even though `gs' has pending generators, since the constraints
1658 // of the polyhedron are not up-to-date, the polyhedron cannot
1659 // have pending generators. By integrating the pending part
1660 // of `gen_sys' we may loose sortedness.
1661 gen_sys.sys.index_first_pending = gen_sys.num_rows();
1662 gen_sys.set_sorted(false);
1663 }
1664 set_generators_up_to_date();
1665 clear_empty();
1666 PPL_ASSERT_HEAVY(OK());
1667 return;
1668 }
1669
1670 if (can_have_something_pending()) {
1671 // Here we do not require `gen_sys' to be sorted.
1672 // also, we _remove_ (instead of copying) the rows of `gs'
1673 // (which is not a const).
1674 for (dimension_type i = 0; i < gs.num_rows(); ++i) {
1675 gs.sys.rows[i].set_topology(topology());
1676 gen_sys.insert_pending(gs.sys.rows[i], Recycle_Input());
1677 }
1678 gs.clear();
1679
1680 set_generators_pending();
1681 }
1682 else {
1683 // Here we do not require `gen_sys' to be sorted.
1684 // also, we _remove_ (instead of copying) the coefficients of `gs'
1685 // (which is not a const).
1686 for (dimension_type i = 0; i < gs.num_rows(); ++i) {
1687 gs.sys.rows[i].set_topology(topology());
1688 gen_sys.insert(gs.sys.rows[i], Recycle_Input());
1689 }
1690 gs.clear();
1691
1692 // Constraints are not up-to-date and generators are not minimized.
1693 clear_constraints_up_to_date();
1694 clear_generators_minimized();
1695 }
1696 PPL_ASSERT_HEAVY(OK(true));
1697 }
1698
1699 void
add_generators(const Generator_System & gs)1700 PPL::Polyhedron::add_generators(const Generator_System& gs) {
1701 // TODO: this is just an executable specification.
1702 Generator_System gs_copy = gs;
1703 add_recycled_generators(gs_copy);
1704 }
1705
1706 void
add_congruences(const Congruence_System & cgs)1707 PPL::Polyhedron::add_congruences(const Congruence_System& cgs) {
1708 // Dimension-compatibility check.
1709 if (space_dim < cgs.space_dimension()) {
1710 throw_dimension_incompatible("add_congruences(cgs)", "cgs", cgs);
1711 }
1712
1713 Constraint_System cs;
1714 bool inserted = false;
1715 for (Congruence_System::const_iterator i = cgs.begin(),
1716 cgs_end = cgs.end(); i != cgs_end; ++i) {
1717 const Congruence& cg = *i;
1718 if (cg.is_equality()) {
1719 Linear_Expression le(cg.expression());
1720 const Constraint c(le, Constraint::EQUALITY, NECESSARILY_CLOSED);
1721
1722 // TODO: Consider stealing the row in c when adding it to cs.
1723 cs.insert(c);
1724 inserted = true;
1725 }
1726 else {
1727 PPL_ASSERT(cg.is_proper_congruence());
1728 if (cg.is_inconsistent()) {
1729 set_empty();
1730 return;
1731 }
1732 if (!cg.is_tautological()) {
1733 throw_invalid_argument("add_congruences(cgs)",
1734 "cgs has a non-trivial, proper congruence");
1735 }
1736 }
1737 }
1738 // Only add cs if it contains something.
1739 if (inserted) {
1740 add_recycled_constraints(cs);
1741 }
1742 }
1743
1744 void
refine_with_constraint(const Constraint & c)1745 PPL::Polyhedron::refine_with_constraint(const Constraint& c) {
1746 // Dimension-compatibility check.
1747 if (space_dim < c.space_dimension()) {
1748 throw_dimension_incompatible("refine_with_constraint(c)", "c", c);
1749 }
1750 // If the polyhedron is known to be empty, do nothing.
1751 if (!marked_empty()) {
1752 refine_no_check(c);
1753 }
1754 }
1755
1756 void
refine_with_congruence(const Congruence & cg)1757 PPL::Polyhedron::refine_with_congruence(const Congruence& cg) {
1758 // Dimension-compatibility check.
1759 if (space_dim < cg.space_dimension()) {
1760 throw_dimension_incompatible("refine_with_congruence(cg)", "cg", cg);
1761 }
1762
1763 // If the polyhedron is known to be empty, do nothing.
1764 if (marked_empty()) {
1765 return;
1766 }
1767
1768 // Dealing with a zero-dimensional space polyhedron first.
1769 if (space_dim == 0) {
1770 if (!cg.is_tautological()) {
1771 set_empty();
1772 }
1773 return;
1774 }
1775
1776 if (cg.is_equality()) {
1777 Linear_Expression le(cg.expression());
1778 const Constraint c(le, Constraint::EQUALITY, NECESSARILY_CLOSED);
1779 refine_no_check(c);
1780 }
1781 }
1782
1783 void
refine_with_constraints(const Constraint_System & cs)1784 PPL::Polyhedron::refine_with_constraints(const Constraint_System& cs) {
1785 // TODO: this is just an executable specification.
1786
1787 // Dimension-compatibility check.
1788 const dimension_type cs_space_dim = cs.space_dimension();
1789 if (space_dim < cs_space_dim) {
1790 throw_dimension_incompatible("refine_with_constraints(cs)a",
1791 "cs", cs);
1792 }
1793 // Adding no constraints is a no-op.
1794 if (cs.has_no_rows()) {
1795 return;
1796 }
1797
1798 if (space_dim == 0) {
1799 // In a 0-dimensional space the constraints are
1800 // tautologies (e.g., 0 == 0 or 1 >= 0 or 1 > 0) or
1801 // inconsistent (e.g., 1 == 0 or -1 >= 0 or 0 > 0).
1802 // In a system of constraints `begin()' and `end()' are equal
1803 // if and only if the system only contains tautologies.
1804 if (cs.begin() != cs.end()) {
1805 // There is a constraint, it must be inconsistent,
1806 // the polyhedron is empty.
1807 status.set_empty();
1808 }
1809 return;
1810 }
1811
1812 if (marked_empty()) {
1813 return;
1814 }
1815
1816 // The constraints (possibly with pending rows) are required.
1817 if (has_pending_generators()) {
1818 process_pending_generators();
1819 }
1820 else if (!constraints_are_up_to_date()) {
1821 update_constraints();
1822 }
1823
1824 const bool adding_pending = can_have_something_pending();
1825 for (dimension_type i = cs.num_rows(); i-- > 0; ) {
1826 const Constraint& c = cs[i];
1827
1828 if (c.is_necessarily_closed() || !is_necessarily_closed()) {
1829 // Since `con_sys' is not empty, the topology and space dimension
1830 // of the inserted constraint are automatically adjusted.
1831 if (adding_pending) {
1832 con_sys.insert_pending(c);
1833 }
1834 else {
1835 con_sys.insert(c);
1836 }
1837 }
1838 else {
1839 // Here we know that *this is necessarily closed so even if c is
1840 // topologically closed, by barely invoking `con_sys.insert(c)' we
1841 // would cause a change in the topology of `con_sys', which is
1842 // wrong. Thus, we insert a topology closed and "topology
1843 // corrected" version of `c'.
1844 const Linear_Expression nc_expr(c.expression());
1845 if (c.is_equality()) {
1846 if (adding_pending) {
1847 con_sys.insert_pending(nc_expr == 0);
1848 }
1849 else {
1850 con_sys.insert(nc_expr == 0);
1851 }
1852 }
1853 else {
1854 if (adding_pending) {
1855 con_sys.insert_pending(nc_expr >= 0);
1856 }
1857 else {
1858 con_sys.insert(nc_expr >= 0);
1859 }
1860 }
1861 }
1862 }
1863
1864 if (adding_pending) {
1865 set_constraints_pending();
1866 }
1867 else {
1868 // Constraints are not minimized and generators are not up-to-date.
1869 clear_constraints_minimized();
1870 clear_generators_up_to_date();
1871 }
1872
1873 // Note: the constraint system may have become unsatisfiable, thus
1874 // we do not check for satisfiability.
1875 PPL_ASSERT_HEAVY(OK());
1876 }
1877
1878 void
refine_with_congruences(const Congruence_System & cgs)1879 PPL::Polyhedron::refine_with_congruences(const Congruence_System& cgs) {
1880 // Dimension-compatibility check.
1881 if (space_dim < cgs.space_dimension()) {
1882 throw_dimension_incompatible("refine_with_congruences(cgs)", "cgs", cgs);
1883 }
1884
1885 Constraint_System cs;
1886 bool inserted = false;
1887 for (Congruence_System::const_iterator i = cgs.begin(),
1888 cgs_end = cgs.end(); i != cgs_end; ++i) {
1889 if (i->is_equality()) {
1890 Linear_Expression le(i->expression());
1891 const Constraint c(le, Constraint::EQUALITY, NECESSARILY_CLOSED);
1892
1893 // TODO: Consider stealing the row in c when adding it to cs.
1894 cs.insert(c);
1895 inserted = true;
1896 }
1897 else if (i->is_inconsistent()) {
1898 set_empty();
1899 return;
1900 }
1901 }
1902 // Only add cgs if congruences were inserted into cgs, as the
1903 // dimension of cs must be at most that of the polyhedron.
1904 if (inserted) {
1905 add_recycled_constraints(cs);
1906 }
1907 }
1908
1909 void
unconstrain(const Variable var)1910 PPL::Polyhedron::unconstrain(const Variable var) {
1911 // Dimension-compatibility check.
1912 if (space_dim < var.space_dimension()) {
1913 throw_dimension_incompatible("unconstrain(var)", var.space_dimension());
1914 }
1915 // Do something only if the polyhedron is non-empty.
1916 if (marked_empty()
1917 || (has_pending_constraints() && !process_pending_constraints())
1918 || (!generators_are_up_to_date() && !update_generators())) {
1919 // Empty: do nothing.
1920 return;
1921 }
1922
1923 PPL_ASSERT(generators_are_up_to_date());
1924 // Since `gen_sys' is not empty, the topology and space dimension
1925 // of the inserted generator are automatically adjusted.
1926 if (can_have_something_pending()) {
1927 gen_sys.insert_pending(Generator::line(var));
1928 set_generators_pending();
1929 }
1930 else {
1931 gen_sys.insert(Generator::line(var));
1932 // After adding the new generator,
1933 // constraints are no longer up-to-date.
1934 clear_generators_minimized();
1935 clear_constraints_up_to_date();
1936 }
1937 PPL_ASSERT_HEAVY(OK(true));
1938 }
1939
1940 void
unconstrain(const Variables_Set & vars)1941 PPL::Polyhedron::unconstrain(const Variables_Set& vars) {
1942 // The cylindrification with respect to no dimensions is a no-op.
1943 // This case also captures the only legal cylindrification
1944 // of a polyhedron in a 0-dim space.
1945 if (vars.empty()) {
1946 return;
1947 }
1948 // Dimension-compatibility check.
1949 const dimension_type min_space_dim = vars.space_dimension();
1950 if (space_dim < min_space_dim) {
1951 throw_dimension_incompatible("unconstrain(vs)", min_space_dim);
1952 }
1953
1954 // Do something only if the polyhedron is non-empty.
1955 if (marked_empty()
1956 || (has_pending_constraints() && !process_pending_constraints())
1957 || (!generators_are_up_to_date() && !update_generators())) {
1958 // Empty: do nothing.
1959 return;
1960 }
1961
1962 PPL_ASSERT(generators_are_up_to_date());
1963 // Since `gen_sys' is not empty, the topology and space dimension
1964 // of the inserted generators are automatically adjusted.
1965 Variables_Set::const_iterator vsi = vars.begin();
1966 Variables_Set::const_iterator vsi_end = vars.end();
1967 if (can_have_something_pending()) {
1968 for ( ; vsi != vsi_end; ++vsi) {
1969 gen_sys.insert_pending(Generator::line(Variable(*vsi)));
1970 }
1971 set_generators_pending();
1972 }
1973 else {
1974 for ( ; vsi != vsi_end; ++vsi) {
1975 gen_sys.insert(Generator::line(Variable(*vsi)));
1976 }
1977 // After adding the new generators,
1978 // constraints are no longer up-to-date.
1979 clear_generators_minimized();
1980 clear_constraints_up_to_date();
1981 }
1982 PPL_ASSERT_HEAVY(OK(true));
1983 }
1984
1985 void
intersection_assign(const Polyhedron & y)1986 PPL::Polyhedron::intersection_assign(const Polyhedron& y) {
1987 Polyhedron& x = *this;
1988 // Topology compatibility check.
1989 if (x.topology() != y.topology()) {
1990 throw_topology_incompatible("intersection_assign(y)", "y", y);
1991 }
1992 // Dimension-compatibility check.
1993 if (x.space_dim != y.space_dim) {
1994 throw_dimension_incompatible("intersection_assign(y)", "y", y);
1995 }
1996 // If one of the two polyhedra is empty, the intersection is empty.
1997 if (x.marked_empty()) {
1998 return;
1999 }
2000 if (y.marked_empty()) {
2001 x.set_empty();
2002 return;
2003 }
2004
2005 // If both polyhedra are zero-dimensional,
2006 // then at this point they are necessarily non-empty,
2007 // so that their intersection is non-empty too.
2008 if (x.space_dim == 0) {
2009 return;
2010 }
2011
2012 // Both systems of constraints have to be up-to-date,
2013 // possibly having pending constraints.
2014 if (x.has_pending_generators()) {
2015 x.process_pending_generators();
2016 }
2017 else if (!x.constraints_are_up_to_date()) {
2018 x.update_constraints();
2019 }
2020
2021 if (y.has_pending_generators()) {
2022 y.process_pending_generators();
2023 }
2024 else if (!y.constraints_are_up_to_date()) {
2025 y.update_constraints();
2026 }
2027
2028 // Here both systems are up-to-date and possibly have pending constraints
2029 // (but they cannot have pending generators).
2030 PPL_ASSERT(!x.has_pending_generators() && x.constraints_are_up_to_date());
2031 PPL_ASSERT(!y.has_pending_generators() && y.constraints_are_up_to_date());
2032
2033 // If `x' can support pending constraints,
2034 // the constraints of `y' are added as pending constraints of `x'.
2035 if (x.can_have_something_pending()) {
2036 x.con_sys.insert_pending(y.con_sys);
2037 x.set_constraints_pending();
2038 }
2039 else {
2040 // `x' cannot support pending constraints.
2041 // If both constraint systems are (fully) sorted, then we can
2042 // merge them; otherwise we simply add the second to the first.
2043 if (x.con_sys.is_sorted()
2044 && y.con_sys.is_sorted() && !y.has_pending_constraints()) {
2045 x.con_sys.merge_rows_assign(y.con_sys);
2046 }
2047 else {
2048 x.con_sys.insert(y.con_sys);
2049 }
2050 // Generators are no longer up-to-date and constraints are no
2051 // longer minimized.
2052 x.clear_generators_up_to_date();
2053 x.clear_constraints_minimized();
2054 }
2055 PPL_ASSERT_HEAVY(x.OK() && y.OK());
2056 }
2057
2058 namespace {
2059
2060 struct Ruled_Out_Pair {
2061 PPL::dimension_type constraint_index;
2062 PPL::dimension_type num_ruled_out;
2063 };
2064
2065 struct Ruled_Out_Less_Than {
operator ()__anon9eeb4caa0111::Ruled_Out_Less_Than2066 bool operator()(const Ruled_Out_Pair& x,
2067 const Ruled_Out_Pair& y) const {
2068 return x.num_ruled_out > y.num_ruled_out;
2069 }
2070 };
2071
2072 /*
2073 Modifies the vector of pointers \p ineqs_p, setting to 0 those entries
2074 that point to redundant inequalities or masked equalities.
2075 The redundancy test is based on saturation matrix \p sat and
2076 on knowing that there exists \p rank non-redundant equalities
2077 (they are implicit, i.e., not explicitly listed in \p ineqs_p).
2078 */
2079 void
drop_redundant_inequalities(std::vector<const PPL::Constraint * > & ineqs_p,const PPL::Topology topology,const PPL::Bit_Matrix & sat,const PPL::dimension_type rank)2080 drop_redundant_inequalities(std::vector<const PPL::Constraint*>& ineqs_p,
2081 const PPL::Topology topology,
2082 const PPL::Bit_Matrix& sat,
2083 const PPL::dimension_type rank) {
2084 using namespace Parma_Polyhedra_Library;
2085 const dimension_type num_rows = ineqs_p.size();
2086 PPL_ASSERT(num_rows > 0);
2087 // `rank' is the rank of the (implicit) system of equalities.
2088 const dimension_type space_dim = ineqs_p[0]->space_dimension();
2089 PPL_ASSERT(space_dim > 0 && space_dim >= rank);
2090 const dimension_type num_coefficients
2091 = space_dim + ((topology == NECESSARILY_CLOSED) ? 0U : 1U);
2092 const dimension_type min_sat = num_coefficients - rank;
2093 const dimension_type num_cols_sat = sat.num_columns();
2094
2095 // Perform quick redundancy test based on the number of saturators.
2096 for (dimension_type i = num_rows; i-- > 0; ) {
2097 if (sat[i].empty()) {
2098 // Masked equalities are redundant.
2099 ineqs_p[i] = 0;
2100 }
2101 else {
2102 const dimension_type num_sat = num_cols_sat - sat[i].count_ones();
2103 if (num_sat < min_sat) {
2104 ineqs_p[i] = 0;
2105 }
2106 }
2107 }
2108
2109 // Re-examine remaining inequalities.
2110 // Iteration index `i' is _intentionally_ increasing.
2111 for (dimension_type i = 0; i < num_rows; ++i) {
2112 if (ineqs_p[i] != 0) {
2113 for (dimension_type j = 0; j < num_rows; ++j) {
2114 bool strict_subset;
2115 if (ineqs_p[j] != 0 && i != j
2116 && subset_or_equal(sat[j], sat[i], strict_subset)) {
2117 if (strict_subset) {
2118 ineqs_p[i] = 0;
2119 break;
2120 }
2121 else {
2122 // Here `sat[j] == sat[i]'.
2123 ineqs_p[j] = 0;
2124 }
2125 }
2126 }
2127 }
2128 }
2129 }
2130
2131 } // namespace
2132
2133 template <typename Linear_System1, typename Row2>
2134 bool
add_to_system_and_check_independence(Linear_System1 & eq_sys,const Row2 & eq)2135 PPL::Polyhedron::add_to_system_and_check_independence(Linear_System1& eq_sys,
2136 const Row2& eq) {
2137 // Check if eq is linearly independent from eq_sys.
2138 PPL_ASSERT(eq.is_line_or_equality());
2139 eq_sys.insert(eq);
2140 const PPL::dimension_type eq_sys_num_rows = eq_sys.num_rows();
2141 const PPL::dimension_type rank = eq_sys.gauss(eq_sys_num_rows);
2142 if (rank == eq_sys_num_rows) {
2143 // eq is linearly independent from eq_sys.
2144 return true;
2145 }
2146 else {
2147 // eq is not linearly independent from eq_sys.
2148 PPL_ASSERT(rank == eq_sys_num_rows - 1);
2149 eq_sys.remove_trailing_rows(1);
2150 return false;
2151 }
2152 }
2153
2154 bool
simplify_using_context_assign(const Polyhedron & y)2155 PPL::Polyhedron::simplify_using_context_assign(const Polyhedron& y) {
2156 Polyhedron& x = *this;
2157 // Topology compatibility check.
2158 if (x.topology() != y.topology()) {
2159 throw_topology_incompatible("simplify_using_context_assign(y)", "y", y);
2160 }
2161 // Dimension-compatibility check.
2162 if (x.space_dim != y.space_dim) {
2163 throw_dimension_incompatible("simplify_using_context_assign(y)", "y", y);
2164 }
2165
2166 // Filter away the zero-dimensional case.
2167 if (x.space_dim == 0) {
2168 if (y.is_empty()) {
2169 x.set_zero_dim_univ();
2170 return false;
2171 }
2172 else {
2173 return !x.is_empty();
2174 }
2175 }
2176
2177 // If `y' is empty, the biggest enlargement for `x' is the universe.
2178 if (!y.minimize()) {
2179 Polyhedron ph(x.topology(), x.space_dim, UNIVERSE);
2180 m_swap(ph);
2181 return false;
2182 }
2183
2184 // If `x' is empty, the intersection is empty.
2185 if (!x.minimize()) {
2186 // Search for a constraint of `y' that is not a tautology.
2187 PPL_ASSERT(!y.has_pending_generators() && y.constraints_are_up_to_date());
2188 for (dimension_type i = y.con_sys.num_rows(); i-- > 0; ) {
2189 const Constraint& y_con_sys_i = y.con_sys[i];
2190 if (!y_con_sys_i.is_tautological()) {
2191 // Found: we obtain a constraint `c' contradicting the one we
2192 // found, and assign to `x' the polyhedron `ph' with `c' as
2193 // the only constraint.
2194 Polyhedron ph(x.topology(), x.space_dim, UNIVERSE);
2195 const Linear_Expression le(y_con_sys_i.expression());
2196 switch (y_con_sys_i.type()) {
2197 case Constraint::EQUALITY:
2198 ph.refine_no_check(le == 1);
2199 break;
2200 case Constraint::NONSTRICT_INEQUALITY:
2201 ph.refine_no_check(le <= -1);
2202 break;
2203 case Constraint::STRICT_INEQUALITY:
2204 ph.refine_no_check(le == 0);
2205 break;
2206 }
2207 m_swap(ph);
2208 PPL_ASSERT_HEAVY(OK());
2209 return false;
2210 }
2211 }
2212 // `y' is the universe: `x' cannot be enlarged.
2213 return false;
2214 }
2215
2216 PPL_ASSERT(x.constraints_are_minimized()
2217 && !x.has_something_pending()
2218 && y.generators_are_minimized()
2219 && !y.has_something_pending());
2220 const Constraint_System& x_cs = x.con_sys;
2221 const dimension_type x_cs_num_rows = x_cs.num_rows();
2222 const Generator_System& y_gs = y.gen_sys;
2223
2224 // Record into `redundant_by_y' the info about which constraints of
2225 // `x' are redundant in the context `y'. Count the number of
2226 // redundancies found.
2227 std::vector<bool> redundant_by_y(x_cs_num_rows, false);
2228 dimension_type num_redundant_by_y = 0;
2229 for (dimension_type i = 0; i < x_cs_num_rows; ++i) {
2230 if (y_gs.satisfied_by_all_generators(x_cs[i])) {
2231 redundant_by_y[i] = true;
2232 ++num_redundant_by_y;
2233 }
2234 }
2235
2236 Constraint_System result_cs;
2237
2238 if (num_redundant_by_y < x_cs_num_rows) {
2239 // Some constraints were not identified as redundant (yet?).
2240 const Constraint_System& y_cs = y.con_sys;
2241 const dimension_type y_cs_num_rows = y_cs.num_rows();
2242 // Compute into `z' the minimized intersection of `x' and `y'.
2243 const bool x_first = (x_cs_num_rows > y_cs_num_rows);
2244 Polyhedron z(x_first ? x : y);
2245 if (x_first) {
2246 z.add_constraints(y_cs);
2247 }
2248 else {
2249 // Only copy (and then recycle) the non-redundant constraints.
2250 Constraint_System tmp_cs;
2251 for (dimension_type i = 0; i < x_cs_num_rows; ++i) {
2252 if (!redundant_by_y[i]) {
2253 tmp_cs.insert(x_cs[i]);
2254 }
2255 }
2256 z.add_recycled_constraints(tmp_cs);
2257 }
2258 if (!z.minimize()) {
2259 // For necessarily closed polyhedra, the objective function is
2260 // the default, zero. Moreover, the default maximization mode is
2261 // OK, since we are only interested in satisfiability.
2262 MIP_Problem lp;
2263 if (x.is_necessarily_closed()) {
2264 lp.add_space_dimensions_and_embed(x.space_dim);
2265 lp.add_constraints(y_cs);
2266 }
2267 else {
2268 // KLUDGE: temporarily mark `y_cs' if it was necessarily
2269 // closed, so that we can interpret the epsilon dimension as a
2270 // standard dimension. Be careful to reset the topology of `cs'
2271 // even on exceptional execution path.
2272 const_cast<Constraint_System&>(y_cs).mark_as_necessarily_closed();
2273 try {
2274 lp.add_space_dimensions_and_embed(x.space_dim+1);
2275 lp.add_constraints(y_cs);
2276 const_cast<Constraint_System&>(y_cs).mark_as_not_necessarily_closed();
2277 }
2278 catch (...) {
2279 const_cast<Constraint_System&>(y_cs).mark_as_not_necessarily_closed();
2280 throw;
2281 }
2282 // For not necessarily closed polyhedra, we maximize the
2283 // epsilon dimension as we want to rule out the possibility
2284 // that all solutions have epsilon = 0. We are in this case
2285 // if and only if the maximal value for epsilon is 0.
2286 lp.set_objective_function(Variable(x.space_dim));
2287 }
2288 // We apply the following heuristics here: constraints of `x' that
2289 // are not made redundant by `y' are added to `lp' depending on
2290 // the number of generators of `y' they rule out (the more generators
2291 // they rule out, the sooner they are added). Of course, as soon
2292 // as `lp' becomes unsatisfiable, we stop adding.
2293 std::vector<Ruled_Out_Pair>
2294 ruled_out_vec(x_cs_num_rows - num_redundant_by_y);
2295 for (dimension_type i = 0, j = 0; i < x_cs_num_rows; ++i) {
2296 if (!redundant_by_y[i]) {
2297 const Constraint& c = x_cs[i];
2298 const Topology_Adjusted_Scalar_Product_Sign sps(c);
2299 dimension_type num_ruled_out_generators = 0;
2300 for (Generator_System::const_iterator k = y_gs.begin(),
2301 y_gs_end = y_gs.end(); k != y_gs_end; ++k) {
2302 const Generator& g = *k;
2303 const int sp_sign = sps(g, c);
2304 if (x.is_necessarily_closed()) {
2305 if (g.is_line()) {
2306 // Lines must saturate the constraint.
2307 if (sp_sign != 0) {
2308 goto ruled_out;
2309 }
2310 }
2311 else {
2312 // `g' is either a ray, a point or a closure point.
2313 if (c.is_inequality()) {
2314 // `c' is a non-strict inequality.
2315 if (sp_sign < 0) {
2316 goto ruled_out;
2317 }
2318 }
2319 else {
2320 // `c' is an equality.
2321 if (sp_sign != 0) {
2322 goto ruled_out;
2323 }
2324 }
2325 }
2326 }
2327 else {
2328 // The topology is not necessarily closed.
2329 switch (g.type()) {
2330 case Generator::LINE:
2331 // Lines must saturate the constraint.
2332 if (sp_sign != 0) {
2333 goto ruled_out;
2334 }
2335 break;
2336 case Generator::POINT:
2337 // Have to perform the special test when dealing with
2338 // a strict inequality.
2339 switch (c.type()) {
2340 case Constraint::EQUALITY:
2341 if (sp_sign != 0) {
2342 goto ruled_out;
2343 }
2344 break;
2345 case Constraint::NONSTRICT_INEQUALITY:
2346 if (sp_sign < 0) {
2347 goto ruled_out;
2348 }
2349 break;
2350 case Constraint::STRICT_INEQUALITY:
2351 if (sp_sign <= 0) {
2352 goto ruled_out;
2353 }
2354 break;
2355 }
2356 break;
2357 case Generator::RAY:
2358 // Intentionally fall through.
2359 case Generator::CLOSURE_POINT:
2360 if (c.is_inequality()) {
2361 // Constraint `c' is either a strict or a non-strict
2362 // inequality.
2363 if (sp_sign < 0) {
2364 goto ruled_out;
2365 }
2366 }
2367 else {
2368 // Constraint `c' is an equality.
2369 if (sp_sign != 0) {
2370 goto ruled_out;
2371 }
2372 }
2373 break;
2374 }
2375 }
2376 // If we reach this point, `g' satisfies `c'.
2377 continue;
2378 ruled_out:
2379 ++num_ruled_out_generators;
2380 }
2381 ruled_out_vec[j].constraint_index = i;
2382 ruled_out_vec[j].num_ruled_out = num_ruled_out_generators;
2383 ++j;
2384 }
2385 }
2386 std::sort(ruled_out_vec.begin(), ruled_out_vec.end(),
2387 Ruled_Out_Less_Than());
2388
2389 for (std::vector<Ruled_Out_Pair>::const_iterator
2390 j = ruled_out_vec.begin(),
2391 ruled_out_vec_end = ruled_out_vec.end();
2392 j != ruled_out_vec_end;
2393 ++j) {
2394 const Constraint& c = x_cs[j->constraint_index];
2395 result_cs.insert(c);
2396 if (x.is_necessarily_closed()) {
2397 lp.add_constraint(c);
2398 }
2399 else {
2400 // KLUDGE: temporarily mark `c' as if it was necessarily
2401 // closed, so that we can interpret the epsilon dimension as a
2402 // standard dimension. Be careful to reset the topology of `c'
2403 // also on exceptional execution paths.
2404 const_cast<Constraint&>(c).mark_as_necessarily_closed();
2405 try {
2406 lp.add_constraint(c);
2407 const_cast<Constraint&>(c).mark_as_not_necessarily_closed();
2408 }
2409 catch (...) {
2410 const_cast<Constraint&>(c).mark_as_not_necessarily_closed();
2411 throw;
2412 }
2413 }
2414 const MIP_Problem_Status status = lp.solve();
2415 switch (status) {
2416 case UNFEASIBLE_MIP_PROBLEM:
2417 break;
2418 case OPTIMIZED_MIP_PROBLEM:
2419 if (x.is_necessarily_closed()) {
2420 continue;
2421 }
2422 else {
2423 Coefficient num;
2424 Coefficient den;
2425 lp.optimal_value(num, den);
2426 if (num != 0) {
2427 continue;
2428 }
2429 }
2430 break;
2431 default:
2432 PPL_UNREACHABLE;
2433 break;
2434 }
2435 Polyhedron result_ph(x.topology(), x.space_dim, UNIVERSE);
2436 result_ph.add_constraints(result_cs);
2437 swap(x, result_ph);
2438 PPL_ASSERT_HEAVY(x.OK());
2439 return false;
2440 }
2441 // Cannot exit from here.
2442 PPL_UNREACHABLE;
2443 }
2444 else {
2445 // Here `z' is not empty and minimized.
2446 PPL_ASSERT(z.constraints_are_minimized()
2447 && z.generators_are_minimized()
2448 && !z.has_something_pending());
2449 const Constraint_System& z_cs = z.con_sys;
2450 const Generator_System& z_gs = z.gen_sys;
2451 const dimension_type z_gs_num_rows = z_gs.num_rows();
2452
2453 // Compute the number of equalities in x_cs, y_cs and z_cs
2454 // (exploiting minimal form knowledge).
2455 dimension_type x_cs_num_eq = 0;
2456 while (x_cs[x_cs_num_eq].is_equality()) {
2457 ++x_cs_num_eq;
2458 }
2459 dimension_type y_cs_num_eq = 0;
2460 while (y_cs[y_cs_num_eq].is_equality()) {
2461 ++y_cs_num_eq;
2462 }
2463 dimension_type z_cs_num_eq = 0;
2464 while (z_cs[z_cs_num_eq].is_equality()) {
2465 ++z_cs_num_eq;
2466 }
2467 PPL_ASSERT(x_cs_num_eq <= z_cs_num_eq && y_cs_num_eq <= z_cs_num_eq);
2468
2469 // Identify non-redundant equalities.
2470 Constraint_System non_redundant_eq;
2471 dimension_type num_non_redundant_eq = 0;
2472 const dimension_type needed_non_redundant_eq = z_cs_num_eq - y_cs_num_eq;
2473 Constraint_System eqs(x.topology());
2474 if (needed_non_redundant_eq > 0) {
2475 // Populate eqs with the equalities from y.
2476 for (dimension_type i = 0; i < y_cs_num_eq; ++i) {
2477 eqs.insert(y_cs[i]);
2478 }
2479 // Try to find another `needed_non_redundant_eq' linear independent
2480 // equalities among those from x.
2481 for (dimension_type i = 0; i < x_cs_num_eq; ++i) {
2482 const Constraint& x_cs_i = x_cs[i];
2483 if (add_to_system_and_check_independence(eqs, x_cs_i)) {
2484 // x_cs_i is linear independent.
2485 non_redundant_eq.insert(x_cs_i);
2486 ++num_non_redundant_eq;
2487 if (num_non_redundant_eq == needed_non_redundant_eq) {
2488 // Already found all the needed equalities.
2489 break;
2490 }
2491 }
2492 }
2493 // NOTE: if num_non_redundant_eq < needed_non_redundant_eq
2494 // then we haven't found all the needed equalities yet:
2495 // this means that some inequalities from x actually holds
2496 // as "masked" equalities in the context of y.
2497 PPL_ASSERT(eqs.num_rows() <= z_cs_num_eq);
2498 PPL_ASSERT(num_non_redundant_eq <= needed_non_redundant_eq);
2499 PPL_ASSERT(z_cs_num_eq - eqs.num_rows()
2500 == needed_non_redundant_eq - num_non_redundant_eq);
2501 }
2502
2503 // Identify non-redundant inequalities.
2504 // Avoid useless copies (no modifications are needed).
2505 std::vector<const Constraint*> non_redundant_ineq_p;
2506 // Fill non_redundant_ineq_p with (pointers to) inequalities
2507 // from y_cs ...
2508 for (dimension_type i = y_cs_num_eq; i < y_cs_num_rows; ++i) {
2509 non_redundant_ineq_p.push_back(&y_cs[i]);
2510 }
2511 // ... and (pointers to) non-redundant inequalities from x_cs.
2512 for (dimension_type i = x_cs_num_eq; i < x_cs_num_rows; ++i) {
2513 if (!redundant_by_y[i]) {
2514 non_redundant_ineq_p.push_back(&x_cs[i]);
2515 }
2516 }
2517 const dimension_type non_redundant_ineq_p_size
2518 = non_redundant_ineq_p.size();
2519 const dimension_type y_cs_num_ineq = y_cs_num_rows - y_cs_num_eq;
2520
2521 // Compute saturation info.
2522 const dimension_type sat_num_rows = non_redundant_ineq_p_size;
2523 Bit_Matrix sat(sat_num_rows, z_gs_num_rows);
2524 for (dimension_type i = sat_num_rows; i-- > 0; ) {
2525 const Constraint& non_redundant_ineq_i = *(non_redundant_ineq_p[i]);
2526 Bit_Row& sat_i = sat[i];
2527 for (dimension_type j = z_gs_num_rows; j-- > 0; ) {
2528 if (Scalar_Products::sign(non_redundant_ineq_i, z_gs[j]) != 0) {
2529 sat_i.set(j);
2530 }
2531 }
2532 if (sat_i.empty() && num_non_redundant_eq < needed_non_redundant_eq) {
2533 // `non_redundant_ineq_i' is actually masking an equality
2534 // and we are still looking for some masked inequalities.
2535 // Iteration goes downwards, so the inequality comes from x_cs.
2536 PPL_ASSERT(i >= y_cs_num_ineq);
2537 // Check if the equality is independent in eqs.
2538 Constraint masked_eq = non_redundant_ineq_i;
2539 masked_eq.set_is_line_or_equality();
2540 masked_eq.sign_normalize();
2541 if (add_to_system_and_check_independence(eqs, masked_eq)) {
2542 // It is independent: add the _inequality_ to non_redundant_eq.
2543 non_redundant_eq.insert(non_redundant_ineq_i);
2544 ++num_non_redundant_eq;
2545 }
2546 }
2547 }
2548 // Here we have already found all the needed (masked) equalities.
2549 PPL_ASSERT(num_non_redundant_eq == needed_non_redundant_eq);
2550
2551 drop_redundant_inequalities(non_redundant_ineq_p, x.topology(),
2552 sat, z_cs_num_eq);
2553
2554 // Place the non-redundant (masked) equalities into result_cs.
2555 swap(result_cs, non_redundant_eq);
2556 // Add to result_cs the non-redundant inequalities from x_cs,
2557 // i.e., those having indices no smaller than y_cs_num_ineq.
2558 for (dimension_type i = y_cs_num_ineq;
2559 i < non_redundant_ineq_p_size;
2560 ++i) {
2561 if (non_redundant_ineq_p[i] != 0) {
2562 result_cs.insert(*non_redundant_ineq_p[i]);
2563 }
2564 }
2565 }
2566 }
2567
2568 Polyhedron result_ph(x.topology(), x.space_dim, UNIVERSE);
2569 result_ph.add_recycled_constraints(result_cs);
2570 swap(x, result_ph);
2571 PPL_ASSERT_HEAVY(x.OK());
2572 return true;
2573 }
2574
2575 void
poly_hull_assign(const Polyhedron & y)2576 PPL::Polyhedron::poly_hull_assign(const Polyhedron& y) {
2577 Polyhedron& x = *this;
2578 // Topology compatibility check.
2579 if (x.topology() != y.topology()) {
2580 throw_topology_incompatible("poly_hull_assign(y)", "y", y);
2581 }
2582 // Dimension-compatibility check.
2583 if (x.space_dim != y.space_dim) {
2584 throw_dimension_incompatible("poly_hull_assign(y)", "y", y);
2585 }
2586
2587 // The poly-hull of a polyhedron `p' with an empty polyhedron is `p'.
2588 if (y.marked_empty()) {
2589 return;
2590 }
2591 if (x.marked_empty()) {
2592 x = y;
2593 return;
2594 }
2595
2596 // If both polyhedra are zero-dimensional,
2597 // then at this point they are necessarily universe polyhedra,
2598 // so that their poly-hull is the universe polyhedron too.
2599 if (x.space_dim == 0) {
2600 return;
2601 }
2602
2603 // Both systems of generators have to be up-to-date,
2604 // possibly having pending generators.
2605 if ((x.has_pending_constraints() && !x.process_pending_constraints())
2606 || (!x.generators_are_up_to_date() && !x.update_generators())) {
2607 // Discovered `x' empty when updating generators.
2608 x = y;
2609 return;
2610 }
2611 if ((y.has_pending_constraints() && !y.process_pending_constraints())
2612 || (!y.generators_are_up_to_date() && !y.update_generators())) {
2613 // Discovered `y' empty when updating generators.
2614 return;
2615 }
2616 // Here both systems are up-to-date and possibly have pending generators
2617 // (but they cannot have pending constraints).
2618 PPL_ASSERT(!x.has_pending_constraints() && x.generators_are_up_to_date());
2619 PPL_ASSERT(!y.has_pending_constraints() && y.generators_are_up_to_date());
2620
2621 // If `x' can support pending generators,
2622 // the generators of `y' are added as pending generators of `x'.
2623 if (x.can_have_something_pending()) {
2624 x.gen_sys.insert_pending(y.gen_sys);
2625 x.set_generators_pending();
2626 }
2627 else {
2628 // `x' cannot support pending generators.
2629 // If both generator systems are (fully) sorted, then we can merge
2630 // them; otherwise we simply add the second to the first.
2631 if (x.gen_sys.is_sorted()
2632 && y.gen_sys.is_sorted() && !y.has_pending_generators()) {
2633 x.gen_sys.merge_rows_assign(y.gen_sys);
2634 }
2635 else {
2636 x.gen_sys.insert(y.gen_sys);
2637 }
2638 // Constraints are no longer up-to-date
2639 // and generators are no longer minimized.
2640 x.clear_constraints_up_to_date();
2641 x.clear_generators_minimized();
2642 }
2643 // At this point both `x' and `y' are not empty.
2644 PPL_ASSERT_HEAVY(x.OK(true) && y.OK(true));
2645 }
2646
2647 void
poly_difference_assign(const Polyhedron & y)2648 PPL::Polyhedron::poly_difference_assign(const Polyhedron& y) {
2649 Polyhedron& x = *this;
2650 // Topology compatibility check.
2651 if (x.topology() != y.topology()) {
2652 throw_topology_incompatible("poly_difference_assign(y)", "y", y);
2653 // Dimension-compatibility check.
2654 }
2655 if (x.space_dim != y.space_dim) {
2656 throw_dimension_incompatible("poly_difference_assign(y)", "y", y);
2657 }
2658
2659 // The difference of a polyhedron `p' and an empty polyhedron is `p'.
2660 if (y.marked_empty()) {
2661 return;
2662 }
2663 // The difference of an empty polyhedron and of a polyhedron `p' is empty.
2664 if (x.marked_empty()) {
2665 return;
2666 }
2667
2668 // If both polyhedra are zero-dimensional,
2669 // then at this point they are necessarily universe polyhedra,
2670 // so that their difference is empty.
2671 if (x.space_dim == 0) {
2672 x.set_empty();
2673 return;
2674 }
2675
2676 // TODO: This is just an executable specification.
2677 // Have to find a more efficient method.
2678
2679 if (y.contains(x)) {
2680 x.set_empty();
2681 return;
2682 }
2683
2684 // Being lazy here is only harmful.
2685 // `minimize()' will process any pending constraints or generators.
2686 if (!y.minimize()) {
2687 return;
2688 }
2689 x.minimize();
2690
2691 Polyhedron new_polyhedron(topology(), x.space_dim, EMPTY);
2692
2693 const Constraint_System& y_cs = y.constraints();
2694 for (Constraint_System::const_iterator i = y_cs.begin(),
2695 y_cs_end = y_cs.end(); i != y_cs_end; ++i) {
2696 const Constraint& c = *i;
2697 PPL_ASSERT(!c.is_tautological());
2698 PPL_ASSERT(!c.is_inconsistent());
2699 // If the polyhedron `x' is included in the polyhedron defined by
2700 // `c', then `c' can be skipped, as adding its complement to `x'
2701 // would result in the empty polyhedron. Moreover, if we operate
2702 // on C-polyhedra and `c' is a non-strict inequality, c _must_ be
2703 // skipped for otherwise we would obtain a result that is less
2704 // precise than the poly-difference.
2705 if (x.relation_with(c).implies(Poly_Con_Relation::is_included())) {
2706 continue;
2707 }
2708 Polyhedron z = x;
2709 const Linear_Expression e(c.expression());
2710 switch (c.type()) {
2711 case Constraint::NONSTRICT_INEQUALITY:
2712 if (is_necessarily_closed()) {
2713 z.refine_no_check(e <= 0);
2714 }
2715 else {
2716 z.refine_no_check(e < 0);
2717 }
2718 break;
2719 case Constraint::STRICT_INEQUALITY:
2720 z.refine_no_check(e <= 0);
2721 break;
2722 case Constraint::EQUALITY:
2723 if (is_necessarily_closed()) {
2724 // We have already filtered out the case
2725 // when `x' is included in `y': the result is `x'.
2726 return;
2727 }
2728 else {
2729 Polyhedron w = x;
2730 w.refine_no_check(e < 0);
2731 new_polyhedron.poly_hull_assign(w);
2732 z.refine_no_check(e > 0);
2733 }
2734 break;
2735 }
2736 new_polyhedron.poly_hull_assign(z);
2737 }
2738 *this = new_polyhedron;
2739
2740 PPL_ASSERT_HEAVY(OK());
2741 }
2742
2743 void
2744 PPL::Polyhedron::
affine_image(const Variable var,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)2745 affine_image(const Variable var,
2746 const Linear_Expression& expr,
2747 Coefficient_traits::const_reference denominator) {
2748 // The denominator cannot be zero.
2749 if (denominator == 0) {
2750 throw_invalid_argument("affine_image(v, e, d)", "d == 0");
2751 }
2752 // Dimension-compatibility checks.
2753 // The dimension of `expr' should not be greater than the dimension
2754 // of `*this'.
2755 if (space_dim < expr.space_dimension()) {
2756 throw_dimension_incompatible("affine_image(v, e, d)", "e", expr);
2757 // `var' should be one of the dimensions of the polyhedron.
2758 }
2759 const dimension_type var_space_dim = var.space_dimension();
2760 if (space_dim < var_space_dim) {
2761 throw_dimension_incompatible("affine_image(v, e, d)", "v", var);
2762 }
2763 if (marked_empty()) {
2764 return;
2765 }
2766
2767 if (expr.coefficient(var) != 0) {
2768 // The transformation is invertible:
2769 // minimality and saturators are preserved, so that
2770 // pending rows, if present, are correctly handled.
2771 if (generators_are_up_to_date()) {
2772 // Generator_System::affine_image() requires the third argument
2773 // to be a positive Coefficient.
2774 if (denominator > 0) {
2775 gen_sys.affine_image(var, expr, denominator);
2776 }
2777 else {
2778 gen_sys.affine_image(var, -expr, -denominator);
2779 }
2780 }
2781 if (constraints_are_up_to_date()) {
2782 // To build the inverse transformation,
2783 // after copying and negating `expr',
2784 // we exchange the roles of `expr[var_space_dim]' and `denominator'.
2785 Linear_Expression inverse;
2786 Coefficient_traits::const_reference c = expr.coefficient(var);
2787 if (c > 0) {
2788 inverse = -expr;
2789 inverse.set_coefficient(var, denominator);
2790 con_sys.affine_preimage(var, inverse, c);
2791 }
2792 else {
2793 // The new denominator is negative: we negate everything once
2794 // more, as Constraint_System::affine_preimage() requires the
2795 // third argument to be positive.
2796 inverse = expr;
2797 inverse.set_coefficient(var, -denominator);
2798 con_sys.affine_preimage(var, inverse, -c);
2799 }
2800 }
2801 }
2802 else {
2803 // The transformation is not invertible.
2804 // We need an up-to-date system of generators.
2805 if (has_something_pending()) {
2806 remove_pending_to_obtain_generators();
2807 }
2808 else if (!generators_are_up_to_date()) {
2809 minimize();
2810 }
2811 if (!marked_empty()) {
2812 // Generator_System::affine_image() requires the third argument
2813 // to be a positive Coefficient.
2814 if (denominator > 0) {
2815 gen_sys.affine_image(var, expr, denominator);
2816 }
2817 else {
2818 gen_sys.affine_image(var, -expr, -denominator);
2819 }
2820
2821 clear_constraints_up_to_date();
2822 clear_generators_minimized();
2823 clear_sat_c_up_to_date();
2824 clear_sat_g_up_to_date();
2825 }
2826 }
2827 PPL_ASSERT_HEAVY(OK());
2828 }
2829
2830
2831 void
2832 PPL::Polyhedron::
affine_preimage(const Variable var,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)2833 affine_preimage(const Variable var,
2834 const Linear_Expression& expr,
2835 Coefficient_traits::const_reference denominator) {
2836 // The denominator cannot be zero.
2837 if (denominator == 0) {
2838 throw_invalid_argument("affine_preimage(v, e, d)", "d == 0");
2839 }
2840
2841 // Dimension-compatibility checks.
2842 // The dimension of `expr' should not be greater than the dimension
2843 // of `*this'.
2844 if (space_dim < expr.space_dimension()) {
2845 throw_dimension_incompatible("affine_preimage(v, e, d)", "e", expr);
2846 // `var' should be one of the dimensions of the polyhedron.
2847 }
2848 const dimension_type var_space_dim = var.space_dimension();
2849 if (space_dim < var_space_dim) {
2850 throw_dimension_incompatible("affine_preimage(v, e, d)", "v", var);
2851 }
2852
2853 if (marked_empty()) {
2854 return;
2855 }
2856
2857 if (expr.coefficient(var) != 0) {
2858 // The transformation is invertible:
2859 // minimality and saturators are preserved.
2860 if (constraints_are_up_to_date()) {
2861 // Constraint_System::affine_preimage() requires the third argument
2862 // to be a positive Coefficient.
2863 if (denominator > 0) {
2864 con_sys.affine_preimage(var, expr, denominator);
2865 }
2866 else {
2867 con_sys.affine_preimage(var, -expr, -denominator);
2868 }
2869 }
2870 if (generators_are_up_to_date()) {
2871 // To build the inverse transformation,
2872 // after copying and negating `expr',
2873 // we exchange the roles of `expr[var_space_dim]' and `denominator'.
2874 Linear_Expression inverse;
2875 Coefficient_traits::const_reference c = expr.coefficient(var);
2876 if (c > 0) {
2877 inverse = -expr;
2878 inverse.set_coefficient(var, denominator);
2879 gen_sys.affine_image(var, inverse, c);
2880 }
2881 else {
2882 // The new denominator is negative:
2883 // we negate everything once more, as Generator_System::affine_image()
2884 // requires the third argument to be positive.
2885 inverse = expr;
2886 inverse.set_coefficient(var, -denominator);
2887 gen_sys.affine_image(var, inverse, -c);
2888 }
2889 }
2890 }
2891 else {
2892 // The transformation is not invertible.
2893 // We need an up-to-date system of constraints.
2894 if (has_something_pending()) {
2895 remove_pending_to_obtain_constraints();
2896 }
2897 else if (!constraints_are_up_to_date()) {
2898 minimize();
2899 }
2900 // Constraint_System::affine_preimage() requires the third argument
2901 // to be a positive Coefficient.
2902 if (denominator > 0) {
2903 con_sys.affine_preimage(var, expr, denominator);
2904 }
2905 else {
2906 con_sys.affine_preimage(var, -expr, -denominator);
2907 }
2908 // Generators, minimality and saturators are no longer valid.
2909 clear_generators_up_to_date();
2910 clear_constraints_minimized();
2911 clear_sat_c_up_to_date();
2912 clear_sat_g_up_to_date();
2913 }
2914 PPL_ASSERT_HEAVY(OK());
2915 }
2916
2917 void
2918 PPL::Polyhedron::
bounded_affine_image(const Variable var,const Linear_Expression & lb_expr,const Linear_Expression & ub_expr,Coefficient_traits::const_reference denominator)2919 bounded_affine_image(const Variable var,
2920 const Linear_Expression& lb_expr,
2921 const Linear_Expression& ub_expr,
2922 Coefficient_traits::const_reference denominator) {
2923 // The denominator cannot be zero.
2924 if (denominator == 0) {
2925 throw_invalid_argument("bounded_affine_image(v, lb, ub, d)", "d == 0");
2926 }
2927
2928 // Dimension-compatibility checks.
2929 // `var' should be one of the dimensions of the polyhedron.
2930 const dimension_type var_space_dim = var.space_dimension();
2931 if (space_dim < var_space_dim) {
2932 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
2933 "v", var);
2934 }
2935 // The dimension of `lb_expr' and `ub_expr' should not be
2936 // greater than the dimension of `*this'.
2937 const dimension_type lb_space_dim = lb_expr.space_dimension();
2938 if (space_dim < lb_space_dim) {
2939 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
2940 "lb", lb_expr);
2941 }
2942 const dimension_type ub_space_dim = ub_expr.space_dimension();
2943 if (space_dim < ub_space_dim) {
2944 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
2945 "ub", ub_expr);
2946 }
2947
2948 // Any image of an empty polyhedron is empty.
2949 if (marked_empty()) {
2950 return;
2951 }
2952
2953 // Check whether `var' occurs in `lb_expr' and/or `ub_expr'.
2954 if (lb_expr.coefficient(var) == 0) {
2955 // Here `var' may only occur in `ub_expr'.
2956 generalized_affine_image(var,
2957 LESS_OR_EQUAL,
2958 ub_expr,
2959 denominator);
2960 if (denominator > 0) {
2961 refine_no_check(lb_expr <= denominator*var);
2962 }
2963 else {
2964 refine_no_check(denominator*var <= lb_expr);
2965 }
2966 }
2967 else if (ub_expr.coefficient(var) == 0) {
2968 // Here `var' only occurs in `lb_expr'.
2969 generalized_affine_image(var,
2970 GREATER_OR_EQUAL,
2971 lb_expr,
2972 denominator);
2973 if (denominator > 0) {
2974 refine_no_check(denominator*var <= ub_expr);
2975 }
2976 else {
2977 refine_no_check(ub_expr <= denominator*var);
2978 }
2979 }
2980 else {
2981 // Here `var' occurs in both `lb_expr' and `ub_expr'.
2982 // To ease the computation, we add an additional dimension.
2983 const Variable new_var(space_dim);
2984 add_space_dimensions_and_embed(1);
2985 // Constrain the new dimension to be equal to `ub_expr'.
2986 refine_no_check(denominator*new_var == ub_expr);
2987 // Apply the affine lower bound.
2988 generalized_affine_image(var,
2989 GREATER_OR_EQUAL,
2990 lb_expr,
2991 denominator);
2992 if (!marked_empty()) {
2993 // Now apply the affine upper bound, as recorded in `new_var'.
2994 refine_no_check(new_var >= var);
2995 }
2996 // Remove the temporarily added dimension.
2997 remove_higher_space_dimensions(space_dim-1);
2998 }
2999 PPL_ASSERT_HEAVY(OK());
3000 }
3001
3002 void
3003 PPL::Polyhedron::
bounded_affine_preimage(const Variable var,const Linear_Expression & lb_expr,const Linear_Expression & ub_expr,Coefficient_traits::const_reference denominator)3004 bounded_affine_preimage(const Variable var,
3005 const Linear_Expression& lb_expr,
3006 const Linear_Expression& ub_expr,
3007 Coefficient_traits::const_reference denominator) {
3008 // The denominator cannot be zero.
3009 if (denominator == 0) {
3010 throw_invalid_argument("bounded_affine_preimage(v, lb, ub, d)", "d == 0");
3011 }
3012
3013 // Dimension-compatibility checks.
3014 // `var' should be one of the dimensions of the polyhedron.
3015 const dimension_type var_space_dim = var.space_dimension();
3016 if (space_dim < var_space_dim) {
3017 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
3018 "v", var);
3019 }
3020 // The dimension of `lb_expr' and `ub_expr' should not be
3021 // greater than the dimension of `*this'.
3022 const dimension_type lb_space_dim = lb_expr.space_dimension();
3023 if (space_dim < lb_space_dim) {
3024 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
3025 "lb", lb_expr);
3026 }
3027 const dimension_type ub_space_dim = ub_expr.space_dimension();
3028 if (space_dim < ub_space_dim) {
3029 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
3030 "ub", ub_expr);
3031 }
3032
3033 // Any preimage of an empty polyhedron is empty.
3034 if (marked_empty()) {
3035 return;
3036 }
3037
3038 // Check whether `var' occurs in neither `lb_expr' nor `ub_expr'.
3039 if (lb_expr.coefficient(var) == 0 && ub_expr.coefficient(var) == 0) {
3040 if (denominator > 0) {
3041 refine_no_check(lb_expr <= denominator*var);
3042 refine_no_check(denominator*var <= ub_expr);
3043 }
3044 else {
3045 refine_no_check(ub_expr <= denominator*var);
3046 refine_no_check(denominator*var <= lb_expr);
3047 }
3048 unconstrain(var);
3049 }
3050 else {
3051 // Here `var' occurs in `lb_expr' or `ub_expr'.
3052 // To ease the computation, add an additional dimension.
3053 const Variable new_var(space_dim);
3054 add_space_dimensions_and_embed(1);
3055
3056 // Swap dimensions `var' and `new_var'.
3057 if (constraints_are_up_to_date()) {
3058 con_sys.swap_space_dimensions(var, new_var);
3059 }
3060 if (generators_are_up_to_date()) {
3061 gen_sys.swap_space_dimensions(var, new_var);
3062 }
3063
3064 // Constrain the new dimension as dictated by `lb_expr' and `ub_expr'.
3065 // (we force minimization because we will need the generators).
3066 if (denominator > 0) {
3067 refine_no_check(lb_expr <= denominator*new_var);
3068 refine_no_check(denominator*new_var <= ub_expr);
3069 }
3070 else {
3071 refine_no_check(ub_expr <= denominator*new_var);
3072 refine_no_check(denominator*new_var <= lb_expr);
3073 }
3074 // Remove the temporarily added dimension.
3075 remove_higher_space_dimensions(space_dim-1);
3076 }
3077 PPL_ASSERT_HEAVY(OK());
3078 }
3079
3080 void
3081 PPL::Polyhedron::
generalized_affine_image(const Variable var,const Relation_Symbol relsym,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)3082 generalized_affine_image(const Variable var,
3083 const Relation_Symbol relsym,
3084 const Linear_Expression& expr,
3085 Coefficient_traits::const_reference denominator) {
3086 // The denominator cannot be zero.
3087 if (denominator == 0) {
3088 throw_invalid_argument("generalized_affine_image(v, r, e, d)", "d == 0");
3089 }
3090
3091 // Dimension-compatibility checks.
3092 // The dimension of `expr' should not be greater than the dimension
3093 // of `*this'.
3094 if (space_dim < expr.space_dimension()) {
3095 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
3096 "e", expr);
3097 }
3098 // `var' should be one of the dimensions of the polyhedron.
3099 const dimension_type var_space_dim = var.space_dimension();
3100 if (space_dim < var_space_dim) {
3101 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
3102 "v", var);
3103 }
3104
3105 // Strict relation symbols are only admitted for NNC polyhedra.
3106 if (is_necessarily_closed()
3107 && (relsym == LESS_THAN || relsym == GREATER_THAN)) {
3108 throw_invalid_argument("generalized_affine_image(v, r, e, d)",
3109 "r is a strict relation symbol");
3110 }
3111 // The relation symbol cannot be a disequality.
3112 if (relsym == NOT_EQUAL) {
3113 throw_invalid_argument("generalized_affine_image(v, r, e, d)",
3114 "r is the disequality relation symbol");
3115 }
3116
3117 // First compute the affine image.
3118 affine_image(var, expr, denominator);
3119
3120 if (relsym == EQUAL) {
3121 // The affine relation is indeed an affine function.
3122 return;
3123 }
3124 // Any image of an empty polyhedron is empty.
3125 // Note: DO check for emptiness here, as we will later add a ray.
3126 if (is_empty()) {
3127 return;
3128 }
3129
3130 switch (relsym) {
3131 case LESS_OR_EQUAL:
3132 add_generator(ray(-var));
3133 break;
3134 case GREATER_OR_EQUAL:
3135 add_generator(ray(var));
3136 break;
3137 case LESS_THAN:
3138 // Intentionally fall through.
3139 case GREATER_THAN:
3140 {
3141 // The relation symbol is strict.
3142 PPL_ASSERT(!is_necessarily_closed());
3143 // While adding the ray, we minimize the generators
3144 // in order to avoid adding too many redundant generators later.
3145 add_generator(ray((relsym == GREATER_THAN) ? var : -var));
3146 minimize();
3147
3148 // We split each point of the generator system into two generators:
3149 // a closure point, having the same coordinates of the given point,
3150 // and another point, having the same coordinates for all but the
3151 // `var' dimension, which is displaced along the direction of the
3152 // newly introduced ray.
3153 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
3154 const Generator& gen_i = gen_sys.sys.rows[i];
3155 if (gen_i.is_point()) {
3156 // Add a copy of `gen_i' at the end of the system.
3157 // Note: copying is really meant, to avoid undefined behavior.
3158 gen_sys.sys.rows.push_back(Generator(gen_i));
3159 // Note: (re-)compute references (invalidated by push_back).
3160 Generator& old_gen = gen_sys.sys.rows[i];
3161 Generator& new_gen = gen_sys.sys.rows.back();
3162 // Transform `old_gen' into a closure point.
3163 old_gen.set_epsilon_coefficient(0);
3164 old_gen.expr.normalize();
3165 PPL_ASSERT(old_gen.OK());
3166 // Displace `new_gen' by `var' (i.e., along the new ray).
3167 if (relsym == GREATER_THAN) {
3168 new_gen.expr += var;
3169 }
3170 else {
3171 new_gen.expr -= var;
3172 }
3173 new_gen.expr.normalize();
3174 PPL_ASSERT(new_gen.OK());
3175 }
3176 }
3177 // Sortedness no longer hold.
3178 gen_sys.set_sorted(false);
3179 gen_sys.unset_pending_rows();
3180 PPL_ASSERT(gen_sys.sys.OK());
3181
3182 clear_constraints_up_to_date();
3183 clear_generators_minimized();
3184 clear_sat_c_up_to_date();
3185 clear_sat_g_up_to_date();
3186 }
3187 break;
3188 default:
3189 // The EQUAL and NOT_EQUAL cases have been already dealt with.
3190 PPL_UNREACHABLE;
3191 break;
3192 }
3193 PPL_ASSERT_HEAVY(OK());
3194 }
3195
3196 void
3197 PPL::Polyhedron::
generalized_affine_preimage(const Variable var,const Relation_Symbol relsym,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)3198 generalized_affine_preimage(const Variable var,
3199 const Relation_Symbol relsym,
3200 const Linear_Expression& expr,
3201 Coefficient_traits::const_reference denominator) {
3202 // The denominator cannot be zero.
3203 if (denominator == 0) {
3204 throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
3205 "d == 0");
3206 }
3207
3208 // Dimension-compatibility checks.
3209 // The dimension of `expr' should not be greater than the dimension
3210 // of `*this'.
3211 if (space_dim < expr.space_dimension()) {
3212 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
3213 "e", expr);
3214 }
3215 // `var' should be one of the dimensions of the polyhedron.
3216 const dimension_type var_space_dim = var.space_dimension();
3217 if (space_dim < var_space_dim) {
3218 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
3219 "v", var);
3220 }
3221
3222 // Strict relation symbols are only admitted for NNC polyhedra.
3223 if (is_necessarily_closed()
3224 && (relsym == LESS_THAN || relsym == GREATER_THAN)) {
3225 throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
3226 "r is a strict relation symbol");
3227 }
3228 // The relation symbol cannot be a disequality.
3229 if (relsym == NOT_EQUAL) {
3230 throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
3231 "r is the disequality relation symbol");
3232 }
3233
3234 // Check whether the affine relation is indeed an affine function.
3235 if (relsym == EQUAL) {
3236 affine_preimage(var, expr, denominator);
3237 return;
3238 }
3239
3240 // Compute the reversed relation symbol to simplify later coding.
3241 Relation_Symbol reversed_relsym;
3242 switch (relsym) {
3243 case LESS_THAN:
3244 reversed_relsym = GREATER_THAN;
3245 break;
3246 case LESS_OR_EQUAL:
3247 reversed_relsym = GREATER_OR_EQUAL;
3248 break;
3249 case GREATER_OR_EQUAL:
3250 reversed_relsym = LESS_OR_EQUAL;
3251 break;
3252 case GREATER_THAN:
3253 reversed_relsym = LESS_THAN;
3254 break;
3255 default:
3256 // The EQUAL and NOT_EQUAL cases have been already dealt with.
3257 PPL_UNREACHABLE;
3258 return;
3259 }
3260
3261 // Check whether the preimage of this affine relation can be easily
3262 // computed as the image of its inverse relation.
3263 const Coefficient& var_coefficient = expr.coefficient(var);
3264 if (var_coefficient != 0) {
3265 const Linear_Expression inverse_expr
3266 = expr - (denominator + var_coefficient) * var;
3267 PPL_DIRTY_TEMP_COEFFICIENT(inverse_denominator);
3268 neg_assign(inverse_denominator, var_coefficient);
3269 const Relation_Symbol inverse_relsym
3270 = (sgn(denominator) == sgn(inverse_denominator))
3271 ? relsym : reversed_relsym;
3272 generalized_affine_image(var, inverse_relsym, inverse_expr,
3273 inverse_denominator);
3274 return;
3275 }
3276
3277 // Here `var_coefficient == 0', so that the preimage cannot
3278 // be easily computed by inverting the affine relation.
3279 // Shrink the polyhedron by adding the constraint induced
3280 // by the affine relation.
3281 const Relation_Symbol corrected_relsym
3282 = (denominator > 0) ? relsym : reversed_relsym;
3283 switch (corrected_relsym) {
3284 case LESS_THAN:
3285 refine_no_check(denominator*var < expr);
3286 break;
3287 case LESS_OR_EQUAL:
3288 refine_no_check(denominator*var <= expr);
3289 break;
3290 case GREATER_OR_EQUAL:
3291 refine_no_check(denominator*var >= expr);
3292 break;
3293 case GREATER_THAN:
3294 refine_no_check(denominator*var > expr);
3295 break;
3296 default:
3297 // The EQUAL and NOT_EQUAL cases have been already dealt with.
3298 PPL_UNREACHABLE;
3299 break;
3300 }
3301 unconstrain(var);
3302 PPL_ASSERT_HEAVY(OK());
3303 }
3304
3305 void
generalized_affine_image(const Linear_Expression & lhs,const Relation_Symbol relsym,const Linear_Expression & rhs)3306 PPL::Polyhedron::generalized_affine_image(const Linear_Expression& lhs,
3307 const Relation_Symbol relsym,
3308 const Linear_Expression& rhs) {
3309 // Dimension-compatibility checks.
3310 // The dimension of `lhs' should not be greater than the dimension
3311 // of `*this'.
3312 dimension_type lhs_space_dim = lhs.space_dimension();
3313 if (space_dim < lhs_space_dim) {
3314 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
3315 "e1", lhs);
3316 }
3317 // The dimension of `rhs' should not be greater than the dimension
3318 // of `*this'.
3319 const dimension_type rhs_space_dim = rhs.space_dimension();
3320 if (space_dim < rhs_space_dim) {
3321 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
3322 "e2", rhs);
3323 }
3324
3325 // Strict relation symbols are only admitted for NNC polyhedra.
3326 if (is_necessarily_closed()
3327 && (relsym == LESS_THAN || relsym == GREATER_THAN)) {
3328 throw_invalid_argument("generalized_affine_image(e1, r, e2)",
3329 "r is a strict relation symbol");
3330 }
3331 // The relation symbol cannot be a disequality.
3332 if (relsym == NOT_EQUAL) {
3333 throw_invalid_argument("generalized_affine_image(e1, r, e2)",
3334 "r is the disequality relation symbol");
3335 }
3336 // Any image of an empty polyhedron is empty.
3337 if (marked_empty()) {
3338 return;
3339 }
3340
3341 // Compute the actual space dimension of `lhs',
3342 // i.e., the highest dimension having a non-zero coefficient in `lhs'.
3343 lhs_space_dim = lhs.last_nonzero();
3344
3345 // If all variables have a zero coefficient, then `lhs' is a constant:
3346 // we can simply add the constraint `lhs relsym rhs'.
3347 if (lhs_space_dim == 0) {
3348 switch (relsym) {
3349 case LESS_THAN:
3350 refine_no_check(lhs < rhs);
3351 break;
3352 case LESS_OR_EQUAL:
3353 refine_no_check(lhs <= rhs);
3354 break;
3355 case EQUAL:
3356 refine_no_check(lhs == rhs);
3357 break;
3358 case GREATER_OR_EQUAL:
3359 refine_no_check(lhs >= rhs);
3360 break;
3361 case GREATER_THAN:
3362 refine_no_check(lhs > rhs);
3363 break;
3364 case NOT_EQUAL:
3365 // The NOT_EQUAL case has been already dealt with.
3366 PPL_UNREACHABLE;
3367 break;
3368 }
3369 return;
3370 }
3371
3372 // Gather in `new_lines' the collections of all the lines having
3373 // the direction of variables occurring in `lhs'.
3374 // While at it, check whether or not there exists a variable
3375 // occurring in both `lhs' and `rhs'.
3376 Generator_System new_lines;
3377 for (Linear_Expression::const_iterator
3378 i = lhs.begin(), i_end = lhs.end(); i != i_end; ++i) {
3379 new_lines.insert(line(i.variable()));
3380 }
3381
3382 const dimension_type num_common_dims
3383 = std::min(lhs.space_dimension(), rhs.space_dimension());
3384 if (lhs.have_a_common_variable(rhs, Variable(0), Variable(num_common_dims))) {
3385 // Some variables in `lhs' also occur in `rhs'.
3386 // To ease the computation, we add an additional dimension.
3387 const Variable new_var(space_dim);
3388 add_space_dimensions_and_embed(1);
3389
3390 // Constrain the new dimension to be equal to the right hand side.
3391 // (check for emptiness because we will add lines).
3392 refine_no_check(new_var == rhs);
3393 if (!is_empty()) {
3394 // Existentially quantify the variables in the left hand side.
3395 add_recycled_generators(new_lines);
3396
3397 // Constrain the new dimension so that it is related to
3398 // the left hand side as dictated by `relsym'
3399 // (we force minimization because we will need the generators).
3400 switch (relsym) {
3401 case LESS_THAN:
3402 refine_no_check(lhs < new_var);
3403 break;
3404 case LESS_OR_EQUAL:
3405 refine_no_check(lhs <= new_var);
3406 break;
3407 case EQUAL:
3408 refine_no_check(lhs == new_var);
3409 break;
3410 case GREATER_OR_EQUAL:
3411 refine_no_check(lhs >= new_var);
3412 break;
3413 case GREATER_THAN:
3414 refine_no_check(lhs > new_var);
3415 break;
3416 case NOT_EQUAL:
3417 // The NOT_EQUAL case has been already dealt with.
3418 PPL_UNREACHABLE;
3419 break;
3420 }
3421 }
3422 // Remove the temporarily added dimension.
3423 remove_higher_space_dimensions(space_dim-1);
3424 }
3425 else {
3426 // `lhs' and `rhs' variables are disjoint:
3427 // there is no need to add a further dimension.
3428
3429 // Any image of an empty polyhedron is empty.
3430 // Note: DO check for emptiness here, as we will add lines.
3431 if (is_empty()) {
3432 return;
3433 }
3434
3435 // Existentially quantify the variables in the left hand side.
3436 add_recycled_generators(new_lines);
3437
3438 // Constrain the left hand side expression so that it is related to
3439 // the right hand side expression as dictated by `relsym'.
3440 switch (relsym) {
3441 case LESS_THAN:
3442 refine_no_check(lhs < rhs);
3443 break;
3444 case LESS_OR_EQUAL:
3445 refine_no_check(lhs <= rhs);
3446 break;
3447 case EQUAL:
3448 refine_no_check(lhs == rhs);
3449 break;
3450 case GREATER_OR_EQUAL:
3451 refine_no_check(lhs >= rhs);
3452 break;
3453 case GREATER_THAN:
3454 refine_no_check(lhs > rhs);
3455 break;
3456 case NOT_EQUAL:
3457 // The NOT_EQUAL case has been already dealt with.
3458 PPL_UNREACHABLE;
3459 break;
3460 }
3461 }
3462 PPL_ASSERT_HEAVY(OK());
3463 }
3464
3465 void
generalized_affine_preimage(const Linear_Expression & lhs,const Relation_Symbol relsym,const Linear_Expression & rhs)3466 PPL::Polyhedron::generalized_affine_preimage(const Linear_Expression& lhs,
3467 const Relation_Symbol relsym,
3468 const Linear_Expression& rhs) {
3469 // Dimension-compatibility checks.
3470 // The dimension of `lhs' should not be greater than the dimension
3471 // of `*this'.
3472 dimension_type lhs_space_dim = lhs.space_dimension();
3473 if (space_dim < lhs_space_dim) {
3474 throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
3475 "e1", lhs);
3476 }
3477 // The dimension of `rhs' should not be greater than the dimension
3478 // of `*this'.
3479 const dimension_type rhs_space_dim = rhs.space_dimension();
3480 if (space_dim < rhs_space_dim) {
3481 throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
3482 "e2", rhs);
3483 }
3484
3485 // Strict relation symbols are only admitted for NNC polyhedra.
3486 if (is_necessarily_closed()
3487 && (relsym == LESS_THAN || relsym == GREATER_THAN)) {
3488 throw_invalid_argument("generalized_affine_preimage(e1, r, e2)",
3489 "r is a strict relation symbol");
3490 }
3491 // The relation symbol cannot be a disequality.
3492 if (relsym == NOT_EQUAL) {
3493 throw_invalid_argument("generalized_affine_preimage(e1, r, e2)",
3494 "r is the disequality relation symbol");
3495 }
3496 // Any preimage of an empty polyhedron is empty.
3497 if (marked_empty()) {
3498 return;
3499 }
3500
3501 // Compute the actual space dimension of `lhs',
3502 // i.e., the highest dimension having a non-zero coefficient in `lhs'.
3503 lhs_space_dim = lhs.last_nonzero();
3504
3505 // If all variables have a zero coefficient, then `lhs' is a constant:
3506 // in this case, preimage and image happen to be the same.
3507 if (lhs_space_dim == 0) {
3508 generalized_affine_image(lhs, relsym, rhs);
3509 return;
3510 }
3511
3512 // Gather in `new_lines' the collections of all the lines having
3513 // the direction of variables occurring in `lhs'.
3514 // While at it, check whether or not there exists a variable
3515 // occurring in both `lhs' and `rhs'.
3516 Generator_System new_lines;
3517 for (Linear_Expression::const_iterator
3518 i = lhs.begin(), i_end = lhs.end(); i != i_end; ++i) {
3519 new_lines.insert(line(i.variable()));
3520 }
3521
3522 const dimension_type num_common_dims
3523 = std::min(lhs.space_dimension(), rhs.space_dimension());
3524 if (lhs.have_a_common_variable(rhs, Variable(0), Variable(num_common_dims))) {
3525 // Some variables in `lhs' also occur in `rhs'.
3526 // To ease the computation, we add an additional dimension.
3527 const Variable new_var(space_dim);
3528 add_space_dimensions_and_embed(1);
3529
3530 // Constrain the new dimension to be equal to `lhs'
3531 // (also check for emptiness because we have to add lines).
3532 refine_no_check(new_var == lhs);
3533 if (!is_empty()) {
3534 // Existentially quantify the variables in the left hand side.
3535 add_recycled_generators(new_lines);
3536
3537 // Constrain the new dimension so that it is related to
3538 // the right hand side as dictated by `relsym'.
3539 switch (relsym) {
3540 case LESS_THAN:
3541 refine_no_check(new_var < rhs);
3542 break;
3543 case LESS_OR_EQUAL:
3544 refine_no_check(new_var <= rhs);
3545 break;
3546 case EQUAL:
3547 refine_no_check(new_var == rhs);
3548 break;
3549 case GREATER_OR_EQUAL:
3550 refine_no_check(new_var >= rhs);
3551 break;
3552 case GREATER_THAN:
3553 refine_no_check(new_var > rhs);
3554 break;
3555 case NOT_EQUAL:
3556 // The NOT_EQUAL case has been already dealt with.
3557 PPL_UNREACHABLE;
3558 break;
3559 }
3560 }
3561 // Remove the temporarily added dimension.
3562 remove_higher_space_dimensions(space_dim-1);
3563 }
3564 else {
3565 // `lhs' and `rhs' variables are disjoint:
3566 // there is no need to add a further dimension.
3567
3568 // Constrain the left hand side expression so that it is related to
3569 // the right hand side expression as dictated by `relsym'.
3570 switch (relsym) {
3571 case LESS_THAN:
3572 refine_no_check(lhs < rhs);
3573 break;
3574 case LESS_OR_EQUAL:
3575 refine_no_check(lhs <= rhs);
3576 break;
3577 case EQUAL:
3578 refine_no_check(lhs == rhs);
3579 break;
3580 case GREATER_OR_EQUAL:
3581 refine_no_check(lhs >= rhs);
3582 break;
3583 case GREATER_THAN:
3584 refine_no_check(lhs > rhs);
3585 break;
3586 case NOT_EQUAL:
3587 // The NOT_EQUAL case has been already dealt with.
3588 PPL_UNREACHABLE;
3589 break;
3590 }
3591 // Any image of an empty polyhedron is empty.
3592 // Note: DO check for emptiness here, as we will add lines.
3593 if (is_empty()) {
3594 return;
3595 }
3596 // Existentially quantify all the variables occurring in `lhs'.
3597 add_recycled_generators(new_lines);
3598 }
3599 PPL_ASSERT_HEAVY(OK());
3600 }
3601
3602 void
time_elapse_assign(const Polyhedron & y)3603 PPL::Polyhedron::time_elapse_assign(const Polyhedron& y) {
3604 Polyhedron& x = *this;
3605 // Topology compatibility check.
3606 if (x.topology() != y.topology()) {
3607 throw_topology_incompatible("time_elapse_assign(y)", "y", y);
3608 }
3609 // Dimension-compatibility checks.
3610 if (x.space_dim != y.space_dim) {
3611 throw_dimension_incompatible("time_elapse_assign(y)", "y", y);
3612 }
3613
3614 // Dealing with the zero-dimensional case.
3615 if (x.space_dim == 0) {
3616 if (y.marked_empty()) {
3617 x.set_empty();
3618 }
3619 return;
3620 }
3621
3622 // If either one of `x' or `y' is empty, the result is empty too.
3623 if (x.marked_empty() || y.marked_empty()
3624 || (x.has_pending_constraints() && !x.process_pending_constraints())
3625 || (!x.generators_are_up_to_date() && !x.update_generators())
3626 || (y.has_pending_constraints() && !y.process_pending_constraints())
3627 || (!y.generators_are_up_to_date() && !y.update_generators())) {
3628 x.set_empty();
3629 return;
3630 }
3631
3632 // At this point both generator systems are up-to-date,
3633 // possibly containing pending generators.
3634 Generator_System gs = y.gen_sys;
3635 const dimension_type old_gs_num_rows = gs.num_rows();
3636 dimension_type gs_num_rows = old_gs_num_rows;
3637
3638 if (!x.is_necessarily_closed()) {
3639 // `x' and `y' are NNC polyhedra.
3640 for (dimension_type i = gs_num_rows; i-- > 0; ) {
3641 Generator& g = gs.sys.rows[i];
3642 switch (g.type()) {
3643 case Generator::POINT:
3644 // The points of `gs' can be erased,
3645 // since their role can be played by closure points.
3646 --gs_num_rows;
3647 swap(g, gs.sys.rows[gs_num_rows]);
3648 break;
3649 case Generator::CLOSURE_POINT:
3650 {
3651 // If it is the origin, erase it.
3652 if (g.expr.all_homogeneous_terms_are_zero()) {
3653 --gs_num_rows;
3654 swap(g, gs.sys.rows[gs_num_rows]);
3655 }
3656 // Otherwise, transform the closure point into a ray.
3657 else {
3658 g.expr.set_inhomogeneous_term(0);
3659 // Enforce normalization.
3660 g.expr.normalize();
3661 PPL_ASSERT(g.OK());
3662 }
3663 }
3664 break;
3665 case Generator::RAY:
3666 case Generator::LINE:
3667 // For rays and lines, nothing to be done.
3668 break;
3669 }
3670 }
3671 }
3672 else {
3673 // `x' and `y' are C polyhedra.
3674 for (dimension_type i = gs_num_rows; i-- > 0; ) {
3675 // For rays and lines, nothing to be done.
3676 if (gs[i].is_point()) {
3677 Generator& p = gs.sys.rows[i];
3678 // If it is the origin, erase it.
3679 if (p.expression().all_homogeneous_terms_are_zero()) {
3680 --gs_num_rows;
3681 swap(p, gs.sys.rows[gs_num_rows]);
3682 }
3683 // Otherwise, transform the point into a ray.
3684 else {
3685 p.expr.set_inhomogeneous_term(0);
3686 // Enforce normalization.
3687 p.expr.normalize();
3688 PPL_ASSERT(p.OK());
3689 }
3690 }
3691 }
3692 }
3693 // If it was present, erase the origin point or closure point,
3694 // which cannot be transformed into a valid ray or line.
3695 // For NNC polyhedra, also erase all the points of `gs',
3696 // whose role can be played by the closure points.
3697 // These have been previously moved to the end of `gs'.
3698 gs.sys.rows.resize(gs_num_rows);
3699
3700 gs.unset_pending_rows();
3701 PPL_ASSERT(gs.sys.OK());
3702
3703 // `gs' may now have no rows.
3704 // Namely, this happens when `y' was the singleton polyhedron
3705 // having the origin as the one and only point.
3706 // In such a case, the resulting polyhedron is equal to `x'.
3707 if (gs_num_rows == 0) {
3708 return;
3709 }
3710 // If the polyhedron can have something pending, we add `gs'
3711 // to `gen_sys' as pending rows
3712 if (x.can_have_something_pending()) {
3713 x.gen_sys.insert_pending(gs);
3714 x.set_generators_pending();
3715 }
3716 // Otherwise, the two systems are merged.
3717 // `Linear_System::merge_rows_assign()' requires both systems to be sorted.
3718 else {
3719 if (!x.gen_sys.is_sorted()) {
3720 x.gen_sys.sort_rows();
3721 }
3722 gs.sort_rows();
3723 x.gen_sys.merge_rows_assign(gs);
3724 // Only the system of generators is up-to-date.
3725 x.clear_constraints_up_to_date();
3726 x.clear_generators_minimized();
3727 }
3728 PPL_ASSERT_HEAVY(x.OK(true) && y.OK(true));
3729 }
3730
3731 bool
frequency(const Linear_Expression & expr,Coefficient & freq_n,Coefficient & freq_d,Coefficient & val_n,Coefficient & val_d) const3732 PPL::Polyhedron::frequency(const Linear_Expression& expr,
3733 Coefficient& freq_n, Coefficient& freq_d,
3734 Coefficient& val_n, Coefficient& val_d) const {
3735 // The dimension of `expr' must be at most the dimension of *this.
3736 if (space_dim < expr.space_dimension()) {
3737 throw_dimension_incompatible("frequency(e, ...)", "e", expr);
3738 }
3739
3740 // If the `expr' has a constant value, then the frequency
3741 // `freq_n' is 0. Otherwise the values for \p expr are not discrete
3742 // and we return false.
3743
3744 // Space dimension is 0: if empty, then return false;
3745 // otherwise the frequency is 1 and the value is 0.
3746 if (space_dim == 0) {
3747 if (is_empty()) {
3748 return false;
3749 }
3750 freq_n = 0;
3751 freq_d = 1;
3752 val_n = expr.inhomogeneous_term();
3753 val_d = 1;
3754 return true;
3755 }
3756
3757 // For an empty polyhedron, we simply return false.
3758 if (marked_empty()
3759 || (has_pending_constraints() && !process_pending_constraints())
3760 || (!generators_are_up_to_date() && !update_generators())) {
3761 return false;
3762 }
3763
3764 // The polyhedron has updated, possibly pending generators.
3765 // The following loop will iterate through the generator
3766 // to see if `expr' has a constant value.
3767 PPL_DIRTY_TEMP(mpq_class, value);
3768
3769 // True if we have no other candidate value to compare with.
3770 bool first_candidate = true;
3771
3772 PPL_DIRTY_TEMP_COEFFICIENT(sp);
3773 PPL_DIRTY_TEMP(mpq_class, candidate);
3774 for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
3775 const Generator& gen_sys_i = gen_sys[i];
3776 Scalar_Products::homogeneous_assign(sp, expr, gen_sys_i);
3777 // Lines and rays in `*this' can cause `expr' to be non-constant.
3778 if (gen_sys_i.is_line_or_ray()) {
3779 const int sp_sign = sgn(sp);
3780 if (sp_sign != 0) {
3781 // `expr' is unbounded in `*this'.
3782 return false;
3783 }
3784 }
3785 else {
3786 // We have a point or a closure point.
3787 PPL_ASSERT(gen_sys_i.is_point() || gen_sys_i.is_closure_point());
3788 // Notice that we are ignoring the constant term in `expr' here.
3789 // We will add it to the value if there is a constant value.
3790 assign_r(candidate.get_num(), sp, ROUND_NOT_NEEDED);
3791 assign_r(candidate.get_den(), gen_sys_i.expr.inhomogeneous_term(), ROUND_NOT_NEEDED);
3792 candidate.canonicalize();
3793 if (first_candidate) {
3794 // We have a (new) candidate value.
3795 first_candidate = false;
3796 value = candidate;
3797 }
3798 else if (candidate != value) {
3799 return false;
3800 }
3801 }
3802 }
3803
3804 // Add in the constant term in `expr'.
3805 PPL_DIRTY_TEMP(mpz_class, n);
3806 assign_r(n, expr.inhomogeneous_term(), ROUND_NOT_NEEDED);
3807 value += n;
3808 val_n = value.get_num();
3809 val_d = value.get_den();
3810
3811 freq_n = 0;
3812 freq_d = 1;
3813 return true;
3814 }
3815
3816 void
topological_closure_assign()3817 PPL::Polyhedron::topological_closure_assign() {
3818 // Necessarily closed polyhedra are trivially closed.
3819 if (is_necessarily_closed()) {
3820 return;
3821 }
3822 // Any empty or zero-dimensional polyhedron is closed.
3823 if (marked_empty() || space_dim == 0) {
3824 return;
3825 }
3826
3827 // The computation can be done using constraints or generators.
3828 // If we use constraints, we will change them, so that having pending
3829 // constraints would be useless. If we use generators, we add generators,
3830 // so that having pending generators still makes sense.
3831
3832 // Process any pending constraints.
3833 if (has_pending_constraints() && !process_pending_constraints()) {
3834 return;
3835 }
3836
3837 // Use constraints only if they are available and
3838 // there are no pending generators.
3839 if (!has_pending_generators() && constraints_are_up_to_date()) {
3840 bool changed = false;
3841
3842 // Transform all strict inequalities into non-strict ones.
3843 for (dimension_type i = con_sys.num_rows(); i-- > 0; ) {
3844 Constraint& c = con_sys.sys.rows[i];
3845 if (c.epsilon_coefficient() < 0 && !c.is_tautological()) {
3846 c.set_epsilon_coefficient(0);
3847 // Enforce normalization.
3848 c.expr.normalize();
3849 PPL_ASSERT(c.OK());
3850 changed = true;
3851 }
3852 }
3853
3854 PPL_ASSERT(con_sys.sys.OK());
3855
3856 if (changed) {
3857 con_sys.insert(Constraint::epsilon_leq_one());
3858 con_sys.set_sorted(false);
3859 // After changing the system of constraints, the generators
3860 // are no longer up-to-date and the constraints are no longer
3861 // minimized.
3862 clear_generators_up_to_date();
3863 clear_constraints_minimized();
3864 }
3865 }
3866 else {
3867 // Here we use generators, possibly keeping constraints.
3868 PPL_ASSERT(generators_are_up_to_date());
3869 // Add the corresponding point to each closure point.
3870 gen_sys.add_corresponding_points();
3871 if (can_have_something_pending()) {
3872 set_generators_pending();
3873 }
3874 else {
3875 // We cannot have pending generators; this also implies
3876 // that generators may have lost their sortedness.
3877 gen_sys.unset_pending_rows();
3878 gen_sys.set_sorted(false);
3879 // Constraints are not up-to-date and generators are not minimized.
3880 clear_constraints_up_to_date();
3881 clear_generators_minimized();
3882 }
3883 }
3884 PPL_ASSERT_HEAVY(OK());
3885 }
3886
3887 /*! \relates Parma_Polyhedra_Library::Polyhedron */
3888 bool
operator ==(const Polyhedron & x,const Polyhedron & y)3889 PPL::operator==(const Polyhedron& x, const Polyhedron& y) {
3890 // If the two polyhedra are topology-incompatible or dimension-incompatible,
3891 // then they cannot be the same polyhedron.
3892 if (x.topology() != y.topology() || x.space_dim != y.space_dim) {
3893 return false;
3894 }
3895
3896 if (x.marked_empty()) {
3897 return y.is_empty();
3898 }
3899 else if (y.marked_empty()) {
3900 return x.is_empty();
3901 }
3902 else if (x.space_dim == 0) {
3903 return true;
3904 }
3905
3906 switch (x.quick_equivalence_test(y)) {
3907 case Polyhedron::TVB_TRUE:
3908 return true;
3909 case Polyhedron::TVB_FALSE:
3910 return false;
3911 default:
3912 if (x.is_included_in(y)) {
3913 if (x.marked_empty()) {
3914 return y.is_empty();
3915 }
3916 else {
3917 return y.is_included_in(x);
3918 }
3919 }
3920 else {
3921 return false;
3922 }
3923 }
3924 }
3925
3926 bool
contains(const Polyhedron & y) const3927 PPL::Polyhedron::contains(const Polyhedron& y) const {
3928 const Polyhedron& x = *this;
3929
3930 // Topology compatibility check.
3931 if (x.topology() != y.topology()) {
3932 throw_topology_incompatible("contains(y)", "y", y);
3933 }
3934
3935 // Dimension-compatibility check.
3936 if (x.space_dim != y.space_dim) {
3937 throw_dimension_incompatible("contains(y)", "y", y);
3938 }
3939
3940 if (y.marked_empty()) {
3941 return true;
3942 }
3943 else if (x.marked_empty()) {
3944 return y.is_empty();
3945 }
3946 else if (y.space_dim == 0) {
3947 return true;
3948 }
3949 else if (x.quick_equivalence_test(y) == Polyhedron::TVB_TRUE) {
3950 return true;
3951 }
3952 else {
3953 return y.is_included_in(x);
3954 }
3955 }
3956
3957 bool
is_disjoint_from(const Polyhedron & y) const3958 PPL::Polyhedron::is_disjoint_from(const Polyhedron& y) const {
3959 Polyhedron z = *this;
3960 z.intersection_assign(y);
3961 return z.is_empty();
3962 }
3963
3964 void
ascii_dump(std::ostream & s) const3965 PPL::Polyhedron::ascii_dump(std::ostream& s) const {
3966 s << "space_dim " << space_dim << "\n";
3967 status.ascii_dump(s);
3968 s << "\ncon_sys ("
3969 << (constraints_are_up_to_date() ? "" : "not_")
3970 << "up-to-date)"
3971 << "\n";
3972 con_sys.ascii_dump(s);
3973 s << "\ngen_sys ("
3974 << (generators_are_up_to_date() ? "" : "not_")
3975 << "up-to-date)"
3976 << "\n";
3977 gen_sys.ascii_dump(s);
3978 s << "\nsat_c\n";
3979 sat_c.ascii_dump(s);
3980 s << "\nsat_g\n";
3981 sat_g.ascii_dump(s);
3982 s << "\n";
3983 }
3984
PPL_OUTPUT_DEFINITIONS(Polyhedron)3985 PPL_OUTPUT_DEFINITIONS(Polyhedron)
3986
3987 bool
3988 PPL::Polyhedron::ascii_load(std::istream& s) {
3989 std::string str;
3990
3991 if (!(s >> str) || str != "space_dim") {
3992 return false;
3993 }
3994
3995 if (!(s >> space_dim)) {
3996 return false;
3997 }
3998
3999 if (!status.ascii_load(s)) {
4000 return false;
4001 }
4002
4003 if (!(s >> str) || str != "con_sys") {
4004 return false;
4005 }
4006
4007 if (!(s >> str) || (str != "(not_up-to-date)" && str != "(up-to-date)")) {
4008 return false;
4009 }
4010
4011 if (!con_sys.ascii_load(s)) {
4012 return false;
4013 }
4014
4015 if (!(s >> str) || str != "gen_sys") {
4016 return false;
4017 }
4018
4019 if (!(s >> str) || (str != "(not_up-to-date)" && str != "(up-to-date)")) {
4020 return false;
4021 }
4022
4023 if (!gen_sys.ascii_load(s)) {
4024 return false;
4025 }
4026
4027 if (!(s >> str) || str != "sat_c") {
4028 return false;
4029 }
4030
4031 if (!sat_c.ascii_load(s)) {
4032 return false;
4033 }
4034
4035 if (!(s >> str) || str != "sat_g") {
4036 return false;
4037 }
4038
4039 if (!sat_g.ascii_load(s)) {
4040 return false;
4041 }
4042
4043 // Check invariants.
4044 PPL_ASSERT_HEAVY(OK());
4045 return true;
4046 }
4047
4048 PPL::memory_size_type
external_memory_in_bytes() const4049 PPL::Polyhedron::external_memory_in_bytes() const {
4050 return
4051 con_sys.external_memory_in_bytes()
4052 + gen_sys.external_memory_in_bytes()
4053 + sat_c.external_memory_in_bytes()
4054 + sat_g.external_memory_in_bytes();
4055 }
4056
4057 void
wrap_assign(const Variables_Set & vars,Bounded_Integer_Type_Width w,Bounded_Integer_Type_Representation r,Bounded_Integer_Type_Overflow o,const Constraint_System * cs_p,unsigned complexity_threshold,bool wrap_individually)4058 PPL::Polyhedron::wrap_assign(const Variables_Set& vars,
4059 Bounded_Integer_Type_Width w,
4060 Bounded_Integer_Type_Representation r,
4061 Bounded_Integer_Type_Overflow o,
4062 const Constraint_System* cs_p,
4063 unsigned complexity_threshold,
4064 bool wrap_individually) {
4065 if (is_necessarily_closed()) {
4066 Implementation::wrap_assign(static_cast<C_Polyhedron&>(*this),
4067 vars, w, r, o, cs_p,
4068 complexity_threshold, wrap_individually,
4069 "C_Polyhedron");
4070 }
4071 else {
4072 Implementation::wrap_assign(static_cast<NNC_Polyhedron&>(*this),
4073 vars, w, r, o, cs_p,
4074 complexity_threshold, wrap_individually,
4075 "NNC_Polyhedron");
4076 }
4077 }
4078
4079 /*! \relates Parma_Polyhedra_Library::Polyhedron */
4080 std::ostream&
operator <<(std::ostream & s,const Polyhedron & ph)4081 PPL::IO_Operators::operator<<(std::ostream& s, const Polyhedron& ph) {
4082 if (ph.is_empty()) {
4083 s << "false";
4084 }
4085 else {
4086 s << ph.minimized_constraints();
4087 }
4088 return s;
4089 }
4090