1 /* Utilities for termination analysis: non-inline, non-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 #include "ppl-config.h"
25 #include "termination_defs.hh"
26 #include "NNC_Polyhedron_defs.hh"
27
28 namespace Parma_Polyhedra_Library {
29
30 namespace Implementation {
31
32 namespace Termination {
33
34 void
assign_all_inequalities_approximation(const Constraint_System & cs_in,Constraint_System & cs_out)35 assign_all_inequalities_approximation(const Constraint_System& cs_in,
36 Constraint_System& cs_out) {
37 if (cs_in.has_strict_inequalities() || cs_in.has_equalities()) {
38 // Here we have some strict inequality and/or equality constraints:
39 // translate them into non-strict inequality constraints.
40 for (Constraint_System::const_iterator i = cs_in.begin(),
41 i_end = cs_in.end(); i != i_end; ++i) {
42 const Constraint& c = *i;
43 if (c.is_equality()) {
44 // Insert the two corresponding opposing inequalities.
45 const Linear_Expression expr(c.expression());
46 cs_out.insert(expr >= 0);
47 cs_out.insert(expr <= 0);
48 }
49 else if (c.is_strict_inequality()) {
50 // Insert the non-strict approximation.
51 const Linear_Expression expr(c.expression());
52 cs_out.insert(expr >= 0);
53 }
54 else {
55 // Insert as is.
56 cs_out.insert(c);
57 }
58 }
59 }
60 else {
61 // No strict inequality and no equality constraints.
62 cs_out = cs_in;
63 }
64 }
65
66 template <>
67 void
assign_all_inequalities_approximation(const C_Polyhedron & ph,Constraint_System & cs)68 assign_all_inequalities_approximation(const C_Polyhedron& ph,
69 Constraint_System& cs) {
70 const Constraint_System& ph_cs = ph.minimized_constraints();
71 if (ph_cs.has_equalities()) {
72 // Translate equalities into inequalities.
73 for (Constraint_System::const_iterator i = ph_cs.begin(),
74 i_end = ph_cs.end(); i != i_end; ++i) {
75 const Constraint& c = *i;
76 if (c.is_equality()) {
77 // Insert the two corresponding opposing inequalities.
78 const Linear_Expression expr(c.expression());
79 cs.insert(expr >= 0);
80 cs.insert(expr <= 0);
81 }
82 else {
83 // Insert as is.
84 cs.insert(c);
85 }
86 }
87 }
88 else {
89 // No equality constraints (and no strict inequalities).
90 cs = ph_cs;
91 }
92 }
93
94 /*! \brief
95 Fill the constraint system(s) for the application of the
96 Mesnard and Serebrenik improved termination tests.
97
98 \param cs
99 The input constraint system, where variables indices are allocated
100 as follows:
101 - \f$ x'_1, \ldots, x'_n \f$ go onto space dimensions
102 \f$ 0, \ldots, n-1 \f$,
103 - \f$ x_1, \ldots, x_n \f$ go onto space dimensions
104 \f$ n, \ldots, 2n-1 \f$.
105 .
106 The system does not contain any equality.
107
108 \param cs_out1
109 The first output constraint system.
110
111 \param cs_out2
112 The second output constraint system, if any: it may be an alias
113 for \p cs_out1.
114
115 The allocation of variable indices in the output constraint
116 systems \p cs_out1 and \p cs_out2 is as follows,
117 where \f$m\f$ is the number of constraints in \p cs:
118 - \f$ \mu_1, \ldots, \mu_n \f$ go onto space dimensions
119 \f$ 0, \ldots, n-1 \f$;
120 - \f$ \mu_0\f$ goes onto space dimension \f$ n \f$;
121 - \f$ y_1, \ldots, y_m \f$ go onto space dimensions
122 \f$ n+1, \ldots, n+m \f$;
123 .
124 if we use the same constraint system, that is
125 <CODE>&cs_out1 == &cs_out2</CODE>, then
126 - \f$ z_1, ..., z_m, z_{m+1}, z_{m+2} \f$ go onto space dimensions
127 \f$ n+m+1, ..., n+2*m+2 \f$;
128 .
129 otherwise
130 - \f$ z_1, ..., z_m, z_{m+1}, z_{m+2} \f$ go onto space dimensions
131 \f$ n+1, ..., n+m+2 \f$.
132 */
133 void
fill_constraint_systems_MS(const Constraint_System & cs,Constraint_System & cs_out1,Constraint_System & cs_out2)134 fill_constraint_systems_MS(const Constraint_System& cs,
135 Constraint_System& cs_out1,
136 Constraint_System& cs_out2) {
137 PPL_ASSERT(cs.space_dimension() % 2 == 0);
138 const dimension_type n = cs.space_dimension() / 2;
139 const dimension_type m = num_constraints(cs);
140
141 #if PRINT_DEBUG_INFO
142 Variable::output_function_type* p_default_output_function
143 = Variable::get_output_function();
144 Variable::set_output_function(output_function_MS);
145
146 output_function_MS_n = n;
147 output_function_MS_m = m;
148
149 std::cout << "*** cs ***" << std::endl;
150 output_function_MS_which = 0;
151 using namespace IO_Operators;
152 std::cout << cs << std::endl;
153 #endif
154
155 const dimension_type y_begin = n+1;
156 const dimension_type z_begin = y_begin + ((&cs_out1 == &cs_out2) ? m : 0);
157
158 // Make sure linear expressions have the correct space dimension.
159 Linear_Expression y_le;
160 y_le.set_space_dimension(y_begin + m);
161 Linear_Expression z_le;
162 z_le.set_space_dimension(z_begin + m + 2);
163 std::vector<Linear_Expression> y_les(2*n, y_le);
164 std::vector<Linear_Expression> z_les(2*n + 1, z_le);
165
166 dimension_type y = y_begin;
167 dimension_type z = z_begin;
168 for (Constraint_System::const_iterator i = cs.begin(),
169 cs_end = cs.end(); i != cs_end; ++i) {
170 const Variable v_y(y);
171 const Variable v_z(z);
172 ++y;
173 ++z;
174 cs_out1.insert(v_y >= 0);
175 cs_out2.insert(v_z >= 0);
176 const Constraint& c_i = *i;
177 Coefficient_traits::const_reference b_i = c_i.inhomogeneous_term();
178 if (b_i != 0) {
179 // Note that b_i is to the left ot the relation sign, hence here
180 // we have -= and not += just to avoid negating b_i.
181 sub_mul_assign(y_le, b_i, v_y);
182 sub_mul_assign(z_le, b_i, v_z);
183 }
184 for (Constraint::expr_type::const_iterator j = c_i.expression().begin(),
185 j_end = c_i.expression().end(); j != j_end; ++j) {
186 Coefficient_traits::const_reference a_i_j = *j;
187 const Variable v = j.variable();
188 add_mul_assign(y_les[v.id()], a_i_j, v_y);
189 add_mul_assign(z_les[v.id()], a_i_j, v_z);
190 }
191 }
192 z_le += Variable(z);
193 z_les[2*n] += Variable(z);
194 cs_out2.insert(Variable(z) >= 0);
195 ++z;
196 z_le -= Variable(z);
197 z_les[2*n] -= Variable(z);
198 cs_out2.insert(Variable(z) >= 0);
199 cs_out1.insert(y_le >= 1);
200 cs_out2.insert(z_le >= 0);
201 dimension_type j = 2*n;
202 while (j-- > n) {
203 cs_out1.insert(y_les[j] == Variable(j-n));
204 cs_out2.insert(z_les[j] == Variable(j-n));
205 }
206 ++j;
207 while (j-- > 0) {
208 cs_out1.insert(y_les[j] == -Variable(j));
209 cs_out2.insert(z_les[j] == 0);
210 }
211 cs_out2.insert(z_les[2*n] == Variable(n));
212
213 #if PRINT_DEBUG_INFO
214 if (&cs_out1 == &cs_out2) {
215 std::cout << "*** cs_mip ***" << std::endl;
216 output_function_MS_which = 3;
217 using namespace IO_Operators;
218 std::cout << cs_mip << std::endl;
219 }
220 else {
221 std::cout << "*** cs_out1 ***" << std::endl;
222 output_function_MS_which = 1;
223 using namespace IO_Operators;
224 std::cout << cs_out1 << std::endl;
225
226 std::cout << "*** cs_out2 ***" << std::endl;
227 output_function_MS_which = 2;
228 using namespace IO_Operators;
229 std::cout << cs_out2 << std::endl;
230 }
231
232 Variable::set_output_function(p_default_output_function);
233 #endif
234 }
235
236 /*! \brief
237 Fill the constraint system(s) for the application of the
238 Podelski and Rybalchenko improved termination tests.
239
240 \param cs_before
241 The input constraint system describing the state <EM>before</EM>
242 the execution of the loop body, where variables indices are allocated
243 as follows:
244 - \f$ x_1, \ldots, x_n \f$ go onto space dimensions
245 \f$ 0, \ldots, n-1 \f$.
246 .
247 The system does not contain any equality.
248
249 \param cs_after
250 The input constraint system describing the state <EM>after</EM>
251 the execution of the loop body, where variables indices are allocated
252 as follows:
253 - \f$ x'_1, \ldots, x'_n \f$ go onto space dimensions
254 \f$ 0, \ldots, n-1 \f$,
255 - \f$ x_1, \ldots, x_n \f$ go onto space dimensions
256 \f$ n, \ldots, 2n-1 \f$.
257 .
258 The system does not contain any equality.
259
260 \param cs_out
261 The output constraint system, where variables indices are allocated
262 as follows:
263 - \f$ u_3 \f$ goes onto space dimensions \f$ 0, \ldots, s-1 \f$;
264 - \f$ u_2 \f$ goes onto space dimensions \f$ s, \ldots, s+r-1 \f$;
265 - \f$ u_1 \f$ goes onto space dimensions \f$ s+r, \ldots, s+2r-1 \f$.
266
267 The improved Podelski-Rybalchenko method described in the paper
268 is based on a loop encoding of the form
269 \f[
270 \begin{pmatrix}
271 A_B & \vect{0} \\
272 A_C & A'_C
273 \end{pmatrix}
274 \begin{pmatrix}
275 \vect{x} \\ \vect{x}'
276 \end{pmatrix}
277 \leq
278 \begin{pmatrix}
279 \vect{b}_B \\ \vect{b}_C
280 \end{pmatrix},
281 \f]
282 where
283 \f$ A_B \in \Qset^{r \times n} \f$,
284 \f$ A_C \in \Qset^{s \times n} \f$,
285 \f$ A'_C \in \Qset^{s \times n} \f$,
286 \f$ \vect{b}_B \in \Qset^r \f$,
287 \f$ \vect{b}_C \in \Qset^s \f$.
288 The corresponding system is:
289 \f[
290 \begin{aligned}
291 (\vect{v}_1-\vect{v}_2)^\transpose A_B
292 - \vect{v}_3^\transpose A_C
293 &= \vect{0}^\transpose, \\
294 \vect{v}_2^\transpose A_B
295 + \vect{v}_3^\transpose (A_C+A_C')
296 &= \vect{0}^\transpose, \\
297 \vect{v}_2 \vect{b}_B + \vect{v}_3 \vect{b}_C
298 &< 0,
299 \end{aligned}
300 \f]
301 where \f$ \vect{v}_1 \in \Qset_+^r \f$, \f$ \vect{v}_2 \in \Qset_+^r \f$,
302 \f$ \vect{v}_3 \in \Qset_+^s \f$.
303 The space of ranking functions is then spanned by
304 \f$ \vect{v}_3^\transpose A_C' \vect x \f$.
305
306 In contrast, our encoding is of the form
307 \f[
308 \begin{pmatrix}
309 \vect{0} & E_B \\
310 E'_C & E_C
311 \end{pmatrix}
312 \begin{pmatrix}
313 \vect{x}' \\ \vect{x}
314 \end{pmatrix}
315 +
316 \begin{pmatrix}
317 \vect{d}_B \\ \vect{d}_C
318 \end{pmatrix}
319 \geq
320 \vect{0},
321 \f]
322 where \f$ {E}_B = -{A}_B \f$, \f$ {E}_C = -{A}_C \f$,
323 \f$ {E}'_C = -{A}'_C \f$, \f$ \vect{d}_B = \vect{b}_B \f$
324 and \f$ \vect{d}_C = \vect{b}_C \f$.
325 The corresponding system is:
326 \f[
327 \begin{aligned}
328 (\vect{u}_1-\vect{u}_2)^\transpose E_B
329 - \vect{u}_3^\transpose E_C
330 &= \vect{0}^\transpose, \\
331 \vect{u}_2^\transpose E_B
332 + \vect{u}_3^\transpose (E_C+E_C')
333 &= \vect{0}^\transpose, \\
334 \vect{u}_2 \vect{d}_B + \vect{u}_3 \vect{d}_C
335 &> 0,
336 \end{aligned}
337 \f]
338 where \f$ \vect{u}_1 \in \Qset_-^r \f$, \f$ \vect{u}_2 \in \Qset_-^r \f$,
339 \f$ \vect{u}_3 \in \Qset_-^s \f$.
340 The space of ranking functions is then spanned by
341 \f$ \vect{u}_3^\transpose E_C' \vect x \f$.
342
343 \param le_out
344 The expression to be minimized in the context of \p cs_out:
345 a value of \f$ -1 \f$ or less entails termination.
346 */
347 void
fill_constraint_system_PR(const Constraint_System & cs_before,const Constraint_System & cs_after,Constraint_System & cs_out,Linear_Expression & le_out)348 fill_constraint_system_PR(const Constraint_System& cs_before,
349 const Constraint_System& cs_after,
350 Constraint_System& cs_out,
351 Linear_Expression& le_out) {
352 PPL_ASSERT(cs_after.space_dimension() % 2 == 0);
353 PPL_ASSERT(2*cs_before.space_dimension() == cs_after.space_dimension());
354 const dimension_type n = cs_before.space_dimension();
355 const dimension_type r = num_constraints(cs_before);
356 const dimension_type s = num_constraints(cs_after);
357 const dimension_type m = r + s;
358
359 // Make sure linear expressions are not reallocated multiple times.
360 if (m > 0) {
361 le_out.set_space_dimension(m + r);
362 }
363 std::vector<Linear_Expression> les_eq(2*n, le_out);
364
365 dimension_type row_index = 0;
366 for (Constraint_System::const_iterator i = cs_before.begin(),
367 cs_before_end = cs_before.end();
368 i != cs_before_end;
369 ++i, ++row_index) {
370 const Variable u1_i(m + row_index);
371 const Variable u2_i(s + row_index);
372 const Constraint::expr_type e_i = i->expression();
373 for (Constraint::expr_type::const_iterator
374 j = e_i.begin(), j_end = e_i.end(); j != j_end; ++j) {
375 Coefficient_traits::const_reference A_ij_B = *j;
376 const Variable v = j.variable();
377 // (u1 - u2) A_B, in the context of j-th constraint.
378 add_mul_assign(les_eq[v.id()], A_ij_B, u1_i);
379 sub_mul_assign(les_eq[v.id()], A_ij_B, u2_i);
380 // u2 A_B, in the context of (j+n)-th constraint.
381 add_mul_assign(les_eq[v.id() + n], A_ij_B, u2_i);
382 }
383 Coefficient_traits::const_reference b_B = e_i.inhomogeneous_term();
384 if (b_B != 0) {
385 // u2 b_B, in the context of the strict inequality constraint.
386 add_mul_assign(le_out, b_B, u2_i);
387 }
388 }
389
390 row_index = 0;
391 for (Constraint_System::const_iterator i = cs_after.begin(),
392 cs_after_end = cs_after.end();
393 i != cs_after_end;
394 ++i, ++row_index) {
395 const Variable u3_i(row_index);
396 const Constraint::expr_type e_i = i->expression();
397 for (Constraint::expr_type::const_iterator
398 j = e_i.lower_bound(Variable(n)),
399 j_end = e_i.end(); j != j_end; ++j) {
400 Coefficient_traits::const_reference A_ij_C = *j;
401 const Variable v = j.variable();
402 // - u3 A_C, in the context of the j-th constraint.
403 sub_mul_assign(les_eq[v.id() - n], A_ij_C, u3_i);
404 // u3 A_C, in the context of the (j+n)-th constraint.
405 add_mul_assign(les_eq[v.id()], A_ij_C, u3_i);
406 }
407 for (Constraint::expr_type::const_iterator j = e_i.begin(),
408 j_end = e_i.lower_bound(Variable(n)); j != j_end; ++j) {
409 Coefficient_traits::const_reference Ap_ij_C = *j;
410 // u3 Ap_C, in the context of the (j+n)-th constraint.
411 add_mul_assign(les_eq[j.variable().id() + n], Ap_ij_C, u3_i);
412 }
413 Coefficient_traits::const_reference b_C = e_i.inhomogeneous_term();
414 if (b_C != 0) {
415 // u3 b_C, in the context of the strict inequality constraint.
416 add_mul_assign(le_out, b_C, u3_i);
417 }
418 }
419
420 // Add the nonnegativity constraints for u_1, u_2 and u_3.
421 for (dimension_type i = s + 2*r; i-- > 0; ) {
422 cs_out.insert(Variable(i) >= 0);
423 }
424
425 for (dimension_type j = 2*n; j-- > 0; ) {
426 cs_out.insert(les_eq[j] == 0);
427 }
428 }
429
430 void
fill_constraint_system_PR_original(const Constraint_System & cs,Constraint_System & cs_out,Linear_Expression & le_out)431 fill_constraint_system_PR_original(const Constraint_System& cs,
432 Constraint_System& cs_out,
433 Linear_Expression& le_out) {
434 PPL_ASSERT(cs.space_dimension() % 2 == 0);
435 const dimension_type n = cs.space_dimension() / 2;
436 const dimension_type m = num_constraints(cs);
437
438 // Make sure linear expressions are not reallocated multiple times.
439 if (m > 0) {
440 le_out.set_space_dimension(2*m);
441 }
442 std::vector<Linear_Expression> les_eq(3*n, le_out);
443
444 dimension_type row_index = 0;
445 for (Constraint_System::const_iterator i = cs.begin(),
446 cs_end = cs.end(); i != cs_end; ++i, ++row_index) {
447 const Constraint::expr_type e_i = i->expression();
448 const Variable lambda1_i(row_index);
449 const Variable lambda2_i(m + row_index);
450 for (Constraint::expr_type::const_iterator j = e_i.begin(),
451 j_end = e_i.lower_bound(Variable(n)); j != j_end; ++j) {
452 Coefficient_traits::const_reference Ap_ij = *j;
453 const Variable v = j.variable();
454 // lambda_1 A'
455 add_mul_assign(les_eq[v.id()], Ap_ij, lambda1_i);
456 // lambda_2 A'
457 add_mul_assign(les_eq[v.id()+n+n], Ap_ij, lambda2_i);
458 }
459 for (Constraint::expr_type::const_iterator
460 j = e_i.lower_bound(Variable(n)),
461 j_end = e_i.end(); j != j_end; ++j) {
462 Coefficient_traits::const_reference A_ij = *j;
463 const Variable v = j.variable();
464 // (lambda_1 - lambda_2) A
465 add_mul_assign(les_eq[v.id()], A_ij, lambda1_i);
466 sub_mul_assign(les_eq[v.id()], A_ij, lambda2_i);
467 // lambda_2 A
468 add_mul_assign(les_eq[v.id()+n], A_ij, lambda2_i);
469 }
470 Coefficient_traits::const_reference b = e_i.inhomogeneous_term();
471 if (b != 0) {
472 // lambda2 b
473 add_mul_assign(le_out, b, lambda2_i);
474 }
475 }
476
477 // Add the non-negativity constraints for lambda_1 and lambda_2.
478 for (dimension_type i = 2*m; i-- > 0; ) {
479 cs_out.insert(Variable(i) >= 0);
480 }
481
482 for (dimension_type j = 3*n; j-- > 0; ) {
483 cs_out.insert(les_eq[j] == 0);
484 }
485 }
486
487 bool
termination_test_MS(const Constraint_System & cs)488 termination_test_MS(const Constraint_System& cs) {
489 Constraint_System cs_mip;
490 fill_constraint_systems_MS(cs, cs_mip, cs_mip);
491
492 const MIP_Problem mip = MIP_Problem(cs_mip.space_dimension(), cs_mip);
493 return mip.is_satisfiable();
494 }
495
496 bool
one_affine_ranking_function_MS(const Constraint_System & cs,Generator & mu)497 one_affine_ranking_function_MS(const Constraint_System& cs, Generator& mu) {
498 Constraint_System cs_mip;
499 fill_constraint_systems_MS(cs, cs_mip, cs_mip);
500
501 const MIP_Problem mip = MIP_Problem(cs_mip.space_dimension(), cs_mip);
502 if (!mip.is_satisfiable()) {
503 return false;
504 }
505 const Generator fp = mip.feasible_point();
506 PPL_ASSERT(fp.is_point());
507 const dimension_type n = cs.space_dimension() / 2;
508 const Linear_Expression le(fp.expression(), n + 1);
509 mu = point(le, fp.divisor());
510 return true;
511 }
512
513 void
all_affine_ranking_functions_MS(const Constraint_System & cs,C_Polyhedron & mu_space)514 all_affine_ranking_functions_MS(const Constraint_System& cs,
515 C_Polyhedron& mu_space) {
516 Constraint_System cs_out1;
517 Constraint_System cs_out2;
518 fill_constraint_systems_MS(cs, cs_out1, cs_out2);
519
520 C_Polyhedron ph1(cs_out1);
521 C_Polyhedron ph2(cs_out2);
522 const dimension_type n = cs.space_dimension() / 2;
523 ph1.remove_higher_space_dimensions(n);
524 ph1.add_space_dimensions_and_embed(1);
525 ph2.remove_higher_space_dimensions(n+1);
526
527 #if PRINT_DEBUG_INFO
528 Variable::output_function_type* p_default_output_function
529 = Variable::get_output_function();
530 Variable::set_output_function(output_function_MS);
531
532 output_function_MS_n = n;
533 output_function_MS_m = num_constraints(cs);
534
535 std::cout << "*** ph1 projected ***" << std::endl;
536 output_function_MS_which = 4;
537 using namespace IO_Operators;
538 std::cout << ph1.minimized_constraints() << std::endl;
539
540 std::cout << "*** ph2 projected ***" << std::endl;
541 std::cout << ph2.minimized_constraints() << std::endl;
542 #endif
543
544 ph1.intersection_assign(ph2);
545
546 #if PRINT_DEBUG_INFO
547 std::cout << "*** intersection ***" << std::endl;
548 using namespace IO_Operators;
549 std::cout << ph1.minimized_constraints() << std::endl;
550
551 Variable::set_output_function(p_default_output_function);
552 #endif
553
554 mu_space.m_swap(ph1);
555 }
556
557 void
all_affine_quasi_ranking_functions_MS(const Constraint_System & cs,C_Polyhedron & decreasing_mu_space,C_Polyhedron & bounded_mu_space)558 all_affine_quasi_ranking_functions_MS(const Constraint_System& cs,
559 C_Polyhedron& decreasing_mu_space,
560 C_Polyhedron& bounded_mu_space) {
561 Constraint_System cs_out1;
562 Constraint_System cs_out2;
563 fill_constraint_systems_MS(cs, cs_out1, cs_out2);
564
565 C_Polyhedron ph1(cs_out1);
566 C_Polyhedron ph2(cs_out2);
567 const dimension_type n = cs.space_dimension() / 2;
568 ph1.remove_higher_space_dimensions(n);
569 ph1.add_space_dimensions_and_embed(1);
570 ph2.remove_higher_space_dimensions(n+1);
571
572 #if PRINT_DEBUG_INFO
573 Variable::output_function_type* p_default_output_function
574 = Variable::get_output_function();
575 Variable::set_output_function(output_function_MS);
576
577 output_function_MS_n = n;
578 output_function_MS_m = num_constraints(cs);
579
580 std::cout << "*** ph1 projected ***" << std::endl;
581 output_function_MS_which = 4;
582 using namespace IO_Operators;
583 std::cout << ph1.minimized_constraints() << std::endl;
584
585 std::cout << "*** ph2 projected ***" << std::endl;
586 std::cout << ph2.minimized_constraints() << std::endl;
587
588 Variable::set_output_function(p_default_output_function);
589 #endif
590
591 decreasing_mu_space.m_swap(ph1);
592 bounded_mu_space.m_swap(ph2);
593 }
594
595 bool
termination_test_PR_original(const Constraint_System & cs)596 termination_test_PR_original(const Constraint_System& cs) {
597 PPL_ASSERT(cs.space_dimension() % 2 == 0);
598
599 Constraint_System cs_mip;
600 Linear_Expression le_ineq;
601 fill_constraint_system_PR_original(cs, cs_mip, le_ineq);
602
603 // Turn minimization problem into satisfiability.
604 cs_mip.insert(le_ineq <= -1);
605
606 const MIP_Problem mip(cs_mip.space_dimension(), cs_mip);
607 return mip.is_satisfiable();
608 }
609
610 bool
termination_test_PR(const Constraint_System & cs_before,const Constraint_System & cs_after)611 termination_test_PR(const Constraint_System& cs_before,
612 const Constraint_System& cs_after) {
613 Constraint_System cs_mip;
614 Linear_Expression le_ineq;
615 fill_constraint_system_PR(cs_before, cs_after, cs_mip, le_ineq);
616
617 #if PRINT_DEBUG_INFO
618 Variable::output_function_type* p_default_output_function
619 = Variable::get_output_function();
620 Variable::set_output_function(output_function_PR);
621
622 output_function_PR_r = num_constraints(cs_before);
623 output_function_PR_s = num_constraints(cs_after);
624
625 std::cout << "*** cs_mip ***" << std::endl;
626 using namespace IO_Operators;
627 std::cout << cs_mip << std::endl;
628 std::cout << "*** le_ineq ***" << std::endl;
629 std::cout << le_ineq << std::endl;
630
631 Variable::set_output_function(p_default_output_function);
632 #endif
633
634 // Turn minimization problem into satisfiability.
635 cs_mip.insert(le_ineq <= -1);
636
637 const MIP_Problem mip(cs_mip.space_dimension(), cs_mip);
638 return mip.is_satisfiable();
639 }
640
641 bool
one_affine_ranking_function_PR(const Constraint_System & cs_before,const Constraint_System & cs_after,Generator & mu)642 one_affine_ranking_function_PR(const Constraint_System& cs_before,
643 const Constraint_System& cs_after,
644 Generator& mu) {
645 return Termination_Helpers
646 ::one_affine_ranking_function_PR(cs_before, cs_after, mu);
647 }
648
649 bool
one_affine_ranking_function_PR_original(const Constraint_System & cs,Generator & mu)650 one_affine_ranking_function_PR_original(const Constraint_System& cs,
651 Generator& mu) {
652 return Termination_Helpers::one_affine_ranking_function_PR_original(cs, mu);
653 }
654
655 void
all_affine_ranking_functions_PR(const Constraint_System & cs_before,const Constraint_System & cs_after,NNC_Polyhedron & mu_space)656 all_affine_ranking_functions_PR(const Constraint_System& cs_before,
657 const Constraint_System& cs_after,
658 NNC_Polyhedron& mu_space) {
659 Termination_Helpers::all_affine_ranking_functions_PR(cs_before, cs_after,
660 mu_space);
661 }
662
663 void
all_affine_ranking_functions_PR_original(const Constraint_System & cs,NNC_Polyhedron & mu_space)664 all_affine_ranking_functions_PR_original(const Constraint_System& cs,
665 NNC_Polyhedron& mu_space) {
666 Termination_Helpers::all_affine_ranking_functions_PR_original(cs, mu_space);
667 }
668
669 } // namespace Termination
670
671 } // namespace Implementation
672
673 bool
674 Termination_Helpers
one_affine_ranking_function_PR(const Constraint_System & cs_before,const Constraint_System & cs_after,Generator & mu)675 ::one_affine_ranking_function_PR(const Constraint_System& cs_before,
676 const Constraint_System& cs_after,
677 Generator& mu) {
678 Constraint_System cs_mip;
679 Linear_Expression le_ineq;
680 Parma_Polyhedra_Library::Implementation::Termination
681 ::fill_constraint_system_PR(cs_before, cs_after, cs_mip, le_ineq);
682
683 #if PRINT_DEBUG_INFO
684 Variable::output_function_type* p_default_output_function
685 = Variable::get_output_function();
686 Variable::set_output_function(output_function_PR);
687
688 output_function_PR_r = num_constraints(cs_before);
689 output_function_PR_s = num_constraints(cs_after);
690
691 std::cout << "*** cs_mip ***" << std::endl;
692 using namespace IO_Operators;
693 std::cout << cs_mip << std::endl;
694 std::cout << "*** le_ineq ***" << std::endl;
695 std::cout << le_ineq << std::endl;
696
697 Variable::set_output_function(p_default_output_function);
698 #endif
699
700 // Turn minimization problem into satisfiability.
701 cs_mip.insert(le_ineq <= -1);
702
703 const MIP_Problem mip(cs_mip.space_dimension(), cs_mip);
704 if (!mip.is_satisfiable()) {
705 return false;
706 }
707
708 const Generator& fp = mip.feasible_point();
709 PPL_ASSERT(fp.is_point());
710
711 // u_3 corresponds to space dimensions 0, ..., s - 1.
712 const dimension_type n = cs_before.space_dimension();
713 // mu_0 is zero: properly set space dimension.
714 Linear_Expression le;
715 le.set_space_dimension(1 + n);
716 // Multiply u_3 by E'_C to obtain mu_1, ..., mu_n.
717 dimension_type row_index = 0;
718 for (Constraint_System::const_iterator i = cs_after.begin(),
719 cs_after_end = cs_after.end();
720 i != cs_after_end;
721 ++i, ++row_index) {
722 Coefficient_traits::const_reference
723 fp_i = fp.coefficient(Variable(row_index));
724 if (fp_i != 0) {
725 le.linear_combine(i->expr, 1, -fp_i, 1, n + 1);
726 }
727 }
728 // Note that we can neglect the divisor of `fp' since it is positive.
729 mu = point(le);
730 return true;
731 }
732
733 bool
734 Termination_Helpers
one_affine_ranking_function_PR_original(const Constraint_System & cs,Generator & mu)735 ::one_affine_ranking_function_PR_original(const Constraint_System& cs,
736 Generator& mu) {
737 PPL_ASSERT(cs.space_dimension() % 2 == 0);
738 const dimension_type n = cs.space_dimension() / 2;
739 const dimension_type m = Implementation::num_constraints(cs);
740
741 Constraint_System cs_mip;
742 Linear_Expression le_ineq;
743 Parma_Polyhedra_Library::Implementation::Termination
744 ::fill_constraint_system_PR_original(cs, cs_mip, le_ineq);
745
746 #if PRINT_DEBUG_INFO
747 std::cout << "*** cs_mip ***" << std::endl;
748 using namespace IO_Operators;
749 std::cout << cs_mip << std::endl;
750 std::cout << "*** le_ineq ***" << std::endl;
751 std::cout << le_ineq << std::endl;
752 #endif
753
754 // Turn minimization problem into satisfiability.
755 cs_mip.insert(le_ineq <= -1);
756
757 const MIP_Problem mip = MIP_Problem(cs_mip.space_dimension(), cs_mip);
758 if (!mip.is_satisfiable()) {
759 return false;
760 }
761 const Generator& fp = mip.feasible_point();
762 PPL_ASSERT(fp.is_point());
763 // mu_0 is zero: properly set space dimension.
764 Linear_Expression le;
765 le.set_space_dimension(1 + n);
766 // Multiply -lambda_2 by A' to obtain mu_1, ..., mu_n.
767 // lambda_2 corresponds to space dimensions m, ..., 2*m - 1.
768 dimension_type row_index = m;
769 for (Constraint_System::const_iterator i = cs.begin(),
770 cs_end = cs.end(); i != cs_end; ++i, ++row_index) {
771 const Variable lambda_2(row_index);
772 Coefficient_traits::const_reference fp_i = fp.coefficient(lambda_2);
773 if (fp_i != 0) {
774 le.linear_combine(i->expr, 1, -fp_i, 1, n + 1);
775 }
776 }
777 // Note that we can neglect the divisor of `fp' since it is positive.
778 mu = point(le);
779 return true;
780 }
781
782 void
783 Termination_Helpers
all_affine_ranking_functions_PR(const Constraint_System & cs_before,const Constraint_System & cs_after,NNC_Polyhedron & mu_space)784 ::all_affine_ranking_functions_PR(const Constraint_System& cs_before,
785 const Constraint_System& cs_after,
786 NNC_Polyhedron& mu_space) {
787 Constraint_System cs_eqs;
788 Linear_Expression le_ineq;
789 Parma_Polyhedra_Library::Implementation::Termination
790 ::fill_constraint_system_PR(cs_before, cs_after, cs_eqs, le_ineq);
791
792 #if PRINT_DEBUG_INFO
793 Variable::output_function_type* p_default_output_function
794 = Variable::get_output_function();
795 Variable::set_output_function(output_function_PR);
796
797 output_function_PR_r = num_constraints(cs_before);
798 output_function_PR_s = num_constraints(cs_after);
799
800 std::cout << "*** cs_eqs ***" << std::endl;
801 using namespace IO_Operators;
802 std::cout << cs_eqs << std::endl;
803 std::cout << "*** le_ineq ***" << std::endl;
804 std::cout << le_ineq << std::endl;
805 #endif
806
807 NNC_Polyhedron ph(cs_eqs);
808 ph.add_constraint(le_ineq < 0);
809 // u_3 corresponds to space dimensions 0, ..., s - 1.
810 const dimension_type s = Implementation::num_constraints(cs_after);
811 ph.remove_higher_space_dimensions(s);
812
813 #if PRINT_DEBUG_INFO
814 std::cout << "*** ph ***" << std::endl;
815 std::cout << ph << std::endl;
816
817 Variable::set_output_function(p_default_output_function);
818 #endif
819
820 const dimension_type n = cs_before.space_dimension();
821
822 const Generator_System& gs_in = ph.generators();
823 Generator_System gs_out;
824 Generator_System::const_iterator gs_in_it = gs_in.begin();
825 Generator_System::const_iterator gs_in_end = gs_in.end();
826 if (gs_in_it == gs_in_end) {
827 // The system is unsatisfiable.
828 mu_space = NNC_Polyhedron(n + 1, EMPTY);
829 }
830 else {
831 for ( ; gs_in_it != gs_in_end; ++gs_in_it) {
832 const Generator& g = *gs_in_it;
833 Linear_Expression le;
834 le.set_space_dimension(n);
835 // Set le to the multiplication of Linear_Expression(g) by E'_C.
836 dimension_type row_index = 0;
837 for (Constraint_System::const_iterator i = cs_after.begin(),
838 cs_after_end = cs_after.end();
839 i != cs_after_end;
840 ++i, ++row_index) {
841 Coefficient_traits::const_reference
842 g_i = g.coefficient(Variable(row_index));
843 if (g_i != 0) {
844 le.linear_combine(i->expr, 1, -g_i, 1, n + 1);
845 }
846 }
847
848 // Add to gs_out the transformed generator.
849 switch (g.type()) {
850 case Generator::LINE:
851 if (!le.all_homogeneous_terms_are_zero()) {
852 gs_out.insert(line(le));
853 }
854 break;
855 case Generator::RAY:
856 if (!le.all_homogeneous_terms_are_zero()) {
857 gs_out.insert(ray(le));
858 }
859 break;
860 case Generator::POINT:
861 gs_out.insert(point(le, g.divisor()));
862 break;
863 case Generator::CLOSURE_POINT:
864 gs_out.insert(closure_point(le, g.divisor()));
865 break;
866 }
867 }
868
869 mu_space = NNC_Polyhedron(gs_out);
870 // mu_0 is zero.
871 mu_space.add_space_dimensions_and_embed(1);
872 }
873 }
874
875 void
876 Termination_Helpers
all_affine_ranking_functions_PR_original(const Constraint_System & cs,NNC_Polyhedron & mu_space)877 ::all_affine_ranking_functions_PR_original(const Constraint_System& cs,
878 NNC_Polyhedron& mu_space) {
879 PPL_ASSERT(cs.space_dimension() % 2 == 0);
880 const dimension_type n = cs.space_dimension() / 2;
881 const dimension_type m = Implementation::num_constraints(cs);
882
883 if (m == 0) {
884 // If there are no constraints at all, we have non-termination,
885 // i.e., no affine ranking function at all.
886 mu_space = NNC_Polyhedron(n + 1, EMPTY);
887 return;
888 }
889
890 Constraint_System cs_eqs;
891 Linear_Expression le_ineq;
892 Parma_Polyhedra_Library::Implementation::Termination
893 ::fill_constraint_system_PR_original(cs, cs_eqs, le_ineq);
894
895 NNC_Polyhedron ph(cs_eqs);
896 ph.add_constraint(le_ineq < 0);
897 // lambda_2 corresponds to space dimensions m, ..., 2*m-1.
898 const Variables_Set lambda1(Variable(0), Variable(m-1));
899 ph.remove_space_dimensions(lambda1);
900
901 #if PRINT_DEBUG_INFO
902 std::cout << "*** ph ***" << std::endl;
903 std::cout << ph << std::endl;
904
905 Variable::set_output_function(p_default_output_function);
906 #endif
907
908 const Generator_System& gs_in = ph.generators();
909 Generator_System::const_iterator gs_in_it = gs_in.begin();
910 Generator_System::const_iterator gs_in_end = gs_in.end();
911 if (gs_in_it == gs_in_end) {
912 // The system is unsatisfiable.
913 mu_space = NNC_Polyhedron(n + 1, EMPTY);
914 }
915 else {
916 Generator_System gs_out;
917 for ( ; gs_in_it != gs_in_end; ++gs_in_it) {
918 const Generator& g = *gs_in_it;
919 Linear_Expression le;
920 le.set_space_dimension(n);
921 // Set le to the multiplication of Linear_Expression(g) by E'_C.
922 dimension_type row_index = 0;
923 for (Constraint_System::const_iterator i = cs.begin(),
924 cs_end = cs.end(); i != cs_end; ++i, ++row_index) {
925 const Variable lambda2_i(row_index);
926 Coefficient_traits::const_reference g_i = g.coefficient(lambda2_i);
927 if (g_i != 0) {
928 le.linear_combine(i->expr, 1, -g_i, 1, n + 1);
929 }
930 }
931
932 // Add to gs_out the transformed generator.
933 switch (g.type()) {
934 case Generator::LINE:
935 if (!le.all_homogeneous_terms_are_zero()) {
936 gs_out.insert(line(le));
937 }
938 break;
939 case Generator::RAY:
940 if (!le.all_homogeneous_terms_are_zero()) {
941 gs_out.insert(ray(le));
942 }
943 break;
944 case Generator::POINT:
945 gs_out.insert(point(le, g.divisor()));
946 break;
947 case Generator::CLOSURE_POINT:
948 gs_out.insert(closure_point(le, g.divisor()));
949 break;
950 }
951 }
952
953 mu_space = NNC_Polyhedron(gs_out);
954 // mu_0 is zero.
955 mu_space.add_space_dimensions_and_embed(1);
956 }
957 }
958
959 } // namespace Parma_Polyhedra_Library
960