00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #ifndef PPL_OR_Matrix_inlines_hh
00025 #define PPL_OR_Matrix_inlines_hh 1
00026
00027 #include "globals.defs.hh"
00028 #include "Checked_Number.defs.hh"
00029 #include "C_Polyhedron.defs.hh"
00030 #include "distances.defs.hh"
00031 #include <cassert>
00032 #include <algorithm>
00033 #include "checked.defs.hh"
00034
00035 namespace Parma_Polyhedra_Library {
00036
00037 template <typename T>
00038 inline dimension_type
00039 OR_Matrix<T>::row_first_element_index(const dimension_type k) {
00040 return ((k+1)*(k+1))/2;
00041 }
00042
00043 template <typename T>
00044 inline dimension_type
00045 OR_Matrix<T>::row_size(const dimension_type k) {
00046 return (k+2) & ~dimension_type(1);
00047 }
00048
00049 #if PPL_OR_MATRIX_EXTRA_DEBUG
00050
00051 template <typename T>
00052 template <typename U>
00053 inline dimension_type
00054 OR_Matrix<T>::Pseudo_Row<U>::size() const {
00055 return size_;
00056 }
00057
00058 #endif // PPL_OR_MATRIX_EXTRA_DEBUG
00059
00060 template <typename T>
00061 template <typename U>
00062 inline
00063 OR_Matrix<T>::Pseudo_Row<U>::Pseudo_Row()
00064 : first(0)
00065 #if PPL_OR_MATRIX_EXTRA_DEBUG
00066 , size_(0)
00067 #endif
00068 {
00069
00070 }
00071
00072 template <typename T>
00073 template <typename U>
00074 inline
00075 OR_Matrix<T>::Pseudo_Row<U>::Pseudo_Row(U& y
00076 #if PPL_OR_MATRIX_EXTRA_DEBUG
00077 , dimension_type s
00078 #endif
00079 )
00080 : first(&y)
00081 #if PPL_OR_MATRIX_EXTRA_DEBUG
00082 , size_(s)
00083 #endif
00084 {
00085 }
00086
00087 template <typename T>
00088 template <typename U>
00089 template <typename V>
00090 inline
00091 OR_Matrix<T>::Pseudo_Row<U>::Pseudo_Row(const Pseudo_Row<V>& y)
00092 : first(y.first)
00093 #if PPL_OR_MATRIX_EXTRA_DEBUG
00094 , size_(y.size_)
00095 #endif
00096 {
00097 }
00098
00099 template <typename T>
00100 template <typename U>
00101 inline OR_Matrix<T>::Pseudo_Row<U>&
00102 OR_Matrix<T>::Pseudo_Row<U>::operator=(const Pseudo_Row& y) {
00103 first = y.first;
00104 #if PPL_OR_MATRIX_EXTRA_DEBUG
00105 size_ = y.size_;
00106 #endif
00107 return *this;
00108 }
00109
00110 template <typename T>
00111 template <typename U>
00112 inline
00113 OR_Matrix<T>::Pseudo_Row<U>::~Pseudo_Row() {
00114 }
00115
00116 template <typename T>
00117 template <typename U>
00118 inline U&
00119 OR_Matrix<T>::Pseudo_Row<U>::operator[](const dimension_type k) const {
00120 #if PPL_OR_MATRIX_EXTRA_DEBUG
00121 assert(k < size_);
00122 #endif
00123 return *(first + k);
00124 }
00125
00126 template <typename T>
00127 template <typename U>
00128 inline
00129 OR_Matrix<T>::any_row_iterator<U>
00130 ::any_row_iterator(const dimension_type n_rows)
00131 : value(),
00132 e(n_rows)
00133
00134 {
00135 }
00136
00137 template <typename T>
00138 template <typename U>
00139 inline
00140 OR_Matrix<T>::any_row_iterator<U>::any_row_iterator(U& base)
00141 : value(base
00142 #if PPL_OR_MATRIX_EXTRA_DEBUG
00143 , OR_Matrix<T>::row_size(0)
00144 #endif
00145 ),
00146 e(0),
00147 i(0) {
00148 }
00149
00150 template <typename T>
00151 template <typename U>
00152 template <typename V>
00153 inline
00154 OR_Matrix<T>::any_row_iterator<U>
00155 ::any_row_iterator(const any_row_iterator<V>& y)
00156 : value(y.value),
00157 e(y.e),
00158 i(y.i) {
00159 }
00160
00161 template <typename T>
00162 template <typename U>
00163 template <typename V>
00164 inline typename OR_Matrix<T>::template any_row_iterator<U>&
00165 OR_Matrix<T>::any_row_iterator<U>::operator=(const any_row_iterator<V>& y) {
00166 value = y.value;
00167 e = y.e;
00168 i = y.i;
00169 return *this;
00170 }
00171
00172 template <typename T>
00173 template <typename U>
00174 inline typename OR_Matrix<T>::template any_row_iterator<U>::reference
00175 OR_Matrix<T>::any_row_iterator<U>::operator*() const {
00176 return value;
00177 }
00178
00179 template <typename T>
00180 template <typename U>
00181 inline typename OR_Matrix<T>::template any_row_iterator<U>::pointer
00182 OR_Matrix<T>::any_row_iterator<U>::operator->() const {
00183 return &value;
00184 }
00185
00186 template <typename T>
00187 template <typename U>
00188 inline typename OR_Matrix<T>::template any_row_iterator<U>&
00189 OR_Matrix<T>::any_row_iterator<U>::operator++() {
00190 ++e;
00191 dimension_type increment = e;
00192 if (e % 2) {
00193 ++increment;
00194 #if PPL_OR_MATRIX_EXTRA_DEBUG
00195 value.size_ += 2;
00196 #endif
00197 }
00198 i += increment;
00199 value.first += increment;
00200 return *this;
00201 }
00202
00203 template <typename T>
00204 template <typename U>
00205 inline typename OR_Matrix<T>::template any_row_iterator<U>
00206 OR_Matrix<T>::any_row_iterator<U>::operator++(int) {
00207 any_row_iterator old = *this;
00208 ++(*this);
00209 return old;
00210 }
00211
00212 template <typename T>
00213 template <typename U>
00214 inline typename OR_Matrix<T>::template any_row_iterator<U>&
00215 OR_Matrix<T>::any_row_iterator<U>::operator--() {
00216 dimension_type decrement = e + 1;
00217 --e;
00218 if (e % 2) {
00219 ++decrement;
00220 #if PPL_OR_MATRIX_EXTRA_DEBUG
00221 value.size_ -= 2;
00222 #endif
00223 }
00224 i -= decrement;
00225 value.first -= decrement;
00226 return *this;
00227 }
00228
00229 template <typename T>
00230 template <typename U>
00231 inline typename OR_Matrix<T>::template any_row_iterator<U>
00232 OR_Matrix<T>::any_row_iterator<U>::operator--(int) {
00233 any_row_iterator old = *this;
00234 --(*this);
00235 return old;
00236 }
00237
00238 template <typename T>
00239 template <typename U>
00240 inline typename OR_Matrix<T>::template any_row_iterator<U>&
00241 OR_Matrix<T>::any_row_iterator<U>::operator+=(difference_type m) {
00242 difference_type increment = m + m*m/2 + m*e;
00243 if (e%2 == 0 && m%2 == 1)
00244 ++increment;
00245 e += m;
00246 i += increment;
00247 value.first += increment;
00248 #if PPL_OR_MATRIX_EXTRA_DEBUG
00249
00250 value.size_ = OR_Matrix::row_size(e);
00251 #endif
00252 return *this;
00253 }
00254
00255 template <typename T>
00256 template <typename U>
00257 inline typename OR_Matrix<T>::template any_row_iterator<U>&
00258 OR_Matrix<T>::any_row_iterator<U>::operator-=(difference_type m) {
00259 return *this += -m;
00260 }
00261
00262 template <typename T>
00263 template <typename U>
00264 inline typename OR_Matrix<T>::template any_row_iterator<U>::difference_type
00265 OR_Matrix<T>::any_row_iterator<U>::operator-(const any_row_iterator& y) const {
00266 return e - y.e;
00267 }
00268
00269 template <typename T>
00270 template <typename U>
00271 inline typename OR_Matrix<T>::template any_row_iterator<U>
00272 OR_Matrix<T>::any_row_iterator<U>::operator+(difference_type m) const {
00273 any_row_iterator r = *this;
00274 r += m;
00275 return r;
00276 }
00277
00278 template <typename T>
00279 template <typename U>
00280 inline typename OR_Matrix<T>::template any_row_iterator<U>
00281 OR_Matrix<T>::any_row_iterator<U>::operator-(const difference_type m) const {
00282 any_row_iterator r = *this;
00283 r -= m;
00284 return r;
00285 }
00286
00287 template <typename T>
00288 template <typename U>
00289 inline bool
00290 OR_Matrix<T>::any_row_iterator<U>
00291 ::operator==(const any_row_iterator& y) const {
00292 return e == y.e;
00293 }
00294
00295 template <typename T>
00296 template <typename U>
00297 inline bool
00298 OR_Matrix<T>::any_row_iterator<U>
00299 ::operator!=(const any_row_iterator& y) const {
00300 return e != y.e;
00301 }
00302
00303 template <typename T>
00304 template <typename U>
00305 inline bool
00306 OR_Matrix<T>::any_row_iterator<U>::operator<(const any_row_iterator& y) const {
00307 return e < y.e;
00308 }
00309
00310 template <typename T>
00311 template <typename U>
00312 inline bool
00313 OR_Matrix<T>::any_row_iterator<U>
00314 ::operator<=(const any_row_iterator& y) const {
00315 return e <= y.e;
00316 }
00317
00318 template <typename T>
00319 template <typename U>
00320 inline bool
00321 OR_Matrix<T>::any_row_iterator<U>::operator>(const any_row_iterator& y) const {
00322 return e > y.e;
00323 }
00324
00325 template <typename T>
00326 template <typename U>
00327 inline bool
00328 OR_Matrix<T>::any_row_iterator<U>
00329 ::operator>=(const any_row_iterator& y) const {
00330 return e >= y.e;
00331 }
00332
00333 template <typename T>
00334 template <typename U>
00335 inline dimension_type
00336 OR_Matrix<T>::any_row_iterator<U>::row_size() const {
00337 return (e+2) & ~dimension_type(1);
00338 }
00339
00340 template <typename T>
00341 template <typename U>
00342 inline dimension_type
00343 OR_Matrix<T>::any_row_iterator<U>::index() const {
00344 return e;
00345 }
00346
00347 template <typename T>
00348 inline typename OR_Matrix<T>::row_iterator
00349 OR_Matrix<T>::row_begin() {
00350 return num_rows() == 0 ? row_iterator(0) : row_iterator(vec[0]);
00351 }
00352
00353 template <typename T>
00354 inline typename OR_Matrix<T>::row_iterator
00355 OR_Matrix<T>::row_end() {
00356 return row_iterator(num_rows());
00357 }
00358
00359 template <typename T>
00360 inline typename OR_Matrix<T>::const_row_iterator
00361 OR_Matrix<T>::row_begin() const {
00362 return num_rows() == 0 ? const_row_iterator(0) : const_row_iterator(vec[0]);
00363 }
00364
00365 template <typename T>
00366 inline typename OR_Matrix<T>::const_row_iterator
00367 OR_Matrix<T>::row_end() const {
00368 return const_row_iterator(num_rows());
00369 }
00370
00371 template <typename T>
00372 inline typename OR_Matrix<T>::element_iterator
00373 OR_Matrix<T>::element_begin() {
00374 return vec.begin();
00375 }
00376
00377 template <typename T>
00378 inline typename OR_Matrix<T>::element_iterator
00379 OR_Matrix<T>::element_end() {
00380 return vec.end();
00381 }
00382
00383 template <typename T>
00384 inline typename OR_Matrix<T>::const_element_iterator
00385 OR_Matrix<T>::element_begin() const {
00386 return vec.begin();
00387 }
00388
00389 template <typename T>
00390 inline typename OR_Matrix<T>::const_element_iterator
00391 OR_Matrix<T>::element_end() const {
00392 return vec.end();
00393 }
00394
00395 template <typename T>
00396 inline void
00397 OR_Matrix<T>::swap(OR_Matrix& y) {
00398 std::swap(vec, y.vec);
00399 std::swap(space_dim, y.space_dim);
00400 std::swap(vec_capacity, y.vec_capacity);
00401 }
00402
00404 inline unsigned long
00405 isqrt(unsigned long x) {
00406 unsigned long r = 0;
00407 for (unsigned long t = 0x40000000; t; t >>= 2) {
00408 unsigned long s = r + t;
00409 if (s <= x) {
00410 x -= s;
00411 r = s + t;
00412 }
00413 r >>= 1;
00414 }
00415 return r;
00416 }
00417
00418 template <typename T>
00419 inline dimension_type
00420 OR_Matrix<T>::max_num_rows() {
00421
00422
00423
00424 dimension_type k = isqrt(2*DB_Row<T>::max_size() + 1);
00425 return (k-1) & ~dimension_type(1);
00426 }
00427
00428 template <typename T>
00429 inline memory_size_type
00430 OR_Matrix<T>::total_memory_in_bytes() const {
00431 return sizeof(*this) + external_memory_in_bytes();
00432 }
00433
00434 template <typename T>
00435 inline
00436 OR_Matrix<T>::OR_Matrix(const dimension_type dim)
00437 : vec(2*dim*(dim+1)),
00438 space_dim(dim),
00439 vec_capacity(vec.size()) {
00440 }
00441
00442 template <typename T>
00443 inline
00444 OR_Matrix<T>::~OR_Matrix() {
00445 }
00446
00447 template <typename T>
00448 inline typename OR_Matrix<T>::row_reference_type
00449 OR_Matrix<T>::operator[](dimension_type k) {
00450 return row_reference_type(vec[row_first_element_index(k)]
00451 #if PPL_OR_MATRIX_EXTRA_DEBUG
00452 , row_size(k)
00453 #endif
00454 );
00455 }
00456
00457 template <typename T>
00458 inline typename OR_Matrix<T>::const_row_reference_type
00459 OR_Matrix<T>::operator[](dimension_type k) const {
00460 return const_row_reference_type(vec[row_first_element_index(k)]
00461 #if PPL_OR_MATRIX_EXTRA_DEBUG
00462 , row_size(k)
00463 #endif
00464 );
00465 }
00466
00467 template <typename T>
00468 inline dimension_type
00469 OR_Matrix<T>::space_dimension() const {
00470 return space_dim;
00471 }
00472
00473 template <typename T>
00474 inline dimension_type
00475 OR_Matrix<T>::num_rows() const {
00476 return 2*space_dimension();
00477 }
00478
00479 template <typename T>
00480 inline void
00481 OR_Matrix<T>::clear() {
00482 OR_Matrix<T>(0).swap(*this);
00483 }
00484
00485 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00486
00487 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
00488 template <typename T>
00489 inline bool
00490 operator==(const OR_Matrix<T>& x, const OR_Matrix<T>& y) {
00491 return x.space_dim == y.space_dim && x.vec == y.vec;
00492 }
00493
00494 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00495
00496 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
00497 template <typename T>
00498 inline bool
00499 operator!=(const OR_Matrix<T>& x, const OR_Matrix<T>& y) {
00500 return !(x == y);
00501 }
00502
00503 template <typename T>
00504 inline
00505 OR_Matrix<T>::OR_Matrix(const OR_Matrix& y)
00506 : vec(y.vec),
00507 space_dim(y.space_dim),
00508 vec_capacity(compute_capacity(y.vec.size())) {
00509 }
00510
00511 template <typename T>
00512 template <typename U>
00513 inline
00514 OR_Matrix<T>::OR_Matrix(const OR_Matrix<U>& y)
00515 : vec(),
00516 space_dim(y.space_dim),
00517 vec_capacity(compute_capacity(y.vec.size())) {
00518 vec.construct_upward_approximation(y.vec, vec_capacity);
00519 assert(OK());
00520 }
00521
00522 template <typename T>
00523 inline OR_Matrix<T>&
00524 OR_Matrix<T>::operator=(const OR_Matrix& y) {
00525 vec = y.vec;
00526 space_dim = y.space_dim;
00527 vec_capacity = compute_capacity(y.vec.size());
00528 return *this;
00529 }
00530
00531 template <typename T>
00532 inline void
00533 OR_Matrix<T>::grow(const dimension_type new_dim) {
00534 assert(new_dim >= space_dim);
00535 if (new_dim > space_dim) {
00536 const dimension_type new_size = 2*new_dim*(new_dim + 1);
00537 if (new_size <= vec_capacity) {
00538
00539 vec.expand_within_capacity(new_size);
00540 space_dim = new_dim;
00541 }
00542 else {
00543
00544 OR_Matrix<T> new_matrix(new_dim);
00545 element_iterator j = new_matrix.element_begin();
00546 for (element_iterator i = element_begin(),
00547 mend = element_end(); i != mend; ++i, ++j)
00548 assign_or_swap(*j, *i);
00549 swap(new_matrix);
00550 }
00551 }
00552 }
00553
00554 template <typename T>
00555 inline void
00556 OR_Matrix<T>::shrink(const dimension_type new_dim) {
00557 assert(new_dim <= space_dim);
00558 const dimension_type new_size = 2*new_dim*(new_dim + 1);
00559 vec.shrink(new_size);
00560 space_dim = new_dim;
00561 }
00562
00563 template <typename T>
00564 inline void
00565 OR_Matrix<T>::resize_no_copy(const dimension_type new_dim) {
00566 if (new_dim > space_dim) {
00567 const dimension_type new_size = 2*new_dim*(new_dim + 1);
00568 if (new_size <= vec_capacity) {
00569
00570 vec.expand_within_capacity(new_size);
00571 space_dim = new_dim;
00572 }
00573 else {
00574
00575 OR_Matrix<T> new_matrix(new_dim);
00576 swap(new_matrix);
00577 }
00578 }
00579 else if (new_dim < space_dim)
00580 shrink(new_dim);
00581 }
00582
00583 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00584
00585 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
00586 template <typename Specialization, typename Temp, typename To, typename T>
00587 inline bool
00588 l_m_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
00589 const OR_Matrix<T>& x,
00590 const OR_Matrix<T>& y,
00591 const Rounding_Dir dir,
00592 Temp& tmp0,
00593 Temp& tmp1,
00594 Temp& tmp2) {
00595 if (x.num_rows() != y.num_rows())
00596 return false;
00597 assign_r(tmp0, 0, ROUND_NOT_NEEDED);
00598 for (typename OR_Matrix<T>::const_element_iterator
00599 i = x.element_begin(), j = y.element_begin(),
00600 mat_end = x.element_end(); i != mat_end; ++i, ++j) {
00601 const T& x_i = *i;
00602 const T& y_i = *j;
00603 if (is_plus_infinity(x_i)) {
00604 if (is_plus_infinity(y_i))
00605 continue;
00606 else {
00607 pinf:
00608 assign_r(r, PLUS_INFINITY, ROUND_NOT_NEEDED);
00609 return true;
00610 }
00611 }
00612 else if (is_plus_infinity(y_i))
00613 goto pinf;
00614
00615 const Temp* tmp1p;
00616 const Temp* tmp2p;
00617 if (x_i > y_i) {
00618 maybe_assign(tmp1p, tmp1, x_i, dir);
00619 maybe_assign(tmp2p, tmp2, y_i, inverse(dir));
00620 }
00621 else {
00622 maybe_assign(tmp1p, tmp1, y_i, dir);
00623 maybe_assign(tmp2p, tmp2, x_i, inverse(dir));
00624 }
00625 sub_assign_r(tmp1, *tmp1p, *tmp2p, dir);
00626 assert(sgn(tmp1) >= 0);
00627 Specialization::combine(tmp0, tmp1, dir);
00628 }
00629
00630 Specialization::finalize(tmp0, dir);
00631 assign_r(r, tmp0, dir);
00632 return true;
00633 }
00634
00635
00636 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00637
00638 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
00639 template <typename Temp, typename To, typename T>
00640 inline bool
00641 rectilinear_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
00642 const OR_Matrix<T>& x,
00643 const OR_Matrix<T>& y,
00644 const Rounding_Dir dir,
00645 Temp& tmp0,
00646 Temp& tmp1,
00647 Temp& tmp2) {
00648 return
00649 l_m_distance_assign<Rectilinear_Distance_Specialization<Temp> >(r, x, y,
00650 dir,
00651 tmp0,
00652 tmp1,
00653 tmp2);
00654 }
00655
00656 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00657
00658 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
00659 template <typename Temp, typename To, typename T>
00660 inline bool
00661 euclidean_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
00662 const OR_Matrix<T>& x,
00663 const OR_Matrix<T>& y,
00664 const Rounding_Dir dir,
00665 Temp& tmp0,
00666 Temp& tmp1,
00667 Temp& tmp2) {
00668 return
00669 l_m_distance_assign<Euclidean_Distance_Specialization<Temp> >(r, x, y,
00670 dir,
00671 tmp0,
00672 tmp1,
00673 tmp2);
00674 }
00675
00676 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00677
00678 #endif // defined(PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS)
00679 template <typename Temp, typename To, typename T>
00680 inline bool
00681 l_infinity_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
00682 const OR_Matrix<T>& x,
00683 const OR_Matrix<T>& y,
00684 const Rounding_Dir dir,
00685 Temp& tmp0,
00686 Temp& tmp1,
00687 Temp& tmp2) {
00688 return
00689 l_m_distance_assign<L_Infinity_Distance_Specialization<Temp> >(r, x, y,
00690 dir,
00691 tmp0,
00692 tmp1,
00693 tmp2);
00694 }
00695
00696 }
00697
00698 namespace std {
00699
00701 template <typename T>
00702 inline void
00703 swap(Parma_Polyhedra_Library::OR_Matrix<T>& x,
00704 Parma_Polyhedra_Library::OR_Matrix<T>& y) {
00705 x.swap(y);
00706 }
00707
00708 }
00709
00710
00711 #endif // !defined(PPL_OR_Matrix_inlines_hh)