1 /* Linear_Form 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_Linear_Form_templates_hh
25 #define PPL_Linear_Form_templates_hh 1
26
27 #include "Linear_Form_defs.hh"
28 #include "Linear_Expression_defs.hh"
29 #include "Box_defs.hh"
30 #include <stdexcept>
31 #include <iostream>
32 #include <cmath>
33
34 namespace Parma_Polyhedra_Library {
35
36 template <typename C>
Linear_Form(const Variable v)37 Linear_Form<C>::Linear_Form(const Variable v)
38 : vec() {
39 const dimension_type space_dim = v.space_dimension();
40 if (space_dim > max_space_dimension()) {
41 throw std::length_error("Linear_Form<C>::"
42 "Linear_Form(v):\n"
43 "v exceeds the maximum allowed "
44 "space dimension.");
45 }
46 vec.reserve(compute_capacity(space_dim+1, vec_type().max_size()));
47 vec.resize(space_dim+1, zero);
48 vec[v.space_dimension()] = C(typename C::boundary_type(1));
49 }
50
51 template <typename C>
Linear_Form(const Variable v,const Variable w)52 Linear_Form<C>::Linear_Form(const Variable v, const Variable w)
53 : vec() {
54 const dimension_type v_space_dim = v.space_dimension();
55 const dimension_type w_space_dim = w.space_dimension();
56 const dimension_type space_dim = std::max(v_space_dim, w_space_dim);
57 if (space_dim > max_space_dimension()) {
58 throw std::length_error("Linear_Form<C>::"
59 "Linear_Form(v, w):\n"
60 "v or w exceed the maximum allowed "
61 "space dimension.");
62 }
63 vec.reserve(compute_capacity(space_dim+1, vec_type().max_size()));
64 vec.resize(space_dim+1, zero);
65 if (v_space_dim != w_space_dim) {
66 vec[v_space_dim] = C(typename C::boundary_type(1));
67 vec[w_space_dim] = C(typename C::boundary_type(-1));
68 }
69 }
70
71 template <typename C>
Linear_Form(const Linear_Expression & e)72 Linear_Form<C>::Linear_Form(const Linear_Expression& e)
73 : vec() {
74 const dimension_type space_dim = e.space_dimension();
75 if (space_dim > max_space_dimension()) {
76 throw std::length_error("Linear_Form<C>::"
77 "Linear_Form(e):\n"
78 "e exceeds the maximum allowed "
79 "space dimension.");
80 }
81 vec.reserve(compute_capacity(space_dim+1, vec_type().max_size()));
82 vec.resize(space_dim+1);
83 for (dimension_type i = space_dim; i-- > 0; )
84 vec[i+1] = e.coefficient(Variable(i));
85 vec[0] = e.inhomogeneous_term();
86 }
87
88 /*! \relates Linear_Form */
89 template <typename C>
90 Linear_Form<C>
operator +(const Linear_Form<C> & f1,const Linear_Form<C> & f2)91 operator+(const Linear_Form<C>& f1, const Linear_Form<C>& f2) {
92 dimension_type f1_size = f1.size();
93 dimension_type f2_size = f2.size();
94 dimension_type min_size;
95 dimension_type max_size;
96 const Linear_Form<C>* p_e_max;
97 if (f1_size > f2_size) {
98 min_size = f2_size;
99 max_size = f1_size;
100 p_e_max = &f1;
101 }
102 else {
103 min_size = f1_size;
104 max_size = f2_size;
105 p_e_max = &f2;
106 }
107
108 Linear_Form<C> r(max_size, false);
109 dimension_type i = max_size;
110 while (i > min_size) {
111 --i;
112 r[i] = p_e_max->vec[i];
113 }
114 while (i > 0) {
115 --i;
116 r[i] = f1[i];
117 r[i] += f2[i];
118 }
119 return r;
120 }
121
122 /*! \relates Linear_Form */
123 template <typename C>
124 Linear_Form<C>
operator +(const Variable v,const Linear_Form<C> & f)125 operator+(const Variable v, const Linear_Form<C>& f) {
126 const dimension_type v_space_dim = v.space_dimension();
127 if (v_space_dim > Linear_Form<C>::max_space_dimension()) {
128 throw std::length_error("Linear_Form "
129 "operator+(v, f):\n"
130 "v exceeds the maximum allowed "
131 "space dimension.");
132 }
133 Linear_Form<C> r(f);
134 if (v_space_dim > f.space_dimension()) {
135 r.extend(v_space_dim+1);
136 }
137 r[v_space_dim] += C(typename C::boundary_type(1));
138 return r;
139 }
140
141 /*! \relates Linear_Form */
142 template <typename C>
143 Linear_Form<C>
operator +(const C & n,const Linear_Form<C> & f)144 operator+(const C& n, const Linear_Form<C>& f) {
145 Linear_Form<C> r(f);
146 r[0] += n;
147 return r;
148 }
149
150 /*! \relates Linear_Form */
151 template <typename C>
152 Linear_Form<C>
operator -(const Linear_Form<C> & f)153 operator-(const Linear_Form<C>& f) {
154 Linear_Form<C> r(f);
155 for (dimension_type i = f.size(); i-- > 0; ) {
156 r[i].neg_assign(r[i]);
157 }
158 return r;
159 }
160
161 /*! \relates Linear_Form */
162 template <typename C>
163 Linear_Form<C>
operator -(const Linear_Form<C> & f1,const Linear_Form<C> & f2)164 operator-(const Linear_Form<C>& f1, const Linear_Form<C>& f2) {
165 dimension_type f1_size = f1.size();
166 dimension_type f2_size = f2.size();
167 if (f1_size > f2_size) {
168 Linear_Form<C> r(f1_size, false);
169 dimension_type i = f1_size;
170 while (i > f2_size) {
171 --i;
172 r[i] = f1[i];
173 }
174 while (i > 0) {
175 --i;
176 r[i] = f1[i];
177 r[i] -= f2[i];
178 }
179 return r;
180 }
181 else {
182 Linear_Form<C> r(f2_size, false);
183 dimension_type i = f2_size;
184 while (i > f1_size) {
185 --i;
186 r[i].neg_assign(f2[i]);
187 }
188 while (i > 0) {
189 --i;
190 r[i] = f1[i];
191 r[i] -= f2[i];
192 }
193 return r;
194 }
195 }
196
197 /*! \relates Linear_Form */
198 template <typename C>
199 Linear_Form<C>
operator -(const Variable v,const Linear_Form<C> & f)200 operator-(const Variable v, const Linear_Form<C>& f) {
201 const dimension_type v_space_dim = v.space_dimension();
202 if (v_space_dim > Linear_Form<C>::max_space_dimension()) {
203 throw std::length_error("Linear_Form "
204 "operator-(v, e):\n"
205 "v exceeds the maximum allowed "
206 "space dimension.");
207 }
208 Linear_Form<C> r(f);
209 if (v_space_dim > f.space_dimension()) {
210 r.extend(v_space_dim+1);
211 }
212 for (dimension_type i = f.size(); i-- > 0; ) {
213 r[i].neg_assign(r[i]);
214 }
215 r[v_space_dim] += C(typename C::boundary_type(1));
216 return r;
217 }
218
219 /*! \relates Linear_Form */
220 template <typename C>
221 Linear_Form<C>
operator -(const Linear_Form<C> & f,const Variable v)222 operator-(const Linear_Form<C>& f, const Variable v) {
223 const dimension_type v_space_dim = v.space_dimension();
224 if (v_space_dim > Linear_Form<C>::max_space_dimension()) {
225 throw std::length_error("Linear_Form "
226 "operator-(e, v):\n"
227 "v exceeds the maximum allowed "
228 "space dimension.");
229 }
230 Linear_Form<C> r(f);
231 if (v_space_dim > f.space_dimension()) {
232 r.extend(v_space_dim+1);
233 }
234 r[v_space_dim] -= C(typename C::boundary_type(1));
235 return r;
236 }
237
238 /*! \relates Linear_Form */
239 template <typename C>
240 Linear_Form<C>
operator -(const C & n,const Linear_Form<C> & f)241 operator-(const C& n, const Linear_Form<C>& f) {
242 Linear_Form<C> r(f);
243 for (dimension_type i = f.size(); i-- > 0; ) {
244 r[i].neg_assign(r[i]);
245 }
246 r[0] += n;
247 return r;
248 }
249
250 /*! \relates Linear_Form */
251 template <typename C>
252 Linear_Form<C>
operator *(const C & n,const Linear_Form<C> & f)253 operator*(const C& n, const Linear_Form<C>& f) {
254 Linear_Form<C> r(f);
255 for (dimension_type i = f.size(); i-- > 0; ) {
256 r[i] *= n;
257 }
258 return r;
259 }
260
261 /*! \relates Linear_Form */
262 template <typename C>
263 Linear_Form<C>&
operator +=(Linear_Form<C> & f1,const Linear_Form<C> & f2)264 operator+=(Linear_Form<C>& f1, const Linear_Form<C>& f2) {
265 dimension_type f1_size = f1.size();
266 dimension_type f2_size = f2.size();
267 if (f1_size < f2_size) {
268 f1.extend(f2_size);
269 }
270 for (dimension_type i = f2_size; i-- > 0; ) {
271 f1[i] += f2[i];
272 }
273 return f1;
274 }
275
276 /*! \relates Linear_Form */
277 template <typename C>
278 Linear_Form<C>&
operator +=(Linear_Form<C> & f,const Variable v)279 operator+=(Linear_Form<C>& f, const Variable v) {
280 const dimension_type v_space_dim = v.space_dimension();
281 if (v_space_dim > Linear_Form<C>::max_space_dimension()) {
282 throw std::length_error("Linear_Form<C>& "
283 "operator+=(e, v):\n"
284 "v exceeds the maximum allowed space dimension.");
285 }
286 if (v_space_dim > f.space_dimension()) {
287 f.extend(v_space_dim+1);
288 }
289 f[v_space_dim] += C(typename C::boundary_type(1));
290 return f;
291 }
292
293 /*! \relates Linear_Form */
294 template <typename C>
295 Linear_Form<C>&
operator -=(Linear_Form<C> & f1,const Linear_Form<C> & f2)296 operator-=(Linear_Form<C>& f1, const Linear_Form<C>& f2) {
297 dimension_type f1_size = f1.size();
298 dimension_type f2_size = f2.size();
299 if (f1_size < f2_size) {
300 f1.extend(f2_size);
301 }
302 for (dimension_type i = f2_size; i-- > 0; ) {
303 f1[i] -= f2[i];
304 }
305 return f1;
306 }
307
308 /*! \relates Linear_Form */
309 template <typename C>
310 Linear_Form<C>&
operator -=(Linear_Form<C> & f,const Variable v)311 operator-=(Linear_Form<C>& f, const Variable v) {
312 const dimension_type v_space_dim = v.space_dimension();
313 if (v_space_dim > Linear_Form<C>::max_space_dimension()) {
314 throw std::length_error("Linear_Form<C>& "
315 "operator-=(e, v):\n"
316 "v exceeds the maximum allowed space dimension.");
317 }
318 if (v_space_dim > f.space_dimension()) {
319 f.extend(v_space_dim+1);
320 }
321 f[v_space_dim] -= C(typename C::boundary_type(1));
322 return f;
323 }
324
325 /*! \relates Linear_Form */
326 template <typename C>
327 Linear_Form<C>&
operator *=(Linear_Form<C> & f,const C & n)328 operator*=(Linear_Form<C>& f, const C& n) {
329 dimension_type f_size = f.size();
330 for (dimension_type i = f_size; i-- > 0; ) {
331 f[i] *= n;
332 }
333 return f;
334 }
335
336 /*! \relates Linear_Form */
337 template <typename C>
338 Linear_Form<C>&
operator /=(Linear_Form<C> & f,const C & n)339 operator/=(Linear_Form<C>& f, const C& n) {
340 dimension_type f_size = f.size();
341 for (dimension_type i = f_size; i-- > 0; ) {
342 f[i] /= n;
343 }
344 return f;
345 }
346
347 /*! \relates Linear_Form */
348 template <typename C>
349 inline bool
operator ==(const Linear_Form<C> & x,const Linear_Form<C> & y)350 operator==(const Linear_Form<C>& x, const Linear_Form<C>& y) {
351 const dimension_type x_size = x.size();
352 const dimension_type y_size = y.size();
353 if (x_size >= y_size) {
354 for (dimension_type i = y_size; i-- > 0; ) {
355 if (x[i] != y[i]) {
356 return false;
357 }
358 }
359 for (dimension_type i = x_size; --i >= y_size; ) {
360 if (x[i] != x.zero) {
361 return false;
362 }
363 }
364 }
365 else {
366 for (dimension_type i = x_size; i-- > 0; ) {
367 if (x[i] != y[i]) {
368 return false;
369 }
370 }
371 for (dimension_type i = y_size; --i >= x_size; ) {
372 if (y[i] != x.zero) {
373 return false;
374 }
375 }
376
377 }
378
379 return true;
380 }
381
382 template <typename C>
383 void
negate()384 Linear_Form<C>::negate() {
385 for (dimension_type i = vec.size(); i-- > 0; ) {
386 vec[i].neg_assign(vec[i]);
387 }
388 return;
389 }
390
391 template <typename C>
392 inline memory_size_type
external_memory_in_bytes() const393 Linear_Form<C>::external_memory_in_bytes() const {
394 memory_size_type n = 0;
395 for (dimension_type i = size(); i-- > 0; ) {
396 n += vec[i].external_memory_in_bytes();
397 }
398 n += vec.capacity()*sizeof(C);
399 return n;
400 }
401
402 template <typename C>
403 bool
OK() const404 Linear_Form<C>::OK() const {
405 for (dimension_type i = size(); i-- > 0; ) {
406 if (!vec[i].OK()) {
407 return false;
408 }
409 }
410 return true;
411 }
412
413 // Floating point analysis related methods.
414 template <typename C>
415 void
relative_error(const Floating_Point_Format analyzed_format,Linear_Form & result) const416 Linear_Form<C>::relative_error(
417 const Floating_Point_Format analyzed_format,
418 Linear_Form& result) const {
419 typedef typename C::boundary_type analyzer_format;
420
421 // Get the necessary information on the analyzed's format.
422 unsigned int f_base;
423 unsigned int f_mantissa_bits;
424 switch (analyzed_format) {
425 case IEEE754_HALF:
426 f_base = float_ieee754_half::BASE;
427 f_mantissa_bits = float_ieee754_half::MANTISSA_BITS;
428 break;
429 case IEEE754_SINGLE:
430 f_base = float_ieee754_single::BASE;
431 f_mantissa_bits = float_ieee754_single::MANTISSA_BITS;
432 break;
433 case IEEE754_DOUBLE:
434 f_base = float_ieee754_double::BASE;
435 f_mantissa_bits = float_ieee754_double::MANTISSA_BITS;
436 break;
437 case IBM_SINGLE:
438 f_base = float_ibm_single::BASE;
439 f_mantissa_bits = float_ibm_single::MANTISSA_BITS;
440 break;
441 case IEEE754_QUAD:
442 f_base = float_ieee754_quad::BASE;
443 f_mantissa_bits = float_ieee754_quad::MANTISSA_BITS;
444 break;
445 case INTEL_DOUBLE_EXTENDED:
446 f_base = float_intel_double_extended::BASE;
447 f_mantissa_bits = float_intel_double_extended::MANTISSA_BITS;
448 break;
449 default:
450 PPL_UNREACHABLE;
451 break;
452 }
453
454 C error_propagator;
455 // We assume that f_base is a power of 2.
456 unsigned int u_power = msb_position(f_base) * f_mantissa_bits;
457 int neg_power = -static_cast<int>(u_power);
458 analyzer_format lb = static_cast<analyzer_format>(ldexp(1.0, neg_power));
459
460 error_propagator.build(i_constraint(GREATER_OR_EQUAL, -lb),
461 i_constraint(LESS_OR_EQUAL, lb));
462
463 // Handle the inhomogeneous term.
464 const C* current_term = &inhomogeneous_term();
465 assert(current_term->is_bounded());
466
467 C current_multiplier(std::max(std::abs(current_term->lower()),
468 std::abs(current_term->upper())));
469 Linear_Form current_result_term(current_multiplier);
470 current_result_term *= error_propagator;
471 result = Linear_Form(current_result_term);
472
473 // Handle the other terms.
474 dimension_type dimension = space_dimension();
475 for (dimension_type i = 0; i < dimension; ++i) {
476 current_term = &coefficient(Variable(i));
477 assert(current_term->is_bounded());
478 current_multiplier = C(std::max(std::abs(current_term->lower()),
479 std::abs(current_term->upper())));
480 current_result_term = Linear_Form(Variable(i));
481 current_result_term *= current_multiplier;
482 current_result_term *= error_propagator;
483 result += current_result_term;
484 }
485
486 return;
487 }
488
489 template <typename C>
490 template <typename Target>
491 bool
intervalize(const FP_Oracle<Target,C> & oracle,C & result) const492 Linear_Form<C>::intervalize(const FP_Oracle<Target,C>& oracle,
493 C& result) const {
494 result = C(inhomogeneous_term());
495 dimension_type dimension = space_dimension();
496 for (dimension_type i = 0; i < dimension; ++i) {
497 C current_addend = coefficient(Variable(i));
498 C curr_int;
499 if (!oracle.get_interval(i, curr_int)) {
500 return false;
501 }
502 current_addend *= curr_int;
503 result += current_addend;
504 }
505
506 return true;
507 }
508
509 /*! \relates Parma_Polyhedra_Library::Linear_Form */
510 template <typename C>
511 std::ostream&
operator <<(std::ostream & s,const Linear_Form<C> & f)512 IO_Operators::operator<<(std::ostream& s, const Linear_Form<C>& f) {
513 const dimension_type num_variables = f.space_dimension();
514 bool first = true;
515 for (dimension_type v = 0; v < num_variables; ++v) {
516 const C& fv = f[v+1];
517 if (fv != typename C::boundary_type(0)) {
518 if (first) {
519 if (fv == typename C::boundary_type(-1)) {
520 s << "-";
521 }
522 else if (fv != typename C::boundary_type(1)) {
523 s << fv << "*";
524 }
525 first = false;
526 }
527 else {
528 if (fv == typename C::boundary_type(-1)) {
529 s << " - ";
530 }
531 else {
532 s << " + ";
533 if (fv != typename C::boundary_type(1)) {
534 s << fv << "*";
535 }
536 }
537 }
538 s << Variable(v);
539 }
540 }
541 // Inhomogeneous term.
542 const C& it = f[0];
543 if (it != 0) {
544 if (!first) {
545 s << " + ";
546 }
547 else {
548 first = false;
549 }
550 s << it;
551 }
552
553 if (first) {
554 // The null linear form.
555 s << Linear_Form<C>::zero;
556 }
557 return s;
558 }
559
560 PPL_OUTPUT_TEMPLATE_DEFINITIONS(C, Linear_Form<C>)
561
562 template <typename C>
563 C Linear_Form<C>::zero(typename C::boundary_type(0));
564
565 } // namespace Parma_Polyhedra_Library
566
567 #endif // !defined(PPL_Linear_Form_templates_hh)
568