1 /* Pointset_Powerset class implementation: non-inline template 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 #ifndef PPL_Pointset_Powerset_templates_hh
25 #define PPL_Pointset_Powerset_templates_hh 1
26
27 #include "Constraint_defs.hh"
28 #include "Constraint_System_defs.hh"
29 #include "Constraint_System_inlines.hh"
30 #include "C_Polyhedron_defs.hh"
31 #include "NNC_Polyhedron_defs.hh"
32 #include "Variables_Set_defs.hh"
33 #include <algorithm>
34 #include <deque>
35 #include <string>
36 #include <iostream>
37 #include <sstream>
38 #include <stdexcept>
39
40 namespace Parma_Polyhedra_Library {
41
42 template <typename PSET>
43 void
add_disjunct(const PSET & ph)44 Pointset_Powerset<PSET>::add_disjunct(const PSET& ph) {
45 Pointset_Powerset& x = *this;
46 if (x.space_dimension() != ph.space_dimension()) {
47 std::ostringstream s;
48 s << "PPL::Pointset_Powerset<PSET>::add_disjunct(ph):\n"
49 << "this->space_dimension() == " << x.space_dimension() << ", "
50 << "ph.space_dimension() == " << ph.space_dimension() << ".";
51 throw std::invalid_argument(s.str());
52 }
53 x.sequence.push_back(Determinate<PSET>(ph));
54 x.reduced = false;
55 PPL_ASSERT_HEAVY(x.OK());
56 }
57
58 template <>
59 template <typename QH>
60 Pointset_Powerset<NNC_Polyhedron>
Pointset_Powerset(const Pointset_Powerset<QH> & y,Complexity_Class complexity)61 ::Pointset_Powerset(const Pointset_Powerset<QH>& y,
62 Complexity_Class complexity)
63 : Base(), space_dim(y.space_dimension()) {
64 Pointset_Powerset& x = *this;
65 for (typename Pointset_Powerset<QH>::const_iterator i = y.begin(),
66 y_end = y.end(); i != y_end; ++i) {
67 x.sequence.push_back(Determinate<NNC_Polyhedron>
68 (NNC_Polyhedron(i->pointset(), complexity)));
69
70 }
71 // FIXME: If the domain elements can be represented _exactly_ as NNC
72 // polyhedra, then having x.reduced = y.reduced is correct. This is
73 // the case if the domains are both linear and convex which holds
74 // for all the currently supported instantiations except for
75 // Grids; for this reason the Grid specialization has a
76 // separate implementation. For any non-linear or non-convex
77 // domains (e.g., a domain of Intervals with restrictions or a
78 // domain of circles) that may be supported in the future, the
79 // assignment x.reduced = y.reduced will be a bug.
80 x.reduced = y.reduced;
81
82 PPL_ASSERT_HEAVY(x.OK());
83 }
84
85 template <typename PSET>
86 template <typename QH>
87 Pointset_Powerset<PSET>
Pointset_Powerset(const Pointset_Powerset<QH> & y,Complexity_Class complexity)88 ::Pointset_Powerset(const Pointset_Powerset<QH>& y,
89 Complexity_Class complexity)
90 : Base(), space_dim(y.space_dimension()) {
91 Pointset_Powerset& x = *this;
92 for (typename Pointset_Powerset<QH>::const_iterator i = y.begin(),
93 y_end = y.end(); i != y_end; ++i) {
94 x.sequence.push_back(Determinate<PSET>(PSET(i->pointset(), complexity)));
95 }
96 // Note: this might be non-reduced even when `y' is known to be
97 // omega-reduced, because the constructor of PSET may have made
98 // different QH elements to become comparable.
99 x.reduced = false;
100 PPL_ASSERT_HEAVY(x.OK());
101 }
102
103 template <typename PSET>
104 void
concatenate_assign(const Pointset_Powerset & y)105 Pointset_Powerset<PSET>::concatenate_assign(const Pointset_Powerset& y) {
106 Pointset_Powerset& x = *this;
107 // Ensure omega-reduction here, since what follows has quadratic complexity.
108 x.omega_reduce();
109 y.omega_reduce();
110 Pointset_Powerset<PSET> new_x(x.space_dim + y.space_dim, EMPTY);
111 for (const_iterator xi = x.begin(), x_end = x.end(),
112 y_begin = y.begin(), y_end = y.end(); xi != x_end; ) {
113 for (const_iterator yi = y_begin; yi != y_end; ++yi) {
114 Det_PSET zi = *xi;
115 zi.concatenate_assign(*yi);
116 PPL_ASSERT_HEAVY(!zi.is_bottom());
117 new_x.sequence.push_back(zi);
118 }
119 ++xi;
120 if ((abandon_expensive_computations != 0)
121 && (xi != x_end) && (y_begin != y_end)) {
122 // Hurry up!
123 PSET x_ph = xi->pointset();
124 for (++xi; xi != x_end; ++xi) {
125 x_ph.upper_bound_assign(xi->pointset());
126 }
127 const_iterator yi = y_begin;
128 PSET y_ph = yi->pointset();
129 for (++yi; yi != y_end; ++yi) {
130 y_ph.upper_bound_assign(yi->pointset());
131 }
132 x_ph.concatenate_assign(y_ph);
133 swap(x, new_x);
134 x.add_disjunct(x_ph);
135 PPL_ASSERT_HEAVY(x.OK());
136 return;
137 }
138 }
139 swap(x, new_x);
140 PPL_ASSERT_HEAVY(x.OK());
141 }
142
143 template <typename PSET>
144 void
add_constraint(const Constraint & c)145 Pointset_Powerset<PSET>::add_constraint(const Constraint& c) {
146 Pointset_Powerset& x = *this;
147 for (Sequence_iterator si = x.sequence.begin(),
148 s_end = x.sequence.end(); si != s_end; ++si) {
149 si->pointset().add_constraint(c);
150 }
151 x.reduced = false;
152 PPL_ASSERT_HEAVY(x.OK());
153 }
154
155 template <typename PSET>
156 void
refine_with_constraint(const Constraint & c)157 Pointset_Powerset<PSET>::refine_with_constraint(const Constraint& c) {
158 Pointset_Powerset& x = *this;
159 for (Sequence_iterator si = x.sequence.begin(),
160 s_end = x.sequence.end(); si != s_end; ++si) {
161 si->pointset().refine_with_constraint(c);
162 }
163 x.reduced = false;
164 PPL_ASSERT_HEAVY(x.OK());
165 }
166
167 template <typename PSET>
168 void
add_constraints(const Constraint_System & cs)169 Pointset_Powerset<PSET>::add_constraints(const Constraint_System& cs) {
170 Pointset_Powerset& x = *this;
171 for (Sequence_iterator si = x.sequence.begin(),
172 s_end = x.sequence.end(); si != s_end; ++si) {
173 si->pointset().add_constraints(cs);
174 }
175 x.reduced = false;
176 PPL_ASSERT_HEAVY(x.OK());
177 }
178
179 template <typename PSET>
180 void
refine_with_constraints(const Constraint_System & cs)181 Pointset_Powerset<PSET>::refine_with_constraints(const Constraint_System& cs) {
182 Pointset_Powerset& x = *this;
183 for (Sequence_iterator si = x.sequence.begin(),
184 s_end = x.sequence.end(); si != s_end; ++si) {
185 si->pointset().refine_with_constraints(cs);
186 }
187 x.reduced = false;
188 PPL_ASSERT_HEAVY(x.OK());
189 }
190
191 template <typename PSET>
192 void
add_congruence(const Congruence & cg)193 Pointset_Powerset<PSET>::add_congruence(const Congruence& cg) {
194 Pointset_Powerset& x = *this;
195 for (Sequence_iterator si = x.sequence.begin(),
196 s_end = x.sequence.end(); si != s_end; ++si) {
197 si->pointset().add_congruence(cg);
198 }
199 x.reduced = false;
200 PPL_ASSERT_HEAVY(x.OK());
201 }
202
203 template <typename PSET>
204 void
refine_with_congruence(const Congruence & cg)205 Pointset_Powerset<PSET>::refine_with_congruence(const Congruence& cg) {
206 Pointset_Powerset& x = *this;
207 for (Sequence_iterator si = x.sequence.begin(),
208 s_end = x.sequence.end(); si != s_end; ++si) {
209 si->pointset().refine_with_congruence(cg);
210 }
211 x.reduced = false;
212 PPL_ASSERT_HEAVY(x.OK());
213 }
214
215 template <typename PSET>
216 void
add_congruences(const Congruence_System & cgs)217 Pointset_Powerset<PSET>::add_congruences(const Congruence_System& cgs) {
218 Pointset_Powerset& x = *this;
219 for (Sequence_iterator si = x.sequence.begin(),
220 s_end = x.sequence.end(); si != s_end; ++si) {
221 si->pointset().add_congruences(cgs);
222 }
223 x.reduced = false;
224 PPL_ASSERT_HEAVY(x.OK());
225 }
226
227 template <typename PSET>
228 void
refine_with_congruences(const Congruence_System & cgs)229 Pointset_Powerset<PSET>::refine_with_congruences(const Congruence_System& cgs) {
230 Pointset_Powerset& x = *this;
231 for (Sequence_iterator si = x.sequence.begin(),
232 s_end = x.sequence.end(); si != s_end; ++si) {
233 si->pointset().refine_with_congruences(cgs);
234 }
235 x.reduced = false;
236 PPL_ASSERT_HEAVY(x.OK());
237 }
238
239 template <typename PSET>
240 void
unconstrain(const Variable var)241 Pointset_Powerset<PSET>::unconstrain(const Variable var) {
242 Pointset_Powerset& x = *this;
243 for (Sequence_iterator si = x.sequence.begin(),
244 s_end = x.sequence.end(); si != s_end; ++si) {
245 si->pointset().unconstrain(var);
246 x.reduced = false;
247 }
248 PPL_ASSERT_HEAVY(x.OK());
249 }
250
251 template <typename PSET>
252 void
unconstrain(const Variables_Set & vars)253 Pointset_Powerset<PSET>::unconstrain(const Variables_Set& vars) {
254 Pointset_Powerset& x = *this;
255 for (Sequence_iterator si = x.sequence.begin(),
256 s_end = x.sequence.end(); si != s_end; ++si) {
257 si->pointset().unconstrain(vars);
258 x.reduced = false;
259 }
260 PPL_ASSERT_HEAVY(x.OK());
261 }
262
263 template <typename PSET>
264 void
add_space_dimensions_and_embed(dimension_type m)265 Pointset_Powerset<PSET>::add_space_dimensions_and_embed(dimension_type m) {
266 Pointset_Powerset& x = *this;
267 for (Sequence_iterator si = x.sequence.begin(),
268 s_end = x.sequence.end(); si != s_end; ++si) {
269 si->pointset().add_space_dimensions_and_embed(m);
270 }
271 x.space_dim += m;
272 PPL_ASSERT_HEAVY(x.OK());
273 }
274
275 template <typename PSET>
276 void
add_space_dimensions_and_project(dimension_type m)277 Pointset_Powerset<PSET>::add_space_dimensions_and_project(dimension_type m) {
278 Pointset_Powerset& x = *this;
279 for (Sequence_iterator si = x.sequence.begin(),
280 s_end = x.sequence.end(); si != s_end; ++si) {
281 si->pointset().add_space_dimensions_and_project(m);
282 }
283 x.space_dim += m;
284 PPL_ASSERT_HEAVY(x.OK());
285 }
286
287 template <typename PSET>
288 void
remove_space_dimensions(const Variables_Set & vars)289 Pointset_Powerset<PSET>::remove_space_dimensions(const Variables_Set& vars) {
290 Pointset_Powerset& x = *this;
291 Variables_Set::size_type num_removed = vars.size();
292 if (num_removed > 0) {
293 for (Sequence_iterator si = x.sequence.begin(),
294 s_end = x.sequence.end(); si != s_end; ++si) {
295 si->pointset().remove_space_dimensions(vars);
296 x.reduced = false;
297 }
298 x.space_dim -= num_removed;
299 PPL_ASSERT_HEAVY(x.OK());
300 }
301 }
302
303 template <typename PSET>
304 void
305 Pointset_Powerset<PSET>
remove_higher_space_dimensions(dimension_type new_dimension)306 ::remove_higher_space_dimensions(dimension_type new_dimension) {
307 Pointset_Powerset& x = *this;
308 if (new_dimension < x.space_dim) {
309 for (Sequence_iterator si = x.sequence.begin(),
310 s_end = x.sequence.end(); si != s_end; ++si) {
311 si->pointset().remove_higher_space_dimensions(new_dimension);
312 x.reduced = false;
313 }
314 x.space_dim = new_dimension;
315 PPL_ASSERT_HEAVY(x.OK());
316 }
317 }
318
319 template <typename PSET>
320 template <typename Partial_Function>
321 void
map_space_dimensions(const Partial_Function & pfunc)322 Pointset_Powerset<PSET>::map_space_dimensions(const Partial_Function& pfunc) {
323 Pointset_Powerset& x = *this;
324 if (x.is_bottom()) {
325 dimension_type n = 0;
326 for (dimension_type i = x.space_dim; i-- > 0; ) {
327 dimension_type new_i;
328 if (pfunc.maps(i, new_i)) {
329 ++n;
330 }
331 }
332 x.space_dim = n;
333 }
334 else {
335 Sequence_iterator s_begin = x.sequence.begin();
336 for (Sequence_iterator si = s_begin,
337 s_end = x.sequence.end(); si != s_end; ++si) {
338 si->pointset().map_space_dimensions(pfunc);
339 }
340 x.space_dim = s_begin->pointset().space_dimension();
341 x.reduced = false;
342 }
343 PPL_ASSERT_HEAVY(x.OK());
344 }
345
346 template <typename PSET>
347 void
expand_space_dimension(Variable var,dimension_type m)348 Pointset_Powerset<PSET>::expand_space_dimension(Variable var,
349 dimension_type m) {
350 Pointset_Powerset& x = *this;
351 for (Sequence_iterator si = x.sequence.begin(),
352 s_end = x.sequence.end(); si != s_end; ++si) {
353 si->pointset().expand_space_dimension(var, m);
354 }
355 x.space_dim += m;
356 PPL_ASSERT_HEAVY(x.OK());
357 }
358
359 template <typename PSET>
360 void
fold_space_dimensions(const Variables_Set & vars,Variable dest)361 Pointset_Powerset<PSET>::fold_space_dimensions(const Variables_Set& vars,
362 Variable dest) {
363 Pointset_Powerset& x = *this;
364 Variables_Set::size_type num_folded = vars.size();
365 if (num_folded > 0) {
366 for (Sequence_iterator si = x.sequence.begin(),
367 s_end = x.sequence.end(); si != s_end; ++si) {
368 si->pointset().fold_space_dimensions(vars, dest);
369 }
370 }
371 x.space_dim -= num_folded;
372 PPL_ASSERT_HEAVY(x.OK());
373 }
374
375 template <typename PSET>
376 void
affine_image(Variable var,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)377 Pointset_Powerset<PSET>::affine_image(Variable var,
378 const Linear_Expression& expr,
379 Coefficient_traits::const_reference
380 denominator) {
381 Pointset_Powerset& x = *this;
382 for (Sequence_iterator si = x.sequence.begin(),
383 s_end = x.sequence.end(); si != s_end; ++si) {
384 si->pointset().affine_image(var, expr, denominator);
385 // Note that the underlying domain can apply conservative approximation:
386 // that is why it would not be correct to make the loss of reduction
387 // conditional on `var' and `expr'.
388 x.reduced = false;
389 }
390 PPL_ASSERT_HEAVY(x.OK());
391 }
392
393 template <typename PSET>
394 void
affine_preimage(Variable var,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)395 Pointset_Powerset<PSET>::affine_preimage(Variable var,
396 const Linear_Expression& expr,
397 Coefficient_traits::const_reference
398 denominator) {
399 Pointset_Powerset& x = *this;
400 for (Sequence_iterator si = x.sequence.begin(),
401 s_end = x.sequence.end(); si != s_end; ++si) {
402 si->pointset().affine_preimage(var, expr, denominator);
403 // Note that the underlying domain can apply conservative approximation:
404 // that is why it would not be correct to make the loss of reduction
405 // conditional on `var' and `expr'.
406 x.reduced = false;
407 }
408 PPL_ASSERT_HEAVY(x.OK());
409 }
410
411
412 template <typename PSET>
413 void
414 Pointset_Powerset<PSET>
generalized_affine_image(const Linear_Expression & lhs,const Relation_Symbol relsym,const Linear_Expression & rhs)415 ::generalized_affine_image(const Linear_Expression& lhs,
416 const Relation_Symbol relsym,
417 const Linear_Expression& rhs) {
418 Pointset_Powerset& x = *this;
419 for (Sequence_iterator si = x.sequence.begin(),
420 s_end = x.sequence.end(); si != s_end; ++si) {
421 si->pointset().generalized_affine_image(lhs, relsym, rhs);
422 x.reduced = false;
423 }
424 PPL_ASSERT_HEAVY(x.OK());
425 }
426
427 template <typename PSET>
428 void
429 Pointset_Powerset<PSET>
generalized_affine_preimage(const Linear_Expression & lhs,const Relation_Symbol relsym,const Linear_Expression & rhs)430 ::generalized_affine_preimage(const Linear_Expression& lhs,
431 const Relation_Symbol relsym,
432 const Linear_Expression& rhs) {
433 Pointset_Powerset& x = *this;
434 for (Sequence_iterator si = x.sequence.begin(),
435 s_end = x.sequence.end(); si != s_end; ++si) {
436 si->pointset().generalized_affine_preimage(lhs, relsym, rhs);
437 x.reduced = false;
438 }
439 PPL_ASSERT_HEAVY(x.OK());
440 }
441
442 template <typename PSET>
443 void
444 Pointset_Powerset<PSET>
generalized_affine_image(Variable var,const Relation_Symbol relsym,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)445 ::generalized_affine_image(Variable var,
446 const Relation_Symbol relsym,
447 const Linear_Expression& expr,
448 Coefficient_traits::const_reference denominator) {
449 Pointset_Powerset& x = *this;
450 for (Sequence_iterator si = x.sequence.begin(),
451 s_end = x.sequence.end(); si != s_end; ++si) {
452 si->pointset().generalized_affine_image(var, relsym, expr, denominator);
453 x.reduced = false;
454 }
455 PPL_ASSERT_HEAVY(x.OK());
456 }
457
458 template <typename PSET>
459 void
460 Pointset_Powerset<PSET>
generalized_affine_preimage(Variable var,const Relation_Symbol relsym,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)461 ::generalized_affine_preimage(Variable var,
462 const Relation_Symbol relsym,
463 const Linear_Expression& expr,
464 Coefficient_traits::const_reference
465 denominator) {
466 Pointset_Powerset& x = *this;
467 for (Sequence_iterator si = x.sequence.begin(),
468 s_end = x.sequence.end(); si != s_end; ++si) {
469 si->pointset().generalized_affine_preimage(var, relsym, expr, denominator);
470 x.reduced = false;
471 }
472 PPL_ASSERT_HEAVY(x.OK());
473 }
474
475
476 template <typename PSET>
477 void
478 Pointset_Powerset<PSET>
bounded_affine_image(Variable var,const Linear_Expression & lb_expr,const Linear_Expression & ub_expr,Coefficient_traits::const_reference denominator)479 ::bounded_affine_image(Variable var,
480 const Linear_Expression& lb_expr,
481 const Linear_Expression& ub_expr,
482 Coefficient_traits::const_reference denominator) {
483 Pointset_Powerset& x = *this;
484 for (Sequence_iterator si = x.sequence.begin(),
485 s_end = x.sequence.end(); si != s_end; ++si) {
486 si->pointset().bounded_affine_image(var, lb_expr, ub_expr, denominator);
487 x.reduced = false;
488 }
489 PPL_ASSERT_HEAVY(x.OK());
490 }
491
492 template <typename PSET>
493 void
494 Pointset_Powerset<PSET>
bounded_affine_preimage(Variable var,const Linear_Expression & lb_expr,const Linear_Expression & ub_expr,Coefficient_traits::const_reference denominator)495 ::bounded_affine_preimage(Variable var,
496 const Linear_Expression& lb_expr,
497 const Linear_Expression& ub_expr,
498 Coefficient_traits::const_reference denominator) {
499 Pointset_Powerset& x = *this;
500 for (Sequence_iterator si = x.sequence.begin(),
501 s_end = x.sequence.end(); si != s_end; ++si) {
502 si->pointset().bounded_affine_preimage(var, lb_expr, ub_expr,
503 denominator);
504 x.reduced = false;
505 }
506 PPL_ASSERT_HEAVY(x.OK());
507 }
508
509 template <typename PSET>
510 dimension_type
affine_dimension() const511 Pointset_Powerset<PSET>::affine_dimension() const {
512 // The affine dimension of the powerset is the affine dimension of
513 // the smallest vector space in which it can be embedded.
514 const Pointset_Powerset& x = *this;
515 C_Polyhedron x_ph(space_dim, EMPTY);
516
517 for (Sequence_const_iterator si = x.sequence.begin(),
518 s_end = x.sequence.end(); si != s_end; ++si) {
519 PSET pi(si->pointset());
520 if (!pi.is_empty()) {
521 C_Polyhedron phi(space_dim);
522 const Constraint_System& cs = pi.minimized_constraints();
523 for (Constraint_System::const_iterator i = cs.begin(),
524 cs_end = cs.end(); i != cs_end; ++i) {
525 const Constraint& c = *i;
526 if (c.is_equality()) {
527 phi.add_constraint(c);
528 }
529 }
530 x_ph.poly_hull_assign(phi);
531 }
532 }
533
534 return x_ph.affine_dimension();
535 }
536
537 template <typename PSET>
538 bool
is_universe() const539 Pointset_Powerset<PSET>::is_universe() const {
540 const Pointset_Powerset& x = *this;
541 // Exploit omega-reduction, if already computed.
542 if (x.is_omega_reduced()) {
543 return x.size() == 1 && x.begin()->pointset().is_universe();
544 }
545
546 // A powerset is universe iff one of its disjuncts is.
547 for (const_iterator x_i = x.begin(), x_end = x.end();
548 x_i != x_end; ++x_i) {
549 if (x_i->pointset().is_universe()) {
550 // Speculative omega-reduction, if it is worth.
551 if (x.size() > 1) {
552 Pointset_Powerset<PSET> universe(x.space_dimension(), UNIVERSE);
553 Pointset_Powerset& xx = const_cast<Pointset_Powerset&>(x);
554 swap(xx, universe);
555 }
556 return true;
557 }
558 }
559 return false;
560 }
561
562 template <typename PSET>
563 bool
is_empty() const564 Pointset_Powerset<PSET>::is_empty() const {
565 const Pointset_Powerset& x = *this;
566 for (Sequence_const_iterator si = x.sequence.begin(),
567 s_end = x.sequence.end(); si != s_end; ++si) {
568 if (!si->pointset().is_empty()) {
569 return false;
570 }
571 }
572 return true;
573 }
574
575 template <typename PSET>
576 bool
is_discrete() const577 Pointset_Powerset<PSET>::is_discrete() const {
578 const Pointset_Powerset& x = *this;
579 for (Sequence_const_iterator si = x.sequence.begin(),
580 s_end = x.sequence.end(); si != s_end; ++si) {
581 if (!si->pointset().is_discrete()) {
582 return false;
583 }
584 }
585 return true;
586 }
587
588 template <typename PSET>
589 bool
is_topologically_closed() const590 Pointset_Powerset<PSET>::is_topologically_closed() const {
591 const Pointset_Powerset& x = *this;
592 // The powerset must be omega-reduced before checking
593 // topological closure.
594 x.omega_reduce();
595 for (Sequence_const_iterator si = x.sequence.begin(),
596 s_end = x.sequence.end(); si != s_end; ++si) {
597 if (!si->pointset().is_topologically_closed()) {
598 return false;
599 }
600 }
601 return true;
602 }
603
604 template <typename PSET>
605 bool
is_bounded() const606 Pointset_Powerset<PSET>::is_bounded() const {
607 const Pointset_Powerset& x = *this;
608 for (Sequence_const_iterator si = x.sequence.begin(),
609 s_end = x.sequence.end(); si != s_end; ++si) {
610 if (!si->pointset().is_bounded()) {
611 return false;
612 }
613 }
614 return true;
615 }
616
617 template <typename PSET>
618 bool
constrains(Variable var) const619 Pointset_Powerset<PSET>::constrains(Variable var) const {
620 const Pointset_Powerset& x = *this;
621 // `var' should be one of the dimensions of the powerset.
622 const dimension_type var_space_dim = var.space_dimension();
623 if (x.space_dimension() < var_space_dim) {
624 std::ostringstream s;
625 s << "PPL::Pointset_Powerset<PSET>::constrains(v):\n"
626 << "this->space_dimension() == " << x.space_dimension() << ", "
627 << "v.space_dimension() == " << var_space_dim << ".";
628 throw std::invalid_argument(s.str());
629 }
630 // omega_reduction needed, since a redundant disjunct may constrain var.
631 x.omega_reduce();
632 // An empty powerset constrains all variables.
633 if (x.is_empty()) {
634 return true;
635 }
636 for (const_iterator x_i = x.begin(), x_end = x.end();
637 x_i != x_end; ++x_i) {
638 if (x_i->pointset().constrains(var)) {
639 return true;
640 }
641 }
642 return false;
643 }
644
645 template <typename PSET>
646 bool
is_disjoint_from(const Pointset_Powerset & y) const647 Pointset_Powerset<PSET>::is_disjoint_from(const Pointset_Powerset& y) const {
648 const Pointset_Powerset& x = *this;
649 for (Sequence_const_iterator si = x.sequence.begin(),
650 x_s_end = x.sequence.end(); si != x_s_end; ++si) {
651 const PSET& pi = si->pointset();
652 for (Sequence_const_iterator sj = y.sequence.begin(),
653 y_s_end = y.sequence.end(); sj != y_s_end; ++sj) {
654 const PSET& pj = sj->pointset();
655 if (!pi.is_disjoint_from(pj)) {
656 return false;
657 }
658 }
659 }
660 return true;
661 }
662
663 template <typename PSET>
664 void
665 Pointset_Powerset<PSET>
drop_some_non_integer_points(const Variables_Set & vars,Complexity_Class complexity)666 ::drop_some_non_integer_points(const Variables_Set& vars,
667 Complexity_Class complexity) {
668 Pointset_Powerset& x = *this;
669 for (Sequence_iterator si = x.sequence.begin(),
670 s_end = x.sequence.end(); si != s_end; ++si) {
671 si->pointset().drop_some_non_integer_points(vars, complexity);
672 }
673 x.reduced = false;
674 PPL_ASSERT_HEAVY(x.OK());
675 }
676
677 template <typename PSET>
678 void
679 Pointset_Powerset<PSET>
drop_some_non_integer_points(Complexity_Class complexity)680 ::drop_some_non_integer_points(Complexity_Class complexity) {
681 Pointset_Powerset& x = *this;
682 for (Sequence_iterator si = x.sequence.begin(),
683 s_end = x.sequence.end(); si != s_end; ++si) {
684 si->pointset().drop_some_non_integer_points(complexity);
685 }
686 x.reduced = false;
687 PPL_ASSERT_HEAVY(x.OK());
688 }
689
690 template <typename PSET>
691 void
topological_closure_assign()692 Pointset_Powerset<PSET>::topological_closure_assign() {
693 Pointset_Powerset& x = *this;
694 for (Sequence_iterator si = x.sequence.begin(),
695 s_end = x.sequence.end(); si != s_end; ++si) {
696 si->pointset().topological_closure_assign();
697 }
698 PPL_ASSERT_HEAVY(x.OK());
699 }
700
701 template <typename PSET>
702 bool
703 Pointset_Powerset<PSET>
intersection_preserving_enlarge_element(PSET & dest) const704 ::intersection_preserving_enlarge_element(PSET& dest) const {
705 // FIXME: this is just an executable specification.
706 const Pointset_Powerset& context = *this;
707 PPL_ASSERT(context.space_dimension() == dest.space_dimension());
708 bool nonempty_intersection = false;
709 // TODO: maybe use a *sorted* constraint system?
710 PSET enlarged(context.space_dimension(), UNIVERSE);
711 for (Sequence_const_iterator si = context.sequence.begin(),
712 s_end = context.sequence.end(); si != s_end; ++si) {
713 PSET context_i(si->pointset());
714 context_i.intersection_assign(enlarged);
715 PSET enlarged_i(dest);
716 if (enlarged_i.simplify_using_context_assign(context_i)) {
717 nonempty_intersection = true;
718 }
719 // TODO: merge the sorted constraints of `enlarged' and `enlarged_i'?
720 enlarged.intersection_assign(enlarged_i);
721 }
722 swap(dest, enlarged);
723 return nonempty_intersection;
724 }
725
726 template <typename PSET>
727 bool
728 Pointset_Powerset<PSET>
simplify_using_context_assign(const Pointset_Powerset & y)729 ::simplify_using_context_assign(const Pointset_Powerset& y) {
730 Pointset_Powerset& x = *this;
731
732 // Omega reduction is required.
733 // TODO: check whether it would be more efficient to Omega-reduce x
734 // during the simplification process: when examining *si, we check
735 // if it has been made redundant by any of the elements preceding it
736 // (which have been already simplified).
737 x.omega_reduce();
738 if (x.is_empty()) {
739 return false;
740 }
741 y.omega_reduce();
742 if (y.is_empty()) {
743 x = y;
744 return false;
745 }
746
747 if (y.size() == 1) {
748 // More efficient, special handling of the singleton context case.
749 const PSET& y_i = y.sequence.begin()->pointset();
750 for (Sequence_iterator si = x.sequence.begin(),
751 s_end = x.sequence.end(); si != s_end; ) {
752 PSET& x_i = si->pointset();
753 if (x_i.simplify_using_context_assign(y_i)) {
754 ++si;
755 }
756 else {
757 // Intersection is empty: drop the disjunct.
758 si = x.sequence.erase(si);
759 }
760 }
761 }
762 else {
763 // The context is not a singleton.
764 for (Sequence_iterator si = x.sequence.begin(),
765 s_end = x.sequence.end(); si != s_end; ) {
766 if (y.intersection_preserving_enlarge_element(si->pointset())) {
767 ++si;
768 }
769 else {
770 // Intersection with `*si' is empty: drop the disjunct.
771 si = x.sequence.erase(si);
772 }
773 }
774 }
775 x.reduced = false;
776 PPL_ASSERT_HEAVY(x.OK());
777 return !x.sequence.empty();
778 }
779
780 template <typename PSET>
781 bool
contains(const Pointset_Powerset & y) const782 Pointset_Powerset<PSET>::contains(const Pointset_Powerset& y) const {
783 const Pointset_Powerset& x = *this;
784 for (Sequence_const_iterator si = y.sequence.begin(),
785 y_s_end = y.sequence.end(); si != y_s_end; ++si) {
786 const PSET& pi = si->pointset();
787 bool pi_is_contained = false;
788 for (Sequence_const_iterator sj = x.sequence.begin(),
789 x_s_end = x.sequence.end();
790 (sj != x_s_end && !pi_is_contained);
791 ++sj) {
792 const PSET& pj = sj->pointset();
793 if (pj.contains(pi)) {
794 pi_is_contained = true;
795 }
796 }
797 if (!pi_is_contained) {
798 return false;
799 }
800 }
801 return true;
802 }
803
804 template <typename PSET>
805 bool
strictly_contains(const Pointset_Powerset & y) const806 Pointset_Powerset<PSET>::strictly_contains(const Pointset_Powerset& y) const {
807 /* omega reduction ensures that a disjunct of y cannot be strictly
808 contained in one disjunct and also contained but not strictly
809 contained in another disjunct of *this */
810 const Pointset_Powerset& x = *this;
811 x.omega_reduce();
812 for (Sequence_const_iterator si = y.sequence.begin(),
813 y_s_end = y.sequence.end(); si != y_s_end; ++si) {
814 const PSET& pi = si->pointset();
815 bool pi_is_strictly_contained = false;
816 for (Sequence_const_iterator sj = x.sequence.begin(),
817 x_s_end = x.sequence.end();
818 (sj != x_s_end && !pi_is_strictly_contained); ++sj) {
819 const PSET& pj = sj->pointset();
820 if (pj.strictly_contains(pi)) {
821 pi_is_strictly_contained = true;
822 }
823 }
824 if (!pi_is_strictly_contained) {
825 return false;
826 }
827 }
828 return true;
829 }
830
831 template <typename PSET>
832 template <typename Cons_or_Congr>
833 Poly_Con_Relation
relation_with_aux(const Cons_or_Congr & c) const834 Pointset_Powerset<PSET>::relation_with_aux(const Cons_or_Congr& c) const {
835 const Pointset_Powerset& x = *this;
836
837 /* *this is included in c if every disjunct is included in c */
838 bool is_included = true;
839 /* *this is disjoint with c if every disjunct is disjoint with c */
840 bool is_disjoint = true;
841 /* *this strictly_intersects with c if:
842 - some disjunct strictly intersects with c
843 or
844 - there exists two disjoints d1 and d2
845 such that d1 is included in c and d2 is disjoint with c
846 */
847 bool is_strictly_intersecting = false;
848 bool included_once = false;
849 bool disjoint_once = false;
850 /* *this saturates c if all disjuncts saturate c */
851 bool saturates = true;
852 for (Sequence_const_iterator si = x.sequence.begin(),
853 s_end = x.sequence.end(); si != s_end; ++si) {
854 Poly_Con_Relation relation_i = si->pointset().relation_with(c);
855 if (relation_i.implies(Poly_Con_Relation::is_included())) {
856 included_once = true;
857 }
858 else {
859 is_included = false;
860 }
861 if (relation_i.implies(Poly_Con_Relation::is_disjoint())) {
862 disjoint_once = true;
863 }
864 else {
865 is_disjoint = false;
866 }
867 if (relation_i.implies(Poly_Con_Relation::strictly_intersects())) {
868 is_strictly_intersecting = true;
869 }
870 if (!relation_i.implies(Poly_Con_Relation::saturates())) {
871 saturates = false;
872 }
873 }
874
875 Poly_Con_Relation result = Poly_Con_Relation::nothing();
876 if (is_included) {
877 result = result && Poly_Con_Relation::is_included();
878 }
879 if (is_disjoint) {
880 result = result && Poly_Con_Relation::is_disjoint();
881 }
882 if (is_strictly_intersecting || (included_once && disjoint_once)) {
883 result = result && Poly_Con_Relation::strictly_intersects();
884 }
885 if (saturates) {
886 result = result && Poly_Con_Relation::saturates();
887 }
888 return result;
889 }
890
891 template <typename PSET>
892 Poly_Gen_Relation
relation_with(const Generator & g) const893 Pointset_Powerset<PSET>::relation_with(const Generator& g) const {
894 const Pointset_Powerset& x = *this;
895
896 for (Sequence_const_iterator si = x.sequence.begin(),
897 s_end = x.sequence.end(); si != s_end; ++si) {
898 Poly_Gen_Relation relation_i = si->pointset().relation_with(g);
899 if (relation_i.implies(Poly_Gen_Relation::subsumes())) {
900 return Poly_Gen_Relation::subsumes();
901 }
902 }
903
904 return Poly_Gen_Relation::nothing();
905 }
906
907 template <typename PSET>
908 bool
909 Pointset_Powerset<PSET>
bounds_from_above(const Linear_Expression & expr) const910 ::bounds_from_above(const Linear_Expression& expr) const {
911 const Pointset_Powerset& x = *this;
912 x.omega_reduce();
913 for (Sequence_const_iterator si = x.sequence.begin(),
914 s_end = x.sequence.end(); si != s_end; ++si) {
915 if (!si->pointset().bounds_from_above(expr)) {
916 return false;
917 }
918 }
919 return true;
920 }
921
922 template <typename PSET>
923 bool
924 Pointset_Powerset<PSET>
bounds_from_below(const Linear_Expression & expr) const925 ::bounds_from_below(const Linear_Expression& expr) const {
926 const Pointset_Powerset& x = *this;
927 x.omega_reduce();
928 for (Sequence_const_iterator si = x.sequence.begin(),
929 s_end = x.sequence.end(); si != s_end; ++si) {
930 if (!si->pointset().bounds_from_below(expr)) {
931 return false;
932 }
933 }
934 return true;
935 }
936
937 template <typename PSET>
938 bool
maximize(const Linear_Expression & expr,Coefficient & sup_n,Coefficient & sup_d,bool & maximum) const939 Pointset_Powerset<PSET>::maximize(const Linear_Expression& expr,
940 Coefficient& sup_n,
941 Coefficient& sup_d,
942 bool& maximum) const {
943 const Pointset_Powerset& x = *this;
944 x.omega_reduce();
945 if (x.is_empty()) {
946 return false;
947 }
948
949 bool first = true;
950
951 PPL_DIRTY_TEMP_COEFFICIENT(best_sup_n);
952 PPL_DIRTY_TEMP_COEFFICIENT(best_sup_d);
953 best_sup_n = 0;
954 best_sup_d = 1;
955 bool best_max = false;
956
957 PPL_DIRTY_TEMP_COEFFICIENT(iter_sup_n);
958 PPL_DIRTY_TEMP_COEFFICIENT(iter_sup_d);
959 iter_sup_n = 0;
960 iter_sup_d = 1;
961 bool iter_max = false;
962
963 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
964
965 for (Sequence_const_iterator si = x.sequence.begin(),
966 s_end = x.sequence.end(); si != s_end; ++si) {
967 if (!si->pointset().maximize(expr, iter_sup_n, iter_sup_d, iter_max)) {
968 return false;
969 }
970 else
971 if (first) {
972 first = false;
973 best_sup_n = iter_sup_n;
974 best_sup_d = iter_sup_d;
975 best_max = iter_max;
976 }
977 else {
978 tmp = (best_sup_n * iter_sup_d) - (iter_sup_n * best_sup_d);
979 if (tmp < 0) {
980 best_sup_n = iter_sup_n;
981 best_sup_d = iter_sup_d;
982 best_max = iter_max;
983 }
984 else if (tmp == 0) {
985 best_max = (best_max || iter_max);
986 }
987 }
988 }
989 sup_n = best_sup_n;
990 sup_d = best_sup_d;
991 maximum = best_max;
992 return true;
993 }
994
995 template <typename PSET>
996 bool
maximize(const Linear_Expression & expr,Coefficient & sup_n,Coefficient & sup_d,bool & maximum,Generator & g) const997 Pointset_Powerset<PSET>::maximize(const Linear_Expression& expr,
998 Coefficient& sup_n,
999 Coefficient& sup_d,
1000 bool& maximum,
1001 Generator& g) const {
1002 const Pointset_Powerset& x = *this;
1003 x.omega_reduce();
1004 if (x.is_empty()) {
1005 return false;
1006 }
1007
1008 bool first = true;
1009
1010 PPL_DIRTY_TEMP_COEFFICIENT(best_sup_n);
1011 PPL_DIRTY_TEMP_COEFFICIENT(best_sup_d);
1012 best_sup_n = 0;
1013 best_sup_d = 1;
1014 bool best_max = false;
1015 Generator best_g = point();
1016
1017 PPL_DIRTY_TEMP_COEFFICIENT(iter_sup_n);
1018 PPL_DIRTY_TEMP_COEFFICIENT(iter_sup_d);
1019 iter_sup_n = 0;
1020 iter_sup_d = 1;
1021 bool iter_max = false;
1022 Generator iter_g = point();
1023
1024 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
1025
1026 for (Sequence_const_iterator si = x.sequence.begin(),
1027 s_end = x.sequence.end(); si != s_end; ++si) {
1028 if (!si->pointset().maximize(expr,
1029 iter_sup_n, iter_sup_d, iter_max, iter_g)) {
1030 return false;
1031 }
1032 else {
1033 if (first) {
1034 first = false;
1035 best_sup_n = iter_sup_n;
1036 best_sup_d = iter_sup_d;
1037 best_max = iter_max;
1038 best_g = iter_g;
1039 }
1040 else {
1041 tmp = (best_sup_n * iter_sup_d) - (iter_sup_n * best_sup_d);
1042 if (tmp < 0) {
1043 best_sup_n = iter_sup_n;
1044 best_sup_d = iter_sup_d;
1045 best_max = iter_max;
1046 best_g = iter_g;
1047 }
1048 else if (tmp == 0) {
1049 best_max = (best_max || iter_max);
1050 best_g = iter_g;
1051 }
1052 }
1053 }
1054 }
1055 sup_n = best_sup_n;
1056 sup_d = best_sup_d;
1057 maximum = best_max;
1058 g = best_g;
1059 return true;
1060 }
1061
1062 template <typename PSET>
1063 bool
minimize(const Linear_Expression & expr,Coefficient & inf_n,Coefficient & inf_d,bool & minimum) const1064 Pointset_Powerset<PSET>::minimize(const Linear_Expression& expr,
1065 Coefficient& inf_n,
1066 Coefficient& inf_d,
1067 bool& minimum) const {
1068 const Pointset_Powerset& x = *this;
1069 x.omega_reduce();
1070 if (x.is_empty()) {
1071 return false;
1072 }
1073
1074 bool first = true;
1075
1076 PPL_DIRTY_TEMP_COEFFICIENT(best_inf_n);
1077 PPL_DIRTY_TEMP_COEFFICIENT(best_inf_d);
1078 best_inf_n = 0;
1079 best_inf_d = 1;
1080 bool best_min = false;
1081
1082 PPL_DIRTY_TEMP_COEFFICIENT(iter_inf_n);
1083 PPL_DIRTY_TEMP_COEFFICIENT(iter_inf_d);
1084 iter_inf_n = 0;
1085 iter_inf_d = 1;
1086 bool iter_min = false;
1087
1088 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
1089
1090 for (Sequence_const_iterator si = x.sequence.begin(),
1091 s_end = x.sequence.end(); si != s_end; ++si) {
1092 if (!si->pointset().minimize(expr, iter_inf_n, iter_inf_d, iter_min)) {
1093 return false;
1094 }
1095 else {
1096 if (first) {
1097 first = false;
1098 best_inf_n = iter_inf_n;
1099 best_inf_d = iter_inf_d;
1100 best_min = iter_min;
1101 }
1102 else {
1103 tmp = (best_inf_n * iter_inf_d) - (iter_inf_n * best_inf_d);
1104 if (tmp > 0) {
1105 best_inf_n = iter_inf_n;
1106 best_inf_d = iter_inf_d;
1107 best_min = iter_min;
1108 }
1109 else if (tmp == 0) {
1110 best_min = (best_min || iter_min);
1111 }
1112 }
1113 }
1114 }
1115 inf_n = best_inf_n;
1116 inf_d = best_inf_d;
1117 minimum = best_min;
1118 return true;
1119 }
1120
1121 template <typename PSET>
1122 bool
minimize(const Linear_Expression & expr,Coefficient & inf_n,Coefficient & inf_d,bool & minimum,Generator & g) const1123 Pointset_Powerset<PSET>::minimize(const Linear_Expression& expr,
1124 Coefficient& inf_n,
1125 Coefficient& inf_d,
1126 bool& minimum,
1127 Generator& g) const {
1128 const Pointset_Powerset& x = *this;
1129 x.omega_reduce();
1130 if (x.is_empty()) {
1131 return false;
1132 }
1133
1134 bool first = true;
1135
1136 PPL_DIRTY_TEMP_COEFFICIENT(best_inf_n);
1137 PPL_DIRTY_TEMP_COEFFICIENT(best_inf_d);
1138 best_inf_n = 0;
1139 best_inf_d = 1;
1140 bool best_min = false;
1141 Generator best_g = point();
1142
1143 PPL_DIRTY_TEMP_COEFFICIENT(iter_inf_n);
1144 PPL_DIRTY_TEMP_COEFFICIENT(iter_inf_d);
1145 iter_inf_n = 0;
1146 iter_inf_d = 1;
1147 bool iter_min = false;
1148 Generator iter_g = point();
1149
1150 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
1151
1152 for (Sequence_const_iterator si = x.sequence.begin(),
1153 s_end = x.sequence.end(); si != s_end; ++si) {
1154 if (!si->pointset().minimize(expr,
1155 iter_inf_n, iter_inf_d, iter_min, iter_g)) {
1156 return false;
1157 }
1158 else {
1159 if (first) {
1160 first = false;
1161 best_inf_n = iter_inf_n;
1162 best_inf_d = iter_inf_d;
1163 best_min = iter_min;
1164 best_g = iter_g;
1165 }
1166 else {
1167 tmp = (best_inf_n * iter_inf_d) - (iter_inf_n * best_inf_d);
1168 if (tmp > 0) {
1169 best_inf_n = iter_inf_n;
1170 best_inf_d = iter_inf_d;
1171 best_min = iter_min;
1172 best_g = iter_g;
1173 }
1174 else if (tmp == 0) {
1175 best_min = (best_min || iter_min);
1176 best_g = iter_g;
1177 }
1178 }
1179 }
1180 }
1181 inf_n = best_inf_n;
1182 inf_d = best_inf_d;
1183 minimum = best_min;
1184 g = best_g;
1185 return true;
1186 }
1187
1188 template <typename PSET>
1189 bool
contains_integer_point() const1190 Pointset_Powerset<PSET>::contains_integer_point() const {
1191 const Pointset_Powerset& x = *this;
1192 for (Sequence_const_iterator si = x.sequence.begin(),
1193 s_end = x.sequence.end(); si != s_end; ++si) {
1194 if (si->pointset().contains_integer_point()) {
1195 return true;
1196 }
1197 }
1198 return false;
1199 }
1200
1201 template <typename PSET>
1202 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)1203 Pointset_Powerset<PSET>::wrap_assign(const Variables_Set& vars,
1204 Bounded_Integer_Type_Width w,
1205 Bounded_Integer_Type_Representation r,
1206 Bounded_Integer_Type_Overflow o,
1207 const Constraint_System* cs_p,
1208 unsigned complexity_threshold,
1209 bool wrap_individually) {
1210 Pointset_Powerset& x = *this;
1211 for (Sequence_iterator si = x.sequence.begin(),
1212 s_end = x.sequence.end(); si != s_end; ++si) {
1213 si->pointset().wrap_assign(vars, w, r, o, cs_p,
1214 complexity_threshold, wrap_individually);
1215 }
1216 x.reduced = false;
1217 PPL_ASSERT_HEAVY(x.OK());
1218 }
1219
1220 template <typename PSET>
1221 void
pairwise_reduce()1222 Pointset_Powerset<PSET>::pairwise_reduce() {
1223 Pointset_Powerset& x = *this;
1224 // It is wise to omega-reduce before pairwise-reducing.
1225 x.omega_reduce();
1226
1227 size_type n = x.size();
1228 size_type deleted;
1229 do {
1230 Pointset_Powerset new_x(x.space_dim, EMPTY);
1231 std::deque<bool> marked(n, false);
1232 deleted = 0;
1233 Sequence_iterator s_begin = x.sequence.begin();
1234 Sequence_iterator s_end = x.sequence.end();
1235 unsigned si_index = 0;
1236 for (Sequence_iterator si = s_begin; si != s_end; ++si, ++si_index) {
1237 if (marked[si_index]) {
1238 continue;
1239 }
1240 PSET& pi = si->pointset();
1241 Sequence_const_iterator sj = si;
1242 unsigned sj_index = si_index;
1243 for (++sj, ++sj_index; sj != s_end; ++sj, ++sj_index) {
1244 if (marked[sj_index]) {
1245 continue;
1246 }
1247 const PSET& pj = sj->pointset();
1248 if (pi.upper_bound_assign_if_exact(pj)) {
1249 marked[si_index] = true;
1250 marked[sj_index] = true;
1251 new_x.add_non_bottom_disjunct_preserve_reduction(pi);
1252 ++deleted;
1253 goto next;
1254 }
1255 }
1256 next:
1257 ;
1258 }
1259 iterator new_x_begin = new_x.begin();
1260 iterator new_x_end = new_x.end();
1261 unsigned xi_index = 0;
1262 for (const_iterator xi = x.begin(),
1263 x_end = x.end(); xi != x_end; ++xi, ++xi_index) {
1264 if (!marked[xi_index]) {
1265 new_x_begin
1266 = new_x.add_non_bottom_disjunct_preserve_reduction(*xi,
1267 new_x_begin,
1268 new_x_end);
1269 }
1270 }
1271 using std::swap;
1272 swap(x.sequence, new_x.sequence);
1273 n -= deleted;
1274 } while (deleted > 0);
1275 PPL_ASSERT_HEAVY(x.OK());
1276 }
1277
1278 template <typename PSET>
1279 template <typename Widening>
1280 void
1281 Pointset_Powerset<PSET>::
BGP99_heuristics_assign(const Pointset_Powerset & y,Widening widen_fun)1282 BGP99_heuristics_assign(const Pointset_Powerset& y, Widening widen_fun) {
1283 // `x' is the current iteration value.
1284 Pointset_Powerset& x = *this;
1285
1286 #ifndef NDEBUG
1287 {
1288 // We assume that `y' entails `x'.
1289 const Pointset_Powerset<PSET> x_copy = x;
1290 const Pointset_Powerset<PSET> y_copy = y;
1291 PPL_ASSERT_HEAVY(y_copy.definitely_entails(x_copy));
1292 }
1293 #endif
1294
1295 size_type n = x.size();
1296 Pointset_Powerset new_x(x.space_dim, EMPTY);
1297 std::deque<bool> marked(n, false);
1298 const_iterator x_begin = x.begin();
1299 const_iterator x_end = x.end();
1300 unsigned i_index = 0;
1301 for (const_iterator i = x_begin,
1302 y_begin = y.begin(), y_end = y.end();
1303 i != x_end; ++i, ++i_index) {
1304 for (const_iterator j = y_begin; j != y_end; ++j) {
1305 const PSET& pi = i->pointset();
1306 const PSET& pj = j->pointset();
1307 if (pi.contains(pj)) {
1308 PSET pi_copy = pi;
1309 widen_fun(pi_copy, pj);
1310 new_x.add_non_bottom_disjunct_preserve_reduction(pi_copy);
1311 marked[i_index] = true;
1312 }
1313 }
1314 }
1315 iterator new_x_begin = new_x.begin();
1316 iterator new_x_end = new_x.end();
1317 i_index = 0;
1318 for (const_iterator i = x_begin; i != x_end; ++i, ++i_index) {
1319 if (!marked[i_index]) {
1320 new_x_begin
1321 = new_x.add_non_bottom_disjunct_preserve_reduction(*i,
1322 new_x_begin,
1323 new_x_end);
1324 }
1325 }
1326 using std::swap;
1327 swap(x.sequence, new_x.sequence);
1328 PPL_ASSERT_HEAVY(x.OK());
1329 PPL_ASSERT(x.is_omega_reduced());
1330 }
1331
1332 template <typename PSET>
1333 template <typename Widening>
1334 void
1335 Pointset_Powerset<PSET>::
BGP99_extrapolation_assign(const Pointset_Powerset & y,Widening widen_fun,unsigned max_disjuncts)1336 BGP99_extrapolation_assign(const Pointset_Powerset& y,
1337 Widening widen_fun,
1338 unsigned max_disjuncts) {
1339 // `x' is the current iteration value.
1340 Pointset_Powerset& x = *this;
1341
1342 #ifndef NDEBUG
1343 {
1344 // We assume that `y' entails `x'.
1345 const Pointset_Powerset<PSET> x_copy = x;
1346 const Pointset_Powerset<PSET> y_copy = y;
1347 PPL_ASSERT_HEAVY(y_copy.definitely_entails(x_copy));
1348 }
1349 #endif
1350
1351 x.pairwise_reduce();
1352 if (max_disjuncts != 0) {
1353 x.collapse(max_disjuncts);
1354 }
1355 x.BGP99_heuristics_assign(y, widen_fun);
1356 }
1357
1358 template <typename PSET>
1359 template <typename Cert>
1360 void
1361 Pointset_Powerset<PSET>::
collect_certificates(std::map<Cert,size_type,typename Cert::Compare> & cert_ms) const1362 collect_certificates(std::map<Cert, size_type,
1363 typename Cert::Compare>& cert_ms) const {
1364 const Pointset_Powerset& x = *this;
1365 PPL_ASSERT(x.is_omega_reduced());
1366 PPL_ASSERT(cert_ms.size() == 0);
1367 for (const_iterator i = x.begin(), end = x.end(); i != end; ++i) {
1368 Cert ph_cert(i->pointset());
1369 ++cert_ms[ph_cert];
1370 }
1371 }
1372
1373 template <typename PSET>
1374 template <typename Cert>
1375 bool
1376 Pointset_Powerset<PSET>::
is_cert_multiset_stabilizing(const std::map<Cert,size_type,typename Cert::Compare> & y_cert_ms) const1377 is_cert_multiset_stabilizing(const std::map<Cert, size_type,
1378 typename Cert::Compare>& y_cert_ms) const {
1379 typedef std::map<Cert, size_type, typename Cert::Compare> Cert_Multiset;
1380 Cert_Multiset x_cert_ms;
1381 collect_certificates(x_cert_ms);
1382 typename Cert_Multiset::const_iterator xi = x_cert_ms.begin();
1383 typename Cert_Multiset::const_iterator x_cert_ms_end = x_cert_ms.end();
1384 typename Cert_Multiset::const_iterator yi = y_cert_ms.begin();
1385 typename Cert_Multiset::const_iterator y_cert_ms_end = y_cert_ms.end();
1386 while (xi != x_cert_ms_end && yi != y_cert_ms_end) {
1387 const Cert& xi_cert = xi->first;
1388 const Cert& yi_cert = yi->first;
1389 switch (xi_cert.compare(yi_cert)) {
1390 case 0:
1391 // xi_cert == yi_cert: check the number of multiset occurrences.
1392 {
1393 const size_type& xi_count = xi->second;
1394 const size_type& yi_count = yi->second;
1395 if (xi_count == yi_count) {
1396 // Same number of occurrences: compare the next pair.
1397 ++xi;
1398 ++yi;
1399 }
1400 else {
1401 // Different number of occurrences: can decide ordering.
1402 return xi_count < yi_count;
1403 }
1404 break;
1405 }
1406 case 1:
1407 // xi_cert > yi_cert: it is not stabilizing.
1408 return false;
1409
1410 case -1:
1411 // xi_cert < yi_cert: it is stabilizing.
1412 return true;
1413 }
1414 }
1415 // Here xi == x_cert_ms_end or yi == y_cert_ms_end.
1416 // Stabilization is achieved if `y_cert_ms' still has other elements.
1417 return yi != y_cert_ms_end;
1418 }
1419
1420 template <typename PSET>
1421 template <typename Cert, typename Widening>
1422 void
BHZ03_widening_assign(const Pointset_Powerset & y,Widening widen_fun)1423 Pointset_Powerset<PSET>::BHZ03_widening_assign(const Pointset_Powerset& y,
1424 Widening widen_fun) {
1425 // `x' is the current iteration value.
1426 Pointset_Powerset& x = *this;
1427
1428 #ifndef NDEBUG
1429 {
1430 // We assume that `y' entails `x'.
1431 const Pointset_Powerset<PSET> x_copy = x;
1432 const Pointset_Powerset<PSET> y_copy = y;
1433 PPL_ASSERT_HEAVY(y_copy.definitely_entails(x_copy));
1434 }
1435 #endif
1436
1437 // First widening technique: do nothing.
1438
1439 // If `y' is the empty collection, do nothing.
1440 PPL_ASSERT(x.size() > 0);
1441 if (y.size() == 0) {
1442 return;
1443 }
1444
1445 // Compute the poly-hull of `x'.
1446 PSET x_hull(x.space_dim, EMPTY);
1447 for (const_iterator i = x.begin(), x_end = x.end(); i != x_end; ++i) {
1448 x_hull.upper_bound_assign(i->pointset());
1449 }
1450
1451 // Compute the poly-hull of `y'.
1452 PSET y_hull(y.space_dim, EMPTY);
1453 for (const_iterator i = y.begin(), y_end = y.end(); i != y_end; ++i) {
1454 y_hull.upper_bound_assign(i->pointset());
1455 }
1456 // Compute the certificate for `y_hull'.
1457 const Cert y_hull_cert(y_hull);
1458
1459 // If the hull is stabilizing, do nothing.
1460 int hull_stabilization = y_hull_cert.compare(x_hull);
1461 if (hull_stabilization == 1) {
1462 return;
1463 }
1464
1465 // Multiset ordering is only useful when `y' is not a singleton.
1466 const bool y_is_not_a_singleton = y.size() > 1;
1467
1468 // The multiset certificate for `y':
1469 // we want to be lazy about its computation.
1470 typedef std::map<Cert, size_type, typename Cert::Compare> Cert_Multiset;
1471 Cert_Multiset y_cert_ms;
1472 bool y_cert_ms_computed = false;
1473
1474 if (hull_stabilization == 0 && y_is_not_a_singleton) {
1475 // Collect the multiset certificate for `y'.
1476 y.collect_certificates(y_cert_ms);
1477 y_cert_ms_computed = true;
1478 // If multiset ordering is stabilizing, do nothing.
1479 if (x.is_cert_multiset_stabilizing(y_cert_ms)) {
1480 return;
1481 }
1482 }
1483
1484 // Second widening technique: try the BGP99 powerset heuristics.
1485 Pointset_Powerset<PSET> bgp99_heuristics = x;
1486 bgp99_heuristics.BGP99_heuristics_assign(y, widen_fun);
1487
1488 // Compute the poly-hull of `bgp99_heuristics'.
1489 PSET bgp99_heuristics_hull(x.space_dim, EMPTY);
1490 for (const_iterator i = bgp99_heuristics.begin(),
1491 b_h_end = bgp99_heuristics.end(); i != b_h_end; ++i) {
1492 bgp99_heuristics_hull.upper_bound_assign(i->pointset());
1493 }
1494
1495 // Check for stabilization and, if successful,
1496 // commit to the result of the extrapolation.
1497 hull_stabilization = y_hull_cert.compare(bgp99_heuristics_hull);
1498 if (hull_stabilization == 1) {
1499 // The poly-hull is stabilizing.
1500 swap(x, bgp99_heuristics);
1501 return;
1502 }
1503 else if (hull_stabilization == 0 && y_is_not_a_singleton) {
1504 // If not already done, compute multiset certificate for `y'.
1505 if (!y_cert_ms_computed) {
1506 y.collect_certificates(y_cert_ms);
1507 y_cert_ms_computed = true;
1508 }
1509 if (bgp99_heuristics.is_cert_multiset_stabilizing(y_cert_ms)) {
1510 swap(x, bgp99_heuristics);
1511 return;
1512 }
1513 // Third widening technique: pairwise-reduction on `bgp99_heuristics'.
1514 // Note that pairwise-reduction does not affect the computation
1515 // of the poly-hulls, so that we only have to check the multiset
1516 // certificate relation.
1517 Pointset_Powerset<PSET> reduced_bgp99_heuristics(bgp99_heuristics);
1518 reduced_bgp99_heuristics.pairwise_reduce();
1519 if (reduced_bgp99_heuristics.is_cert_multiset_stabilizing(y_cert_ms)) {
1520 swap(x, reduced_bgp99_heuristics);
1521 return;
1522 }
1523 }
1524
1525 // Fourth widening technique: this is applicable only when
1526 // `y_hull' is a proper subset of `bgp99_heuristics_hull'.
1527 if (bgp99_heuristics_hull.strictly_contains(y_hull)) {
1528 // Compute (y_hull \widen bgp99_heuristics_hull).
1529 PSET ph = bgp99_heuristics_hull;
1530 widen_fun(ph, y_hull);
1531 // Compute the difference between `ph' and `bgp99_heuristics_hull'.
1532 ph.difference_assign(bgp99_heuristics_hull);
1533 x.add_disjunct(ph);
1534 return;
1535 }
1536
1537 // Fall back to the computation of the poly-hull.
1538 Pointset_Powerset<PSET> x_hull_singleton(x.space_dim, EMPTY);
1539 x_hull_singleton.add_disjunct(x_hull);
1540 swap(x, x_hull_singleton);
1541 }
1542
1543 template <typename PSET>
1544 void
ascii_dump(std::ostream & s) const1545 Pointset_Powerset<PSET>::ascii_dump(std::ostream& s) const {
1546 const Pointset_Powerset& x = *this;
1547 s << "size " << x.size()
1548 << "\nspace_dim " << x.space_dim
1549 << "\n";
1550 for (const_iterator xi = x.begin(), x_end = x.end();
1551 xi != x_end; ++xi) {
1552 xi->pointset().ascii_dump(s);
1553 }
1554 }
1555
PPL_OUTPUT_TEMPLATE_DEFINITIONS(PSET,Pointset_Powerset<PSET>)1556 PPL_OUTPUT_TEMPLATE_DEFINITIONS(PSET, Pointset_Powerset<PSET>)
1557
1558 template <typename PSET>
1559 bool
1560 Pointset_Powerset<PSET>::ascii_load(std::istream& s) {
1561 Pointset_Powerset& x = *this;
1562 std::string str;
1563
1564 if (!(s >> str) || str != "size") {
1565 return false;
1566 }
1567
1568 size_type sz;
1569
1570 if (!(s >> sz)) {
1571 return false;
1572 }
1573
1574 if (!(s >> str) || str != "space_dim") {
1575 return false;
1576 }
1577
1578 if (!(s >> x.space_dim)) {
1579 return false;
1580 }
1581
1582 Pointset_Powerset new_x(x.space_dim, EMPTY);
1583 while (sz-- > 0) {
1584 PSET ph;
1585 if (!ph.ascii_load(s)) {
1586 return false;
1587 }
1588 new_x.add_disjunct(ph);
1589 }
1590 swap(x, new_x);
1591
1592 // Check invariants.
1593 PPL_ASSERT_HEAVY(x.OK());
1594 return true;
1595 }
1596
1597 template <typename PSET>
1598 bool
OK() const1599 Pointset_Powerset<PSET>::OK() const {
1600 const Pointset_Powerset& x = *this;
1601 for (const_iterator xi = x.begin(), x_end = x.end(); xi != x_end; ++xi) {
1602 const PSET& pi = xi->pointset();
1603 if (pi.space_dimension() != x.space_dim) {
1604 #ifndef NDEBUG
1605 std::cerr << "Space dimension mismatch: is " << pi.space_dimension()
1606 << " in an element of the sequence,\nshould be "
1607 << x.space_dim << "."
1608 << std::endl;
1609 #endif
1610 return false;
1611 }
1612 }
1613 return x.Base::OK();
1614 }
1615
1616 namespace Implementation {
1617
1618 namespace Pointset_Powersets {
1619
1620 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
1621 //! Partitions polyhedron \p pset according to constraint \p c.
1622 /*! \relates Parma_Polyhedra_Library::Pointset_Powerset
1623 On exit, the intersection of \p pset and constraint \p c is stored
1624 in \p pset, whereas the intersection of \p pset with the negation of \p c
1625 is added as a new disjunct of the powerset \p r.
1626 */
1627 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
1628 template <typename PSET>
1629 void
linear_partition_aux(const Constraint & c,PSET & pset,Pointset_Powerset<NNC_Polyhedron> & r)1630 linear_partition_aux(const Constraint& c,
1631 PSET& pset,
1632 Pointset_Powerset<NNC_Polyhedron>& r) {
1633 const Linear_Expression le(c.expression());
1634 const Constraint& neg_c = c.is_strict_inequality() ? (le <= 0) : (le < 0);
1635 NNC_Polyhedron nnc_ph_pset(pset);
1636 nnc_ph_pset.add_constraint(neg_c);
1637 if (!nnc_ph_pset.is_empty()) {
1638 r.add_disjunct(nnc_ph_pset);
1639 }
1640 pset.add_constraint(c);
1641 }
1642
1643 } // namespace Pointset_Powersets
1644
1645 } // namespace Implementation
1646
1647
1648 /*! \relates Pointset_Powerset */
1649 template <typename PSET>
1650 std::pair<PSET, Pointset_Powerset<NNC_Polyhedron> >
linear_partition(const PSET & p,const PSET & q)1651 linear_partition(const PSET& p, const PSET& q) {
1652 using Implementation::Pointset_Powersets::linear_partition_aux;
1653
1654 Pointset_Powerset<NNC_Polyhedron> r(p.space_dimension(), EMPTY);
1655 PSET pset = q;
1656 const Constraint_System& p_constraints = p.constraints();
1657 for (Constraint_System::const_iterator i = p_constraints.begin(),
1658 p_constraints_end = p_constraints.end();
1659 i != p_constraints_end;
1660 ++i) {
1661 const Constraint& c = *i;
1662 if (c.is_equality()) {
1663 const Linear_Expression le(c.expression());
1664 linear_partition_aux(le <= 0, pset, r);
1665 linear_partition_aux(le >= 0, pset, r);
1666 }
1667 else {
1668 linear_partition_aux(c, pset, r);
1669 }
1670 }
1671 return std::make_pair(pset, r);
1672 }
1673
1674 } // namespace Parma_Polyhedra_Library
1675
1676 #endif // !defined(PPL_Pointset_Powerset_templates_hh)
1677