// Copyright (c) 2005 Stanford University (USA). // All rights reserved. // // This file is part of CGAL (www.cgal.org); you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public License as // published by the Free Software Foundation; either version 3 of the License, // or (at your option) any later version. // // Licensees holding a valid commercial license may use this file in // accordance with the commercial license agreement provided with the software. // // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // // $URL$ // $Id$ // // // Author(s) : Daniel Russel #ifndef CGAL_POLYNOMIAL_NUMERIC_ROOT_STACK_H #define CGAL_POLYNOMIAL_NUMERIC_ROOT_STACK_H #include #include #include #include #include #ifdef CGAL_USE_GSL #include #endif namespace CGAL { namespace POLYNOMIAL { namespace internal { template ::has_infinity> struct Numeric_pick_root { typedef double Root; //typedef internal::Double_with_infinity Root; }; template <> struct Numeric_pick_root { typedef internal::Double_with_infinity Root; }; } } } //namespace CGAL::POLYNOMIAL::internal namespace CGAL { namespace POLYNOMIAL { #ifdef CGAL_USE_GSL #define CGAL_DEFAULT_NUMERIC_SOLVER CGAL::POLYNOMIAL::internal::GSL_numeric_solver #define CGAL_DEFAULT_CLEANED_NUMERIC_SOLVER CGAL::POLYNOMIAL::internal::GSL_cleaned_numeric_solver #else #define CGAL_DEFAULT_NUMERIC_SOLVER CGAL::POLYNOMIAL::internal::Turkowski_numeric_solver #define CGAL_DEFAULT_CLEANED_NUMERIC_SOLVER CGAL::POLYNOMIAL::internal::Turkowski_cleaned_numeric_solver #endif template class Numeric_root_stack { public: typedef internal::Numeric_pick_root<>::Root Root; protected: /* All this mess is to handle when roots are not doubles and when the coefficients are not doubles */ template void initialize_2(const double *b, const double *e, double lb, double ub, std::vector &roots) { std::vector lroots; ns_(b, e, lb, ub, lroots); roots.insert(roots.end(), lroots.begin(), lroots.end()); } void initialize_2(const double *b, const double *e, double lb, double ub, std::vector &roots) { ns_(b, e, lb, ub, roots); } template void initialize(const Fn &f, Root lb, Root ub) { std::vector c(f.degree()+1); for (unsigned int i=0; i<= f.degree(); ++i) { c[i]= to_double(f[i]); } initialize_2(&*c.begin(), &*c.begin()+ f.degree()+1, static_cast(lb), static_cast(ub), roots_); /*if (CLEAN) { polynomial_compute_cleaned_roots(&*c.begin(), &*c.begin()+ f.degree()+1, lb, ub, roots_); } else { polynomial_compute_roots(&*c.begin(), &*c.begin()+ f.degree()+1, lb, ub, roots_); }*/ } void initialize(const Polynomial &f, Root lb, Root ub) { const double *p0= &*f.begin(); initialize_2(p0, p0+ f.degree()+1, static_cast(lb), static_cast(ub), roots_); /*if (CLEAN) { polynomial_compute_cleaned_roots(&*f.begin(), &*f.begin()+ f.degree()+1, lb, ub, roots_); } else { polynomial_compute_roots(&*f.begin(), &*f.begin()+ f.degree()+1, lb, ub, roots_); }*/ } public: typedef typename Solver_traits::Function Function; typedef Solver_traits Traits; Numeric_root_stack(const typename Solver_traits::Function &f, Root lb, Root ub, const Solver_traits&) { //std::cout << "Solving " << f << " from " << lb << " to " << ub; initialize(f, lb, ub); /*if (!roots_.empty()) std::cout << " got " << roots_.back() << std::endl; else std::cout << std::endl;*/ #if 0 for (unsigned int i=1; i < roots_.size(); ++i) { if (roots_[i]> roots_[i-1]){ std::cerr << "ERROR: roots out of order "; std::copy(roots_.begin(), roots_.end(), std::ostream_iterator(std::cerr, " ")); std::cerr << " for " << f << " from " << lb << " to " << ub << std::endl; roots_.clear(); initialize(f,lb, ub); } } #endif //CGAL_Polynomial_postcondition(roots_[i] <= roots_[i-1]); } Numeric_root_stack(){}; void pop() { CGAL_Polynomial_precondition(!roots_.empty()); roots_.pop_back(); } const Root& top() const { CGAL_Polynomial_precondition(!roots_.empty()); return roots_.back(); } bool empty() const { return roots_.empty(); } protected: Numeric_solver ns_; std::vector roots_; }; } } //namespace CGAL::POLYNOMIAL #endif