// Copyright (c) 2009-2014 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org) // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Kernel_d/include/CGAL/Kernel_d/Interval_linear_algebra.h $ // $Id: Interval_linear_algebra.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Samuel Hornus, Olivier Devillers #ifndef CGAL_INTERVAL_LINEAR_ALGEBRA_H #define CGAL_INTERVAL_LINEAR_ALGEBRA_H #include #include namespace CGAL { namespace internal { template Sign // The matrix is row major: M[i] represents row i. /* FIXME : the function DOES MODIFY the matrix M, but calling functions assume M is not modified --> BUGS. (e.g. Side_of_oriented_subsphere) */ sign_of_determinantDxD_with_interval_arithmetic(Matrix & M) // attempts to compute the determinant using interval arithmetic { typedef Interval_nt NT; // Field number type int sign = 1; int curdim = M.dimension().first; for(int col = 0; col < curdim; ++col) { int pivot_line = -1; NT pivot = 0.0; for(int l = col; l < curdim; ++l) { NT candidate = CGAL::abs(M[l][col]); if( CGAL::certainly(0.0 < candidate) ) if( pivot.inf() < candidate.inf() ) { pivot_line = l; pivot = candidate; } } if( -1 == pivot_line ) { throw CGAL::Uncertain_conversion_exception("undecidable interval computation of determinant"); } // if the pivot value is negative, invert the sign if( M[pivot_line][col] < 0.0 ) sign = - sign; // put the pivot cell on the diagonal if( pivot_line > col ) { std::swap(M[col], M[pivot_line]); // this takes constant time with std::vector<> sign = - sign; } // using matrix-line operations, set zero in all cells below the pivot cell. // This costs roughly (curdim-col-1)^2 mults and adds, because we don't actually // compute the zeros below the pivot cell NT denom = NT(1.0) / M[col][col]; for(int line = col + 1; line < curdim; ++line) { NT numer = M[line][col] * denom; for (int c = col + 1; c < curdim; ++c) M[line][c] -= numer * M[col][c]; } } return Sign(sign); } } // end of namespace internal template<> inline Sign Linear_algebraCd::sign_of_determinant(const Matrix & M) { switch( M.dimension().first ) { case 2: return CGAL::sign_of_determinant( M(0,0), M(0,1), M(1,0), M(1,1)); break; case 3: return CGAL::sign_of_determinant( M(0,0), M(0,1), M(0,2), M(1,0), M(1,1), M(1,2), M(2,0), M(2,1), M(2,2)); break; case 4: return CGAL::sign_of_determinant( M(0,0), M(0,1), M(0,2), M(0,3), M(1,0), M(1,1), M(1,2), M(1,3), M(2,0), M(2,1), M(2,2), M(2,3), M(3,0), M(3,1), M(3,2), M(3,3)); break; default: return internal::sign_of_determinantDxD_with_interval_arithmetic(const_cast(M)); break; } } } //namespace CGAL #endif // CGAL_INTERVAL_LINEAR_ALGEBRA_H