// 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_INTERNAL_GSL_NUMERIC_SOLVER_H #define CGAL_POLYNOMIAL_INTERNAL_GSL_NUMERIC_SOLVER_H #include #include #include #include #include namespace CGAL { namespace POLYNOMIAL { namespace internal { /*template inline double gsl_max_error() { if (CLEAN) return .0000005; else return 0; }*/ template inline void gsl_compute_roots(const double *begin, const double *end, double lb, double ub, std::vector &roots_) { //double max_error=gsl_max_error(); gsl_poly_complex_workspace *workspace; int degree= end-begin-1; workspace= gsl_poly_complex_workspace_alloc(degree+1); //! I can't use auto_ptr because it is an array double *roots= new double[2*degree]; /*for ( int i=0; i< 2*fn.degree(); ++i){ roots[i]=0; }*/ int ret = gsl_poly_complex_solve(begin, degree+1, workspace, roots); roots_.reserve(ret); if (ret!=GSL_SUCCESS) { std::cerr << "GSL solver did not converge on "; std::copy(begin, end, std::ostream_iterator(std::cerr, " ")); std::cerr << std::endl; } double last= -std::numeric_limits::infinity(); for ( int i=degree-1; i>=0; --i) { double r= roots[2*i]; double c= roots[2*i+1]; if (root_is_good(r, c, lb, ub)) { roots_.push_back(r); } else if (CLEAN && root_is_good(r,c, last, ub)) { last= r; } } std::sort(roots_.begin(), roots_.end(), std::greater() ); delete roots; gsl_poly_complex_workspace_free(workspace); if (CLEAN) filter_solver_roots(begin,end, lb, ub, last, roots_); return; } template inline void gsl_compute_cubic_roots(const double *begin, const double *end, double lb, double ub, std::vector &roots_) { CGAL_Polynomial_precondition(begin[3] != 0); //double max_error=gsl_max_error(); double r[3]; int num_roots= gsl_poly_solve_cubic(begin[2]/begin[3], begin[1]/begin[3], begin[0]/begin[3], &r[0],&r[1],&r[2]); roots_.reserve(num_roots); double last= -std::numeric_limits::infinity(); // I want reverse sorted roots for (int i=num_roots-1; i>=0; --i) { if (r[i]> lb && r[i] < ub) roots_.push_back(r[i]); else if (CLEAN && r[i] last){ last= r[i]; } } if (CLEAN) filter_solver_roots(begin, end, lb, ub, last, roots_); } inline void gsl_polynomial_compute_roots(const double *begin, const double *end, double lb, double ub, std::vector &roots) { int degree= end-begin-1; switch( degree) { case -1: case 0: break; case 1: compute_linear_roots(begin,end, lb, ub, roots); break; case 2: compute_quadratic_roots(begin, end, lb, ub, roots); break; case 3: gsl_compute_cubic_roots(begin, end, lb, ub, roots); break; default: gsl_compute_roots(begin, end, lb, ub, roots); } } inline void gsl_polynomial_compute_cleaned_roots(const double *begin, const double *end, double lb, double ub, std::vector &roots) { int degree= end-begin-1; switch( degree) { case -1: case 0: break; case 1: compute_linear_cleaned_roots(begin,end, lb, ub, roots); break; case 2: compute_quadratic_cleaned_roots(begin, end, lb, ub, roots); break; case 3: gsl_compute_cubic_roots(begin, end, lb, ub, roots); break; default: gsl_compute_roots(begin, end, lb, ub, roots); } } inline double gsl_evaluate_polynomial(const double *b, const double *e, double t) { if (b==e) return 0; /*double *d= new double[coefs.size()]; for (unsigned int i=0; i< coefs.size(); ++i){ d[i]= coefs[i]; }*/ double v= gsl_poly_eval(b, e-b, t); return v; } /* inline void gsl_polynomial_compute_roots(const double *begin, const double *end, double lb, double ub, std::vector &roots){ CGAL_assertion(0&& no_gsl_support_built); } inline void gsl_polynomial_compute_cleaned_roots(const double *begin, const double *end, double lb, double ub, std::vector &roots){ CGAL_assertion(0&& no_gsl_support_built); } inline double gsl_evaluate_polynomial(const double *b, const double *e, double t) { CGAL_assertion(0&& no_gsl_support_built); return std::numeric_limits::nan(); } */ /*// GSL void gsl_polynomial_compute_roots(const double *begin, const double *end, double lb, double ub, std::vector &roots); void gsl_polynomial_compute_cleaned_roots(const double *begin, const double *end, double lb, double ub, std::vector &roots);*/ struct GSL_numeric_solver { void operator()(const double *begin, const double *end, double lb, double ub, std::vector &roots) const { gsl_polynomial_compute_roots(begin, end, lb, ub, roots); } }; struct GSL_cleaned_numeric_solver { void operator()(const double *begin, const double *end, double lb, double ub, std::vector &roots) const { gsl_polynomial_compute_cleaned_roots(begin, end, lb, ub, roots); } }; } } } //namespace CGAL::POLYNOMIAL::internal #endif