1 /* Dense_Row class implementation (non-inline 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 "Dense_Row_defs.hh"
26 #include "Coefficient_defs.hh"
27 #include "assertions.hh"
28 #include "Sparse_Row_defs.hh"
29 #include <iostream>
30 #include <iomanip>
31
32 namespace PPL = Parma_Polyhedra_Library;
33
Dense_Row(const Sparse_Row & y,dimension_type sz,dimension_type capacity)34 PPL::Dense_Row::Dense_Row(const Sparse_Row& y,
35 dimension_type sz, dimension_type capacity) {
36 resize(sz, capacity);
37 for (Sparse_Row::const_iterator i = y.begin(),
38 i_end = y.lower_bound(std::min(y.size(), sz)); i != i_end; ++i) {
39 (*this)[i.index()] = *i;
40 }
41 PPL_ASSERT(OK());
42 }
43
44 void
resize(dimension_type new_size)45 PPL::Dense_Row::resize(dimension_type new_size) {
46 if (new_size <= size()) {
47 shrink(new_size);
48 }
49 else {
50 if (new_size > capacity()) {
51 // Reallocation is required.
52 // TODO: Consider using realloc() here.
53 // TODO: Consider using a smarter allocation strategy.
54 const dimension_type new_capacity = new_size;
55 Coefficient* const new_vec = impl.coeff_allocator.allocate(new_capacity);
56
57 if (impl.vec != 0) {
58 memcpy(new_vec, impl.vec, sizeof(Coefficient) * impl.size);
59 impl.coeff_allocator.deallocate(impl.vec, impl.capacity);
60 }
61
62 impl.vec = new_vec;
63 impl.capacity = new_capacity;
64 }
65 PPL_ASSERT(new_size <= impl.capacity);
66 // Construct the additional elements.
67 while (impl.size != new_size) {
68 new(&impl.vec[impl.size]) Coefficient();
69 ++impl.size;
70 }
71 }
72 PPL_ASSERT(size() == new_size);
73 PPL_ASSERT(OK());
74 }
75
76 void
resize(dimension_type new_size,dimension_type new_capacity)77 PPL::Dense_Row::resize(dimension_type new_size, dimension_type new_capacity) {
78 PPL_ASSERT(new_size <= new_capacity);
79
80 if (new_capacity == 0) {
81 destroy();
82 impl.vec = 0;
83 impl.size = 0;
84 impl.capacity = 0;
85
86 PPL_ASSERT(size() == new_size);
87 PPL_ASSERT(capacity() == new_capacity);
88 PPL_ASSERT(OK());
89
90 return;
91 }
92
93 if (new_capacity < capacity()) {
94
95 shrink(new_size);
96
97 PPL_ASSERT(impl.size == new_size);
98
99 Coefficient* const new_vec = impl.coeff_allocator.allocate(new_capacity);
100
101 PPL_ASSERT(impl.vec != 0);
102
103 memcpy(new_vec, impl.vec, sizeof(Coefficient) * impl.size);
104
105 impl.coeff_allocator.deallocate(impl.vec, impl.capacity);
106
107 impl.vec = new_vec;
108 impl.capacity = new_capacity;
109 }
110 else {
111 if (new_capacity > capacity()) {
112
113 Coefficient* const new_vec = impl.coeff_allocator.allocate(new_capacity);
114
115 if (impl.vec != 0) {
116 memcpy(new_vec, impl.vec, sizeof(Coefficient) * impl.size);
117 impl.coeff_allocator.deallocate(impl.vec, impl.capacity);
118 }
119
120 impl.vec = new_vec;
121 impl.capacity = new_capacity;
122
123 resize(new_size);
124 }
125 }
126
127 PPL_ASSERT(size() == new_size);
128 PPL_ASSERT(capacity() == new_capacity);
129 PPL_ASSERT(OK());
130 }
131
132 void
clear()133 PPL::Dense_Row::clear() {
134 for (iterator i = begin(), i_end = end(); i != i_end; ++i) {
135 *i = 0;
136 }
137 }
138
139 void
add_zeroes_and_shift(dimension_type n,dimension_type i)140 PPL::Dense_Row::add_zeroes_and_shift(dimension_type n, dimension_type i) {
141 PPL_ASSERT(i <= size());
142 const dimension_type new_size = size() + n;
143 if (new_size > capacity()) {
144 Dense_Row new_row;
145 const dimension_type new_capacity = compute_capacity(new_size, max_size());
146 // This may throw.
147 new_row.impl.vec = new_row.impl.coeff_allocator.allocate(new_capacity);
148 new_row.impl.capacity = new_capacity;
149
150 dimension_type j = i;
151 try {
152 // Construct coefficients with value 0 in
153 // new_row.impl.vec[i ... i + n - 1]
154 for ( ; j < i + n; ++j) {
155 new(&(new_row.impl.vec[j])) Coefficient(0);
156 }
157 } catch (...) {
158 // Destroy the zeroes constructed so far.
159 while (j != i) {
160 --j;
161 new_row.impl.vec[j].~Coefficient();
162 }
163 // The new_row's destructor will de-allocate the memory.
164 throw;
165 }
166
167 // Raw-copy the coefficients.
168 memcpy(new_row.impl.vec, impl.vec, sizeof(Coefficient) * i);
169 memcpy(&(new_row.impl.vec[i + n]), &impl.vec[i],
170 sizeof(Coefficient) * (impl.size - i));
171
172 using std::swap;
173 swap(impl.vec, new_row.impl.vec);
174 swap(impl.capacity, new_row.impl.capacity);
175
176 // *this now owns all coefficients, including the newly-added zeroes.
177 impl.size = new_size;
178
179 // The old vec will be de-allocated at the end of this block.
180
181 }
182 else {
183 memmove(&impl.vec[n + i], &impl.vec[i], sizeof(Coefficient)
184 * (impl.size - i));
185 impl.size = i;
186 const dimension_type target_size = impl.size + n;
187 PPL_ASSERT(target_size == i + n);
188 try {
189 // Construct n zeroes where the moved elements resided.
190 while (impl.size != target_size) {
191 new(&impl.vec[impl.size]) Coefficient(0);
192 ++impl.size;
193 }
194 impl.size = new_size;
195 } catch (...) {
196 // impl.vec[impl.size]..impl.vec[target_size-1] are still unconstructed,
197 // but impl.vec[target_size]..impl.vec[new_size] are constructed,
198 // because the memmove() moved already-constructed objects.
199
200 // NOTE: This loop can't throw, because destructors must not throw.
201 for (dimension_type j = target_size; j < new_size; ++j) {
202 impl.vec[j].~Coefficient();
203 }
204
205 throw;
206 }
207 }
208
209 PPL_ASSERT(OK());
210 }
211
212 void
expand_within_capacity(const dimension_type new_size)213 PPL::Dense_Row::expand_within_capacity(const dimension_type new_size) {
214 PPL_ASSERT(new_size <= impl.capacity);
215 PPL_ASSERT(size() <= new_size && new_size <= max_size());
216 while (impl.size != new_size) {
217 new(&impl.vec[impl.size]) Coefficient();
218 ++impl.size;
219 }
220 PPL_ASSERT(size() == new_size);
221 PPL_ASSERT(OK());
222 }
223
224 void
shrink(dimension_type new_size)225 PPL::Dense_Row::shrink(dimension_type new_size) {
226 PPL_ASSERT(new_size <= size());
227 // Since ~Coefficient() does not throw exceptions, nothing here does.
228
229 // We assume construction was done "forward".
230 // We thus perform destruction "backward".
231 while (impl.size != new_size) {
232 --impl.size;
233 impl.vec[impl.size].~Coefficient();
234 }
235
236 PPL_ASSERT(size() == new_size);
237 PPL_ASSERT(OK());
238 }
239
Dense_Row(const Sparse_Row & row)240 PPL::Dense_Row::Dense_Row(const Sparse_Row& row)
241 : impl() {
242
243 init(row);
244
245 PPL_ASSERT(size() == row.size());
246 PPL_ASSERT(OK());
247 }
248
249 void
init(const Sparse_Row & row)250 PPL::Dense_Row::init(const Sparse_Row& row) {
251 impl.capacity = row.size();
252 impl.vec = impl.coeff_allocator.allocate(impl.capacity);
253 Sparse_Row::const_iterator itr = row.begin();
254 Sparse_Row::const_iterator itr_end = row.end();
255 while (impl.size != impl.capacity) {
256 // Constructs (*this)[impl.size] with row[impl.size].
257 if (itr != itr_end && itr.index() == impl.size) {
258 new(&impl.vec[impl.size]) Coefficient(*itr);
259 ++itr;
260 }
261 else {
262 new(&impl.vec[impl.size]) Coefficient();
263 }
264 ++impl.size;
265 }
266 PPL_ASSERT(size() == row.size());
267 PPL_ASSERT(OK());
268 }
269
270 PPL::Dense_Row&
operator =(const Sparse_Row & row)271 PPL::Dense_Row::operator=(const Sparse_Row& row) {
272 if (size() > row.size()) {
273 // TODO: If the shrink() is modified to reallocate a smaller chunk,
274 // this can be optimized.
275 shrink(row.size());
276 Sparse_Row::const_iterator itr = row.begin();
277 Sparse_Row::const_iterator itr_end = row.end();
278 for (dimension_type i = 0; i < impl.size; ++i) {
279 // Computes (*this)[impl.size] = row[impl.size].
280 if (itr != itr_end && itr.index() == i) {
281 impl.vec[impl.size] = *itr;
282 ++itr;
283 }
284 else {
285 impl.vec[impl.size] = Coefficient_zero();
286 }
287 }
288 }
289 else {
290 if (capacity() >= row.size()) {
291 // size() <= row.size() <= capacity().
292 Sparse_Row::const_iterator itr = row.begin();
293 Sparse_Row::const_iterator itr_end = row.end();
294 for (dimension_type i = 0; i < impl.size; ++i) {
295 // The following code is equivalent to (*this)[i] = row[i].
296 if (itr != itr_end && itr.index() == impl.size) {
297 new(&impl.vec[impl.size]) Coefficient(*itr);
298 ++itr;
299 }
300 else {
301 new(&impl.vec[impl.size]) Coefficient();
302 }
303 }
304 // Construct the additional elements.
305 for ( ; impl.size != row.size(); ++impl.size) {
306 // Constructs (*this)[impl.size] with row[impl.size].
307 if (itr != itr_end && itr.index() == impl.size) {
308 new(&impl.vec[impl.size]) Coefficient(*itr);
309 ++itr;
310 }
311 else {
312 new(&impl.vec[impl.size]) Coefficient();
313 }
314 }
315 }
316 else {
317 // Reallocation is required.
318 destroy();
319 init(row);
320 }
321 }
322 PPL_ASSERT(size() == row.size());
323 PPL_ASSERT(OK());
324
325 return *this;
326 }
327
328 void
normalize()329 PPL::Dense_Row::normalize() {
330 Dense_Row& x = *this;
331 // Compute the GCD of all the coefficients.
332 const dimension_type sz = size();
333 dimension_type i = sz;
334 PPL_DIRTY_TEMP_COEFFICIENT(gcd);
335 while (i > 0) {
336 Coefficient_traits::const_reference x_i = x[--i];
337 if (const int x_i_sign = sgn(x_i)) {
338 gcd = x_i;
339 if (x_i_sign < 0) {
340 neg_assign(gcd);
341 }
342 goto compute_gcd;
343 }
344 }
345 // We reach this point only if all the coefficients were zero.
346 return;
347
348 compute_gcd:
349 if (gcd == 1) {
350 return;
351 }
352 while (i > 0) {
353 Coefficient_traits::const_reference x_i = x[--i];
354 if (x_i != 0) {
355 // Note: we use the ternary version instead of a more concise
356 // gcd_assign(gcd, x_i) to take advantage of the fact that
357 // `gcd' will decrease very rapidly (see D. Knuth, The Art of
358 // Computer Programming, second edition, Section 4.5.2,
359 // Algorithm C, and the discussion following it). Our
360 // implementation of gcd_assign(x, y, z) for checked numbers is
361 // optimized for the case where `z' is smaller than `y', so that
362 // on checked numbers we gain. On the other hand, for the
363 // implementation of gcd_assign(x, y, z) on GMP's unbounded
364 // integers we cannot make any assumption, so here we draw.
365 // Overall, we win.
366 gcd_assign(gcd, x_i, gcd);
367 if (gcd == 1) {
368 return;
369 }
370 }
371 }
372 // Divide the coefficients by the GCD.
373 for (dimension_type j = sz; j-- > 0; ) {
374 Coefficient& x_j = x[j];
375 exact_div_assign(x_j, x_j, gcd);
376 }
377 }
378
379 void
reset(dimension_type first,dimension_type last)380 PPL::Dense_Row::reset(dimension_type first, dimension_type last) {
381 PPL_ASSERT(first <= last);
382 PPL_ASSERT(last <= size());
383 for (dimension_type i = first; i < last; ++i) {
384 (*this)[i] = 0;
385 }
386 }
387
388 void
linear_combine(const Dense_Row & y,Coefficient_traits::const_reference coeff1,Coefficient_traits::const_reference coeff2)389 PPL::Dense_Row::linear_combine(const Dense_Row& y,
390 Coefficient_traits::const_reference coeff1,
391 Coefficient_traits::const_reference coeff2) {
392 Dense_Row& x = *this;
393 PPL_ASSERT(x.size() == y.size());
394
395 x.linear_combine(y, coeff1, coeff2, 0, x.size());
396 }
397
398 void
linear_combine(const Dense_Row & y,Coefficient_traits::const_reference coeff1,Coefficient_traits::const_reference coeff2,dimension_type start,dimension_type end)399 PPL::Dense_Row::linear_combine(const Dense_Row& y,
400 Coefficient_traits::const_reference coeff1,
401 Coefficient_traits::const_reference coeff2,
402 dimension_type start, dimension_type end) {
403 Dense_Row& x = *this;
404 PPL_ASSERT(start <= end);
405 PPL_ASSERT(end <= x.size());
406 PPL_ASSERT(end <= y.size());
407 PPL_ASSERT(coeff1 != 0);
408 PPL_ASSERT(coeff2 != 0);
409
410 // If coeff1 is 1 and/or coeff2 is 1 or -1, we use an optimized
411 // implementation.
412
413 if (coeff1 == 1) {
414 if (coeff2 == 1) {
415 // Optimized implementation for coeff1==1, coeff2==1.
416 for (dimension_type i = start; i < end; ++i) {
417 if (y[i] != 0) {
418 x[i] += y[i];
419 }
420 }
421 return;
422 }
423 if (coeff2 == -1) {
424 // Optimized implementation for coeff1==1, coeff2==-1.
425 for (dimension_type i = start; i < end; ++i) {
426 if (y[i] != 0) {
427 x[i] -= y[i];
428 }
429 }
430 return;
431 }
432 // Optimized implementation for coeff1==1.
433 for (dimension_type i = start; i < end; ++i) {
434 Coefficient& x_i = x[i];
435 // The test against 0 gives rise to a consistent speed up: see
436 // http://www.cs.unipr.it/pipermail/ppl-devel/2009-February/014000.html
437 Coefficient_traits::const_reference y_i = y[i];
438 if (y_i != 0) {
439 add_mul_assign(x_i, y_i, coeff2);
440 }
441 }
442 return;
443 }
444
445 if (coeff2 == 1) {
446 // Optimized implementation for coeff2==1.
447 for (dimension_type i = start; i < end; ++i) {
448 x[i] *= coeff1;
449 if (y[i] != 0) {
450 x[i] += y[i];
451 }
452 }
453 return;
454 }
455 if (coeff2 == -1) {
456 // Optimized implementation for coeff2==-1.
457 for (dimension_type i = start; i < end; ++i) {
458 x[i] *= coeff1;
459 if (y[i] != 0) {
460 x[i] -= y[i];
461 }
462 }
463 return;
464 }
465 // General case.
466 for (dimension_type i = start; i < end; ++i) {
467 Coefficient& x_i = x[i];
468 x[i] *= coeff1;
469 // The test against 0 gives rise to a consistent speed up: see
470 // http://www.cs.unipr.it/pipermail/ppl-devel/2009-February/014000.html
471 Coefficient_traits::const_reference y_i = y[i];
472 if (y_i != 0) {
473 add_mul_assign(x_i, y_i, coeff2);
474 }
475 }
476 }
477
478 void
ascii_dump(std::ostream & s) const479 PPL::Dense_Row::ascii_dump(std::ostream& s) const {
480 const Dense_Row& x = *this;
481 const dimension_type x_size = x.size();
482 s << "size " << x_size << " ";
483 for (dimension_type i = 0; i < x_size; ++i) {
484 s << x[i] << ' ';
485 }
486 s << "\n";
487 }
488
PPL_OUTPUT_DEFINITIONS_ASCII_ONLY(Dense_Row)489 PPL_OUTPUT_DEFINITIONS_ASCII_ONLY(Dense_Row)
490
491 bool
492 PPL::Dense_Row::ascii_load(std::istream& s) {
493 std::string str;
494 if (!(s >> str) || str != "size") {
495 return false;
496 }
497 dimension_type new_size;
498 if (!(s >> new_size)) {
499 return false;
500 }
501
502 resize(new_size);
503
504 for (dimension_type col = 0; col < new_size; ++col) {
505 if (!(s >> (*this)[col])) {
506 return false;
507 }
508 }
509
510 return true;
511 }
512
513 PPL::memory_size_type
external_memory_in_bytes(dimension_type) const514 PPL::Dense_Row::external_memory_in_bytes(dimension_type /* capacity */) const {
515 return external_memory_in_bytes();
516 }
517
518 PPL::memory_size_type
external_memory_in_bytes() const519 PPL::Dense_Row::external_memory_in_bytes() const {
520 memory_size_type n = impl.capacity * sizeof(Coefficient);
521 for (dimension_type i = size(); i-- > 0; ) {
522 n += PPL::external_memory_in_bytes(impl.vec[i]);
523 }
524 return n;
525 }
526
527 bool
OK() const528 PPL::Dense_Row::OK() const {
529 #ifndef NDEBUG
530 using std::endl;
531 using std::cerr;
532 #endif
533
534 bool is_broken = false;
535
536 if (impl.capacity > max_size()) {
537 #ifndef NDEBUG
538 cerr << "Dense_Row capacity exceeds the maximum allowed size:" << endl
539 << "is " << impl.capacity
540 << ", should be less than or equal to " << max_size() << "."
541 << endl;
542 #endif
543 is_broken = true;
544 }
545
546 if (size() > max_size()) {
547 #ifndef NDEBUG
548 cerr << "Dense_Row size exceeds the maximum allowed size:" << endl
549 << "is " << size()
550 << ", should be less than or equal to " << max_size() << "." << endl;
551 #endif
552 is_broken = true;
553 }
554
555 if (impl.capacity < size()) {
556 #ifndef NDEBUG
557 cerr << "Dense_Row is completely broken: capacity is " << impl.capacity
558 << ", size is " << size() << "." << endl;
559 #endif
560 is_broken = true;
561 }
562
563 if (capacity() == 0) {
564 if (impl.vec != 0) {
565 is_broken = true;
566 }
567 }
568 else {
569 if (impl.vec == 0) {
570 is_broken = true;
571 }
572 }
573
574 return !is_broken;
575 }
576
577 bool
OK(const dimension_type row_size) const578 PPL::Dense_Row::OK(const dimension_type row_size) const {
579 #ifndef NDEBUG
580 using std::endl;
581 using std::cerr;
582 #endif
583
584 bool is_broken = !OK();
585
586 // Check the declared size.
587 if (size() != row_size) {
588 #ifndef NDEBUG
589 cerr << "Dense_Row size mismatch: is " << size()
590 << ", should be " << row_size << "." << endl;
591 #endif
592 is_broken = true;
593 }
594 return !is_broken;
595 }
596
597 /*! \relates Parma_Polyhedra_Library::Dense_Row */
598 bool
operator ==(const Dense_Row & x,const Dense_Row & y)599 PPL::operator==(const Dense_Row& x, const Dense_Row& y) {
600 const dimension_type x_size = x.size();
601 const dimension_type y_size = y.size();
602
603 if (x_size != y_size) {
604 return false;
605 }
606
607 for (dimension_type i = x_size; i-- > 0; ) {
608 if (x[i] != y[i]) {
609 return false;
610 }
611 }
612
613 return true;
614 }
615