1 /* Box 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_Box_templates_hh
25 #define PPL_Box_templates_hh 1
26
27 #include "Variables_Set_defs.hh"
28 #include "Constraint_System_defs.hh"
29 #include "Constraint_System_inlines.hh"
30 #include "Generator_System_defs.hh"
31 #include "Generator_System_inlines.hh"
32 #include "Poly_Con_Relation_defs.hh"
33 #include "Poly_Gen_Relation_defs.hh"
34 #include "Polyhedron_defs.hh"
35 #include "Grid_defs.hh"
36 #include "Interval_defs.hh"
37 #include "Linear_Form_defs.hh"
38 #include "BD_Shape_defs.hh"
39 #include "Octagonal_Shape_defs.hh"
40 #include "MIP_Problem_defs.hh"
41 #include "Rational_Interval.hh"
42 #include <vector>
43 #include <map>
44 #include <iostream>
45
46 namespace Parma_Polyhedra_Library {
47
48 template <typename ITV>
49 inline
Box(dimension_type num_dimensions,Degenerate_Element kind)50 Box<ITV>::Box(dimension_type num_dimensions, Degenerate_Element kind)
51 : seq(check_space_dimension_overflow(num_dimensions,
52 max_space_dimension(),
53 "PPL::Box::",
54 "Box(n, k)",
55 "n exceeds the maximum "
56 "allowed space dimension")),
57 status() {
58 // In a box that is marked empty the intervals are completely
59 // meaningless: we exploit this by avoiding their initialization.
60 if (kind == UNIVERSE) {
61 for (dimension_type i = num_dimensions; i-- > 0; ) {
62 seq[i].assign(UNIVERSE);
63 }
64 set_empty_up_to_date();
65 }
66 else {
67 set_empty();
68 }
69 PPL_ASSERT(OK());
70 }
71
72 template <typename ITV>
73 inline
Box(const Constraint_System & cs)74 Box<ITV>::Box(const Constraint_System& cs)
75 : seq(check_space_dimension_overflow(cs.space_dimension(),
76 max_space_dimension(),
77 "PPL::Box::",
78 "Box(cs)",
79 "cs exceeds the maximum "
80 "allowed space dimension")),
81 status() {
82 // FIXME: check whether we can avoid the double initialization.
83 for (dimension_type i = cs.space_dimension(); i-- > 0; ) {
84 seq[i].assign(UNIVERSE);
85 }
86 add_constraints_no_check(cs);
87 }
88
89 template <typename ITV>
90 inline
Box(const Congruence_System & cgs)91 Box<ITV>::Box(const Congruence_System& cgs)
92 : seq(check_space_dimension_overflow(cgs.space_dimension(),
93 max_space_dimension(),
94 "PPL::Box::",
95 "Box(cgs)",
96 "cgs exceeds the maximum "
97 "allowed space dimension")),
98 status() {
99 // FIXME: check whether we can avoid the double initialization.
100 for (dimension_type i = cgs.space_dimension(); i-- > 0; ) {
101 seq[i].assign(UNIVERSE);
102 }
103 add_congruences_no_check(cgs);
104 }
105
106 template <typename ITV>
107 template <typename Other_ITV>
108 inline
Box(const Box<Other_ITV> & y,Complexity_Class)109 Box<ITV>::Box(const Box<Other_ITV>& y, Complexity_Class)
110 : seq(y.space_dimension()),
111 // FIXME: why the following does not work?
112 // status(y.status) {
113 status() {
114 // FIXME: remove when the above is fixed.
115 if (y.marked_empty()) {
116 set_empty();
117 }
118
119 if (!y.marked_empty()) {
120 for (dimension_type k = y.space_dimension(); k-- > 0; ) {
121 seq[k].assign(y.seq[k]);
122 }
123 }
124 PPL_ASSERT(OK());
125 }
126
127 template <typename ITV>
Box(const Generator_System & gs)128 Box<ITV>::Box(const Generator_System& gs)
129 : seq(check_space_dimension_overflow(gs.space_dimension(),
130 max_space_dimension(),
131 "PPL::Box::",
132 "Box(gs)",
133 "gs exceeds the maximum "
134 "allowed space dimension")),
135 status() {
136 const Generator_System::const_iterator gs_begin = gs.begin();
137 const Generator_System::const_iterator gs_end = gs.end();
138 if (gs_begin == gs_end) {
139 // An empty generator system defines the empty box.
140 set_empty();
141 return;
142 }
143
144 // The empty flag will be meaningful, whatever happens from now on.
145 set_empty_up_to_date();
146
147 const dimension_type space_dim = space_dimension();
148 PPL_DIRTY_TEMP(mpq_class, q);
149 bool point_seen = false;
150 // Going through all the points.
151 for (Generator_System::const_iterator
152 gs_i = gs_begin; gs_i != gs_end; ++gs_i) {
153 const Generator& g = *gs_i;
154 if (g.is_point()) {
155 const Coefficient& d = g.divisor();
156 if (point_seen) {
157 // This is not the first point: `seq' already contains valid values.
158 // TODO: If the variables in the expression that have coefficient 0
159 // have no effect on seq[i], this loop can be optimized using
160 // Generator::expr_type::const_iterator.
161 for (dimension_type i = space_dim; i-- > 0; ) {
162 assign_r(q.get_num(), g.coefficient(Variable(i)), ROUND_NOT_NEEDED);
163 assign_r(q.get_den(), d, ROUND_NOT_NEEDED);
164 q.canonicalize();
165 PPL_DIRTY_TEMP(ITV, iq);
166 iq.build(i_constraint(EQUAL, q));
167 seq[i].join_assign(iq);
168 }
169 }
170 else {
171 // This is the first point seen: initialize `seq'.
172 point_seen = true;
173 // TODO: If the variables in the expression that have coefficient 0
174 // have no effect on seq[i], this loop can be optimized using
175 // Generator::expr_type::const_iterator.
176 for (dimension_type i = space_dim; i-- > 0; ) {
177 assign_r(q.get_num(), g.coefficient(Variable(i)), ROUND_NOT_NEEDED);
178 assign_r(q.get_den(), d, ROUND_NOT_NEEDED);
179 q.canonicalize();
180 seq[i].build(i_constraint(EQUAL, q));
181 }
182 }
183 }
184 }
185
186 if (!point_seen) {
187 // The generator system is not empty, but contains no points.
188 throw std::invalid_argument("PPL::Box<ITV>::Box(gs):\n"
189 "the non-empty generator system gs "
190 "contains no points.");
191 }
192
193 // Going through all the lines, rays and closure points.
194 for (Generator_System::const_iterator gs_i = gs_begin;
195 gs_i != gs_end; ++gs_i) {
196 const Generator& g = *gs_i;
197 switch (g.type()) {
198 case Generator::LINE:
199 for (Generator::expr_type::const_iterator i = g.expression().begin(),
200 i_end = g.expression().end();
201 i != i_end; ++i) {
202 seq[i.variable().id()].assign(UNIVERSE);
203 }
204 break;
205 case Generator::RAY:
206 for (Generator::expr_type::const_iterator i = g.expression().begin(),
207 i_end = g.expression().end();
208 i != i_end; ++i) {
209 switch (sgn(*i)) {
210 case 1:
211 seq[i.variable().id()].upper_extend();
212 break;
213 case -1:
214 seq[i.variable().id()].lower_extend();
215 break;
216 default:
217 PPL_UNREACHABLE;
218 break;
219 }
220 }
221 break;
222 case Generator::CLOSURE_POINT:
223 {
224 const Coefficient& d = g.divisor();
225 // TODO: If the variables in the expression that have coefficient 0
226 // have no effect on seq[i], this loop can be optimized using
227 // Generator::expr_type::const_iterator.
228 for (dimension_type i = space_dim; i-- > 0; ) {
229 assign_r(q.get_num(), g.coefficient(Variable(i)), ROUND_NOT_NEEDED);
230 assign_r(q.get_den(), d, ROUND_NOT_NEEDED);
231 q.canonicalize();
232 ITV& seq_i = seq[i];
233 seq_i.lower_extend(i_constraint(GREATER_THAN, q));
234 seq_i.upper_extend(i_constraint(LESS_THAN, q));
235 }
236 }
237 break;
238 default:
239 // Points already dealt with.
240 break;
241 }
242 }
243 PPL_ASSERT(OK());
244 }
245
246 template <typename ITV>
247 template <typename T>
Box(const BD_Shape<T> & bds,Complexity_Class)248 Box<ITV>::Box(const BD_Shape<T>& bds, Complexity_Class)
249 : seq(check_space_dimension_overflow(bds.space_dimension(),
250 max_space_dimension(),
251 "PPL::Box::",
252 "Box(bds)",
253 "bds exceeds the maximum "
254 "allowed space dimension")),
255 status() {
256 // Expose all the interval constraints.
257 bds.shortest_path_closure_assign();
258 if (bds.marked_empty()) {
259 set_empty();
260 PPL_ASSERT(OK());
261 return;
262 }
263
264 // The empty flag will be meaningful, whatever happens from now on.
265 set_empty_up_to_date();
266
267 const dimension_type space_dim = space_dimension();
268 if (space_dim == 0) {
269 PPL_ASSERT(OK());
270 return;
271 }
272
273 typedef typename BD_Shape<T>::coefficient_type Coeff;
274 PPL_DIRTY_TEMP(Coeff, tmp);
275 const DB_Row<Coeff>& dbm_0 = bds.dbm[0];
276 for (dimension_type i = space_dim; i-- > 0; ) {
277 I_Constraint<Coeff> lower;
278 I_Constraint<Coeff> upper;
279 ITV& seq_i = seq[i];
280
281 // Set the upper bound.
282 const Coeff& u = dbm_0[i+1];
283 if (!is_plus_infinity(u)) {
284 upper.set(LESS_OR_EQUAL, u, true);
285 }
286
287 // Set the lower bound.
288 const Coeff& negated_l = bds.dbm[i+1][0];
289 if (!is_plus_infinity(negated_l)) {
290 neg_assign_r(tmp, negated_l, ROUND_DOWN);
291 lower.set(GREATER_OR_EQUAL, tmp);
292 }
293
294 seq_i.build(lower, upper);
295 }
296 PPL_ASSERT(OK());
297 }
298
299 template <typename ITV>
300 template <typename T>
Box(const Octagonal_Shape<T> & oct,Complexity_Class)301 Box<ITV>::Box(const Octagonal_Shape<T>& oct, Complexity_Class)
302 : seq(check_space_dimension_overflow(oct.space_dimension(),
303 max_space_dimension(),
304 "PPL::Box::",
305 "Box(oct)",
306 "oct exceeds the maximum "
307 "allowed space dimension")),
308 status() {
309 // Expose all the interval constraints.
310 oct.strong_closure_assign();
311 if (oct.marked_empty()) {
312 set_empty();
313 return;
314 }
315
316 // The empty flag will be meaningful, whatever happens from now on.
317 set_empty_up_to_date();
318
319 const dimension_type space_dim = space_dimension();
320 if (space_dim == 0) {
321 return;
322 }
323
324 PPL_DIRTY_TEMP(mpq_class, lower_bound);
325 PPL_DIRTY_TEMP(mpq_class, upper_bound);
326 for (dimension_type i = space_dim; i-- > 0; ) {
327 typedef typename Octagonal_Shape<T>::coefficient_type Coeff;
328 I_Constraint<mpq_class> lower;
329 I_Constraint<mpq_class> upper;
330 ITV& seq_i = seq[i];
331 const dimension_type ii = 2*i;
332 const dimension_type cii = ii + 1;
333
334 // Set the upper bound.
335 const Coeff& twice_ub = oct.matrix[cii][ii];
336 if (!is_plus_infinity(twice_ub)) {
337 assign_r(upper_bound, twice_ub, ROUND_NOT_NEEDED);
338 div_2exp_assign_r(upper_bound, upper_bound, 1, ROUND_NOT_NEEDED);
339 upper.set(LESS_OR_EQUAL, upper_bound);
340 }
341
342 // Set the lower bound.
343 const Coeff& twice_lb = oct.matrix[ii][cii];
344 if (!is_plus_infinity(twice_lb)) {
345 assign_r(lower_bound, twice_lb, ROUND_NOT_NEEDED);
346 neg_assign_r(lower_bound, lower_bound, ROUND_NOT_NEEDED);
347 div_2exp_assign_r(lower_bound, lower_bound, 1, ROUND_NOT_NEEDED);
348 lower.set(GREATER_OR_EQUAL, lower_bound);
349 }
350 seq_i.build(lower, upper);
351 }
352 }
353
354 template <typename ITV>
Box(const Polyhedron & ph,Complexity_Class complexity)355 Box<ITV>::Box(const Polyhedron& ph, Complexity_Class complexity)
356 : seq(check_space_dimension_overflow(ph.space_dimension(),
357 max_space_dimension(),
358 "PPL::Box::",
359 "Box(ph)",
360 "ph exceeds the maximum "
361 "allowed space dimension")),
362 status() {
363 // The empty flag will be meaningful, whatever happens from now on.
364 set_empty_up_to_date();
365
366 // We do not need to bother about `complexity' if:
367 // a) the polyhedron is already marked empty; or ...
368 if (ph.marked_empty()) {
369 set_empty();
370 return;
371 }
372
373 // b) the polyhedron is zero-dimensional; or ...
374 const dimension_type space_dim = ph.space_dimension();
375 if (space_dim == 0) {
376 return;
377 }
378
379 // c) the polyhedron is already described by a generator system.
380 if (ph.generators_are_up_to_date() && !ph.has_pending_constraints()) {
381 Box tmp(ph.generators());
382 m_swap(tmp);
383 return;
384 }
385
386 // Here generators are not up-to-date or there are pending constraints.
387 PPL_ASSERT(ph.constraints_are_up_to_date());
388
389 if (complexity == POLYNOMIAL_COMPLEXITY) {
390 // FIXME: is there a way to avoid this initialization?
391 for (dimension_type i = space_dim; i-- > 0; ) {
392 seq[i].assign(UNIVERSE);
393 }
394 // Get a simplified version of the constraints.
395 const Constraint_System cs = ph.simplified_constraints();
396 // Propagate easy-to-find bounds from the constraints,
397 // allowing for a limited number of iterations.
398 // FIXME: 20 is just a wild guess.
399 const dimension_type max_iterations = 20;
400 propagate_constraints_no_check(cs, max_iterations);
401 }
402 else if (complexity == SIMPLEX_COMPLEXITY) {
403 MIP_Problem lp(space_dim);
404 const Constraint_System& ph_cs = ph.constraints();
405 if (!ph_cs.has_strict_inequalities()) {
406 lp.add_constraints(ph_cs);
407 }
408 else {
409 // Adding to `lp' a topologically closed version of `ph_cs'.
410 for (Constraint_System::const_iterator i = ph_cs.begin(),
411 ph_cs_end = ph_cs.end(); i != ph_cs_end; ++i) {
412 const Constraint& c = *i;
413 if (c.is_strict_inequality()) {
414 const Linear_Expression expr(c.expression());
415 lp.add_constraint(expr >= 0);
416 }
417 else {
418 lp.add_constraint(c);
419 }
420 }
421 }
422 // Check for unsatisfiability.
423 if (!lp.is_satisfiable()) {
424 set_empty();
425 return;
426 }
427 // Get all the bounds for the space dimensions.
428 Generator g(point());
429 PPL_DIRTY_TEMP(mpq_class, lower_bound);
430 PPL_DIRTY_TEMP(mpq_class, upper_bound);
431 PPL_DIRTY_TEMP_COEFFICIENT(bound_numer);
432 PPL_DIRTY_TEMP_COEFFICIENT(bound_denom);
433 for (dimension_type i = space_dim; i-- > 0; ) {
434 I_Constraint<mpq_class> lower;
435 I_Constraint<mpq_class> upper;
436 ITV& seq_i = seq[i];
437 lp.set_objective_function(Variable(i));
438 // Evaluate upper bound.
439 lp.set_optimization_mode(MAXIMIZATION);
440 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
441 g = lp.optimizing_point();
442 lp.evaluate_objective_function(g, bound_numer, bound_denom);
443 assign_r(upper_bound.get_num(), bound_numer, ROUND_NOT_NEEDED);
444 assign_r(upper_bound.get_den(), bound_denom, ROUND_NOT_NEEDED);
445 PPL_ASSERT(is_canonical(upper_bound));
446 upper.set(LESS_OR_EQUAL, upper_bound);
447 }
448 // Evaluate optimal lower bound.
449 lp.set_optimization_mode(MINIMIZATION);
450 if (lp.solve() == OPTIMIZED_MIP_PROBLEM) {
451 g = lp.optimizing_point();
452 lp.evaluate_objective_function(g, bound_numer, bound_denom);
453 assign_r(lower_bound.get_num(), bound_numer, ROUND_NOT_NEEDED);
454 assign_r(lower_bound.get_den(), bound_denom, ROUND_NOT_NEEDED);
455 PPL_ASSERT(is_canonical(lower_bound));
456 lower.set(GREATER_OR_EQUAL, lower_bound);
457 }
458 seq_i.build(lower, upper);
459 }
460 }
461 else {
462 PPL_ASSERT(complexity == ANY_COMPLEXITY);
463 if (ph.is_empty()) {
464 set_empty();
465 }
466 else {
467 Box tmp(ph.generators());
468 m_swap(tmp);
469 }
470 }
471 }
472
473 template <typename ITV>
Box(const Grid & gr,Complexity_Class)474 Box<ITV>::Box(const Grid& gr, Complexity_Class)
475 : seq(check_space_dimension_overflow(gr.space_dimension(),
476 max_space_dimension(),
477 "PPL::Box::",
478 "Box(gr)",
479 "gr exceeds the maximum "
480 "allowed space dimension")),
481 status() {
482
483 if (gr.marked_empty()) {
484 set_empty();
485 return;
486 }
487
488 // The empty flag will be meaningful, whatever happens from now on.
489 set_empty_up_to_date();
490
491 const dimension_type space_dim = gr.space_dimension();
492
493 if (space_dim == 0) {
494 return;
495 }
496
497 if (!gr.generators_are_up_to_date() && !gr.update_generators()) {
498 // Updating found the grid empty.
499 set_empty();
500 return;
501 }
502
503 PPL_ASSERT(!gr.gen_sys.empty());
504
505 // For each dimension that is bounded by the grid, set both bounds
506 // of the interval to the value of the associated coefficient in a
507 // generator point.
508 PPL_DIRTY_TEMP(mpq_class, bound);
509 PPL_DIRTY_TEMP_COEFFICIENT(bound_numer);
510 PPL_DIRTY_TEMP_COEFFICIENT(bound_denom);
511 for (dimension_type i = space_dim; i-- > 0; ) {
512 ITV& seq_i = seq[i];
513 Variable var(i);
514 bool max;
515 if (gr.maximize(var, bound_numer, bound_denom, max)) {
516 assign_r(bound.get_num(), bound_numer, ROUND_NOT_NEEDED);
517 assign_r(bound.get_den(), bound_denom, ROUND_NOT_NEEDED);
518 bound.canonicalize();
519 seq_i.build(i_constraint(EQUAL, bound));
520 }
521 else {
522 seq_i.assign(UNIVERSE);
523 }
524 }
525 }
526
527 template <typename ITV>
528 template <typename D1, typename D2, typename R>
Box(const Partially_Reduced_Product<D1,D2,R> & dp,Complexity_Class complexity)529 Box<ITV>::Box(const Partially_Reduced_Product<D1, D2, R>& dp,
530 Complexity_Class complexity)
531 : seq(), status() {
532 check_space_dimension_overflow(dp.space_dimension(),
533 max_space_dimension(),
534 "PPL::Box::",
535 "Box(dp)",
536 "dp exceeds the maximum "
537 "allowed space dimension");
538 Box tmp1(dp.domain1(), complexity);
539 Box tmp2(dp.domain2(), complexity);
540 tmp1.intersection_assign(tmp2);
541 m_swap(tmp1);
542 }
543
544 template <typename ITV>
545 inline void
add_space_dimensions_and_embed(const dimension_type m)546 Box<ITV>::add_space_dimensions_and_embed(const dimension_type m) {
547 // Adding no dimensions is a no-op.
548 if (m == 0) {
549 return;
550 }
551 check_space_dimension_overflow(m, max_space_dimension() - space_dimension(),
552 "PPL::Box::",
553 "add_space_dimensions_and_embed(m)",
554 "adding m new space dimensions exceeds "
555 "the maximum allowed space dimension");
556 // To embed an n-dimension space box in a (n+m)-dimension space,
557 // we just add `m' new universe elements to the sequence.
558 seq.insert(seq.end(), m, ITV(UNIVERSE));
559 PPL_ASSERT(OK());
560 }
561
562 template <typename ITV>
563 inline void
add_space_dimensions_and_project(const dimension_type m)564 Box<ITV>::add_space_dimensions_and_project(const dimension_type m) {
565 // Adding no dimensions is a no-op.
566 if (m == 0) {
567 return;
568 }
569 check_space_dimension_overflow(m, max_space_dimension() - space_dimension(),
570 "PPL::Box::",
571 "add_space_dimensions_and_project(m)",
572 "adding m new space dimensions exceeds "
573 "the maximum allowed space dimension");
574 // Add `m' new zero elements to the sequence.
575 seq.insert(seq.end(), m, ITV(0));
576 PPL_ASSERT(OK());
577 }
578
579 template <typename ITV>
580 bool
operator ==(const Box<ITV> & x,const Box<ITV> & y)581 operator==(const Box<ITV>& x, const Box<ITV>& y) {
582 const dimension_type x_space_dim = x.space_dimension();
583 if (x_space_dim != y.space_dimension()) {
584 return false;
585 }
586
587 if (x.is_empty()) {
588 return y.is_empty();
589 }
590
591 if (y.is_empty()) {
592 return x.is_empty();
593 }
594
595 for (dimension_type k = x_space_dim; k-- > 0; ) {
596 if (x.seq[k] != y.seq[k]) {
597 return false;
598 }
599 }
600 return true;
601 }
602
603 template <typename ITV>
604 bool
bounds(const Linear_Expression & expr,const bool from_above) const605 Box<ITV>::bounds(const Linear_Expression& expr, const bool from_above) const {
606 // `expr' should be dimension-compatible with `*this'.
607 const dimension_type expr_space_dim = expr.space_dimension();
608 const dimension_type space_dim = space_dimension();
609 if (space_dim < expr_space_dim) {
610 throw_dimension_incompatible((from_above
611 ? "bounds_from_above(e)"
612 : "bounds_from_below(e)"), "e", expr);
613 }
614 // A zero-dimensional or empty Box bounds everything.
615 if (space_dim == 0 || is_empty()) {
616 return true;
617 }
618 const int from_above_sign = from_above ? 1 : -1;
619 // TODO: This loop can be optimized more, if needed, exploiting the
620 // (possible) sparseness of expr.
621 for (Linear_Expression::const_iterator i = expr.begin(),
622 i_end = expr.end(); i != i_end; ++i) {
623 const Variable v = i.variable();
624 switch (sgn(*i) * from_above_sign) {
625 case 1:
626 if (seq[v.id()].upper_is_boundary_infinity()) {
627 return false;
628 }
629 break;
630 case 0:
631 PPL_UNREACHABLE;
632 break;
633 case -1:
634 if (seq[v.id()].lower_is_boundary_infinity()) {
635 return false;
636 }
637 break;
638 }
639 }
640 return true;
641 }
642
643 template <typename ITV>
644 Poly_Con_Relation
interval_relation(const ITV & i,const Constraint::Type constraint_type,Coefficient_traits::const_reference numer,Coefficient_traits::const_reference denom)645 interval_relation(const ITV& i,
646 const Constraint::Type constraint_type,
647 Coefficient_traits::const_reference numer,
648 Coefficient_traits::const_reference denom) {
649
650 if (i.is_universe()) {
651 return Poly_Con_Relation::strictly_intersects();
652 }
653
654 PPL_DIRTY_TEMP(mpq_class, bound);
655 assign_r(bound.get_num(), numer, ROUND_NOT_NEEDED);
656 assign_r(bound.get_den(), denom, ROUND_NOT_NEEDED);
657 bound.canonicalize();
658 neg_assign_r(bound, bound, ROUND_NOT_NEEDED);
659 const bool is_lower_bound = (denom > 0);
660
661 PPL_DIRTY_TEMP(mpq_class, bound_diff);
662 if (constraint_type == Constraint::EQUALITY) {
663 if (i.lower_is_boundary_infinity()) {
664 PPL_ASSERT(!i.upper_is_boundary_infinity());
665 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
666 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
667 switch (sgn(bound_diff)) {
668 case 1:
669 return Poly_Con_Relation::strictly_intersects();
670 case 0:
671 return i.upper_is_open()
672 ? Poly_Con_Relation::is_disjoint()
673 : Poly_Con_Relation::strictly_intersects();
674 case -1:
675 return Poly_Con_Relation::is_disjoint();
676 }
677 }
678 else {
679 assign_r(bound_diff, i.lower(), ROUND_NOT_NEEDED);
680 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
681 switch (sgn(bound_diff)) {
682 case 1:
683 return Poly_Con_Relation::is_disjoint();
684 case 0:
685 if (i.lower_is_open()) {
686 return Poly_Con_Relation::is_disjoint();
687 }
688 if (i.is_singleton()) {
689 return Poly_Con_Relation::is_included()
690 && Poly_Con_Relation::saturates();
691 }
692 return Poly_Con_Relation::strictly_intersects();
693 case -1:
694 if (i.upper_is_boundary_infinity()) {
695 return Poly_Con_Relation::strictly_intersects();
696 }
697 else {
698 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
699 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
700 switch (sgn(bound_diff)) {
701 case 1:
702 return Poly_Con_Relation::strictly_intersects();
703 case 0:
704 if (i.upper_is_open()) {
705 return Poly_Con_Relation::is_disjoint();
706 }
707 else {
708 return Poly_Con_Relation::strictly_intersects();
709 }
710 case -1:
711 return Poly_Con_Relation::is_disjoint();
712 }
713 }
714 }
715 }
716 }
717
718 PPL_ASSERT(constraint_type != Constraint::EQUALITY);
719 if (is_lower_bound) {
720 if (i.lower_is_boundary_infinity()) {
721 PPL_ASSERT(!i.upper_is_boundary_infinity());
722 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
723 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
724 switch (sgn(bound_diff)) {
725 case 1:
726 return Poly_Con_Relation::strictly_intersects();
727 case 0:
728 if (constraint_type == Constraint::STRICT_INEQUALITY
729 || i.upper_is_open()) {
730 return Poly_Con_Relation::is_disjoint();
731 }
732 else {
733 return Poly_Con_Relation::strictly_intersects();
734 }
735 case -1:
736 return Poly_Con_Relation::is_disjoint();
737 }
738 }
739 else {
740 assign_r(bound_diff, i.lower(), ROUND_NOT_NEEDED);
741 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
742 switch (sgn(bound_diff)) {
743 case 1:
744 return Poly_Con_Relation::is_included();
745 case 0:
746 if (constraint_type == Constraint::NONSTRICT_INEQUALITY
747 || i.lower_is_open()) {
748 Poly_Con_Relation result = Poly_Con_Relation::is_included();
749 if (i.is_singleton()) {
750 result = result && Poly_Con_Relation::saturates();
751 }
752 return result;
753 }
754 else {
755 PPL_ASSERT(constraint_type == Constraint::STRICT_INEQUALITY
756 && !i.lower_is_open());
757 if (i.is_singleton()) {
758 return Poly_Con_Relation::is_disjoint()
759 && Poly_Con_Relation::saturates();
760 }
761 else {
762 return Poly_Con_Relation::strictly_intersects();
763 }
764 }
765 case -1:
766 if (i.upper_is_boundary_infinity()) {
767 return Poly_Con_Relation::strictly_intersects();
768 }
769 else {
770 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
771 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
772 switch (sgn(bound_diff)) {
773 case 1:
774 return Poly_Con_Relation::strictly_intersects();
775 case 0:
776 if (constraint_type == Constraint::STRICT_INEQUALITY
777 || i.upper_is_open()) {
778 return Poly_Con_Relation::is_disjoint();
779 }
780 else {
781 return Poly_Con_Relation::strictly_intersects();
782 }
783 case -1:
784 return Poly_Con_Relation::is_disjoint();
785 }
786 }
787 }
788 }
789 }
790 else {
791 // `c' is an upper bound.
792 if (i.upper_is_boundary_infinity()) {
793 return Poly_Con_Relation::strictly_intersects();
794 }
795 else {
796 assign_r(bound_diff, i.upper(), ROUND_NOT_NEEDED);
797 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
798 switch (sgn(bound_diff)) {
799 case -1:
800 return Poly_Con_Relation::is_included();
801 case 0:
802 if (constraint_type == Constraint::NONSTRICT_INEQUALITY
803 || i.upper_is_open()) {
804 Poly_Con_Relation result = Poly_Con_Relation::is_included();
805 if (i.is_singleton()) {
806 result = result && Poly_Con_Relation::saturates();
807 }
808 return result;
809 }
810 else {
811 PPL_ASSERT(constraint_type == Constraint::STRICT_INEQUALITY
812 && !i.upper_is_open());
813 if (i.is_singleton()) {
814 return Poly_Con_Relation::is_disjoint()
815 && Poly_Con_Relation::saturates();
816 }
817 else {
818 return Poly_Con_Relation::strictly_intersects();
819 }
820 }
821 case 1:
822 if (i.lower_is_boundary_infinity()) {
823 return Poly_Con_Relation::strictly_intersects();
824 }
825 else {
826 assign_r(bound_diff, i.lower(), ROUND_NOT_NEEDED);
827 sub_assign_r(bound_diff, bound_diff, bound, ROUND_NOT_NEEDED);
828 switch (sgn(bound_diff)) {
829 case -1:
830 return Poly_Con_Relation::strictly_intersects();
831 case 0:
832 if (constraint_type == Constraint::STRICT_INEQUALITY
833 || i.lower_is_open()) {
834 return Poly_Con_Relation::is_disjoint();
835 }
836 else {
837 return Poly_Con_Relation::strictly_intersects();
838 }
839 case 1:
840 return Poly_Con_Relation::is_disjoint();
841 }
842 }
843 }
844 }
845 }
846
847 // Quiet a compiler warning: this program point is unreachable.
848 PPL_UNREACHABLE;
849 return Poly_Con_Relation::nothing();
850 }
851
852 template <typename ITV>
853 Poly_Con_Relation
relation_with(const Congruence & cg) const854 Box<ITV>::relation_with(const Congruence& cg) const {
855 const dimension_type cg_space_dim = cg.space_dimension();
856 const dimension_type space_dim = space_dimension();
857
858 // Dimension-compatibility check.
859 if (cg_space_dim > space_dim) {
860 throw_dimension_incompatible("relation_with(cg)", cg);
861 }
862 if (is_empty()) {
863 return Poly_Con_Relation::saturates()
864 && Poly_Con_Relation::is_included()
865 && Poly_Con_Relation::is_disjoint();
866 }
867
868 if (space_dim == 0) {
869 if (cg.is_inconsistent()) {
870 return Poly_Con_Relation::is_disjoint();
871 }
872 else {
873 return Poly_Con_Relation::saturates()
874 && Poly_Con_Relation::is_included();
875 }
876 }
877
878 if (cg.is_equality()) {
879 const Constraint c(cg);
880 return relation_with(c);
881 }
882
883 PPL_DIRTY_TEMP(Rational_Interval, r);
884 PPL_DIRTY_TEMP(Rational_Interval, t);
885 PPL_DIRTY_TEMP(mpq_class, m);
886 r = 0;
887 for (Congruence::expr_type::const_iterator i = cg.expression().begin(),
888 i_end = cg.expression().end(); i != i_end; ++i) {
889 const Coefficient& cg_i = *i;
890 const Variable v = i.variable();
891 assign_r(m, cg_i, ROUND_NOT_NEEDED);
892 // FIXME: an add_mul_assign() method would come handy here.
893 t.build(seq[v.id()].lower_constraint(), seq[v.id()].upper_constraint());
894 t *= m;
895 r += t;
896 }
897
898 if (r.lower_is_boundary_infinity() || r.upper_is_boundary_infinity()) {
899 return Poly_Con_Relation::strictly_intersects();
900 }
901
902 // Find the value that satisfies the congruence and is
903 // nearest to the lower bound such that the point lies on or above it.
904
905 PPL_DIRTY_TEMP_COEFFICIENT(lower);
906 PPL_DIRTY_TEMP_COEFFICIENT(mod);
907 PPL_DIRTY_TEMP_COEFFICIENT(v);
908 mod = cg.modulus();
909 v = cg.inhomogeneous_term() % mod;
910 assign_r(lower, r.lower(), ROUND_DOWN);
911 v -= ((lower / mod) * mod);
912 if (v + lower > 0) {
913 v -= mod;
914 }
915 return interval_relation(r, Constraint::EQUALITY, v);
916 }
917
918 template <typename ITV>
919 Poly_Con_Relation
relation_with(const Constraint & c) const920 Box<ITV>::relation_with(const Constraint& c) const {
921 const dimension_type c_space_dim = c.space_dimension();
922 const dimension_type space_dim = space_dimension();
923
924 // Dimension-compatibility check.
925 if (c_space_dim > space_dim) {
926 throw_dimension_incompatible("relation_with(c)", c);
927 }
928
929 if (is_empty()) {
930 return Poly_Con_Relation::saturates()
931 && Poly_Con_Relation::is_included()
932 && Poly_Con_Relation::is_disjoint();
933 }
934
935 if (space_dim == 0) {
936 if ((c.is_equality() && c.inhomogeneous_term() != 0)
937 || (c.is_inequality() && c.inhomogeneous_term() < 0)) {
938 return Poly_Con_Relation::is_disjoint();
939 }
940 else if (c.is_strict_inequality() && c.inhomogeneous_term() == 0) {
941 // The constraint 0 > 0 implicitly defines the hyperplane 0 = 0;
942 // thus, the zero-dimensional point also saturates it.
943 return Poly_Con_Relation::saturates()
944 && Poly_Con_Relation::is_disjoint();
945 }
946 else if (c.is_equality() || c.inhomogeneous_term() == 0) {
947 return Poly_Con_Relation::saturates()
948 && Poly_Con_Relation::is_included();
949 }
950 else {
951 // The zero-dimensional point saturates
952 // neither the positivity constraint 1 >= 0,
953 // nor the strict positivity constraint 1 > 0.
954 return Poly_Con_Relation::is_included();
955 }
956 }
957
958 dimension_type c_num_vars = 0;
959 dimension_type c_only_var = 0;
960
961 if (Box_Helpers::extract_interval_constraint(c, c_num_vars, c_only_var)) {
962 if (c_num_vars == 0) {
963 // c is a trivial constraint.
964 switch (sgn(c.inhomogeneous_term())) {
965 case -1:
966 return Poly_Con_Relation::is_disjoint();
967 case 0:
968 if (c.is_strict_inequality()) {
969 return Poly_Con_Relation::saturates()
970 && Poly_Con_Relation::is_disjoint();
971 }
972 else {
973 return Poly_Con_Relation::saturates()
974 && Poly_Con_Relation::is_included();
975 }
976 case 1:
977 return Poly_Con_Relation::is_included();
978 }
979 }
980 else {
981 // c is an interval constraint.
982 return interval_relation(seq[c_only_var],
983 c.type(),
984 c.inhomogeneous_term(),
985 c.coefficient(Variable(c_only_var)));
986 }
987 }
988 else {
989 // Deal with a non-trivial and non-interval constraint.
990 PPL_DIRTY_TEMP(Rational_Interval, r);
991 PPL_DIRTY_TEMP(Rational_Interval, t);
992 PPL_DIRTY_TEMP(mpq_class, m);
993 r = 0;
994 const Constraint::expr_type& e = c.expression();
995 for (Constraint::expr_type::const_iterator i = e.begin(), i_end = e.end();
996 i != i_end; ++i) {
997 assign_r(m, *i, ROUND_NOT_NEEDED);
998 const Variable v = i.variable();
999 // FIXME: an add_mul_assign() method would come handy here.
1000 t.build(seq[v.id()].lower_constraint(), seq[v.id()].upper_constraint());
1001 t *= m;
1002 r += t;
1003 }
1004 return interval_relation(r,
1005 c.type(),
1006 c.inhomogeneous_term());
1007 }
1008
1009 // Quiet a compiler warning: this program point is unreachable.
1010 PPL_UNREACHABLE;
1011 return Poly_Con_Relation::nothing();
1012 }
1013
1014 template <typename ITV>
1015 Poly_Gen_Relation
relation_with(const Generator & g) const1016 Box<ITV>::relation_with(const Generator& g) const {
1017 const dimension_type space_dim = space_dimension();
1018 const dimension_type g_space_dim = g.space_dimension();
1019
1020 // Dimension-compatibility check.
1021 if (space_dim < g_space_dim) {
1022 throw_dimension_incompatible("relation_with(g)", g);
1023 }
1024
1025 // The empty box cannot subsume a generator.
1026 if (is_empty()) {
1027 return Poly_Gen_Relation::nothing();
1028 }
1029
1030 // A universe box in a zero-dimensional space subsumes
1031 // all the generators of a zero-dimensional space.
1032 if (space_dim == 0) {
1033 return Poly_Gen_Relation::subsumes();
1034 }
1035
1036 if (g.is_line_or_ray()) {
1037 if (g.is_line()) {
1038 const Generator::expr_type& e = g.expression();
1039 for (Generator::expr_type::const_iterator i = e.begin(), i_end = e.end();
1040 i != i_end; ++i) {
1041 if (!seq[i.variable().id()].is_universe()) {
1042 return Poly_Gen_Relation::nothing();
1043 }
1044 }
1045 return Poly_Gen_Relation::subsumes();
1046 }
1047 else {
1048 PPL_ASSERT(g.is_ray());
1049 const Generator::expr_type& e = g.expression();
1050 for (Generator::expr_type::const_iterator i = e.begin(), i_end = e.end();
1051 i != i_end; ++i) {
1052 const Variable v = i.variable();
1053 switch (sgn(*i)) {
1054 case 1:
1055 if (!seq[v.id()].upper_is_boundary_infinity()) {
1056 return Poly_Gen_Relation::nothing();
1057 }
1058 break;
1059 case 0:
1060 PPL_UNREACHABLE;
1061 break;
1062 case -1:
1063 if (!seq[v.id()].lower_is_boundary_infinity()) {
1064 return Poly_Gen_Relation::nothing();
1065 }
1066 break;
1067 }
1068 }
1069 return Poly_Gen_Relation::subsumes();
1070 }
1071 }
1072
1073 // Here `g' is a point or closure point.
1074 const Coefficient& g_divisor = g.divisor();
1075 PPL_DIRTY_TEMP(mpq_class, g_coord);
1076 PPL_DIRTY_TEMP(mpq_class, bound);
1077 // TODO: If the variables in the expression that have coefficient 0
1078 // have no effect on seq[i], this loop can be optimized using
1079 // Generator::expr_type::const_iterator.
1080 for (dimension_type i = g_space_dim; i-- > 0; ) {
1081 const ITV& seq_i = seq[i];
1082 if (seq_i.is_universe()) {
1083 continue;
1084 }
1085 assign_r(g_coord.get_num(), g.coefficient(Variable(i)), ROUND_NOT_NEEDED);
1086 assign_r(g_coord.get_den(), g_divisor, ROUND_NOT_NEEDED);
1087 g_coord.canonicalize();
1088 // Check lower bound.
1089 if (!seq_i.lower_is_boundary_infinity()) {
1090 assign_r(bound, seq_i.lower(), ROUND_NOT_NEEDED);
1091 if (g_coord <= bound) {
1092 if (seq_i.lower_is_open()) {
1093 if (g.is_point() || g_coord != bound) {
1094 return Poly_Gen_Relation::nothing();
1095 }
1096 }
1097 else if (g_coord != bound) {
1098 return Poly_Gen_Relation::nothing();
1099 }
1100 }
1101 }
1102 // Check upper bound.
1103 if (!seq_i.upper_is_boundary_infinity()) {
1104 assign_r(bound, seq_i.upper(), ROUND_NOT_NEEDED);
1105 if (g_coord >= bound) {
1106 if (seq_i.upper_is_open()) {
1107 if (g.is_point() || g_coord != bound) {
1108 return Poly_Gen_Relation::nothing();
1109 }
1110 }
1111 else if (g_coord != bound) {
1112 return Poly_Gen_Relation::nothing();
1113 }
1114 }
1115 }
1116 }
1117 return Poly_Gen_Relation::subsumes();
1118 }
1119
1120
1121 template <typename ITV>
1122 bool
max_min(const Linear_Expression & expr,const bool maximize,Coefficient & ext_n,Coefficient & ext_d,bool & included) const1123 Box<ITV>::max_min(const Linear_Expression& expr,
1124 const bool maximize,
1125 Coefficient& ext_n, Coefficient& ext_d,
1126 bool& included) const {
1127 // `expr' should be dimension-compatible with `*this'.
1128 const dimension_type space_dim = space_dimension();
1129 const dimension_type expr_space_dim = expr.space_dimension();
1130 if (space_dim < expr_space_dim) {
1131 throw_dimension_incompatible((maximize
1132 ? "maximize(e, ...)"
1133 : "minimize(e, ...)"), "e", expr);
1134 }
1135
1136 // Deal with zero-dim Box first.
1137 if (space_dim == 0) {
1138 if (marked_empty()) {
1139 return false;
1140 }
1141 else {
1142 ext_n = expr.inhomogeneous_term();
1143 ext_d = 1;
1144 included = true;
1145 return true;
1146 }
1147 }
1148
1149 // For an empty Box we simply return false.
1150 if (is_empty()) {
1151 return false;
1152 }
1153
1154 PPL_DIRTY_TEMP(mpq_class, result);
1155 assign_r(result, expr.inhomogeneous_term(), ROUND_NOT_NEEDED);
1156 bool is_included = true;
1157 const int maximize_sign = maximize ? 1 : -1;
1158 PPL_DIRTY_TEMP(mpq_class, bound_i);
1159 PPL_DIRTY_TEMP(mpq_class, expr_i);
1160 for (Linear_Expression::const_iterator i = expr.begin(),
1161 i_end = expr.end(); i != i_end; ++i) {
1162 const ITV& seq_i = seq[i.variable().id()];
1163 assign_r(expr_i, *i, ROUND_NOT_NEEDED);
1164 switch (sgn(expr_i) * maximize_sign) {
1165 case 1:
1166 if (seq_i.upper_is_boundary_infinity()) {
1167 return false;
1168 }
1169 assign_r(bound_i, seq_i.upper(), ROUND_NOT_NEEDED);
1170 add_mul_assign_r(result, bound_i, expr_i, ROUND_NOT_NEEDED);
1171 if (seq_i.upper_is_open()) {
1172 is_included = false;
1173 }
1174 break;
1175 case 0:
1176 PPL_UNREACHABLE;
1177 break;
1178 case -1:
1179 if (seq_i.lower_is_boundary_infinity()) {
1180 return false;
1181 }
1182 assign_r(bound_i, seq_i.lower(), ROUND_NOT_NEEDED);
1183 add_mul_assign_r(result, bound_i, expr_i, ROUND_NOT_NEEDED);
1184 if (seq_i.lower_is_open()) {
1185 is_included = false;
1186 }
1187 break;
1188 }
1189 }
1190 // Extract output info.
1191 PPL_ASSERT(is_canonical(result));
1192 ext_n = result.get_num();
1193 ext_d = result.get_den();
1194 included = is_included;
1195 return true;
1196 }
1197
1198 template <typename ITV>
1199 bool
max_min(const Linear_Expression & expr,const bool maximize,Coefficient & ext_n,Coefficient & ext_d,bool & included,Generator & g) const1200 Box<ITV>::max_min(const Linear_Expression& expr,
1201 const bool maximize,
1202 Coefficient& ext_n, Coefficient& ext_d,
1203 bool& included,
1204 Generator& g) const {
1205 if (!max_min(expr, maximize, ext_n, ext_d, included)) {
1206 return false;
1207 }
1208 // Compute generator `g'.
1209 Linear_Expression g_expr;
1210 PPL_DIRTY_TEMP_COEFFICIENT(g_divisor);
1211 g_divisor = 1;
1212 const int maximize_sign = maximize ? 1 : -1;
1213 PPL_DIRTY_TEMP(mpq_class, g_coord);
1214 PPL_DIRTY_TEMP_COEFFICIENT(numer);
1215 PPL_DIRTY_TEMP_COEFFICIENT(denom);
1216 PPL_DIRTY_TEMP_COEFFICIENT(lcm);
1217 PPL_DIRTY_TEMP_COEFFICIENT(factor);
1218 // TODO: Check if the following loop can be optimized to exploit the
1219 // (possible) sparseness of expr.
1220 for (dimension_type i = space_dimension(); i-- > 0; ) {
1221 const ITV& seq_i = seq[i];
1222 switch (sgn(expr.coefficient(Variable(i))) * maximize_sign) {
1223 case 1:
1224 assign_r(g_coord, seq_i.upper(), ROUND_NOT_NEEDED);
1225 break;
1226 case 0:
1227 // If 0 belongs to the interval, choose it
1228 // (and directly proceed to the next iteration).
1229 // FIXME: name qualification issue.
1230 if (seq_i.contains(0)) {
1231 continue;
1232 }
1233 if (!seq_i.lower_is_boundary_infinity()) {
1234 if (seq_i.lower_is_open()) {
1235 if (!seq_i.upper_is_boundary_infinity()) {
1236 if (seq_i.upper_is_open()) {
1237 // Bounded and open interval: compute middle point.
1238 assign_r(g_coord, seq_i.lower(), ROUND_NOT_NEEDED);
1239 PPL_DIRTY_TEMP(mpq_class, q_seq_i_upper);
1240 assign_r(q_seq_i_upper, seq_i.upper(), ROUND_NOT_NEEDED);
1241 g_coord += q_seq_i_upper;
1242 g_coord /= 2;
1243 }
1244 else {
1245 // The upper bound is in the interval.
1246 assign_r(g_coord, seq_i.upper(), ROUND_NOT_NEEDED);
1247 }
1248 }
1249 else {
1250 // Lower is open, upper is unbounded.
1251 assign_r(g_coord, seq_i.lower(), ROUND_NOT_NEEDED);
1252 ++g_coord;
1253 }
1254 }
1255 else {
1256 // The lower bound is in the interval.
1257 assign_r(g_coord, seq_i.lower(), ROUND_NOT_NEEDED);
1258 }
1259 }
1260 else {
1261 // Lower is unbounded, hence upper is bounded
1262 // (since we know that 0 does not belong to the interval).
1263 PPL_ASSERT(!seq_i.upper_is_boundary_infinity());
1264 assign_r(g_coord, seq_i.upper(), ROUND_NOT_NEEDED);
1265 if (seq_i.upper_is_open()) {
1266 --g_coord;
1267 }
1268 }
1269 break;
1270 case -1:
1271 assign_r(g_coord, seq_i.lower(), ROUND_NOT_NEEDED);
1272 break;
1273 }
1274 // Add g_coord * Variable(i) to the generator.
1275 assign_r(denom, g_coord.get_den(), ROUND_NOT_NEEDED);
1276 lcm_assign(lcm, g_divisor, denom);
1277 exact_div_assign(factor, lcm, g_divisor);
1278 g_expr *= factor;
1279 exact_div_assign(factor, lcm, denom);
1280 assign_r(numer, g_coord.get_num(), ROUND_NOT_NEEDED);
1281 numer *= factor;
1282 g_expr += numer * Variable(i);
1283 g_divisor = lcm;
1284 }
1285 g = Generator::point(g_expr, g_divisor);
1286 return true;
1287 }
1288
1289 template <typename ITV>
1290 bool
contains(const Box & y) const1291 Box<ITV>::contains(const Box& y) const {
1292 const Box& x = *this;
1293 // Dimension-compatibility check.
1294 if (x.space_dimension() != y.space_dimension()) {
1295 x.throw_dimension_incompatible("contains(y)", y);
1296 }
1297
1298 // If `y' is empty, then `x' contains `y'.
1299 if (y.is_empty()) {
1300 return true;
1301 }
1302
1303 // If `x' is empty, then `x' cannot contain `y'.
1304 if (x.is_empty()) {
1305 return false;
1306 }
1307
1308 for (dimension_type k = x.seq.size(); k-- > 0; ) {
1309 // FIXME: fix this name qualification issue.
1310 if (!x.seq[k].contains(y.seq[k])) {
1311 return false;
1312 }
1313 }
1314 return true;
1315 }
1316
1317 template <typename ITV>
1318 bool
is_disjoint_from(const Box & y) const1319 Box<ITV>::is_disjoint_from(const Box& y) const {
1320 const Box& x = *this;
1321 // Dimension-compatibility check.
1322 if (x.space_dimension() != y.space_dimension()) {
1323 x.throw_dimension_incompatible("is_disjoint_from(y)", y);
1324 }
1325
1326 // If any of `x' or `y' is marked empty, then they are disjoint.
1327 // Note: no need to use `is_empty', as the following loop is anyway correct.
1328 if (x.marked_empty() || y.marked_empty()) {
1329 return true;
1330 }
1331
1332 for (dimension_type k = x.seq.size(); k-- > 0; ) {
1333 // FIXME: fix this name qualification issue.
1334 if (x.seq[k].is_disjoint_from(y.seq[k])) {
1335 return true;
1336 }
1337 }
1338 return false;
1339 }
1340
1341 template <typename ITV>
1342 inline bool
upper_bound_assign_if_exact(const Box & y)1343 Box<ITV>::upper_bound_assign_if_exact(const Box& y) {
1344 Box& x = *this;
1345
1346 // Dimension-compatibility check.
1347 if (x.space_dimension() != y.space_dimension()) {
1348 x.throw_dimension_incompatible("upper_bound_assign_if_exact(y)", y);
1349 }
1350
1351 // The lub of a box with an empty box is equal to the first box.
1352 if (y.is_empty()) {
1353 return true;
1354 }
1355 if (x.is_empty()) {
1356 x = y;
1357 return true;
1358 }
1359
1360 bool x_j_does_not_contain_y_j = false;
1361 bool y_j_does_not_contain_x_j = false;
1362
1363 for (dimension_type i = x.seq.size(); i-- > 0; ) {
1364 const ITV& x_seq_i = x.seq[i];
1365 const ITV& y_seq_i = y.seq[i];
1366
1367 if (!x_seq_i.can_be_exactly_joined_to(y_seq_i)) {
1368 return false;
1369 }
1370
1371 // Note: the use of `y_i_does_not_contain_x_i' is needed
1372 // because we want to temporarily preserve the old value
1373 // of `y_j_does_not_contain_x_j'.
1374 bool y_i_does_not_contain_x_i = !y_seq_i.contains(x_seq_i);
1375 if (y_i_does_not_contain_x_i && x_j_does_not_contain_y_j) {
1376 return false;
1377 }
1378 if (!x_seq_i.contains(y_seq_i)) {
1379 if (y_j_does_not_contain_x_j) {
1380 return false;
1381 }
1382 else {
1383 x_j_does_not_contain_y_j = true;
1384 }
1385 }
1386 if (y_i_does_not_contain_x_i) {
1387 y_j_does_not_contain_x_j = true;
1388 }
1389 }
1390
1391 // The upper bound is exact: compute it into *this.
1392 for (dimension_type k = x.seq.size(); k-- > 0; ) {
1393 x.seq[k].join_assign(y.seq[k]);
1394 }
1395 return true;
1396 }
1397
1398 template <typename ITV>
1399 bool
OK() const1400 Box<ITV>::OK() const {
1401 if (status.test_empty_up_to_date() && !status.test_empty()) {
1402 Box tmp = *this;
1403 tmp.reset_empty_up_to_date();
1404 if (tmp.check_empty()) {
1405 #ifndef NDEBUG
1406 std::cerr << "The box is empty, but it is marked as non-empty."
1407 << std::endl;
1408 #endif // NDEBUG
1409 return false;
1410 }
1411 }
1412
1413 // A box that is not marked empty must have meaningful intervals.
1414 if (!marked_empty()) {
1415 for (dimension_type k = seq.size(); k-- > 0; ) {
1416 if (!seq[k].OK()) {
1417 return false;
1418 }
1419 }
1420 }
1421
1422 return true;
1423 }
1424
1425 template <typename ITV>
1426 dimension_type
affine_dimension() const1427 Box<ITV>::affine_dimension() const {
1428 dimension_type d = space_dimension();
1429 // A zero-space-dim box always has affine dimension zero.
1430 if (d == 0) {
1431 return 0;
1432 }
1433
1434 // An empty box has affine dimension zero.
1435 if (is_empty()) {
1436 return 0;
1437 }
1438
1439 for (dimension_type k = d; k-- > 0; ) {
1440 if (seq[k].is_singleton()) {
1441 --d;
1442 }
1443 }
1444
1445 return d;
1446 }
1447
1448 template <typename ITV>
1449 bool
check_empty() const1450 Box<ITV>::check_empty() const {
1451 PPL_ASSERT(!marked_empty());
1452 Box<ITV>& x = const_cast<Box<ITV>&>(*this);
1453 for (dimension_type k = seq.size(); k-- > 0; ) {
1454 if (seq[k].is_empty()) {
1455 x.set_empty();
1456 return true;
1457 }
1458 }
1459 x.set_nonempty();
1460 return false;
1461 }
1462
1463 template <typename ITV>
1464 bool
is_universe() const1465 Box<ITV>::is_universe() const {
1466 if (marked_empty()) {
1467 return false;
1468 }
1469 for (dimension_type k = seq.size(); k-- > 0; ) {
1470 if (!seq[k].is_universe()) {
1471 return false;
1472 }
1473 }
1474 return true;
1475 }
1476
1477 template <typename ITV>
1478 bool
is_topologically_closed() const1479 Box<ITV>::is_topologically_closed() const {
1480 if (ITV::is_always_topologically_closed() || is_empty()) {
1481 return true;
1482 }
1483
1484 for (dimension_type k = seq.size(); k-- > 0; ) {
1485 if (!seq[k].is_topologically_closed()) {
1486 return false;
1487 }
1488 }
1489 return true;
1490 }
1491
1492 template <typename ITV>
1493 bool
is_discrete() const1494 Box<ITV>::is_discrete() const {
1495 if (is_empty()) {
1496 return true;
1497 }
1498 for (dimension_type k = seq.size(); k-- > 0; ) {
1499 if (!seq[k].is_singleton()) {
1500 return false;
1501 }
1502 }
1503 return true;
1504 }
1505
1506 template <typename ITV>
1507 bool
is_bounded() const1508 Box<ITV>::is_bounded() const {
1509 if (is_empty()) {
1510 return true;
1511 }
1512 for (dimension_type k = seq.size(); k-- > 0; ) {
1513 if (!seq[k].is_bounded()) {
1514 return false;
1515 }
1516 }
1517 return true;
1518 }
1519
1520 template <typename ITV>
1521 bool
contains_integer_point() const1522 Box<ITV>::contains_integer_point() const {
1523 if (marked_empty()) {
1524 return false;
1525 }
1526 for (dimension_type k = seq.size(); k-- > 0; ) {
1527 if (!seq[k].contains_integer_point()) {
1528 return false;
1529 }
1530 }
1531 return true;
1532 }
1533
1534 template <typename ITV>
1535 bool
frequency(const Linear_Expression & expr,Coefficient & freq_n,Coefficient & freq_d,Coefficient & val_n,Coefficient & val_d) const1536 Box<ITV>::frequency(const Linear_Expression& expr,
1537 Coefficient& freq_n, Coefficient& freq_d,
1538 Coefficient& val_n, Coefficient& val_d) const {
1539 dimension_type space_dim = space_dimension();
1540 // The dimension of `expr' must be at most the dimension of *this.
1541 if (space_dim < expr.space_dimension()) {
1542 throw_dimension_incompatible("frequency(e, ...)", "e", expr);
1543 }
1544
1545 // Check if `expr' has a constant value.
1546 // If it is constant, set the frequency `freq_n' to 0
1547 // and return true. Otherwise the values for \p expr
1548 // are not discrete so return false.
1549
1550 // Space dimension is 0: if empty, then return false;
1551 // otherwise the frequency is 0 and the value is the inhomogeneous term.
1552 if (space_dim == 0) {
1553 if (is_empty()) {
1554 return false;
1555 }
1556 freq_n = 0;
1557 freq_d = 1;
1558 val_n = expr.inhomogeneous_term();
1559 val_d = 1;
1560 return true;
1561 }
1562
1563 // For an empty Box, we simply return false.
1564 if (is_empty()) {
1565 return false;
1566 }
1567
1568 // The Box has at least 1 dimension and is not empty.
1569 PPL_DIRTY_TEMP_COEFFICIENT(numer);
1570 PPL_DIRTY_TEMP_COEFFICIENT(denom);
1571 PPL_DIRTY_TEMP(mpq_class, tmp);
1572 Coefficient c = expr.inhomogeneous_term();
1573
1574 PPL_DIRTY_TEMP_COEFFICIENT(val_denom);
1575 val_denom = 1;
1576
1577 for (Linear_Expression::const_iterator i = expr.begin(), i_end = expr.end();
1578 i != i_end; ++i) {
1579 const ITV& seq_i = seq[i.variable().id()];
1580 // Check if `v' is constant in the BD shape.
1581 if (seq_i.is_singleton()) {
1582 // If `v' is constant, replace it in `le' by the value.
1583 assign_r(tmp, seq_i.lower(), ROUND_NOT_NEEDED);
1584 numer = tmp.get_num();
1585 denom = tmp.get_den();
1586 c *= denom;
1587 c += numer * val_denom * (*i);
1588 val_denom *= denom;
1589 continue;
1590 }
1591 // The expression `expr' is not constant.
1592 return false;
1593 }
1594
1595 // The expression `expr' is constant.
1596 freq_n = 0;
1597 freq_d = 1;
1598
1599 // Reduce `val_n' and `val_d'.
1600 normalize2(c, val_denom, val_n, val_d);
1601 return true;
1602 }
1603
1604 template <typename ITV>
1605 bool
constrains(Variable var) const1606 Box<ITV>::constrains(Variable var) const {
1607 // `var' should be one of the dimensions of the polyhedron.
1608 const dimension_type var_space_dim = var.space_dimension();
1609 if (space_dimension() < var_space_dim) {
1610 throw_dimension_incompatible("constrains(v)", "v", var);
1611 }
1612
1613 if (marked_empty() || !seq[var_space_dim-1].is_universe()) {
1614 return true;
1615 }
1616 // Now force an emptiness check.
1617 return is_empty();
1618 }
1619
1620 template <typename ITV>
1621 void
unconstrain(const Variables_Set & vars)1622 Box<ITV>::unconstrain(const Variables_Set& vars) {
1623 // The cylindrification with respect to no dimensions is a no-op.
1624 // This case also captures the only legal cylindrification
1625 // of a box in a 0-dim space.
1626 if (vars.empty()) {
1627 return;
1628 }
1629
1630 // Dimension-compatibility check.
1631 const dimension_type min_space_dim = vars.space_dimension();
1632 if (space_dimension() < min_space_dim) {
1633 throw_dimension_incompatible("unconstrain(vs)", min_space_dim);
1634 }
1635 // If the box is already empty, there is nothing left to do.
1636 if (marked_empty()) {
1637 return;
1638 }
1639
1640 // Here the box might still be empty (but we haven't detected it yet):
1641 // check emptiness of the interval for each of the variables in
1642 // `vars' before cylindrification.
1643 for (Variables_Set::const_iterator vsi = vars.begin(),
1644 vsi_end = vars.end(); vsi != vsi_end; ++vsi) {
1645 ITV& seq_vsi = seq[*vsi];
1646 if (!seq_vsi.is_empty()) {
1647 seq_vsi.assign(UNIVERSE);
1648 }
1649 else {
1650 set_empty();
1651 break;
1652 }
1653 }
1654 PPL_ASSERT(OK());
1655 }
1656
1657 template <typename ITV>
1658 void
topological_closure_assign()1659 Box<ITV>::topological_closure_assign() {
1660 if (ITV::is_always_topologically_closed() || is_empty()) {
1661 return;
1662 }
1663
1664 for (dimension_type k = seq.size(); k-- > 0; ) {
1665 seq[k].topological_closure_assign();
1666 }
1667 }
1668
1669 template <typename ITV>
1670 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)1671 Box<ITV>::wrap_assign(const Variables_Set& vars,
1672 Bounded_Integer_Type_Width w,
1673 Bounded_Integer_Type_Representation r,
1674 Bounded_Integer_Type_Overflow o,
1675 const Constraint_System* cs_p,
1676 unsigned complexity_threshold,
1677 bool wrap_individually) {
1678 #if 0 // Generic implementation commented out.
1679 Implementation::wrap_assign(*this,
1680 vars, w, r, o, cs_p,
1681 complexity_threshold, wrap_individually,
1682 "Box");
1683 #else // Specialized implementation.
1684 PPL_USED(wrap_individually);
1685 PPL_USED(complexity_threshold);
1686 Box& x = *this;
1687
1688 // Dimension-compatibility check for `*cs_p', if any.
1689 const dimension_type vars_space_dim = vars.space_dimension();
1690 if (cs_p != 0 && cs_p->space_dimension() > vars_space_dim) {
1691 std::ostringstream s;
1692 s << "PPL::Box<ITV>::wrap_assign(vars, w, r, o, cs_p, ...):"
1693 << std::endl
1694 << "vars.space_dimension() == " << vars_space_dim
1695 << ", cs_p->space_dimension() == " << cs_p->space_dimension() << ".";
1696 throw std::invalid_argument(s.str());
1697 }
1698
1699 // Wrapping no variable only requires refining with *cs_p, if any.
1700 if (vars.empty()) {
1701 if (cs_p != 0) {
1702 refine_with_constraints(*cs_p);
1703 }
1704 return;
1705 }
1706
1707 // Dimension-compatibility check for `vars'.
1708 const dimension_type space_dim = x.space_dimension();
1709 if (space_dim < vars_space_dim) {
1710 std::ostringstream s;
1711 s << "PPL::Box<ITV>::wrap_assign(vars, ...):"
1712 << std::endl
1713 << "this->space_dimension() == " << space_dim
1714 << ", required space dimension == " << vars_space_dim << ".";
1715 throw std::invalid_argument(s.str());
1716 }
1717
1718 // Wrapping an empty polyhedron is a no-op.
1719 if (x.is_empty()) {
1720 return;
1721 }
1722
1723 // FIXME: temporarily (ab-) using Coefficient.
1724 // Set `min_value' and `max_value' to the minimum and maximum values
1725 // a variable of width `w' and signedness `s' can take.
1726 PPL_DIRTY_TEMP_COEFFICIENT(min_value);
1727 PPL_DIRTY_TEMP_COEFFICIENT(max_value);
1728 if (r == UNSIGNED) {
1729 min_value = 0;
1730 mul_2exp_assign(max_value, Coefficient_one(), w);
1731 --max_value;
1732 }
1733 else {
1734 PPL_ASSERT(r == SIGNED_2_COMPLEMENT);
1735 mul_2exp_assign(max_value, Coefficient_one(), w-1);
1736 neg_assign(min_value, max_value);
1737 --max_value;
1738 }
1739
1740 // FIXME: Build the (integer) quadrant interval.
1741 PPL_DIRTY_TEMP(ITV, integer_quadrant_itv);
1742 PPL_DIRTY_TEMP(ITV, rational_quadrant_itv);
1743 {
1744 I_Constraint<Coefficient> lower = i_constraint(GREATER_OR_EQUAL, min_value);
1745 I_Constraint<Coefficient> upper = i_constraint(LESS_OR_EQUAL, max_value);
1746 integer_quadrant_itv.build(lower, upper);
1747 // The rational quadrant is only needed if overflow is undefined.
1748 if (o == OVERFLOW_UNDEFINED) {
1749 ++max_value;
1750 upper = i_constraint(LESS_THAN, max_value);
1751 rational_quadrant_itv.build(lower, upper);
1752 }
1753 }
1754
1755 const Variables_Set::const_iterator vs_end = vars.end();
1756
1757 if (cs_p == 0) {
1758 // No constraint refinement is needed here.
1759 switch (o) {
1760 case OVERFLOW_WRAPS:
1761 for (Variables_Set::const_iterator i = vars.begin(); i != vs_end; ++i) {
1762 x.seq[*i].wrap_assign(w, r, integer_quadrant_itv);
1763 }
1764 reset_empty_up_to_date();
1765 break;
1766 case OVERFLOW_UNDEFINED:
1767 for (Variables_Set::const_iterator i = vars.begin(); i != vs_end; ++i) {
1768 ITV& x_seq_v = x.seq[*i];
1769 if (!rational_quadrant_itv.contains(x_seq_v)) {
1770 x_seq_v.assign(integer_quadrant_itv);
1771 }
1772 }
1773 break;
1774 case OVERFLOW_IMPOSSIBLE:
1775 for (Variables_Set::const_iterator i = vars.begin(); i != vs_end; ++i) {
1776 x.seq[*i].intersect_assign(integer_quadrant_itv);
1777 }
1778 reset_empty_up_to_date();
1779 break;
1780 }
1781 PPL_ASSERT(x.OK());
1782 return;
1783 }
1784
1785 PPL_ASSERT(cs_p != 0);
1786 const Constraint_System& cs = *cs_p;
1787 // A map associating interval constraints to variable indexes.
1788 typedef std::map<dimension_type, std::vector<const Constraint*> > map_type;
1789 map_type var_cs_map;
1790 for (Constraint_System::const_iterator i = cs.begin(),
1791 i_end = cs.end(); i != i_end; ++i) {
1792 const Constraint& c = *i;
1793 dimension_type c_num_vars = 0;
1794 dimension_type c_only_var = 0;
1795 if (Box_Helpers::extract_interval_constraint(c, c_num_vars, c_only_var)) {
1796 if (c_num_vars == 1) {
1797 // An interval constraint on variable index `c_only_var'.
1798 PPL_ASSERT(c_only_var < space_dim);
1799 // We do care about c if c_only_var is going to be wrapped.
1800 if (vars.find(c_only_var) != vs_end) {
1801 var_cs_map[c_only_var].push_back(&c);
1802 }
1803 }
1804 else {
1805 PPL_ASSERT(c_num_vars == 0);
1806 // Note: tautologies have been filtered out by iterators.
1807 PPL_ASSERT(c.is_inconsistent());
1808 x.set_empty();
1809 return;
1810 }
1811 }
1812 }
1813
1814 PPL_DIRTY_TEMP(ITV, refinement_itv);
1815 const map_type::const_iterator var_cs_map_end = var_cs_map.end();
1816 // Loop through the variable indexes in `vars'.
1817 for (Variables_Set::const_iterator i = vars.begin(); i != vs_end; ++i) {
1818 const dimension_type v = *i;
1819 refinement_itv = integer_quadrant_itv;
1820 // Look for the refinement constraints for space dimension index `v'.
1821 map_type::const_iterator var_cs_map_iter = var_cs_map.find(v);
1822 if (var_cs_map_iter != var_cs_map_end) {
1823 // Refine interval for variable `v'.
1824 const map_type::mapped_type& var_cs = var_cs_map_iter->second;
1825 for (dimension_type j = var_cs.size(); j-- > 0; ) {
1826 const Constraint& c = *var_cs[j];
1827 refine_interval_no_check(refinement_itv,
1828 c.type(),
1829 c.inhomogeneous_term(),
1830 c.coefficient(Variable(v)));
1831 }
1832 }
1833 // Wrap space dimension index `v'.
1834 ITV& x_seq_v = x.seq[v];
1835 switch (o) {
1836 case OVERFLOW_WRAPS:
1837 x_seq_v.wrap_assign(w, r, refinement_itv);
1838 break;
1839 case OVERFLOW_UNDEFINED:
1840 if (!rational_quadrant_itv.contains(x_seq_v)) {
1841 x_seq_v.assign(UNIVERSE);
1842 }
1843 break;
1844 case OVERFLOW_IMPOSSIBLE:
1845 x_seq_v.intersect_assign(refinement_itv);
1846 break;
1847 }
1848 }
1849 PPL_ASSERT(x.OK());
1850 #endif
1851 }
1852
1853 template <typename ITV>
1854 void
drop_some_non_integer_points(Complexity_Class)1855 Box<ITV>::drop_some_non_integer_points(Complexity_Class) {
1856 if (std::numeric_limits<typename ITV::boundary_type>::is_integer
1857 && !ITV::info_type::store_open) {
1858 return;
1859 }
1860
1861 if (marked_empty()) {
1862 return;
1863 }
1864
1865 for (dimension_type k = seq.size(); k-- > 0; ) {
1866 seq[k].drop_some_non_integer_points();
1867 }
1868
1869 PPL_ASSERT(OK());
1870 }
1871
1872 template <typename ITV>
1873 void
drop_some_non_integer_points(const Variables_Set & vars,Complexity_Class)1874 Box<ITV>::drop_some_non_integer_points(const Variables_Set& vars,
1875 Complexity_Class) {
1876 // Dimension-compatibility check.
1877 const dimension_type min_space_dim = vars.space_dimension();
1878 if (space_dimension() < min_space_dim) {
1879 throw_dimension_incompatible("drop_some_non_integer_points(vs, cmpl)",
1880 min_space_dim);
1881 }
1882 if (std::numeric_limits<typename ITV::boundary_type>::is_integer
1883 && !ITV::info_type::store_open) {
1884 return;
1885 }
1886
1887 if (marked_empty()) {
1888 return;
1889 }
1890
1891 for (Variables_Set::const_iterator v_i = vars.begin(),
1892 v_end = vars.end(); v_i != v_end; ++v_i) {
1893 seq[*v_i].drop_some_non_integer_points();
1894 }
1895
1896 PPL_ASSERT(OK());
1897 }
1898
1899 template <typename ITV>
1900 void
intersection_assign(const Box & y)1901 Box<ITV>::intersection_assign(const Box& y) {
1902 Box& x = *this;
1903 const dimension_type space_dim = space_dimension();
1904
1905 // Dimension-compatibility check.
1906 if (space_dim != y.space_dimension()) {
1907 x.throw_dimension_incompatible("intersection_assign(y)", y);
1908 }
1909
1910 // If one of the two boxes is empty, the intersection is empty.
1911 if (x.marked_empty()) {
1912 return;
1913 }
1914 if (y.marked_empty()) {
1915 x.set_empty();
1916 return;
1917 }
1918
1919 // If both boxes are zero-dimensional, then at this point they are
1920 // necessarily non-empty, so that their intersection is non-empty too.
1921 if (space_dim == 0) {
1922 return;
1923 }
1924
1925 // FIXME: here we may conditionally exploit a capability of the
1926 // underlying interval to eagerly detect empty results.
1927 reset_empty_up_to_date();
1928
1929 for (dimension_type k = space_dim; k-- > 0; ) {
1930 x.seq[k].intersect_assign(y.seq[k]);
1931 }
1932
1933 PPL_ASSERT(x.OK());
1934 }
1935
1936 template <typename ITV>
1937 void
upper_bound_assign(const Box & y)1938 Box<ITV>::upper_bound_assign(const Box& y) {
1939 Box& x = *this;
1940
1941 // Dimension-compatibility check.
1942 if (x.space_dimension() != y.space_dimension()) {
1943 x.throw_dimension_incompatible("upper_bound_assign(y)", y);
1944 }
1945
1946 // The lub of a box with an empty box is equal to the first box.
1947 if (y.is_empty()) {
1948 return;
1949 }
1950 if (x.is_empty()) {
1951 x = y;
1952 return;
1953 }
1954
1955 for (dimension_type k = x.seq.size(); k-- > 0; ) {
1956 x.seq[k].join_assign(y.seq[k]);
1957 }
1958
1959 PPL_ASSERT(x.OK());
1960 }
1961
1962 template <typename ITV>
1963 void
concatenate_assign(const Box & y)1964 Box<ITV>::concatenate_assign(const Box& y) {
1965 Box& x = *this;
1966 const dimension_type x_space_dim = x.space_dimension();
1967 const dimension_type y_space_dim = y.space_dimension();
1968
1969 // If `y' is marked empty, the result will be empty too.
1970 if (y.marked_empty()) {
1971 x.set_empty();
1972 }
1973
1974 // If `y' is a 0-dim space box, there is nothing left to do.
1975 if (y_space_dim == 0) {
1976 return;
1977 }
1978 // The resulting space dimension must be at most the maximum.
1979 check_space_dimension_overflow(y.space_dimension(),
1980 max_space_dimension() - space_dimension(),
1981 "PPL::Box::",
1982 "concatenate_assign(y)",
1983 "concatenation exceeds the maximum "
1984 "allowed space dimension");
1985 // Here `y_space_dim > 0', so that a non-trivial concatenation will occur:
1986 // make sure that reallocation will occur once at most.
1987 x.seq.reserve(x_space_dim + y_space_dim);
1988
1989 // If `x' is marked empty, then it is sufficient to adjust
1990 // the dimension of the vector space.
1991 if (x.marked_empty()) {
1992 x.seq.insert(x.seq.end(), y_space_dim, ITV(EMPTY));
1993 PPL_ASSERT(x.OK());
1994 return;
1995 }
1996
1997 // Here neither `x' nor `y' are marked empty: concatenate them.
1998 std::copy(y.seq.begin(), y.seq.end(),
1999 std::back_insert_iterator<Sequence>(x.seq));
2000 // Update the `empty_up_to_date' flag.
2001 if (!y.status.test_empty_up_to_date()) {
2002 reset_empty_up_to_date();
2003 }
2004
2005 PPL_ASSERT(x.OK());
2006 }
2007
2008 template <typename ITV>
2009 void
difference_assign(const Box & y)2010 Box<ITV>::difference_assign(const Box& y) {
2011 const dimension_type space_dim = space_dimension();
2012
2013 // Dimension-compatibility check.
2014 if (space_dim != y.space_dimension()) {
2015 throw_dimension_incompatible("difference_assign(y)", y);
2016 }
2017
2018 Box& x = *this;
2019 if (x.is_empty() || y.is_empty()) {
2020 return;
2021 }
2022 switch (space_dim) {
2023 case 0:
2024 // If `x' is zero-dimensional, then at this point both `x' and `y'
2025 // are the universe box, so that their difference is empty.
2026 x.set_empty();
2027 break;
2028
2029 case 1:
2030 x.seq[0].difference_assign(y.seq[0]);
2031 if (x.seq[0].is_empty()) {
2032 x.set_empty();
2033 }
2034 break;
2035
2036 default:
2037 {
2038 dimension_type index_non_contained = space_dim;
2039 dimension_type number_non_contained = 0;
2040 for (dimension_type i = space_dim; i-- > 0; ) {
2041 if (!y.seq[i].contains(x.seq[i])) {
2042 if (++number_non_contained == 1) {
2043 index_non_contained = i;
2044 }
2045 else {
2046 break;
2047 }
2048 }
2049 }
2050
2051 switch (number_non_contained) {
2052 case 0:
2053 // `y' covers `x': the difference is empty.
2054 x.set_empty();
2055 break;
2056 case 1:
2057 x.seq[index_non_contained]
2058 .difference_assign(y.seq[index_non_contained]);
2059 if (x.seq[index_non_contained].is_empty()) {
2060 x.set_empty();
2061 }
2062 break;
2063 default:
2064 // Nothing to do: the difference is `x'.
2065 break;
2066 }
2067 }
2068 break;
2069 }
2070 PPL_ASSERT(OK());
2071 }
2072
2073 template <typename ITV>
2074 bool
simplify_using_context_assign(const Box & y)2075 Box<ITV>::simplify_using_context_assign(const Box& y) {
2076 Box& x = *this;
2077 const dimension_type num_dims = x.space_dimension();
2078 // Dimension-compatibility check.
2079 if (num_dims != y.space_dimension()) {
2080 x.throw_dimension_incompatible("simplify_using_context_assign(y)", y);
2081 }
2082
2083 // Filter away the zero-dimensional case.
2084 if (num_dims == 0) {
2085 if (y.marked_empty()) {
2086 x.set_nonempty();
2087 return false;
2088 }
2089 else {
2090 return !x.marked_empty();
2091 }
2092 }
2093
2094 // Filter away the case when `y' is empty.
2095 if (y.is_empty()) {
2096 for (dimension_type i = num_dims; i-- > 0; ) {
2097 x.seq[i].assign(UNIVERSE);
2098 }
2099 x.set_nonempty();
2100 return false;
2101 }
2102
2103 if (x.is_empty()) {
2104 // Find in `y' a non-universe interval, if any.
2105 for (dimension_type i = 0; i < num_dims; ++i) {
2106 if (y.seq[i].is_universe()) {
2107 x.seq[i].assign(UNIVERSE);
2108 }
2109 else {
2110 // Set x.seq[i] so as to contradict y.seq[i], if possible.
2111 ITV& seq_i = x.seq[i];
2112 seq_i.empty_intersection_assign(y.seq[i]);
2113 if (seq_i.is_empty()) {
2114 // We were not able to assign to `seq_i' a non-empty interval:
2115 // reset `seq_i' to the universe interval and keep searching.
2116 seq_i.assign(UNIVERSE);
2117 continue;
2118 }
2119 // We assigned to `seq_i' a non-empty interval:
2120 // set the other intervals to universe and return.
2121 for (++i; i < num_dims; ++i) {
2122 x.seq[i].assign(UNIVERSE);
2123 }
2124 x.set_nonempty();
2125 PPL_ASSERT(x.OK());
2126 return false;
2127 }
2128 }
2129 // All intervals in `y' are universe or could not be contradicted:
2130 // simplification can leave the empty box `x' as is.
2131 PPL_ASSERT(x.OK() && x.is_empty());
2132 return false;
2133 }
2134
2135 // Loop index `i' is intentionally going upwards.
2136 for (dimension_type i = 0; i < num_dims; ++i) {
2137 if (!x.seq[i].simplify_using_context_assign(y.seq[i])) {
2138 PPL_ASSERT(!x.seq[i].is_empty());
2139 // The intersection of `x' and `y' is empty due to the i-th interval:
2140 // reset other intervals to UNIVERSE.
2141 for (dimension_type j = num_dims; j-- > i; ) {
2142 x.seq[j].assign(UNIVERSE);
2143 }
2144 for (dimension_type j = i; j-- > 0; ) {
2145 x.seq[j].assign(UNIVERSE);
2146 }
2147 PPL_ASSERT(x.OK());
2148 return false;
2149 }
2150 }
2151 PPL_ASSERT(x.OK());
2152 return true;
2153 }
2154
2155 template <typename ITV>
2156 void
time_elapse_assign(const Box & y)2157 Box<ITV>::time_elapse_assign(const Box& y) {
2158 Box& x = *this;
2159 const dimension_type x_space_dim = x.space_dimension();
2160
2161 // Dimension-compatibility check.
2162 if (x_space_dim != y.space_dimension()) {
2163 x.throw_dimension_incompatible("time_elapse_assign(y)", y);
2164 }
2165
2166 // Dealing with the zero-dimensional case.
2167 if (x_space_dim == 0) {
2168 if (y.marked_empty()) {
2169 x.set_empty();
2170 }
2171 return;
2172 }
2173
2174 // If either one of `x' or `y' is empty, the result is empty too.
2175 // Note: if possible, avoid cost of checking for emptiness.
2176 if (x.marked_empty() || y.marked_empty()
2177 || x.is_empty() || y.is_empty()) {
2178 x.set_empty();
2179 return;
2180 }
2181
2182 for (dimension_type i = x_space_dim; i-- > 0; ) {
2183 ITV& x_seq_i = x.seq[i];
2184 const ITV& y_seq_i = y.seq[i];
2185 if (!x_seq_i.lower_is_boundary_infinity()) {
2186 if (y_seq_i.lower_is_boundary_infinity() || y_seq_i.lower() < 0) {
2187 x_seq_i.lower_extend();
2188 }
2189 }
2190 if (!x_seq_i.upper_is_boundary_infinity()) {
2191 if (y_seq_i.upper_is_boundary_infinity() || y_seq_i.upper() > 0) {
2192 x_seq_i.upper_extend();
2193 }
2194 }
2195 }
2196 PPL_ASSERT(x.OK());
2197 }
2198
2199 template <typename ITV>
2200 inline void
remove_space_dimensions(const Variables_Set & vars)2201 Box<ITV>::remove_space_dimensions(const Variables_Set& vars) {
2202 // The removal of no dimensions from any box is a no-op.
2203 // Note that this case also captures the only legal removal of
2204 // space dimensions from a box in a zero-dimensional space.
2205 if (vars.empty()) {
2206 PPL_ASSERT(OK());
2207 return;
2208 }
2209
2210 const dimension_type old_space_dim = space_dimension();
2211
2212 // Dimension-compatibility check.
2213 const dimension_type vsi_space_dim = vars.space_dimension();
2214 if (old_space_dim < vsi_space_dim) {
2215 throw_dimension_incompatible("remove_space_dimensions(vs)",
2216 vsi_space_dim);
2217 }
2218
2219 const dimension_type new_space_dim = old_space_dim - vars.size();
2220
2221 // If the box is empty (this must be detected), then resizing is all
2222 // what is needed. If it is not empty and we are removing _all_ the
2223 // dimensions then, again, resizing suffices.
2224 if (is_empty() || new_space_dim == 0) {
2225 seq.resize(new_space_dim);
2226 PPL_ASSERT(OK());
2227 return;
2228 }
2229
2230 // For each variable to be removed, we fill the corresponding interval
2231 // by shifting left those intervals that will not be removed.
2232 Variables_Set::const_iterator vsi = vars.begin();
2233 Variables_Set::const_iterator vsi_end = vars.end();
2234 dimension_type dst = *vsi;
2235 dimension_type src = dst + 1;
2236 for (++vsi; vsi != vsi_end; ++vsi) {
2237 const dimension_type vsi_next = *vsi;
2238 // All intervals in between are moved to the left.
2239 while (src < vsi_next) {
2240 swap(seq[dst++], seq[src++]);
2241 }
2242 ++src;
2243 }
2244
2245 // Moving the remaining intervals.
2246 while (src < old_space_dim) {
2247 swap(seq[dst++], seq[src++]);
2248 }
2249
2250 PPL_ASSERT(dst == new_space_dim);
2251 seq.resize(new_space_dim);
2252
2253 PPL_ASSERT(OK());
2254 }
2255
2256 template <typename ITV>
2257 void
remove_higher_space_dimensions(const dimension_type new_dimension)2258 Box<ITV>::remove_higher_space_dimensions(const dimension_type new_dimension) {
2259 // Dimension-compatibility check: the variable having
2260 // maximum index is the one occurring last in the set.
2261 const dimension_type space_dim = space_dimension();
2262 if (new_dimension > space_dim) {
2263 throw_dimension_incompatible("remove_higher_space_dimensions(nd)",
2264 new_dimension);
2265 }
2266
2267 // The removal of no dimensions from any box is a no-op.
2268 // Note that this case also captures the only legal removal of
2269 // dimensions from a zero-dim space box.
2270 if (new_dimension == space_dim) {
2271 PPL_ASSERT(OK());
2272 return;
2273 }
2274
2275 seq.resize(new_dimension);
2276 PPL_ASSERT(OK());
2277 }
2278
2279 template <typename ITV>
2280 template <typename Partial_Function>
2281 void
map_space_dimensions(const Partial_Function & pfunc)2282 Box<ITV>::map_space_dimensions(const Partial_Function& pfunc) {
2283 const dimension_type space_dim = space_dimension();
2284 if (space_dim == 0) {
2285 return;
2286 }
2287
2288 if (pfunc.has_empty_codomain()) {
2289 // All dimensions vanish: the box becomes zero_dimensional.
2290 remove_higher_space_dimensions(0);
2291 return;
2292 }
2293
2294 const dimension_type new_space_dim = pfunc.max_in_codomain() + 1;
2295 // If the box is empty, then simply adjust the space dimension.
2296 if (is_empty()) {
2297 remove_higher_space_dimensions(new_space_dim);
2298 return;
2299 }
2300
2301 // We create a new Box with the new space dimension.
2302 Box<ITV> tmp(new_space_dim);
2303 // Map the intervals, exchanging the indexes.
2304 for (dimension_type i = 0; i < space_dim; ++i) {
2305 dimension_type new_i;
2306 if (pfunc.maps(i, new_i)) {
2307 swap(seq[i], tmp.seq[new_i]);
2308 }
2309 }
2310 m_swap(tmp);
2311 PPL_ASSERT(OK());
2312 }
2313
2314 template <typename ITV>
2315 void
fold_space_dimensions(const Variables_Set & vars,const Variable dest)2316 Box<ITV>::fold_space_dimensions(const Variables_Set& vars,
2317 const Variable dest) {
2318 const dimension_type space_dim = space_dimension();
2319 // `dest' should be one of the dimensions of the box.
2320 if (dest.space_dimension() > space_dim) {
2321 throw_dimension_incompatible("fold_space_dimensions(vs, v)", "v", dest);
2322 }
2323
2324 // The folding of no dimensions is a no-op.
2325 if (vars.empty()) {
2326 return;
2327 }
2328
2329 // All variables in `vars' should be dimensions of the box.
2330 if (vars.space_dimension() > space_dim) {
2331 throw_dimension_incompatible("fold_space_dimensions(vs, v)",
2332 vars.space_dimension());
2333 }
2334 // Moreover, `dest.id()' should not occur in `vars'.
2335 if (vars.find(dest.id()) != vars.end()) {
2336 throw_invalid_argument("fold_space_dimensions(vs, v)",
2337 "v should not occur in vs");
2338 }
2339
2340 // Note: the check for emptiness is needed for correctness.
2341 if (!is_empty()) {
2342 // Join the interval corresponding to variable `dest' with the intervals
2343 // corresponding to the variables in `vars'.
2344 ITV& seq_v = seq[dest.id()];
2345 for (Variables_Set::const_iterator i = vars.begin(),
2346 vs_end = vars.end(); i != vs_end; ++i) {
2347 seq_v.join_assign(seq[*i]);
2348 }
2349 }
2350 remove_space_dimensions(vars);
2351 }
2352
2353 template <typename ITV>
2354 void
add_constraint_no_check(const Constraint & c)2355 Box<ITV>::add_constraint_no_check(const Constraint& c) {
2356 PPL_ASSERT(c.space_dimension() <= space_dimension());
2357
2358 dimension_type c_num_vars = 0;
2359 dimension_type c_only_var = 0;
2360 // Throw an exception if c is not an interval constraints.
2361 if (!Box_Helpers::extract_interval_constraint(c, c_num_vars, c_only_var)) {
2362 throw_invalid_argument("add_constraint(c)",
2363 "c is not an interval constraint");
2364 }
2365
2366 // Throw an exception if c is a nontrivial strict constraint
2367 // and ITV does not support open boundaries.
2368 if (c.is_strict_inequality() && c_num_vars != 0
2369 && ITV::is_always_topologically_closed()) {
2370 throw_invalid_argument("add_constraint(c)",
2371 "c is a nontrivial strict constraint");
2372 }
2373
2374 // Avoid doing useless work if the box is known to be empty.
2375 if (marked_empty()) {
2376 return;
2377 }
2378
2379 const Coefficient& n = c.inhomogeneous_term();
2380 if (c_num_vars == 0) {
2381 // Dealing with a trivial constraint.
2382 if (n < 0
2383 || (c.is_equality() && n != 0)
2384 || (c.is_strict_inequality() && n == 0)) {
2385 set_empty();
2386 }
2387 return;
2388 }
2389
2390 PPL_ASSERT(c_num_vars == 1);
2391 const Coefficient& d = c.coefficient(Variable(c_only_var));
2392 add_interval_constraint_no_check(c_only_var, c.type(), n, d);
2393 }
2394
2395 template <typename ITV>
2396 void
add_constraints_no_check(const Constraint_System & cs)2397 Box<ITV>::add_constraints_no_check(const Constraint_System& cs) {
2398 PPL_ASSERT(cs.space_dimension() <= space_dimension());
2399 // Note: even when the box is known to be empty, we need to go
2400 // through all the constraints to fulfill the method's contract
2401 // for what concerns exception throwing.
2402 for (Constraint_System::const_iterator i = cs.begin(),
2403 cs_end = cs.end(); i != cs_end; ++i) {
2404 add_constraint_no_check(*i);
2405 }
2406 PPL_ASSERT(OK());
2407 }
2408
2409 template <typename ITV>
2410 void
add_congruence_no_check(const Congruence & cg)2411 Box<ITV>::add_congruence_no_check(const Congruence& cg) {
2412 PPL_ASSERT(cg.space_dimension() <= space_dimension());
2413
2414 // Set aside the case of proper congruences.
2415 if (cg.is_proper_congruence()) {
2416 if (cg.is_inconsistent()) {
2417 set_empty();
2418 return;
2419 }
2420 else if (cg.is_tautological()) {
2421 return;
2422 }
2423 else {
2424 throw_invalid_argument("add_congruence(cg)",
2425 "cg is a nontrivial proper congruence");
2426 }
2427 }
2428
2429 PPL_ASSERT(cg.is_equality());
2430 dimension_type cg_num_vars = 0;
2431 dimension_type cg_only_var = 0;
2432 // Throw an exception if c is not an interval congruence.
2433 if (!Box_Helpers::extract_interval_congruence(cg, cg_num_vars, cg_only_var)) {
2434 throw_invalid_argument("add_congruence(cg)",
2435 "cg is not an interval congruence");
2436 }
2437
2438 // Avoid doing useless work if the box is known to be empty.
2439 if (marked_empty()) {
2440 return;
2441 }
2442
2443 const Coefficient& n = cg.inhomogeneous_term();
2444 if (cg_num_vars == 0) {
2445 // Dealing with a trivial equality congruence.
2446 if (n != 0) {
2447 set_empty();
2448 }
2449 return;
2450 }
2451
2452 PPL_ASSERT(cg_num_vars == 1);
2453 const Coefficient& d = cg.coefficient(Variable(cg_only_var));
2454 add_interval_constraint_no_check(cg_only_var, Constraint::EQUALITY, n, d);
2455 }
2456
2457 template <typename ITV>
2458 void
add_congruences_no_check(const Congruence_System & cgs)2459 Box<ITV>::add_congruences_no_check(const Congruence_System& cgs) {
2460 PPL_ASSERT(cgs.space_dimension() <= space_dimension());
2461 // Note: even when the box is known to be empty, we need to go
2462 // through all the congruences to fulfill the method's contract
2463 // for what concerns exception throwing.
2464 for (Congruence_System::const_iterator i = cgs.begin(),
2465 cgs_end = cgs.end(); i != cgs_end; ++i) {
2466 add_congruence_no_check(*i);
2467 }
2468 PPL_ASSERT(OK());
2469 }
2470
2471 template <typename ITV>
2472 void
refine_no_check(const Constraint & c)2473 Box<ITV>::refine_no_check(const Constraint& c) {
2474 PPL_ASSERT(c.space_dimension() <= space_dimension());
2475 PPL_ASSERT(!marked_empty());
2476
2477 dimension_type c_num_vars = 0;
2478 dimension_type c_only_var = 0;
2479 // Non-interval constraints are approximated.
2480 if (!Box_Helpers::extract_interval_constraint(c, c_num_vars, c_only_var)) {
2481 propagate_constraint_no_check(c);
2482 return;
2483 }
2484
2485 const Coefficient& n = c.inhomogeneous_term();
2486 if (c_num_vars == 0) {
2487 // Dealing with a trivial constraint.
2488 if (n < 0
2489 || (c.is_equality() && n != 0)
2490 || (c.is_strict_inequality() && n == 0)) {
2491 set_empty();
2492 }
2493 return;
2494 }
2495
2496 PPL_ASSERT(c_num_vars == 1);
2497 const Coefficient& d = c.coefficient(Variable(c_only_var));
2498 add_interval_constraint_no_check(c_only_var, c.type(), n, d);
2499 }
2500
2501 template <typename ITV>
2502 void
refine_no_check(const Constraint_System & cs)2503 Box<ITV>::refine_no_check(const Constraint_System& cs) {
2504 PPL_ASSERT(cs.space_dimension() <= space_dimension());
2505 for (Constraint_System::const_iterator i = cs.begin(),
2506 cs_end = cs.end(); !marked_empty() && i != cs_end; ++i) {
2507 refine_no_check(*i);
2508 }
2509 PPL_ASSERT(OK());
2510 }
2511
2512 template <typename ITV>
2513 void
refine_no_check(const Congruence & cg)2514 Box<ITV>::refine_no_check(const Congruence& cg) {
2515 PPL_ASSERT(!marked_empty());
2516
2517 PPL_ASSERT(cg.space_dimension() <= space_dimension());
2518
2519 if (cg.is_proper_congruence()) {
2520 // A proper congruences is also an interval constraint
2521 // if and only if it is trivial.
2522 if (cg.is_inconsistent()) {
2523 set_empty();
2524 }
2525 return;
2526 }
2527
2528 PPL_ASSERT(cg.is_equality());
2529 Constraint c(cg);
2530 refine_no_check(c);
2531 }
2532
2533 template <typename ITV>
2534 void
refine_no_check(const Congruence_System & cgs)2535 Box<ITV>::refine_no_check(const Congruence_System& cgs) {
2536 PPL_ASSERT(cgs.space_dimension() <= space_dimension());
2537 for (Congruence_System::const_iterator i = cgs.begin(),
2538 cgs_end = cgs.end(); !marked_empty() && i != cgs_end; ++i) {
2539 refine_no_check(*i);
2540 }
2541 PPL_ASSERT(OK());
2542 }
2543
2544 #if 1 // Alternative implementations for propagate_constraint_no_check.
2545 namespace Implementation {
2546
2547 namespace Boxes {
2548
2549 inline bool
propagate_constraint_check_result(Result r,Ternary & open)2550 propagate_constraint_check_result(Result r, Ternary& open) {
2551 r = result_relation_class(r);
2552 switch (r) {
2553 case V_GT_MINUS_INFINITY:
2554 case V_LT_PLUS_INFINITY:
2555 return true;
2556 case V_LT:
2557 case V_GT:
2558 open = T_YES;
2559 return false;
2560 case V_LE:
2561 case V_GE:
2562 if (open == T_NO) {
2563 open = T_MAYBE;
2564 }
2565 return false;
2566 case V_EQ:
2567 return false;
2568 default:
2569 PPL_UNREACHABLE;
2570 return true;
2571 }
2572 }
2573
2574 } // namespace Boxes
2575
2576 } // namespace Implementation
2577
2578
2579 template <typename ITV>
2580 void
propagate_constraint_no_check(const Constraint & c)2581 Box<ITV>::propagate_constraint_no_check(const Constraint& c) {
2582 using namespace Implementation::Boxes;
2583
2584 PPL_ASSERT(c.space_dimension() <= space_dimension());
2585
2586 typedef
2587 typename Select_Temp_Boundary_Type<typename ITV::boundary_type>::type
2588 Temp_Boundary_Type;
2589
2590 const dimension_type c_space_dim = c.space_dimension();
2591 const Constraint::Type c_type = c.type();
2592 const Coefficient& c_inhomogeneous_term = c.inhomogeneous_term();
2593
2594 // Find a space dimension having a non-zero coefficient (if any).
2595 const dimension_type last_k
2596 = c.expression().last_nonzero(1, c_space_dim + 1);
2597 if (last_k == c_space_dim + 1) {
2598 // Constraint c is trivial: check if it is inconsistent.
2599 if (c_inhomogeneous_term < 0
2600 || (c_inhomogeneous_term == 0
2601 && c_type != Constraint::NONSTRICT_INEQUALITY)) {
2602 set_empty();
2603 }
2604 return;
2605 }
2606
2607 // Here constraint c is non-trivial.
2608 PPL_ASSERT(last_k <= c_space_dim);
2609 Temp_Boundary_Type t_bound;
2610 Temp_Boundary_Type t_a;
2611 Temp_Boundary_Type t_x;
2612 Ternary open;
2613 const Constraint::expr_type c_e = c.expression();
2614 for (Constraint::expr_type::const_iterator k = c_e.begin(),
2615 k_end = c_e.lower_bound(Variable(last_k)); k != k_end; ++k) {
2616 const Coefficient& a_k = *k;
2617 const Variable k_var = k.variable();
2618 const int sgn_a_k = sgn(a_k);
2619 if (sgn_a_k == 0) {
2620 continue;
2621 }
2622 Result r;
2623 if (sgn_a_k > 0) {
2624 open = (c_type == Constraint::STRICT_INEQUALITY) ? T_YES : T_NO;
2625 if (open == T_NO) {
2626 maybe_reset_fpu_inexact<Temp_Boundary_Type>();
2627 }
2628 r = assign_r(t_bound, c_inhomogeneous_term, ROUND_UP);
2629 if (propagate_constraint_check_result(r, open)) {
2630 goto maybe_refine_upper_1;
2631 }
2632 r = neg_assign_r(t_bound, t_bound, ROUND_DOWN);
2633 if (propagate_constraint_check_result(r, open)) {
2634 goto maybe_refine_upper_1;
2635 }
2636 for (Constraint::expr_type::const_iterator i = c_e.begin(),
2637 i_end = c_e.lower_bound(Variable(last_k)); i != i_end; ++i) {
2638 const Variable i_var = i.variable();
2639 if (i_var.id() == k_var.id()) {
2640 continue;
2641 }
2642 const Coefficient& a_i = *i;
2643 const int sgn_a_i = sgn(a_i);
2644 ITV& x_i = seq[i_var.id()];
2645 if (sgn_a_i < 0) {
2646 if (x_i.lower_is_boundary_infinity()) {
2647 goto maybe_refine_upper_1;
2648 }
2649 r = assign_r(t_a, a_i, ROUND_DOWN);
2650 if (propagate_constraint_check_result(r, open)) {
2651 goto maybe_refine_upper_1;
2652 }
2653 r = assign_r(t_x, x_i.lower(), ROUND_DOWN);
2654 if (propagate_constraint_check_result(r, open)) {
2655 goto maybe_refine_upper_1;
2656 }
2657 if (x_i.lower_is_open()) {
2658 open = T_YES;
2659 }
2660 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_DOWN);
2661 if (propagate_constraint_check_result(r, open)) {
2662 goto maybe_refine_upper_1;
2663 }
2664 }
2665 else {
2666 PPL_ASSERT(sgn_a_i > 0);
2667 if (x_i.upper_is_boundary_infinity()) {
2668 goto maybe_refine_upper_1;
2669 }
2670 r = assign_r(t_a, a_i, ROUND_UP);
2671 if (propagate_constraint_check_result(r, open)) {
2672 goto maybe_refine_upper_1;
2673 }
2674 r = assign_r(t_x, x_i.upper(), ROUND_UP);
2675 if (propagate_constraint_check_result(r, open)) {
2676 goto maybe_refine_upper_1;
2677 }
2678 if (x_i.upper_is_open()) {
2679 open = T_YES;
2680 }
2681 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_DOWN);
2682 if (propagate_constraint_check_result(r, open)) {
2683 goto maybe_refine_upper_1;
2684 }
2685 }
2686 }
2687 r = assign_r(t_a, a_k, ROUND_UP);
2688 if (propagate_constraint_check_result(r, open)) {
2689 goto maybe_refine_upper_1;
2690 }
2691 r = div_assign_r(t_bound, t_bound, t_a, ROUND_DOWN);
2692 if (propagate_constraint_check_result(r, open)) {
2693 goto maybe_refine_upper_1;
2694 }
2695
2696 // Refine the lower bound of `seq[k]' with `t_bound'.
2697 if (open == T_MAYBE
2698 && maybe_check_fpu_inexact<Temp_Boundary_Type>() == 1) {
2699 open = T_YES;
2700 }
2701 {
2702 const Relation_Symbol rel
2703 = (open == T_YES) ? GREATER_THAN : GREATER_OR_EQUAL;
2704 seq[k_var.id()].add_constraint(i_constraint(rel, t_bound));
2705 }
2706 reset_empty_up_to_date();
2707 maybe_refine_upper_1:
2708 if (c_type != Constraint::EQUALITY) {
2709 continue;
2710 }
2711 open = T_NO;
2712 maybe_reset_fpu_inexact<Temp_Boundary_Type>();
2713 r = assign_r(t_bound, c_inhomogeneous_term, ROUND_DOWN);
2714 if (propagate_constraint_check_result(r, open)) {
2715 goto next_k;
2716 }
2717 r = neg_assign_r(t_bound, t_bound, ROUND_UP);
2718 if (propagate_constraint_check_result(r, open)) {
2719 goto next_k;
2720 }
2721 for (Constraint::expr_type::const_iterator i = c_e.begin(),
2722 i_end = c_e.lower_bound(Variable(c_space_dim)); i != i_end; ++i) {
2723 const Variable i_var = i.variable();
2724 if (i_var.id() == k_var.id()) {
2725 continue;
2726 }
2727 const Coefficient& a_i = *i;
2728 const int sgn_a_i = sgn(a_i);
2729 ITV& x_i = seq[i_var.id()];
2730 if (sgn_a_i < 0) {
2731 if (x_i.upper_is_boundary_infinity()) {
2732 goto next_k;
2733 }
2734 r = assign_r(t_a, a_i, ROUND_UP);
2735 if (propagate_constraint_check_result(r, open)) {
2736 goto next_k;
2737 }
2738 r = assign_r(t_x, x_i.upper(), ROUND_UP);
2739 if (propagate_constraint_check_result(r, open)) {
2740 goto next_k;
2741 }
2742 if (x_i.upper_is_open()) {
2743 open = T_YES;
2744 }
2745 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
2746 if (propagate_constraint_check_result(r, open)) {
2747 goto next_k;
2748 }
2749 }
2750 else {
2751 PPL_ASSERT(sgn_a_i > 0);
2752 if (x_i.lower_is_boundary_infinity()) {
2753 goto next_k;
2754 }
2755 r = assign_r(t_a, a_i, ROUND_DOWN);
2756 if (propagate_constraint_check_result(r, open)) {
2757 goto next_k;
2758 }
2759 r = assign_r(t_x, x_i.lower(), ROUND_DOWN);
2760 if (propagate_constraint_check_result(r, open)) {
2761 goto next_k;
2762 }
2763 if (x_i.lower_is_open()) {
2764 open = T_YES;
2765 }
2766 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
2767 if (propagate_constraint_check_result(r, open)) {
2768 goto next_k;
2769 }
2770 }
2771 }
2772 r = assign_r(t_a, a_k, ROUND_DOWN);
2773 if (propagate_constraint_check_result(r, open)) {
2774 goto next_k;
2775 }
2776 r = div_assign_r(t_bound, t_bound, t_a, ROUND_UP);
2777 if (propagate_constraint_check_result(r, open)) {
2778 goto next_k;
2779 }
2780
2781 // Refine the upper bound of seq[k] with t_bound.
2782 if (open == T_MAYBE
2783 && maybe_check_fpu_inexact<Temp_Boundary_Type>() == 1) {
2784 open = T_YES;
2785 }
2786 const Relation_Symbol rel
2787 = (open == T_YES) ? LESS_THAN : LESS_OR_EQUAL;
2788 seq[k_var.id()].add_constraint(i_constraint(rel, t_bound));
2789 reset_empty_up_to_date();
2790 }
2791 else {
2792 PPL_ASSERT(sgn_a_k < 0);
2793 open = (c_type == Constraint::STRICT_INEQUALITY) ? T_YES : T_NO;
2794 if (open == T_NO) {
2795 maybe_reset_fpu_inexact<Temp_Boundary_Type>();
2796 }
2797 r = assign_r(t_bound, c_inhomogeneous_term, ROUND_UP);
2798 if (propagate_constraint_check_result(r, open)) {
2799 goto maybe_refine_upper_2;
2800 }
2801 r = neg_assign_r(t_bound, t_bound, ROUND_DOWN);
2802 if (propagate_constraint_check_result(r, open)) {
2803 goto maybe_refine_upper_2;
2804 }
2805 for (Constraint::expr_type::const_iterator i = c_e.begin(),
2806 i_end = c_e.lower_bound(Variable(c_space_dim)); i != i_end; ++i) {
2807 const Variable i_var = i.variable();
2808 if (i_var.id() == k_var.id()) {
2809 continue;
2810 }
2811 const Coefficient& a_i = *i;
2812 const int sgn_a_i = sgn(a_i);
2813 ITV& x_i = seq[i_var.id()];
2814 if (sgn_a_i < 0) {
2815 if (x_i.lower_is_boundary_infinity()) {
2816 goto maybe_refine_upper_2;
2817 }
2818 r = assign_r(t_a, a_i, ROUND_DOWN);
2819 if (propagate_constraint_check_result(r, open)) {
2820 goto maybe_refine_upper_2;
2821 }
2822 r = assign_r(t_x, x_i.lower(), ROUND_DOWN);
2823 if (propagate_constraint_check_result(r, open)) {
2824 goto maybe_refine_upper_2;
2825 }
2826 if (x_i.lower_is_open()) {
2827 open = T_YES;
2828 }
2829 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
2830 if (propagate_constraint_check_result(r, open)) {
2831 goto maybe_refine_upper_2;
2832 }
2833 }
2834 else {
2835 PPL_ASSERT(sgn_a_i > 0);
2836 if (x_i.upper_is_boundary_infinity()) {
2837 goto maybe_refine_upper_2;
2838 }
2839 r = assign_r(t_a, a_i, ROUND_UP);
2840 if (propagate_constraint_check_result(r, open)) {
2841 goto maybe_refine_upper_2;
2842 }
2843 r = assign_r(t_x, x_i.upper(), ROUND_UP);
2844 if (propagate_constraint_check_result(r, open)) {
2845 goto maybe_refine_upper_2;
2846 }
2847 if (x_i.upper_is_open()) {
2848 open = T_YES;
2849 }
2850 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_DOWN);
2851 if (propagate_constraint_check_result(r, open)) {
2852 goto maybe_refine_upper_2;
2853 }
2854 }
2855 }
2856 r = assign_r(t_a, a_k, ROUND_UP);
2857 if (propagate_constraint_check_result(r, open)) {
2858 goto maybe_refine_upper_2;
2859 }
2860 r = div_assign_r(t_bound, t_bound, t_a, ROUND_UP);
2861 if (propagate_constraint_check_result(r, open)) {
2862 goto maybe_refine_upper_2;
2863 }
2864 // Refine the upper bound of seq[k] with t_bound.
2865 if (open == T_MAYBE
2866 && maybe_check_fpu_inexact<Temp_Boundary_Type>() == 1) {
2867 open = T_YES;
2868 }
2869 {
2870 const Relation_Symbol rel
2871 = (open == T_YES) ? LESS_THAN : LESS_OR_EQUAL;
2872 seq[k_var.id()].add_constraint(i_constraint(rel, t_bound));
2873 }
2874 reset_empty_up_to_date();
2875 maybe_refine_upper_2:
2876 if (c_type != Constraint::EQUALITY) {
2877 continue;
2878 }
2879 open = T_NO;
2880 maybe_reset_fpu_inexact<Temp_Boundary_Type>();
2881 r = assign_r(t_bound, c_inhomogeneous_term, ROUND_DOWN);
2882 if (propagate_constraint_check_result(r, open)) {
2883 goto next_k;
2884 }
2885 r = neg_assign_r(t_bound, t_bound, ROUND_UP);
2886 if (propagate_constraint_check_result(r, open)) {
2887 goto next_k;
2888 }
2889 for (Constraint::expr_type::const_iterator i = c_e.begin(),
2890 i_end = c_e.lower_bound(Variable(c_space_dim)); i != i_end; ++i) {
2891 const Variable i_var = i.variable();
2892 if (i_var.id() == k_var.id()) {
2893 continue;
2894 }
2895 const Coefficient& a_i = *i;
2896 const int sgn_a_i = sgn(a_i);
2897 ITV& x_i = seq[i_var.id()];
2898 if (sgn_a_i < 0) {
2899 if (x_i.upper_is_boundary_infinity()) {
2900 goto next_k;
2901 }
2902 r = assign_r(t_a, a_i, ROUND_UP);
2903 if (propagate_constraint_check_result(r, open)) {
2904 goto next_k;
2905 }
2906 r = assign_r(t_x, x_i.upper(), ROUND_UP);
2907 if (propagate_constraint_check_result(r, open)) {
2908 goto next_k;
2909 }
2910 if (x_i.upper_is_open()) {
2911 open = T_YES;
2912 }
2913 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
2914 if (propagate_constraint_check_result(r, open)) {
2915 goto next_k;
2916 }
2917 }
2918 else {
2919 PPL_ASSERT(sgn_a_i > 0);
2920 if (x_i.lower_is_boundary_infinity()) {
2921 goto next_k;
2922 }
2923 r = assign_r(t_a, a_i, ROUND_DOWN);
2924 if (propagate_constraint_check_result(r, open)) {
2925 goto next_k;
2926 }
2927 r = assign_r(t_x, x_i.lower(), ROUND_DOWN);
2928 if (propagate_constraint_check_result(r, open)) {
2929 goto next_k;
2930 }
2931 if (x_i.lower_is_open()) {
2932 open = T_YES;
2933 }
2934 r = sub_mul_assign_r(t_bound, t_a, t_x, ROUND_UP);
2935 if (propagate_constraint_check_result(r, open)) {
2936 goto next_k;
2937 }
2938 }
2939 }
2940 r = assign_r(t_a, a_k, ROUND_DOWN);
2941 if (propagate_constraint_check_result(r, open)) {
2942 goto next_k;
2943 }
2944 r = div_assign_r(t_bound, t_bound, t_a, ROUND_DOWN);
2945 if (propagate_constraint_check_result(r, open)) {
2946 goto next_k;
2947 }
2948
2949 // Refine the lower bound of seq[k] with t_bound.
2950 if (open == T_MAYBE
2951 && maybe_check_fpu_inexact<Temp_Boundary_Type>() == 1) {
2952 open = T_YES;
2953 }
2954 const Relation_Symbol rel
2955 = (open == T_YES) ? GREATER_THAN : GREATER_OR_EQUAL;
2956 seq[k_var.id()].add_constraint(i_constraint(rel, t_bound));
2957 reset_empty_up_to_date();
2958 }
2959 next_k:
2960 ;
2961 }
2962 }
2963
2964 #else // Alternative implementations for propagate_constraint_no_check.
2965
2966 template <typename ITV>
2967 void
propagate_constraint_no_check(const Constraint & c)2968 Box<ITV>::propagate_constraint_no_check(const Constraint& c) {
2969 PPL_ASSERT(c.space_dimension() <= space_dimension());
2970
2971 dimension_type c_space_dim = c.space_dimension();
2972 ITV k[c_space_dim];
2973 ITV p[c_space_dim];
2974 for (Constraint::expr_type::const_iterator i = c_e.begin(),
2975 i_end = c_e.lower_bound(Variable(c_space_dim)); i != i_end; ++i) {
2976 const Variable i_var = i.variable();
2977 k[i_var.id()] = *i;
2978 ITV& p_i = p[i_var.id()];
2979 p_i = seq[i_var.id()];
2980 p_i.mul_assign(p_i, k[i_var.id()]);
2981 }
2982 const Coefficient& inhomogeneous_term = c.inhomogeneous_term();
2983 for (Constraint::expr_type::const_iterator i = c_e.begin(),
2984 i_end = c_e.lower_bound(Variable(c_space_dim)); i != i_end; ++i) {
2985 const Variable i_var = i.variable();
2986 int sgn_coefficient_i = sgn(*i);
2987 ITV q(inhomogeneous_term);
2988 for (Constraint::expr_type::const_iterator j = c_e.begin(),
2989 j_end = c_e.lower_bound(Variable(c_space_dim)); j != j_end; ++j) {
2990 const Variable j_var = j.variable();
2991 if (i_var == j_var) {
2992 continue;
2993 }
2994 q.add_assign(q, p[j_var.id()]);
2995 }
2996 q.div_assign(q, k[i_var.id()]);
2997 q.neg_assign(q);
2998 Relation_Symbol rel;
2999 switch (c.type()) {
3000 case Constraint::EQUALITY:
3001 rel = EQUAL;
3002 break;
3003 case Constraint::NONSTRICT_INEQUALITY:
3004 rel = (sgn_coefficient_i > 0) ? GREATER_OR_EQUAL : LESS_OR_EQUAL;
3005 break;
3006 case Constraint::STRICT_INEQUALITY:
3007 rel = (sgn_coefficient_i > 0) ? GREATER_THAN : LESS_THAN;
3008 break;
3009 }
3010 seq[i_var.id()].add_constraint(i_constraint(rel, q));
3011 // FIXME: could/should we exploit the return value of add_constraint
3012 // in case it is available?
3013 // FIXME: should we instead be lazy and do not even bother about
3014 // the possibility the interval becomes empty apart from setting
3015 // empty_up_to_date = false?
3016 if (seq[i_var.id()].is_empty()) {
3017 set_empty();
3018 break;
3019 }
3020 }
3021
3022 PPL_ASSERT(OK());
3023 }
3024
3025 #endif // Alternative implementations for propagate_constraint_no_check.
3026
3027 template <typename ITV>
3028 void
3029 Box<ITV>
propagate_constraints_no_check(const Constraint_System & cs,const dimension_type max_iterations)3030 ::propagate_constraints_no_check(const Constraint_System& cs,
3031 const dimension_type max_iterations) {
3032 const dimension_type space_dim = space_dimension();
3033 PPL_ASSERT(cs.space_dimension() <= space_dim);
3034
3035 const Constraint_System::const_iterator cs_begin = cs.begin();
3036 const Constraint_System::const_iterator cs_end = cs.end();
3037 const dimension_type propagation_weight
3038 = Implementation::num_constraints(cs) * space_dim;
3039
3040 Sequence copy;
3041 bool changed;
3042 dimension_type num_iterations = 0;
3043 do {
3044 WEIGHT_BEGIN();
3045 ++num_iterations;
3046 copy = seq;
3047 for (Constraint_System::const_iterator i = cs_begin; i != cs_end; ++i) {
3048 propagate_constraint_no_check(*i);
3049 }
3050
3051 WEIGHT_ADD_MUL(40, propagation_weight);
3052 // Check if the client has requested abandoning all expensive
3053 // computations. If so, the exception specified by the client
3054 // is thrown now.
3055 maybe_abandon();
3056
3057 // NOTE: if max_iterations == 0 (i.e., no iteration limit is set)
3058 // the following test will anyway trigger on wrap around.
3059 if (num_iterations == max_iterations) {
3060 break;
3061 }
3062
3063 changed = (copy != seq);
3064 } while (changed);
3065 }
3066
3067 template <typename ITV>
3068 void
affine_image(const Variable var,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)3069 Box<ITV>::affine_image(const Variable var,
3070 const Linear_Expression& expr,
3071 Coefficient_traits::const_reference denominator) {
3072 // The denominator cannot be zero.
3073 if (denominator == 0) {
3074 throw_invalid_argument("affine_image(v, e, d)", "d == 0");
3075 }
3076
3077 // Dimension-compatibility checks.
3078 const dimension_type space_dim = space_dimension();
3079 const dimension_type expr_space_dim = expr.space_dimension();
3080 if (space_dim < expr_space_dim) {
3081 throw_dimension_incompatible("affine_image(v, e, d)", "e", expr);
3082 }
3083 // `var' should be one of the dimensions of the polyhedron.
3084 const dimension_type var_space_dim = var.space_dimension();
3085 if (space_dim < var_space_dim) {
3086 throw_dimension_incompatible("affine_image(v, e, d)", "v", var);
3087 }
3088
3089 if (is_empty()) {
3090 return;
3091 }
3092
3093 Tmp_Interval_Type expr_value;
3094 Tmp_Interval_Type temp0;
3095 Tmp_Interval_Type temp1;
3096 expr_value.assign(expr.inhomogeneous_term());
3097 for (Linear_Expression::const_iterator i = expr.begin(),
3098 i_end = expr.end(); i != i_end; ++i) {
3099 temp0.assign(*i);
3100 temp1.assign(seq[i.variable().id()]);
3101 temp0.mul_assign(temp0, temp1);
3102 expr_value.add_assign(expr_value, temp0);
3103 }
3104 if (denominator != 1) {
3105 temp0.assign(denominator);
3106 expr_value.div_assign(expr_value, temp0);
3107 }
3108 seq[var.id()].assign(expr_value);
3109
3110 PPL_ASSERT(OK());
3111 }
3112
3113 template <typename ITV>
3114 void
affine_form_image(const Variable var,const Linear_Form<ITV> & lf)3115 Box<ITV>::affine_form_image(const Variable var,
3116 const Linear_Form<ITV>& lf) {
3117
3118 // Check that ITV has a floating point boundary type.
3119 PPL_COMPILE_TIME_CHECK(!std::numeric_limits<typename ITV::boundary_type>
3120 ::is_exact, "Box<ITV>::affine_form_image(Variable, Linear_Form):"
3121 "ITV has not a floating point boundary type.");
3122
3123 // Dimension-compatibility checks.
3124 const dimension_type space_dim = space_dimension();
3125 const dimension_type lf_space_dim = lf.space_dimension();
3126 if (space_dim < lf_space_dim) {
3127 throw_dimension_incompatible("affine_form_image(var, lf)", "lf", lf);
3128 }
3129
3130 // `var' should be one of the dimensions of the polyhedron.
3131 const dimension_type var_space_dim = var.space_dimension();
3132 if (space_dim < var_space_dim) {
3133 throw_dimension_incompatible("affine_form_image(var, lf)", "var", var);
3134 }
3135
3136 if (is_empty()) {
3137 return;
3138 }
3139
3140 // Intervalization of 'lf'.
3141 ITV result = lf.inhomogeneous_term();
3142 for (dimension_type i = 0; i < lf_space_dim; ++i) {
3143 ITV current_addend = lf.coefficient(Variable(i));
3144 const ITV& curr_int = seq[i];
3145 current_addend *= curr_int;
3146 result += current_addend;
3147 }
3148
3149 seq[var.id()].assign(result);
3150 PPL_ASSERT(OK());
3151 }
3152
3153 template <typename ITV>
3154 void
affine_preimage(const Variable var,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)3155 Box<ITV>::affine_preimage(const Variable var,
3156 const Linear_Expression& expr,
3157 Coefficient_traits::const_reference
3158 denominator) {
3159 // The denominator cannot be zero.
3160 if (denominator == 0) {
3161 throw_invalid_argument("affine_preimage(v, e, d)", "d == 0");
3162 }
3163
3164 // Dimension-compatibility checks.
3165 const dimension_type x_space_dim = space_dimension();
3166 const dimension_type expr_space_dim = expr.space_dimension();
3167 if (x_space_dim < expr_space_dim) {
3168 throw_dimension_incompatible("affine_preimage(v, e, d)", "e", expr);
3169 }
3170 // `var' should be one of the dimensions of the polyhedron.
3171 const dimension_type var_space_dim = var.space_dimension();
3172 if (x_space_dim < var_space_dim) {
3173 throw_dimension_incompatible("affine_preimage(v, e, d)", "v", var);
3174 }
3175
3176 if (is_empty()) {
3177 return;
3178 }
3179
3180 const Coefficient& expr_v = expr.coefficient(var);
3181 const bool invertible = (expr_v != 0);
3182 if (!invertible) {
3183 Tmp_Interval_Type expr_value;
3184 Tmp_Interval_Type temp0;
3185 Tmp_Interval_Type temp1;
3186 expr_value.assign(expr.inhomogeneous_term());
3187 for (Linear_Expression::const_iterator i = expr.begin(),
3188 i_end = expr.end(); i != i_end; ++i) {
3189 temp0.assign(*i);
3190 temp1.assign(seq[i.variable().id()]);
3191 temp0.mul_assign(temp0, temp1);
3192 expr_value.add_assign(expr_value, temp0);
3193 }
3194 if (denominator != 1) {
3195 temp0.assign(denominator);
3196 expr_value.div_assign(expr_value, temp0);
3197 }
3198 ITV& x_seq_v = seq[var.id()];
3199 expr_value.intersect_assign(x_seq_v);
3200 if (expr_value.is_empty()) {
3201 set_empty();
3202 }
3203 else {
3204 x_seq_v.assign(UNIVERSE);
3205 }
3206 }
3207 else {
3208 // The affine transformation is invertible.
3209 // CHECKME: for efficiency, would it be meaningful to avoid
3210 // the computation of inverse by partially evaluating the call
3211 // to affine_image?
3212 Linear_Expression inverse;
3213 inverse -= expr;
3214 inverse += (expr_v + denominator) * var;
3215 affine_image(var, inverse, expr_v);
3216 }
3217 PPL_ASSERT(OK());
3218 }
3219
3220 template <typename ITV>
3221 void
3222 Box<ITV>
bounded_affine_image(const Variable var,const Linear_Expression & lb_expr,const Linear_Expression & ub_expr,Coefficient_traits::const_reference denominator)3223 ::bounded_affine_image(const Variable var,
3224 const Linear_Expression& lb_expr,
3225 const Linear_Expression& ub_expr,
3226 Coefficient_traits::const_reference denominator) {
3227 // The denominator cannot be zero.
3228 if (denominator == 0) {
3229 throw_invalid_argument("bounded_affine_image(v, lb, ub, d)", "d == 0");
3230 }
3231
3232 // Dimension-compatibility checks.
3233 const dimension_type space_dim = space_dimension();
3234 // The dimension of `lb_expr' and `ub_expr' should not be
3235 // greater than the dimension of `*this'.
3236 const dimension_type lb_space_dim = lb_expr.space_dimension();
3237 if (space_dim < lb_space_dim) {
3238 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
3239 "lb", lb_expr);
3240 }
3241 const dimension_type ub_space_dim = ub_expr.space_dimension();
3242 if (space_dim < ub_space_dim) {
3243 throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
3244 "ub", ub_expr);
3245 }
3246 // `var' should be one of the dimensions of the box.
3247 const dimension_type var_space_dim = var.space_dimension();
3248 if (space_dim < var_space_dim) {
3249 throw_dimension_incompatible("affine_image(v, e, d)", "v", var);
3250 }
3251 // Any image of an empty box is empty.
3252 if (is_empty()) {
3253 return;
3254 }
3255 // Add the constraint implied by the `lb_expr' and `ub_expr'.
3256 if (denominator > 0) {
3257 refine_with_constraint(lb_expr <= ub_expr);
3258 }
3259 else {
3260 refine_with_constraint(lb_expr >= ub_expr);
3261 }
3262
3263 // Check whether `var' occurs in `lb_expr' and/or `ub_expr'.
3264 if (lb_expr.coefficient(var) == 0) {
3265 // Here `var' can only occur in `ub_expr'.
3266 generalized_affine_image(var,
3267 LESS_OR_EQUAL,
3268 ub_expr,
3269 denominator);
3270 if (denominator > 0) {
3271 refine_with_constraint(lb_expr <= denominator*var);
3272 }
3273 else {
3274 refine_with_constraint(denominator*var <= lb_expr);
3275 }
3276 }
3277 else if (ub_expr.coefficient(var) == 0) {
3278 // Here `var' can only occur in `lb_expr'.
3279 generalized_affine_image(var,
3280 GREATER_OR_EQUAL,
3281 lb_expr,
3282 denominator);
3283 if (denominator > 0) {
3284 refine_with_constraint(denominator*var <= ub_expr);
3285 }
3286 else {
3287 refine_with_constraint(ub_expr <= denominator*var);
3288 }
3289 }
3290 else {
3291 // Here `var' occurs in both `lb_expr' and `ub_expr'. As boxes
3292 // can only use the non-relational constraints, we find the
3293 // maximum/minimum values `ub_expr' and `lb_expr' obtain with the
3294 // box and use these instead of the `ub-expr' and `lb-expr'.
3295 PPL_DIRTY_TEMP_COEFFICIENT(max_numer);
3296 PPL_DIRTY_TEMP_COEFFICIENT(max_denom);
3297 bool max_included;
3298 PPL_DIRTY_TEMP_COEFFICIENT(min_numer);
3299 PPL_DIRTY_TEMP_COEFFICIENT(min_denom);
3300 bool min_included;
3301 ITV& seq_v = seq[var.id()];
3302 if (maximize(ub_expr, max_numer, max_denom, max_included)) {
3303 if (minimize(lb_expr, min_numer, min_denom, min_included)) {
3304 // The `ub_expr' has a maximum value and the `lb_expr'
3305 // has a minimum value for the box.
3306 // Set the bounds for `var' using the minimum for `lb_expr'.
3307 min_denom *= denominator;
3308 PPL_DIRTY_TEMP(mpq_class, q1);
3309 PPL_DIRTY_TEMP(mpq_class, q2);
3310 assign_r(q1.get_num(), min_numer, ROUND_NOT_NEEDED);
3311 assign_r(q1.get_den(), min_denom, ROUND_NOT_NEEDED);
3312 q1.canonicalize();
3313 // Now make the maximum of lb_expr the upper bound. If the
3314 // maximum is not at a box point, then inequality is strict.
3315 max_denom *= denominator;
3316 assign_r(q2.get_num(), max_numer, ROUND_NOT_NEEDED);
3317 assign_r(q2.get_den(), max_denom, ROUND_NOT_NEEDED);
3318 q2.canonicalize();
3319
3320 if (denominator > 0) {
3321 Relation_Symbol gr = min_included ? GREATER_OR_EQUAL : GREATER_THAN;
3322 Relation_Symbol lr = max_included ? LESS_OR_EQUAL : LESS_THAN;
3323 seq_v.build(i_constraint(gr, q1), i_constraint(lr, q2));
3324 }
3325 else {
3326 Relation_Symbol gr = max_included ? GREATER_OR_EQUAL : GREATER_THAN;
3327 Relation_Symbol lr = min_included ? LESS_OR_EQUAL : LESS_THAN;
3328 seq_v.build(i_constraint(gr, q2), i_constraint(lr, q1));
3329 }
3330 }
3331 else {
3332 // The `ub_expr' has a maximum value but the `lb_expr'
3333 // has no minimum value for the box.
3334 // Set the bounds for `var' using the maximum for `lb_expr'.
3335 PPL_DIRTY_TEMP(mpq_class, q);
3336 max_denom *= denominator;
3337 assign_r(q.get_num(), max_numer, ROUND_NOT_NEEDED);
3338 assign_r(q.get_den(), max_denom, ROUND_NOT_NEEDED);
3339 q.canonicalize();
3340 Relation_Symbol rel = (denominator > 0)
3341 ? (max_included ? LESS_OR_EQUAL : LESS_THAN)
3342 : (max_included ? GREATER_OR_EQUAL : GREATER_THAN);
3343 seq_v.build(i_constraint(rel, q));
3344 }
3345 }
3346 else if (minimize(lb_expr, min_numer, min_denom, min_included)) {
3347 // The `ub_expr' has no maximum value but the `lb_expr'
3348 // has a minimum value for the box.
3349 // Set the bounds for `var' using the minimum for `lb_expr'.
3350 min_denom *= denominator;
3351 PPL_DIRTY_TEMP(mpq_class, q);
3352 assign_r(q.get_num(), min_numer, ROUND_NOT_NEEDED);
3353 assign_r(q.get_den(), min_denom, ROUND_NOT_NEEDED);
3354 q.canonicalize();
3355
3356 Relation_Symbol rel = (denominator > 0)
3357 ? (min_included ? GREATER_OR_EQUAL : GREATER_THAN)
3358 : (min_included ? LESS_OR_EQUAL : LESS_THAN);
3359 seq_v.build(i_constraint(rel, q));
3360 }
3361 else {
3362 // The `ub_expr' has no maximum value and the `lb_expr'
3363 // has no minimum value for the box.
3364 // So we set the bounds to be unbounded.
3365 seq_v.assign(UNIVERSE);
3366 }
3367 }
3368 PPL_ASSERT(OK());
3369 }
3370
3371 template <typename ITV>
3372 void
3373 Box<ITV>
bounded_affine_preimage(const Variable var,const Linear_Expression & lb_expr,const Linear_Expression & ub_expr,Coefficient_traits::const_reference denominator)3374 ::bounded_affine_preimage(const Variable var,
3375 const Linear_Expression& lb_expr,
3376 const Linear_Expression& ub_expr,
3377 Coefficient_traits::const_reference denominator) {
3378 // The denominator cannot be zero.
3379 const dimension_type space_dim = space_dimension();
3380 if (denominator == 0) {
3381 throw_invalid_argument("bounded_affine_preimage(v, lb, ub, d)", "d == 0");
3382 }
3383
3384 // Dimension-compatibility checks.
3385 // `var' should be one of the dimensions of the polyhedron.
3386 const dimension_type var_space_dim = var.space_dimension();
3387 if (space_dim < var_space_dim) {
3388 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
3389 "v", var);
3390 }
3391 // The dimension of `lb_expr' and `ub_expr' should not be
3392 // greater than the dimension of `*this'.
3393 const dimension_type lb_space_dim = lb_expr.space_dimension();
3394 if (space_dim < lb_space_dim) {
3395 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
3396 "lb", lb_expr);
3397 }
3398 const dimension_type ub_space_dim = ub_expr.space_dimension();
3399 if (space_dim < ub_space_dim) {
3400 throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
3401 "ub", ub_expr);
3402 }
3403
3404 // Any preimage of an empty polyhedron is empty.
3405 if (marked_empty()) {
3406 return;
3407 }
3408
3409 const bool negative_denom = (denominator < 0);
3410 const Coefficient& lb_var_coeff = lb_expr.coefficient(var);
3411 const Coefficient& ub_var_coeff = ub_expr.coefficient(var);
3412
3413 // If the implied constraint between `ub_expr and `lb_expr' is
3414 // independent of `var', then impose it now.
3415 if (lb_var_coeff == ub_var_coeff) {
3416 if (negative_denom) {
3417 refine_with_constraint(lb_expr >= ub_expr);
3418 }
3419 else {
3420 refine_with_constraint(lb_expr <= ub_expr);
3421 }
3422 }
3423
3424 ITV& seq_var = seq[var.id()];
3425 if (!seq_var.is_universe()) {
3426 // We want to work with a positive denominator,
3427 // so the sign and its (unsigned) value are separated.
3428 PPL_DIRTY_TEMP_COEFFICIENT(pos_denominator);
3429 pos_denominator = denominator;
3430 if (negative_denom) {
3431 neg_assign(pos_denominator, pos_denominator);
3432 }
3433 // Store all the information about the upper and lower bounds
3434 // for `var' before making this interval unbounded.
3435 bool open_lower = seq_var.lower_is_open();
3436 bool unbounded_lower = seq_var.lower_is_boundary_infinity();
3437 PPL_DIRTY_TEMP(mpq_class, q_seq_var_lower);
3438 PPL_DIRTY_TEMP_COEFFICIENT(numer_lower);
3439 PPL_DIRTY_TEMP_COEFFICIENT(denom_lower);
3440 if (!unbounded_lower) {
3441 assign_r(q_seq_var_lower, seq_var.lower(), ROUND_NOT_NEEDED);
3442 assign_r(numer_lower, q_seq_var_lower.get_num(), ROUND_NOT_NEEDED);
3443 assign_r(denom_lower, q_seq_var_lower.get_den(), ROUND_NOT_NEEDED);
3444 if (negative_denom) {
3445 neg_assign(denom_lower, denom_lower);
3446 }
3447 numer_lower *= pos_denominator;
3448 seq_var.lower_extend();
3449 }
3450 bool open_upper = seq_var.upper_is_open();
3451 bool unbounded_upper = seq_var.upper_is_boundary_infinity();
3452 PPL_DIRTY_TEMP(mpq_class, q_seq_var_upper);
3453 PPL_DIRTY_TEMP_COEFFICIENT(numer_upper);
3454 PPL_DIRTY_TEMP_COEFFICIENT(denom_upper);
3455 if (!unbounded_upper) {
3456 assign_r(q_seq_var_upper, seq_var.upper(), ROUND_NOT_NEEDED);
3457 assign_r(numer_upper, q_seq_var_upper.get_num(), ROUND_NOT_NEEDED);
3458 assign_r(denom_upper, q_seq_var_upper.get_den(), ROUND_NOT_NEEDED);
3459 if (negative_denom) {
3460 neg_assign(denom_upper, denom_upper);
3461 }
3462 numer_upper *= pos_denominator;
3463 seq_var.upper_extend();
3464 }
3465
3466 if (!unbounded_lower) {
3467 // `lb_expr' is revised by removing the `var' component,
3468 // multiplying by `-' denominator of the lower bound for `var',
3469 // and adding the lower bound for `var' to the inhomogeneous term.
3470 Linear_Expression revised_lb_expr(ub_expr);
3471 revised_lb_expr -= ub_var_coeff * var;
3472 PPL_DIRTY_TEMP_COEFFICIENT(d);
3473 neg_assign(d, denom_lower);
3474 revised_lb_expr *= d;
3475 revised_lb_expr += numer_lower;
3476
3477 // Find the minimum value for the revised lower bound expression
3478 // and use this to refine the appropriate bound.
3479 bool included;
3480 PPL_DIRTY_TEMP_COEFFICIENT(denom);
3481 if (minimize(revised_lb_expr, numer_lower, denom, included)) {
3482 denom_lower *= (denom * ub_var_coeff);
3483 PPL_DIRTY_TEMP(mpq_class, q);
3484 assign_r(q.get_num(), numer_lower, ROUND_NOT_NEEDED);
3485 assign_r(q.get_den(), denom_lower, ROUND_NOT_NEEDED);
3486 q.canonicalize();
3487 if (!included) {
3488 open_lower = true;
3489 }
3490 Relation_Symbol rel;
3491 if ((ub_var_coeff >= 0) ? !negative_denom : negative_denom) {
3492 rel = open_lower ? GREATER_THAN : GREATER_OR_EQUAL;
3493 }
3494 else {
3495 rel = open_lower ? LESS_THAN : LESS_OR_EQUAL;
3496 }
3497 seq_var.add_constraint(i_constraint(rel, q));
3498 if (seq_var.is_empty()) {
3499 set_empty();
3500 return;
3501 }
3502 }
3503 }
3504
3505 if (!unbounded_upper) {
3506 // `ub_expr' is revised by removing the `var' component,
3507 // multiplying by `-' denominator of the upper bound for `var',
3508 // and adding the upper bound for `var' to the inhomogeneous term.
3509 Linear_Expression revised_ub_expr(lb_expr);
3510 revised_ub_expr -= lb_var_coeff * var;
3511 PPL_DIRTY_TEMP_COEFFICIENT(d);
3512 neg_assign(d, denom_upper);
3513 revised_ub_expr *= d;
3514 revised_ub_expr += numer_upper;
3515
3516 // Find the maximum value for the revised upper bound expression
3517 // and use this to refine the appropriate bound.
3518 bool included;
3519 PPL_DIRTY_TEMP_COEFFICIENT(denom);
3520 if (maximize(revised_ub_expr, numer_upper, denom, included)) {
3521 denom_upper *= (denom * lb_var_coeff);
3522 PPL_DIRTY_TEMP(mpq_class, q);
3523 assign_r(q.get_num(), numer_upper, ROUND_NOT_NEEDED);
3524 assign_r(q.get_den(), denom_upper, ROUND_NOT_NEEDED);
3525 q.canonicalize();
3526 if (!included) {
3527 open_upper = true;
3528 }
3529 Relation_Symbol rel;
3530 if ((lb_var_coeff >= 0) ? !negative_denom : negative_denom) {
3531 rel = open_upper ? LESS_THAN : LESS_OR_EQUAL;
3532 }
3533 else {
3534 rel = open_upper ? GREATER_THAN : GREATER_OR_EQUAL;
3535 }
3536 seq_var.add_constraint(i_constraint(rel, q));
3537 if (seq_var.is_empty()) {
3538 set_empty();
3539 return;
3540 }
3541 }
3542 }
3543 }
3544
3545 // If the implied constraint between `ub_expr and `lb_expr' is
3546 // dependent on `var', then impose on the new box.
3547 if (lb_var_coeff != ub_var_coeff) {
3548 if (denominator > 0) {
3549 refine_with_constraint(lb_expr <= ub_expr);
3550 }
3551 else {
3552 refine_with_constraint(lb_expr >= ub_expr);
3553 }
3554 }
3555
3556 PPL_ASSERT(OK());
3557 }
3558
3559 template <typename ITV>
3560 void
3561 Box<ITV>
generalized_affine_image(const Variable var,const Relation_Symbol relsym,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)3562 ::generalized_affine_image(const Variable var,
3563 const Relation_Symbol relsym,
3564 const Linear_Expression& expr,
3565 Coefficient_traits::const_reference denominator) {
3566 // The denominator cannot be zero.
3567 if (denominator == 0) {
3568 throw_invalid_argument("generalized_affine_image(v, r, e, d)", "d == 0");
3569 }
3570
3571 // Dimension-compatibility checks.
3572 const dimension_type space_dim = space_dimension();
3573 // The dimension of `expr' should not be greater than the dimension
3574 // of `*this'.
3575 if (space_dim < expr.space_dimension()) {
3576 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
3577 "e", expr);
3578 }
3579 // `var' should be one of the dimensions of the box.
3580 const dimension_type var_space_dim = var.space_dimension();
3581 if (space_dim < var_space_dim) {
3582 throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
3583 "v", var);
3584 }
3585
3586 // The relation symbol cannot be a disequality.
3587 if (relsym == NOT_EQUAL) {
3588 throw_invalid_argument("generalized_affine_image(v, r, e, d)",
3589 "r is the disequality relation symbol");
3590 }
3591
3592 // First compute the affine image.
3593 affine_image(var, expr, denominator);
3594
3595 if (relsym == EQUAL) {
3596 // The affine relation is indeed an affine function.
3597 return;
3598 }
3599 // Any image of an empty box is empty.
3600 if (is_empty()) {
3601 return;
3602 }
3603
3604 ITV& seq_var = seq[var.id()];
3605 switch (relsym) {
3606 case LESS_OR_EQUAL:
3607 seq_var.lower_extend();
3608 break;
3609 case LESS_THAN:
3610 seq_var.lower_extend();
3611 if (!seq_var.upper_is_boundary_infinity()) {
3612 seq_var.remove_sup();
3613 }
3614 break;
3615 case GREATER_OR_EQUAL:
3616 seq_var.upper_extend();
3617 break;
3618 case GREATER_THAN:
3619 seq_var.upper_extend();
3620 if (!seq_var.lower_is_boundary_infinity()) {
3621 seq_var.remove_inf();
3622 }
3623 break;
3624 default:
3625 // The EQUAL and NOT_EQUAL cases have been already dealt with.
3626 PPL_UNREACHABLE;
3627 break;
3628 }
3629 PPL_ASSERT(OK());
3630 }
3631
3632 template <typename ITV>
3633 void
3634 Box<ITV>
generalized_affine_preimage(const Variable var,const Relation_Symbol relsym,const Linear_Expression & expr,Coefficient_traits::const_reference denominator)3635 ::generalized_affine_preimage(const Variable var,
3636 const Relation_Symbol relsym,
3637 const Linear_Expression& expr,
3638 Coefficient_traits::const_reference denominator)
3639 {
3640 // The denominator cannot be zero.
3641 if (denominator == 0) {
3642 throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
3643 "d == 0");
3644 }
3645 // Dimension-compatibility checks.
3646 const dimension_type space_dim = space_dimension();
3647 // The dimension of `expr' should not be greater than the dimension
3648 // of `*this'.
3649 if (space_dim < expr.space_dimension()) {
3650 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
3651 "e", expr);
3652 }
3653 // `var' should be one of the dimensions of the box.
3654 const dimension_type var_space_dim = var.space_dimension();
3655 if (space_dim < var_space_dim) {
3656 throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
3657 "v", var);
3658 }
3659 // The relation symbol cannot be a disequality.
3660 if (relsym == NOT_EQUAL) {
3661 throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
3662 "r is the disequality relation symbol");
3663 }
3664
3665 // Check whether the affine relation is indeed an affine function.
3666 if (relsym == EQUAL) {
3667 affine_preimage(var, expr, denominator);
3668 return;
3669 }
3670
3671 // Compute the reversed relation symbol to simplify later coding.
3672 Relation_Symbol reversed_relsym;
3673 switch (relsym) {
3674 case LESS_THAN:
3675 reversed_relsym = GREATER_THAN;
3676 break;
3677 case LESS_OR_EQUAL:
3678 reversed_relsym = GREATER_OR_EQUAL;
3679 break;
3680 case GREATER_OR_EQUAL:
3681 reversed_relsym = LESS_OR_EQUAL;
3682 break;
3683 case GREATER_THAN:
3684 reversed_relsym = LESS_THAN;
3685 break;
3686 default:
3687 // The EQUAL and NOT_EQUAL cases have been already dealt with.
3688 PPL_UNREACHABLE;
3689 break;
3690 }
3691
3692 // Check whether the preimage of this affine relation can be easily
3693 // computed as the image of its inverse relation.
3694 const Coefficient& var_coefficient = expr.coefficient(var);
3695 if (var_coefficient != 0) {
3696 Linear_Expression inverse_expr
3697 = expr - (denominator + var_coefficient) * var;
3698 PPL_DIRTY_TEMP_COEFFICIENT(inverse_denominator);
3699 neg_assign(inverse_denominator, var_coefficient);
3700 Relation_Symbol inverse_relsym
3701 = (sgn(denominator) == sgn(inverse_denominator))
3702 ? relsym
3703 : reversed_relsym;
3704 generalized_affine_image(var, inverse_relsym, inverse_expr,
3705 inverse_denominator);
3706 return;
3707 }
3708
3709 // Here `var_coefficient == 0', so that the preimage cannot
3710 // be easily computed by inverting the affine relation.
3711 // Shrink the box by adding the constraint induced
3712 // by the affine relation.
3713 // First, compute the maximum and minimum value reached by
3714 // `denominator*var' on the box as we need to use non-relational
3715 // expressions.
3716 PPL_DIRTY_TEMP_COEFFICIENT(max_numer);
3717 PPL_DIRTY_TEMP_COEFFICIENT(max_denom);
3718 bool max_included;
3719 bool bound_above = maximize(denominator*var, max_numer, max_denom, max_included);
3720 PPL_DIRTY_TEMP_COEFFICIENT(min_numer);
3721 PPL_DIRTY_TEMP_COEFFICIENT(min_denom);
3722 bool min_included;
3723 bool bound_below = minimize(denominator*var, min_numer, min_denom, min_included);
3724 // Use the correct relation symbol
3725 const Relation_Symbol corrected_relsym
3726 = (denominator > 0) ? relsym : reversed_relsym;
3727 // Revise the expression to take into account the denominator of the
3728 // maximum/minimum value for `var'.
3729 Linear_Expression revised_expr;
3730 PPL_DIRTY_TEMP_COEFFICIENT(d);
3731 if (corrected_relsym == LESS_THAN || corrected_relsym == LESS_OR_EQUAL) {
3732 if (bound_below) {
3733 revised_expr = expr;
3734 revised_expr.set_inhomogeneous_term(Coefficient_zero());
3735 revised_expr *= d;
3736 }
3737 }
3738 else {
3739 if (bound_above) {
3740 revised_expr = expr;
3741 revised_expr.set_inhomogeneous_term(Coefficient_zero());
3742 revised_expr *= max_denom;
3743 }
3744 }
3745
3746 switch (corrected_relsym) {
3747 case LESS_THAN:
3748 if (bound_below) {
3749 refine_with_constraint(min_numer < revised_expr);
3750 }
3751 break;
3752 case LESS_OR_EQUAL:
3753 if (bound_below) {
3754 (min_included)
3755 ? refine_with_constraint(min_numer <= revised_expr)
3756 : refine_with_constraint(min_numer < revised_expr);
3757 }
3758 break;
3759 case GREATER_OR_EQUAL:
3760 if (bound_above) {
3761 (max_included)
3762 ? refine_with_constraint(max_numer >= revised_expr)
3763 : refine_with_constraint(max_numer > revised_expr);
3764 }
3765 break;
3766 case GREATER_THAN:
3767 if (bound_above) {
3768 refine_with_constraint(max_numer > revised_expr);
3769 }
3770 break;
3771 default:
3772 // The EQUAL and NOT_EQUAL cases have been already dealt with.
3773 PPL_UNREACHABLE;
3774 break;
3775 }
3776 // If the shrunk box is empty, its preimage is empty too.
3777 if (is_empty()) {
3778 return;
3779 }
3780 ITV& seq_v = seq[var.id()];
3781 seq_v.assign(UNIVERSE);
3782 PPL_ASSERT(OK());
3783 }
3784
3785 template <typename ITV>
3786 void
3787 Box<ITV>
generalized_affine_image(const Linear_Expression & lhs,const Relation_Symbol relsym,const Linear_Expression & rhs)3788 ::generalized_affine_image(const Linear_Expression& lhs,
3789 const Relation_Symbol relsym,
3790 const Linear_Expression& rhs) {
3791 // Dimension-compatibility checks.
3792 // The dimension of `lhs' should not be greater than the dimension
3793 // of `*this'.
3794 dimension_type lhs_space_dim = lhs.space_dimension();
3795 const dimension_type space_dim = space_dimension();
3796 if (space_dim < lhs_space_dim) {
3797 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
3798 "e1", lhs);
3799 }
3800 // The dimension of `rhs' should not be greater than the dimension
3801 // of `*this'.
3802 const dimension_type rhs_space_dim = rhs.space_dimension();
3803 if (space_dim < rhs_space_dim) {
3804 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
3805 "e2", rhs);
3806 }
3807
3808 // The relation symbol cannot be a disequality.
3809 if (relsym == NOT_EQUAL) {
3810 throw_invalid_argument("generalized_affine_image(e1, r, e2)",
3811 "r is the disequality relation symbol");
3812 }
3813
3814 // Any image of an empty box is empty.
3815 if (marked_empty()) {
3816 return;
3817 }
3818
3819 // Compute the maximum and minimum value reached by the rhs on the box.
3820 PPL_DIRTY_TEMP_COEFFICIENT(max_numer);
3821 PPL_DIRTY_TEMP_COEFFICIENT(max_denom);
3822 bool max_included;
3823 bool max_rhs = maximize(rhs, max_numer, max_denom, max_included);
3824 PPL_DIRTY_TEMP_COEFFICIENT(min_numer);
3825 PPL_DIRTY_TEMP_COEFFICIENT(min_denom);
3826 bool min_included;
3827 bool min_rhs = minimize(rhs, min_numer, min_denom, min_included);
3828
3829 // Check whether there is 0, 1 or more than one variable in the lhs
3830 // and record the variable with the highest dimension; set the box
3831 // intervals to be unbounded for all other dimensions with non-zero
3832 // coefficients in the lhs.
3833 bool has_var = false;
3834 dimension_type has_var_id = lhs.last_nonzero();
3835
3836 if (has_var_id != 0) {
3837 has_var = true;
3838 --has_var_id;
3839 dimension_type other_var = lhs.first_nonzero(1, has_var_id + 1);
3840 --other_var;
3841 if (other_var != has_var_id) {
3842 // There is more than one dimension with non-zero coefficient, so
3843 // we cannot have any information about the dimensions in the lhs.
3844 ITV& seq_var = seq[has_var_id];
3845 seq_var.assign(UNIVERSE);
3846 // Since all but the highest dimension with non-zero coefficient
3847 // in the lhs have been set unbounded, it remains to set the
3848 // highest dimension in the lhs unbounded.
3849 ITV& seq_i = seq[other_var];
3850 seq_i.assign(UNIVERSE);
3851 PPL_ASSERT(OK());
3852 return;
3853 }
3854 }
3855
3856 if (has_var) {
3857 // There is exactly one dimension with non-zero coefficient.
3858 ITV& seq_var = seq[has_var_id];
3859
3860 // Compute the new bounds for this dimension defined by the rhs
3861 // expression.
3862 const Coefficient& inhomo = lhs.inhomogeneous_term();
3863 const Coefficient& coeff = lhs.coefficient(Variable(has_var_id));
3864 PPL_DIRTY_TEMP(mpq_class, q_max);
3865 PPL_DIRTY_TEMP(mpq_class, q_min);
3866 if (max_rhs) {
3867 max_numer -= inhomo * max_denom;
3868 max_denom *= coeff;
3869 assign_r(q_max.get_num(), max_numer, ROUND_NOT_NEEDED);
3870 assign_r(q_max.get_den(), max_denom, ROUND_NOT_NEEDED);
3871 q_max.canonicalize();
3872 }
3873 if (min_rhs) {
3874 min_numer -= inhomo * min_denom;
3875 min_denom *= coeff;
3876 assign_r(q_min.get_num(), min_numer, ROUND_NOT_NEEDED);
3877 assign_r(q_min.get_den(), min_denom, ROUND_NOT_NEEDED);
3878 q_min.canonicalize();
3879 }
3880
3881 // The choice as to which bounds should be set depends on the sign
3882 // of the coefficient of the dimension `has_var_id' in the lhs.
3883 if (coeff > 0) {
3884 // The coefficient of the dimension in the lhs is positive.
3885 switch (relsym) {
3886 case LESS_OR_EQUAL:
3887 if (max_rhs) {
3888 Relation_Symbol rel = max_included ? LESS_OR_EQUAL : LESS_THAN;
3889 seq_var.build(i_constraint(rel, q_max));
3890 }
3891 else {
3892 seq_var.assign(UNIVERSE);
3893 }
3894 break;
3895 case LESS_THAN:
3896 if (max_rhs) {
3897 seq_var.build(i_constraint(LESS_THAN, q_max));
3898 }
3899 else {
3900 seq_var.assign(UNIVERSE);
3901 }
3902 break;
3903 case EQUAL:
3904 {
3905 I_Constraint<mpq_class> l;
3906 I_Constraint<mpq_class> u;
3907 if (max_rhs) {
3908 u.set(max_included ? LESS_OR_EQUAL : LESS_THAN, q_max);
3909 }
3910 if (min_rhs) {
3911 l.set(min_included ? GREATER_OR_EQUAL : GREATER_THAN, q_min);
3912 }
3913 seq_var.build(l, u);
3914 break;
3915 }
3916 case GREATER_OR_EQUAL:
3917 if (min_rhs) {
3918 Relation_Symbol rel = min_included ? GREATER_OR_EQUAL : GREATER_THAN;
3919 seq_var.build(i_constraint(rel, q_min));
3920 }
3921 else {
3922 seq_var.assign(UNIVERSE);
3923 }
3924 break;
3925 case GREATER_THAN:
3926 if (min_rhs) {
3927 seq_var.build(i_constraint(GREATER_THAN, q_min));
3928 }
3929 else {
3930 seq_var.assign(UNIVERSE);
3931 }
3932 break;
3933 default:
3934 // The NOT_EQUAL case has been already dealt with.
3935 PPL_UNREACHABLE;
3936 break;
3937 }
3938 }
3939 else {
3940 // The coefficient of the dimension in the lhs is negative.
3941 switch (relsym) {
3942 case GREATER_OR_EQUAL:
3943 if (min_rhs) {
3944 Relation_Symbol rel = min_included ? LESS_OR_EQUAL : LESS_THAN;
3945 seq_var.build(i_constraint(rel, q_min));
3946 }
3947 else {
3948 seq_var.assign(UNIVERSE);
3949 }
3950 break;
3951 case GREATER_THAN:
3952 if (min_rhs) {
3953 seq_var.build(i_constraint(LESS_THAN, q_min));
3954 }
3955 else {
3956 seq_var.assign(UNIVERSE);
3957 }
3958 break;
3959 case EQUAL:
3960 {
3961 I_Constraint<mpq_class> l;
3962 I_Constraint<mpq_class> u;
3963 if (max_rhs) {
3964 l.set(max_included ? GREATER_OR_EQUAL : GREATER_THAN, q_max);
3965 }
3966 if (min_rhs) {
3967 u.set(min_included ? LESS_OR_EQUAL : LESS_THAN, q_min);
3968 }
3969 seq_var.build(l, u);
3970 break;
3971 }
3972 case LESS_OR_EQUAL:
3973 if (max_rhs) {
3974 Relation_Symbol rel = max_included ? GREATER_OR_EQUAL : GREATER_THAN;
3975 seq_var.build(i_constraint(rel, q_max));
3976 }
3977 else {
3978 seq_var.assign(UNIVERSE);
3979 }
3980 break;
3981 case LESS_THAN:
3982 if (max_rhs) {
3983 seq_var.build(i_constraint(GREATER_THAN, q_max));
3984 }
3985 else {
3986 seq_var.assign(UNIVERSE);
3987 }
3988 break;
3989 default:
3990 // The NOT_EQUAL case has been already dealt with.
3991 PPL_UNREACHABLE;
3992 break;
3993 }
3994 }
3995 }
3996
3997 else {
3998 // The lhs is a constant value, so we just need to add the
3999 // appropriate constraint.
4000 const Coefficient& inhomo = lhs.inhomogeneous_term();
4001 switch (relsym) {
4002 case LESS_THAN:
4003 refine_with_constraint(inhomo < rhs);
4004 break;
4005 case LESS_OR_EQUAL:
4006 refine_with_constraint(inhomo <= rhs);
4007 break;
4008 case EQUAL:
4009 refine_with_constraint(inhomo == rhs);
4010 break;
4011 case GREATER_OR_EQUAL:
4012 refine_with_constraint(inhomo >= rhs);
4013 break;
4014 case GREATER_THAN:
4015 refine_with_constraint(inhomo > rhs);
4016 break;
4017 default:
4018 // The NOT_EQUAL case has been already dealt with.
4019 PPL_UNREACHABLE;
4020 break;
4021 }
4022 }
4023 PPL_ASSERT(OK());
4024 }
4025
4026 template <typename ITV>
4027 void
generalized_affine_preimage(const Linear_Expression & lhs,const Relation_Symbol relsym,const Linear_Expression & rhs)4028 Box<ITV>::generalized_affine_preimage(const Linear_Expression& lhs,
4029 const Relation_Symbol relsym,
4030 const Linear_Expression& rhs) {
4031 // Dimension-compatibility checks.
4032 // The dimension of `lhs' should not be greater than the dimension
4033 // of `*this'.
4034 dimension_type lhs_space_dim = lhs.space_dimension();
4035 const dimension_type space_dim = space_dimension();
4036 if (space_dim < lhs_space_dim) {
4037 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
4038 "e1", lhs);
4039 }
4040 // The dimension of `rhs' should not be greater than the dimension
4041 // of `*this'.
4042 const dimension_type rhs_space_dim = rhs.space_dimension();
4043 if (space_dim < rhs_space_dim) {
4044 throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
4045 "e2", rhs);
4046 }
4047
4048 // The relation symbol cannot be a disequality.
4049 if (relsym == NOT_EQUAL) {
4050 throw_invalid_argument("generalized_affine_image(e1, r, e2)",
4051 "r is the disequality relation symbol");
4052 }
4053 // Any image of an empty box is empty.
4054 if (marked_empty()) {
4055 return;
4056 }
4057 // For any dimension occurring in the lhs, swap and change the sign
4058 // of this component for the rhs and lhs. Then use these in a call
4059 // to generalized_affine_image/3.
4060 Linear_Expression revised_lhs = lhs;
4061 Linear_Expression revised_rhs = rhs;
4062 for (Linear_Expression::const_iterator i = lhs.begin(),
4063 i_end = lhs.end(); i != i_end; ++i) {
4064 const Variable var = i.variable();
4065 PPL_DIRTY_TEMP_COEFFICIENT(tmp);
4066 tmp = *i;
4067 tmp += rhs.coefficient(var);
4068 sub_mul_assign(revised_rhs, tmp, var);
4069 sub_mul_assign(revised_lhs, tmp, var);
4070 }
4071 generalized_affine_image(revised_lhs, relsym, revised_rhs);
4072 PPL_ASSERT(OK());
4073 }
4074
4075 template <typename ITV>
4076 template <typename T, typename Iterator>
4077 typename Enable_If<Is_Same<T, Box<ITV> >::value
4078 && Is_Same_Or_Derived<Interval_Base, ITV>::value,
4079 void>::type
CC76_widening_assign(const T & y,Iterator first,Iterator last)4080 Box<ITV>::CC76_widening_assign(const T& y, Iterator first, Iterator last) {
4081 if (y.is_empty()) {
4082 return;
4083 }
4084 for (dimension_type i = seq.size(); i-- > 0; ) {
4085 seq[i].CC76_widening_assign(y.seq[i], first, last);
4086 }
4087
4088 PPL_ASSERT(OK());
4089 }
4090
4091 template <typename ITV>
4092 template <typename T>
4093 typename Enable_If<Is_Same<T, Box<ITV> >::value
4094 && Is_Same_Or_Derived<Interval_Base, ITV>::value,
4095 void>::type
CC76_widening_assign(const T & y,unsigned * tp)4096 Box<ITV>::CC76_widening_assign(const T& y, unsigned* tp) {
4097 static typename ITV::boundary_type stop_points[] = {
4098 typename ITV::boundary_type(-2),
4099 typename ITV::boundary_type(-1),
4100 typename ITV::boundary_type(0),
4101 typename ITV::boundary_type(1),
4102 typename ITV::boundary_type(2)
4103 };
4104
4105 Box& x = *this;
4106 // If there are tokens available, work on a temporary copy.
4107 if (tp != 0 && *tp > 0) {
4108 Box<ITV> x_tmp(x);
4109 x_tmp.CC76_widening_assign(y, 0);
4110 // If the widening was not precise, use one of the available tokens.
4111 if (!x.contains(x_tmp)) {
4112 --(*tp);
4113 }
4114 return;
4115 }
4116 x.CC76_widening_assign(y,
4117 stop_points,
4118 stop_points
4119 + sizeof(stop_points)/sizeof(stop_points[0]));
4120 }
4121
4122 template <typename ITV>
4123 void
get_limiting_box(const Constraint_System & cs,Box & limiting_box) const4124 Box<ITV>::get_limiting_box(const Constraint_System& cs,
4125 Box& limiting_box) const {
4126 // Private method: the caller has to ensure the following.
4127 PPL_ASSERT(cs.space_dimension() <= space_dimension());
4128
4129 for (Constraint_System::const_iterator cs_i = cs.begin(),
4130 cs_end = cs.end(); cs_i != cs_end; ++cs_i) {
4131 const Constraint& c = *cs_i;
4132 dimension_type c_num_vars = 0;
4133 dimension_type c_only_var = 0;
4134 // Constraints that are not interval constraints are ignored.
4135 if (!Box_Helpers::extract_interval_constraint(c, c_num_vars, c_only_var)) {
4136 continue;
4137 }
4138 // Trivial constraints are ignored.
4139 if (c_num_vars != 0) {
4140 // c is a non-trivial interval constraint.
4141 // add interval constraint to limiting box
4142 const Coefficient& n = c.inhomogeneous_term();
4143 const Coefficient& d = c.coefficient(Variable(c_only_var));
4144 if (interval_relation(seq[c_only_var], c.type(), n, d)
4145 == Poly_Con_Relation::is_included()) {
4146 limiting_box.add_interval_constraint_no_check(c_only_var, c.type(),
4147 n, d);
4148 }
4149 }
4150 }
4151 }
4152
4153 template <typename ITV>
4154 void
limited_CC76_extrapolation_assign(const Box & y,const Constraint_System & cs,unsigned * tp)4155 Box<ITV>::limited_CC76_extrapolation_assign(const Box& y,
4156 const Constraint_System& cs,
4157 unsigned* tp) {
4158 Box& x = *this;
4159 const dimension_type space_dim = x.space_dimension();
4160
4161 // Dimension-compatibility check.
4162 if (space_dim != y.space_dimension()) {
4163 throw_dimension_incompatible("limited_CC76_extrapolation_assign(y, cs)",
4164 y);
4165 }
4166 // `cs' must be dimension-compatible with the two boxes.
4167 const dimension_type cs_space_dim = cs.space_dimension();
4168 if (space_dim < cs_space_dim) {
4169 throw_constraint_incompatible("limited_CC76_extrapolation_assign(y, cs)");
4170 }
4171 // The limited CC76-extrapolation between two boxes in a
4172 // zero-dimensional space is also a zero-dimensional box
4173 if (space_dim == 0) {
4174 return;
4175 }
4176 // Assume `y' is contained in or equal to `*this'.
4177 PPL_EXPECT_HEAVY(copy_contains(*this, y));
4178
4179 // If `*this' is empty, since `*this' contains `y', `y' is empty too.
4180 if (marked_empty()) {
4181 return;
4182 }
4183 // If `y' is empty, we return.
4184 if (y.marked_empty()) {
4185 return;
4186 }
4187 // Build a limiting box using all the constraints in cs
4188 // that are satisfied by *this.
4189 Box limiting_box(space_dim, UNIVERSE);
4190 get_limiting_box(cs, limiting_box);
4191
4192 x.CC76_widening_assign(y, tp);
4193
4194 // Intersect the widened box with the limiting box.
4195 intersection_assign(limiting_box);
4196 }
4197
4198 template <typename ITV>
4199 template <typename T>
4200 typename Enable_If<Is_Same<T, Box<ITV> >::value
4201 && Is_Same_Or_Derived<Interval_Base, ITV>::value,
4202 void>::type
CC76_narrowing_assign(const T & y)4203 Box<ITV>::CC76_narrowing_assign(const T& y) {
4204 const dimension_type space_dim = space_dimension();
4205
4206 // Dimension-compatibility check.
4207 if (space_dim != y.space_dimension()) {
4208 throw_dimension_incompatible("CC76_narrowing_assign(y)", y);
4209 }
4210
4211 // Assume `*this' is contained in or equal to `y'.
4212 PPL_EXPECT_HEAVY(copy_contains(y, *this));
4213
4214 // If both boxes are zero-dimensional,
4215 // since `y' contains `*this', we simply return `*this'.
4216 if (space_dim == 0) {
4217 return;
4218 }
4219 // If `y' is empty, since `y' contains `this', `*this' is empty too.
4220 if (y.is_empty()) {
4221 return;
4222 }
4223 // If `*this' is empty, we return.
4224 if (is_empty()) {
4225 return;
4226 }
4227 // Replace each constraint in `*this' by the corresponding constraint
4228 // in `y' if the corresponding inhomogeneous terms are both finite.
4229 for (dimension_type i = space_dim; i-- > 0; ) {
4230 ITV& x_i = seq[i];
4231 const ITV& y_i = y.seq[i];
4232 if (!x_i.lower_is_boundary_infinity()
4233 && !y_i.lower_is_boundary_infinity()
4234 && x_i.lower() != y_i.lower()) {
4235 x_i.lower() = y_i.lower();
4236 }
4237 if (!x_i.upper_is_boundary_infinity()
4238 && !y_i.upper_is_boundary_infinity()
4239 && x_i.upper() != y_i.upper()) {
4240 x_i.upper() = y_i.upper();
4241 }
4242 }
4243 PPL_ASSERT(OK());
4244 }
4245
4246 template <typename ITV>
4247 Constraint_System
constraints() const4248 Box<ITV>::constraints() const {
4249 const dimension_type space_dim = space_dimension();
4250 Constraint_System cs;
4251 cs.set_space_dimension(space_dim);
4252
4253 if (space_dim == 0) {
4254 if (marked_empty()) {
4255 cs = Constraint_System::zero_dim_empty();
4256 }
4257 return cs;
4258 }
4259
4260 if (marked_empty()) {
4261 cs.insert(Constraint::zero_dim_false());
4262 return cs;
4263 }
4264
4265 for (dimension_type k = 0; k < space_dim; ++k) {
4266 const Variable v_k = Variable(k);
4267 PPL_DIRTY_TEMP_COEFFICIENT(n);
4268 PPL_DIRTY_TEMP_COEFFICIENT(d);
4269 bool closed = false;
4270 if (has_lower_bound(v_k, n, d, closed)) {
4271 if (closed) {
4272 cs.insert(d * v_k >= n);
4273 }
4274 else {
4275 cs.insert(d * v_k > n);
4276 }
4277 }
4278 if (has_upper_bound(v_k, n, d, closed)) {
4279 if (closed) {
4280 cs.insert(d * v_k <= n);
4281 }
4282 else {
4283 cs.insert(d * v_k < n);
4284 }
4285 }
4286 }
4287 return cs;
4288 }
4289
4290 template <typename ITV>
4291 Constraint_System
minimized_constraints() const4292 Box<ITV>::minimized_constraints() const {
4293 const dimension_type space_dim = space_dimension();
4294 Constraint_System cs;
4295 cs.set_space_dimension(space_dim);
4296
4297 if (space_dim == 0) {
4298 if (marked_empty()) {
4299 cs = Constraint_System::zero_dim_empty();
4300 }
4301 return cs;
4302 }
4303
4304 // Make sure emptiness is detected.
4305 if (is_empty()) {
4306 cs.insert(Constraint::zero_dim_false());
4307 return cs;
4308 }
4309
4310 for (dimension_type k = 0; k < space_dim; ++k) {
4311 const Variable v_k = Variable(k);
4312 PPL_DIRTY_TEMP_COEFFICIENT(n);
4313 PPL_DIRTY_TEMP_COEFFICIENT(d);
4314 bool closed = false;
4315 if (has_lower_bound(v_k, n, d, closed)) {
4316 if (closed) {
4317 // Make sure equality constraints are detected.
4318 if (seq[k].is_singleton()) {
4319 cs.insert(d * v_k == n);
4320 continue;
4321 }
4322 else {
4323 cs.insert(d * v_k >= n);
4324 }
4325 }
4326 else {
4327 cs.insert(d * v_k > n);
4328 }
4329 }
4330 if (has_upper_bound(v_k, n, d, closed)) {
4331 if (closed) {
4332 cs.insert(d * v_k <= n);
4333 }
4334 else {
4335 cs.insert(d * v_k < n);
4336 }
4337 }
4338 }
4339 return cs;
4340 }
4341
4342 template <typename ITV>
4343 Congruence_System
congruences() const4344 Box<ITV>::congruences() const {
4345 const dimension_type space_dim = space_dimension();
4346 Congruence_System cgs(space_dim);
4347
4348 if (space_dim == 0) {
4349 if (marked_empty()) {
4350 cgs = Congruence_System::zero_dim_empty();
4351 }
4352 return cgs;
4353 }
4354
4355 // Make sure emptiness is detected.
4356 if (is_empty()) {
4357 cgs.insert(Congruence::zero_dim_false());
4358 return cgs;
4359 }
4360
4361 for (dimension_type k = 0; k < space_dim; ++k) {
4362 const Variable v_k = Variable(k);
4363 PPL_DIRTY_TEMP_COEFFICIENT(n);
4364 PPL_DIRTY_TEMP_COEFFICIENT(d);
4365 bool closed = false;
4366 if (has_lower_bound(v_k, n, d, closed) && closed) {
4367 // Make sure equality congruences are detected.
4368 if (seq[k].is_singleton()) {
4369 cgs.insert((d * v_k %= n) / 0);
4370 }
4371 }
4372 }
4373 return cgs;
4374 }
4375
4376 template <typename ITV>
4377 memory_size_type
external_memory_in_bytes() const4378 Box<ITV>::external_memory_in_bytes() const {
4379 memory_size_type n = seq.capacity() * sizeof(ITV);
4380 for (dimension_type k = seq.size(); k-- > 0; ) {
4381 n += seq[k].external_memory_in_bytes();
4382 }
4383 return n;
4384 }
4385
4386 /*! \relates Parma_Polyhedra_Library::Box */
4387 template <typename ITV>
4388 std::ostream&
operator <<(std::ostream & s,const Box<ITV> & box)4389 IO_Operators::operator<<(std::ostream& s, const Box<ITV>& box) {
4390 if (box.is_empty()) {
4391 s << "false";
4392 }
4393 else if (box.is_universe()) {
4394 s << "true";
4395 }
4396 else {
4397 for (dimension_type k = 0,
4398 space_dim = box.space_dimension(); k < space_dim; ) {
4399 s << Variable(k) << " in " << box[k];
4400 ++k;
4401 if (k < space_dim) {
4402 s << ", ";
4403 }
4404 else {
4405 break;
4406 }
4407 }
4408 }
4409 return s;
4410 }
4411
4412 template <typename ITV>
4413 void
ascii_dump(std::ostream & s) const4414 Box<ITV>::ascii_dump(std::ostream& s) const {
4415 const char separator = ' ';
4416 status.ascii_dump(s);
4417 const dimension_type space_dim = space_dimension();
4418 s << "space_dim" << separator << space_dim;
4419 s << "\n";
4420 for (dimension_type i = 0; i < space_dim; ++i) {
4421 seq[i].ascii_dump(s);
4422 }
4423 }
4424
PPL_OUTPUT_TEMPLATE_DEFINITIONS(ITV,Box<ITV>)4425 PPL_OUTPUT_TEMPLATE_DEFINITIONS(ITV, Box<ITV>)
4426
4427 template <typename ITV>
4428 bool
4429 Box<ITV>::ascii_load(std::istream& s) {
4430 if (!status.ascii_load(s)) {
4431 return false;
4432 }
4433 std::string str;
4434 dimension_type space_dim;
4435 if (!(s >> str) || str != "space_dim") {
4436 return false;
4437 }
4438 if (!(s >> space_dim)) {
4439 return false;
4440 }
4441 seq.clear();
4442 ITV seq_i;
4443 for (dimension_type i = 0; i < space_dim; ++i) {
4444 if (seq_i.ascii_load(s)) {
4445 seq.push_back(seq_i);
4446 }
4447 else {
4448 return false;
4449 }
4450 }
4451
4452 // Check invariants.
4453 PPL_ASSERT(OK());
4454 return true;
4455 }
4456
4457 template <typename ITV>
4458 void
throw_dimension_incompatible(const char * method,const Box & y) const4459 Box<ITV>::throw_dimension_incompatible(const char* method,
4460 const Box& y) const {
4461 std::ostringstream s;
4462 s << "PPL::Box::" << method << ":" << std::endl
4463 << "this->space_dimension() == " << this->space_dimension()
4464 << ", y->space_dimension() == " << y.space_dimension() << ".";
4465 throw std::invalid_argument(s.str());
4466 }
4467
4468 template <typename ITV>
4469 void
4470 Box<ITV>
throw_dimension_incompatible(const char * method,dimension_type required_dim) const4471 ::throw_dimension_incompatible(const char* method,
4472 dimension_type required_dim) const {
4473 std::ostringstream s;
4474 s << "PPL::Box::" << method << ":" << std::endl
4475 << "this->space_dimension() == " << space_dimension()
4476 << ", required dimension == " << required_dim << ".";
4477 throw std::invalid_argument(s.str());
4478 }
4479
4480 template <typename ITV>
4481 void
throw_dimension_incompatible(const char * method,const Constraint & c) const4482 Box<ITV>::throw_dimension_incompatible(const char* method,
4483 const Constraint& c) const {
4484 std::ostringstream s;
4485 s << "PPL::Box::" << method << ":" << std::endl
4486 << "this->space_dimension() == " << space_dimension()
4487 << ", c->space_dimension == " << c.space_dimension() << ".";
4488 throw std::invalid_argument(s.str());
4489 }
4490
4491 template <typename ITV>
4492 void
throw_dimension_incompatible(const char * method,const Congruence & cg) const4493 Box<ITV>::throw_dimension_incompatible(const char* method,
4494 const Congruence& cg) const {
4495 std::ostringstream s;
4496 s << "PPL::Box::" << method << ":" << std::endl
4497 << "this->space_dimension() == " << space_dimension()
4498 << ", cg->space_dimension == " << cg.space_dimension() << ".";
4499 throw std::invalid_argument(s.str());
4500 }
4501
4502 template <typename ITV>
4503 void
throw_dimension_incompatible(const char * method,const Constraint_System & cs) const4504 Box<ITV>::throw_dimension_incompatible(const char* method,
4505 const Constraint_System& cs) const {
4506 std::ostringstream s;
4507 s << "PPL::Box::" << method << ":" << std::endl
4508 << "this->space_dimension() == " << space_dimension()
4509 << ", cs->space_dimension == " << cs.space_dimension() << ".";
4510 throw std::invalid_argument(s.str());
4511 }
4512
4513 template <typename ITV>
4514 void
throw_dimension_incompatible(const char * method,const Congruence_System & cgs) const4515 Box<ITV>::throw_dimension_incompatible(const char* method,
4516 const Congruence_System& cgs) const {
4517 std::ostringstream s;
4518 s << "PPL::Box::" << method << ":" << std::endl
4519 << "this->space_dimension() == " << space_dimension()
4520 << ", cgs->space_dimension == " << cgs.space_dimension() << ".";
4521 throw std::invalid_argument(s.str());
4522 }
4523
4524 template <typename ITV>
4525 void
throw_dimension_incompatible(const char * method,const Generator & g) const4526 Box<ITV>::throw_dimension_incompatible(const char* method,
4527 const Generator& g) const {
4528 std::ostringstream s;
4529 s << "PPL::Box::" << method << ":" << std::endl
4530 << "this->space_dimension() == " << space_dimension()
4531 << ", g->space_dimension == " << g.space_dimension() << ".";
4532 throw std::invalid_argument(s.str());
4533 }
4534
4535 template <typename ITV>
4536 void
throw_constraint_incompatible(const char * method)4537 Box<ITV>::throw_constraint_incompatible(const char* method) {
4538 std::ostringstream s;
4539 s << "PPL::Box::" << method << ":" << std::endl
4540 << "the constraint is incompatible.";
4541 throw std::invalid_argument(s.str());
4542 }
4543
4544 template <typename ITV>
4545 void
throw_expression_too_complex(const char * method,const Linear_Expression & le)4546 Box<ITV>::throw_expression_too_complex(const char* method,
4547 const Linear_Expression& le) {
4548 using namespace IO_Operators;
4549 std::ostringstream s;
4550 s << "PPL::Box::" << method << ":" << std::endl
4551 << le << " is too complex.";
4552 throw std::invalid_argument(s.str());
4553 }
4554
4555 template <typename ITV>
4556 void
throw_dimension_incompatible(const char * method,const char * le_name,const Linear_Expression & le) const4557 Box<ITV>::throw_dimension_incompatible(const char* method,
4558 const char* le_name,
4559 const Linear_Expression& le) const {
4560 std::ostringstream s;
4561 s << "PPL::Box::" << method << ":" << std::endl
4562 << "this->space_dimension() == " << space_dimension()
4563 << ", " << le_name << "->space_dimension() == "
4564 << le.space_dimension() << ".";
4565 throw std::invalid_argument(s.str());
4566 }
4567
4568 template <typename ITV>
4569 template <typename C>
4570 void
throw_dimension_incompatible(const char * method,const char * lf_name,const Linear_Form<C> & lf) const4571 Box<ITV>::throw_dimension_incompatible(const char* method,
4572 const char* lf_name,
4573 const Linear_Form<C>& lf) const {
4574 std::ostringstream s;
4575 s << "PPL::Box::" << method << ":\n"
4576 << "this->space_dimension() == " << space_dimension()
4577 << ", " << lf_name << "->space_dimension() == "
4578 << lf.space_dimension() << ".";
4579 throw std::invalid_argument(s.str());
4580 }
4581
4582 template <typename ITV>
4583 void
throw_invalid_argument(const char * method,const char * reason)4584 Box<ITV>::throw_invalid_argument(const char* method, const char* reason) {
4585 std::ostringstream s;
4586 s << "PPL::Box::" << method << ":" << std::endl
4587 << reason;
4588 throw std::invalid_argument(s.str());
4589 }
4590
4591 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
4592 /*! \relates Box */
4593 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
4594 template <typename Specialization,
4595 typename Temp, typename To, typename ITV>
4596 bool
l_m_distance_assign(Checked_Number<To,Extended_Number_Policy> & r,const Box<ITV> & x,const Box<ITV> & y,const Rounding_Dir dir,Temp & tmp0,Temp & tmp1,Temp & tmp2)4597 l_m_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
4598 const Box<ITV>& x, const Box<ITV>& y,
4599 const Rounding_Dir dir,
4600 Temp& tmp0, Temp& tmp1, Temp& tmp2) {
4601 const dimension_type x_space_dim = x.space_dimension();
4602 // Dimension-compatibility check.
4603 if (x_space_dim != y.space_dimension()) {
4604 return false;
4605 }
4606 // Zero-dim boxes are equal if and only if they are both empty or universe.
4607 if (x_space_dim == 0) {
4608 if (x.marked_empty() == y.marked_empty()) {
4609 assign_r(r, 0, ROUND_NOT_NEEDED);
4610 }
4611 else {
4612 assign_r(r, PLUS_INFINITY, ROUND_NOT_NEEDED);
4613 }
4614 return true;
4615 }
4616
4617 // The distance computation requires a check for emptiness.
4618 (void) x.is_empty();
4619 (void) y.is_empty();
4620 // If one of two boxes is empty, then they are equal if and only if
4621 // the other box is empty too.
4622 if (x.marked_empty() || y.marked_empty()) {
4623 if (x.marked_empty() == y.marked_empty()) {
4624 assign_r(r, 0, ROUND_NOT_NEEDED);
4625 return true;
4626 }
4627 else {
4628 goto pinf;
4629 }
4630 }
4631
4632 assign_r(tmp0, 0, ROUND_NOT_NEEDED);
4633 for (dimension_type i = x_space_dim; i-- > 0; ) {
4634 const ITV& x_i = x.seq[i];
4635 const ITV& y_i = y.seq[i];
4636 // Dealing with the lower bounds.
4637 if (x_i.lower_is_boundary_infinity()) {
4638 if (!y_i.lower_is_boundary_infinity()) {
4639 goto pinf;
4640 }
4641 }
4642 else if (y_i.lower_is_boundary_infinity()) {
4643 goto pinf;
4644 }
4645 else {
4646 const Temp* tmp1p;
4647 const Temp* tmp2p;
4648 if (x_i.lower() > y_i.lower()) {
4649 maybe_assign(tmp1p, tmp1, x_i.lower(), dir);
4650 maybe_assign(tmp2p, tmp2, y_i.lower(), inverse(dir));
4651 }
4652 else {
4653 maybe_assign(tmp1p, tmp1, y_i.lower(), dir);
4654 maybe_assign(tmp2p, tmp2, x_i.lower(), inverse(dir));
4655 }
4656 sub_assign_r(tmp1, *tmp1p, *tmp2p, dir);
4657 PPL_ASSERT(sgn(tmp1) >= 0);
4658 Specialization::combine(tmp0, tmp1, dir);
4659 }
4660 // Dealing with the lower bounds.
4661 if (x_i.upper_is_boundary_infinity()) {
4662 if (y_i.upper_is_boundary_infinity()) {
4663 continue;
4664 }
4665 else {
4666 goto pinf;
4667 }
4668 }
4669 else if (y_i.upper_is_boundary_infinity()) {
4670 goto pinf;
4671 }
4672 else {
4673 const Temp* tmp1p;
4674 const Temp* tmp2p;
4675 if (x_i.upper() > y_i.upper()) {
4676 maybe_assign(tmp1p, tmp1, x_i.upper(), dir);
4677 maybe_assign(tmp2p, tmp2, y_i.upper(), inverse(dir));
4678 }
4679 else {
4680 maybe_assign(tmp1p, tmp1, y_i.upper(), dir);
4681 maybe_assign(tmp2p, tmp2, x_i.upper(), inverse(dir));
4682 }
4683 sub_assign_r(tmp1, *tmp1p, *tmp2p, dir);
4684 PPL_ASSERT(sgn(tmp1) >= 0);
4685 Specialization::combine(tmp0, tmp1, dir);
4686 }
4687 }
4688 Specialization::finalize(tmp0, dir);
4689 assign_r(r, tmp0, dir);
4690 return true;
4691
4692 pinf:
4693 assign_r(r, PLUS_INFINITY, ROUND_NOT_NEEDED);
4694 return true;
4695 }
4696
4697 } // namespace Parma_Polyhedra_Library
4698
4699 #endif // !defined(PPL_Box_templates_hh)
4700