diff --git a/.gitattributes b/.gitattributes index 8d992a37c03..eadc8b11a58 100644 --- a/.gitattributes +++ b/.gitattributes @@ -258,6 +258,13 @@ Algebraic_kernel_d/doc_tex/Algebraic_kernel_d_ref/intro.tex -text Algebraic_kernel_d/dont_submit -text Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2.h -text Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Algebraic_real_traits.h -text +Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Approximate_arithmetic_controller.h -text +Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Best_approximation_cache.h -text +Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_bfs.h -text +Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_traits_on_vert_line.h -text +Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Non_generic_position_exception.h -text +Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/alg_real_utils.h -text +Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/enums.h -text Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Curve_analysis_2.h -text Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Curve_pair_analysis_2.h -text Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/LRU_hashed_map.h -text diff --git a/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Approximate_arithmetic_controller.h b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Approximate_arithmetic_controller.h new file mode 100644 index 00000000000..9a8a05b7467 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Approximate_arithmetic_controller.h @@ -0,0 +1,229 @@ +// TODO: Add licence +// +// 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) : Michael Kerber +// Ralf Schindlbeck +// +// ============================================================================ + + +#ifndef CGAL_APPROXIMATE_ARITHMETIC_CONTROLLER +#define CGAL_APPROXIMATE_ARITHMETIC_CONTROLLER 1 + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + +template +class Approximate_arithmetic_controller { + +public: + + typedef Poly_1_ Poly_1; + typedef AlgebraicReal Algebraic_real; + + typedef typename Algebraic_real::Rational Boundary; + + typedef typename Poly_1::NT Coefficient; + + typedef typename boost::numeric::interval Boundary_interval; + + typedef typename CGAL::Coercion_traits::Type + Eval_result; + + typedef typename boost::numeric::interval Eval_interval; + + //Default constructor + Approximate_arithmetic_controller() {} + + // Constrcutor with Polynomial and Algebraic real + Approximate_arithmetic_controller(Poly_1 f, Algebraic_real alpha) + : _m_f(f), _m_alpha(alpha) + { +#if AcX_USE_DERIVATIVE_OPTION || AcX_USE_BISECTION_OPTION + f_x = CGAL::diff(_m_f); + flag = false; +#endif + +#if AcX_USE_BISECTION_OPTION + flag_r = false; +#endif + iter_refine_steps = 0; + } + + //! Obains an interval approximation of f(alpha) + Eval_interval interval_approximation() const { + Boundary_interval interval(_m_alpha.low(),_m_alpha.high()); + #if AcX_USE_DERIVATIVE_OPTION || AcX_USE_BISECTION_OPTION + //! AcX_USE_DERIVATIVE_OPTION On or AcX_USE_BISECTION_OPTION On; + //! flag request, because to save runtime + if(!flag) + { + interval_f_x = evaluate_iv(f_x, interval); + CGAL::Sign sign_f_x = CGAL::sign(interval_f_x.lower()) * CGAL::sign(interval_f_x.upper()); + //! find the sign of the derived function _m_f on their bounds + sign_low = f_x.sign_at(_m_alpha.low()); + + if(sign_f_x > 0) + { + //! f_x have no horizontal tangent => interval good ;-) + flag = true; + //std::cout << "function f is strictly monotone" << std::endl; + + #if AcX_USE_BISECTION_OPTION + // help-variables, needed to save the old/new state of the intervals + low = _m_alpha.low(); + high = _m_alpha.high(); + low_f = _m_f.evaluate(_m_alpha.low()); + high_f = _m_f.evaluate(_m_alpha.high()); + #endif + } + else + { + //! f_x have a horizontal tangent => interval was overestimated + //std::cout << "function f is not monotone" << std::endl; + interval_f = evaluate_iv( _m_f, interval ); + + } + } + if(flag) + { + //! computing direct, without IA => because the interval isn't + //! overestimated + #if !AcX_USE_BISECTION_OPTION + if(sign_low == CGAL::POSITIVE) + { + //! slope of _m_f is positive => interval OK + interval_f = Eval_interval( _m_f.evaluate(_m_alpha.low()), _m_f.evaluate(_m_alpha.high()) ); + } + else + { + //! slope of _m_f is negative => switch interval-boundaries + interval_f = Eval_interval( _m_f.evaluate(_m_alpha.high()), _m_f.evaluate(_m_alpha.low()) ); + } + #else + if(!flag_r) + { + if(sign_low == CGAL::POSITIVE) + { + //! slope of _m_f is positive => interval OK + interval_f = Eval_interval( low_f, high_f ); + } + else + { + //! slope of _m_f is negative => switch interval-boundaries + interval_f = Eval_interval( high_f, low_f ); + } + } + else + { + //! request if the low boundary hasn't change after the refinement + if(low == low_r) + { + high = high_r; + high_f = _m_f.evaluate(_m_alpha.high()); + if(sign_low == CGAL::POSITIVE) + { + interval_f = Eval_interval( low_f , high_f ); + } + else + { + interval_f = Eval_interval( high_f, low_f ); + } + } + //! request if the high boundary hasn't change after the refinement + else if(high == high_r) + { + low = low_r; + low_f = _m_f.evaluate(_m_alpha.low()); + if(sign_low == CGAL::POSITIVE) + { + interval_f = Eval_interval( low_f , high_f ); + } + else + { + interval_f = Eval_interval( high_f, low_f ); + } + } + else // for quadratic refinement method + { + //std::cout << "Refinement method without modified option, step: " << i << std::endl; + low = low_r; + high = high_r; + low_f = _m_f.evaluate(_m_alpha.low()); + high_f = _m_f.evaluate(_m_alpha.high()); + if(sign_low == CGAL::POSITIVE) + { + interval_f = Eval_interval( low_f , high_f ); + } + else + { + interval_f = Eval_interval( high_f, low_f ); + } + } + } + #endif + } + return interval_f; + #else + //! AcX_USE_DERIVATIVE_OPTION Off or AcX_USE_BISECTION_OPTION Off; + return evaluate_iv( _m_f, interval ); + #endif + } + + //! Refines the algebraic real + int refine_value() const { + _m_alpha.refine(); + #if AcX_USE_BISECTION_OPTION + // help-variables, needed to save the new intervals after refinement + flag_r = true; + low_r = _m_alpha.low(); + high_r = _m_alpha.high(); + #endif + iter_refine_steps++; + return iter_refine_steps; + } + + +private: + + Poly_1 _m_f; + Algebraic_real _m_alpha; + //! helpvariables for the "derivate" and "bisection" option +#if AcX_USE_DERIVATIVE_OPTION || AcX_USE_BISECTION_OPTION + // variable: indicate the first dervative of polynome f + Poly_1 f_x; + // variable: denotes the sign of the low boundary of the Algebraic_real alpha + mutable CGAL::Sign sign_low; + // variable: terms the intervals for the poylnome f and its first derivative + mutable Eval_interval interval_f, interval_f_x; + // helpvariable: flag for indication, if f is strictly monotone in the boundary + mutable bool flag; +#endif + //! helpvariables only for the "bisection" option +#if AcX_USE_BISECTION_OPTION + // variable: denotes the old state of the maximum borders of alpha + mutable Boundary low, high; + // helpvariable: flag_r for indication, if one refinement step is made + mutable bool flag_r; + // variable: terms the new state of the maximum borders of alpha + mutable Boundary low_r, high_r; + // variable: denotes the state of the maximum borders of the polynome f + mutable Boundary low_f, high_f; +#endif + // variable, which counts the refinement steps + mutable int iter_refine_steps; + +}; //class Approximate_arithmetic_controller + +} // namespace CGALi + +CGAL_END_NAMESPACE + +#endif diff --git a/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Best_approximation_cache.h b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Best_approximation_cache.h new file mode 100644 index 00000000000..6db896e2b4b --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Best_approximation_cache.h @@ -0,0 +1,122 @@ +// TODO: Add licence +// +// 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) : Michael Kerber +// +// ============================================================================ + +#ifndef CGAL_BEST_APPROXIMATION_CACHE +#define CGAL_BEST_APPROXIMATION_CACHE 1 + + +#include +#include +#include +#include + +CGAL_BEGIN_NAMESPACE + + + namespace CGALi { + + /*! + * \brief Stores computed approximations for coefficients. + * + * During the execution of the Bitstream Descartes method, approximations + * of the coefficients are computed independently by different branches of + * the Descartes tree. To avoid computing them several times, the best + * konwn approximation for a special coefficient is stored in a \c std::map + * structure. The data contains the used precision and the integer that + * represents the coefficient for this precision, according to the + * \c CGAL::CGALi::Bitstream_descartes_rndl_tree_traits definition. + */ + template + class Best_approximation_cache + : public ::CGAL::Handle_with_policy > > { + + public: + + //! The Coefficient type + typedef Coefficient_ Coefficient; + + //! The integer type + typedef Integer_ Integer; + + typedef std::pair Data_pair; + + typedef std::map Map; + + typedef ::CGAL::Handle_with_policy Base; + + typedef Best_approximation_cache Self; + + //! Default constructor + Best_approximation_cache() { + } + + //! Copy constructor + Best_approximation_cache(const Self& s) + : Base(static_cast(s)) + {} + + /*! + * \brief Checks whether the coefficient \c c already has already been + * approximated. + */ + bool approximation_exists(Coefficient c) const { + typename Map::const_iterator it = this->ptr()->find(c); + return it!=this->ptr()->end(); + } + + /*! + * \brief The best approximation known for the coefficient c + * is returned. + * + * It is necessary to check whether an approximation of c exists + * using the \c approximation_exists method beforehand! + */ + void get_best_approximation(Coefficient c,long& prec, Integer& val) const{ + typename Map::const_iterator it = this->ptr()->find(c); + CGAL_assertion(it!=this->ptr()->end()); + prec = it->second.first; + val = it->second.second; + } + + /*! + * \brief Updates the approximation memory + * + * If an approximation for c is known, this method can be used + * to store it for later usage. The update is only performed if the + * precision is indeed higher than the current precision in the map. + */ + void update_approximation(Coefficient c,long prec, const Integer& val) { + typename Map::iterator it = this->ptr()->find(c); + if(it!=this->ptr()->end()) { + if(it->second.first >= prec) { + return; + } + else { + Data_pair better_pair = std::make_pair(prec,val); + it->second=better_pair; + } + } + else { + Data_pair new_pair = std::make_pair(prec,val); + this->ptr()->operator[](c)=new_pair; + } + + } + }; + + + }// namespace CGALi + +CGAL_END_NAMESPACE + +#endif // AcX_BEST_APPROXIMATION_CACHE diff --git a/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_bfs.h b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_bfs.h new file mode 100644 index 00000000000..7195cf6ade7 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_bfs.h @@ -0,0 +1,2038 @@ +// TODO: Add licence +// +// 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) : Michael Kerber +// +// ============================================================================ + +#ifndef CGAL_GENERIC_DESCARTES_BFS +#define CGAL_GENERIC_DESCARTES_BFS 1 + + +#include +#include + +#if BITSTREAM_USES_E08_TREE +#include +#else +#include +#endif + +#include + +CGAL_BEGIN_NAMESPACE + + namespace CGALi { + + //! enum to distinguish between different descartes instances + enum Bitstream_descartes_bfs_type { + GENERIC_DESCARTES = 0, + SQUARE_FREE_DESCARTES = 1, //!< uses Square_free_descartes_tag constructor + M_K_DESCARTES = 2, //!< uses M_k_descartes_tag constructor + BACKSHEAR_DESCARTES = 3, //!< uses Backshear_descartes_tag constructor + M_KI_DESCARTES = 4, //!< uses M_ki_descartes_tag constructor + EXCHANGE_DESCARTES = 5, //!< uses Exchange_descartes_tag constructor + INVERSE_TRANSFORM_DESCARTES = 6, //! < uses Inverse_transform_descartes + VERT_LINE_ADAPTER_DESCARTES = 7 // ! < uses Vert_line_adapter_descartes + }; + +// pre-declaration + template + class Bitstream_descartes_bfs; + +//! Tag for the square free Bitstream Descartes method + struct Square_free_descartes_tag {}; + +//! Tag for the Bitstream m-k-Descartes method + struct M_k_descartes_tag {}; + +//! Tag for the Exchange-Descartes method + struct Exchange_descartes_tag {}; + +//! Tag for the Backshear Descartes method + struct Backshear_descartes_tag {}; + +//! Tag for the Bitstream m-ki-Descartes method + struct M_ki_descartes_tag {}; + +//! Tag for the inverse transformation Descartes method + struct Inverse_transform_descartes_tag {}; + +//! Tag fot the vert-line adapter Descartes method + struct Vert_line_adapter_descartes_tag {}; + +//! forward declaration + template + class Bitstream_descartes_bfs; + + +/* + * \brief Thrown whenever a non-specialised virtual member function is called + */ + class Virtual_method_exception {}; + + +/* + * \brief The base class for all variants of the Bitstream Descartes method. + * + */ + template + class Generic_descartes_bfs_rep + : public Policy::template Hierarchy_base::Type { + + public: + + //! The traits class for approximations + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + //! The Handle class + typedef + Bitstream_descartes_bfs + Handle; + + //! The Coeeficient type of the input polynomial + typedef typename Bitstream_descartes_rndl_tree_traits::Coefficient + Coefficient; + + //! The polynomial type + typedef CGAL::Polynomial Polynomial; + + typedef Generic_descartes_bfs_rep Self; + + //! The type of the used Bitstream Descartes tree +#if BITSTREAM_USES_E08_TREE + typedef NiX::Bitstream_descartes_E08_tree + + Bitstream_tree; +#else + typedef CGAL::CGALi::Bitstream_descartes_rndl_tree + + Bitstream_tree; +#endif + + + //! The used integer type + typedef typename Bitstream_descartes_rndl_tree_traits::Integer Integer; + + //! The type for the iterator of the nodes of the bitstream tree + typedef typename Bitstream_tree::Node_iterator + Node_iterator; + + //! The same as constant iterator + typedef typename Bitstream_tree::Node_const_iterator + Node_const_iterator; + + //! How the boundaries of the isolating intervals are represented + typedef typename Bitstream_descartes_rndl_tree_traits::Boundary + Boundary; + + //! Default constructor (does nothing) + Generic_descartes_bfs_rep(Bitstream_descartes_bfs_type type + =GENERIC_DESCARTES) : + type_(type) { + }; + + /*! + * Constructor computing an interval containing all real roots of \c f, + * and initialising the Bitstream Descartes tree + */ + Generic_descartes_bfs_rep(Bitstream_descartes_bfs_type type, + Polynomial f, + Bitstream_descartes_rndl_tree_traits traits) : + type_(type), + f_(f), + traits_(traits), + is_isolated_(false) { + + Integer lower,upper; + long log_div; + this->get_interval(f,lower,upper,log_div,traits); + //AcX_DSTREAM("f: " << f << std::endl); + if(f.degree()>0) { + bitstream_tree +#if BITSTREAM_USES_E08_TREE + = Bitstream_tree(-log_div, + f.begin(), + f.end(), + typename Bitstream_tree::Monomial_basis_tag(), + traits); +#else + = Bitstream_tree(lower,upper,log_div, + f.begin(), + f.end(), + typename Bitstream_tree::Monomial_basis_tag(), + traits); +#endif + + if(bitstream_tree.begin()==bitstream_tree.end()) { + number_of_intervals=0; + } + else { + number_of_intervals=1; + } + } + else { + number_of_intervals=0; + } + } + + /*! + * Constructor that copies the Bitstream tree given from outside + * and initialising the Bitstream Descartes tree + * The tree must "fit" to the polynomial + */ + Generic_descartes_bfs_rep(Bitstream_descartes_bfs_type type, + Polynomial f, + Bitstream_tree tree, + Bitstream_descartes_rndl_tree_traits traits) : + type_(type), + f_(f), + traits_(traits), + bitstream_tree(tree), + is_isolated_(false) { + + tree.set_traits(traits); + + number_of_intervals=0; + + for(Node_iterator curr=bitstream_tree.begin(); + curr != bitstream_tree.end(); + curr++) { + number_of_intervals++; + } + } + + + + //! Destructor (does nothing) + virtual ~Generic_descartes_bfs_rep() { + } + + //! Needed for the referencing counting mechanism + virtual CGAL::Reference_counted_hierarchy<>* clone() { + return new Generic_descartes_bfs_rep(*this); + } + + /*! + * \brief Computes a better approximation of the \c i th root of the + * polynomial + */ + virtual void refine_interval(int i) const { + CGAL_assertion(i>=0 && i < number_of_intervals); + Node_iterator curr = bitstream_tree.begin(), begin, end, + new_begin, helper; + std::advance(curr,i); + int intervals=1; + end=curr; + ++end; + begin=curr; + do { + //AcX_DSTREAM(bitstream_tree.lower(begin) << " " << bitstream_tree.upper(begin) << std::endl); + //AcX_DSTREAM(bitstream_tree.min_var(begin) << " " << bitstream_tree.max_var(begin) << std::endl); + int new_intervals = bitstream_tree.subdivide(begin,new_begin,helper); + intervals+=new_intervals-1; + begin=new_begin; + curr=helper; + + // Fixes the bug when a interval splits, and the leftmost subintervals + // has no children with sign variation >=1 + if(intervals==1) { + break; + } + if(new_intervals==0) { + continue; + } + + while(curr!=end) { + intervals+=bitstream_tree.subdivide(curr,new_begin,helper)-1; + curr=helper; + } + } + while(intervals!=1); + //AcX_DSTREAM("Refined " << left_boundary(i) << " " << right_boundary(i) << std::endl); + + } + + /*! + * \brief isolates the root of \c f + * + * The mechanism is the following: The \c bitstream_tree member of the + * object is transformed via subdivision until the + * \c termination_condition routine of the object returns true. When this + * happens, the \c process_nodes routine of the object is called. + */ + virtual void isolate() { + + //AcX_DSTREAM("Starting isolation" << std::endl); + + Node_iterator curr = bitstream_tree.begin(),dummy,new_curr; + + if(curr==bitstream_tree.end()) { + is_isolated_=true; + return; + } + + int newly_created; + + while(! this->termination_condition()) { + + + if(curr==bitstream_tree.end()) { + curr=bitstream_tree.begin(); + } + + if(bitstream_tree.max_var(curr)==1) { + ++curr; + } + else { + //AcX_DSTREAM("Subdivision at " + //<< CGAL::to_double(bitstream_tree.lower(curr)) << " " + //<< CGAL::to_double(bitstream_tree.upper(curr)) << std::flush); + newly_created = bitstream_tree.subdivide(curr,dummy,new_curr); + number_of_intervals+=newly_created-1; + curr=new_curr; + //AcX_DSTREAM("done" << std::endl); + } + + } + this->process_nodes(); + is_isolated_ = true; + } + + /*! + * \brief Computes an interval containing all real roots of \c p, + * using the Fujiwara root bound. + * + * So far, the \c log_div variable is always set to zero, this means + * that [lower,upper] is the interval containing all real roots + */ + virtual void get_interval(const Polynomial& p, Integer& lower, + Integer& upper, long& log_div, + Bitstream_descartes_rndl_tree_traits traits) { + + + typename Bitstream_descartes_rndl_tree_traits::Lower_bound_log2_abs + lower_bound_log2_abs = traits.lower_bound_log2_abs_object(); + typename + Bitstream_descartes_rndl_tree_traits::Upper_bound_log2_abs_approximator + upper_bound_log2_abs_approximator + = traits.upper_bound_log2_abs_approximator_object(); + //AcX_DSTREAM("Fujiwara bound.." << p << std::endl); +#if BITSTREAM_USES_E08_TREE + log_div = -NiX::Fujiwara_root_bound_log(p.begin(), + p.end(), + lower_bound_log2_abs, + upper_bound_log2_abs_approximator + ); +#else + + log_div = -CGAL::CGALi + ::Fujiwara_root_bound_log(p.begin(), + p.end(), + lower_bound_log2_abs, + upper_bound_log2_abs_approximator + ); +#endif + + //AcX_DSTREAM("Fujiwara returns " << log_div << std::endl); + // To be sure + log_div--; + lower=Integer(-1); + upper=Integer(1); + return; + + /* + + typename Bitstream_descartes_rndl_tree_traits::Approximator + approx = traits.approximator_object(); + long approx_power = 1; + Integer denom = CGAL::abs(approx(p.lcoeff(),approx_power)); + while(CGAL::compare(denom,Integer(1))!=CGAL::POSITIVE) { + approx_power*=2; + denom = CGAL::abs(approx(p.lcoeff(),approx_power)); + } + denom-=Integer(1); + CGAL_assertion(CGAL::compare(denom,Integer(0))==CGAL::POSITIVE); + Integer num(0); + + for(typename Polynomial::const_iterator curr=p.begin();curr!=p.end();curr++) { + Integer curr_approx = CGAL::abs(approx(*curr,approx_power))+Integer(1); + //AcX_DSTREAM("Coeffs aprox: " << curr_approx << std::endl); + num+=curr_approx; + } + //AcX_DSTREAM("Num: " << num << " Denom: " << denom << " Div: " << int_div(num,denom)+Integer(1) << "log_div: " << log_div << std::endl); + upper=CGAL::div(num,denom)+Integer(1); + lower=-upper; + log_div=0; + return; + */ + } + + //! returns the number of detected isolating intervals + virtual int number_of_real_roots() const { + return number_of_intervals; + } + + //! The lower boundary of the \c i th root + virtual Boundary left_boundary(int i) const { + CGAL_assertion(i>=0 && i < number_of_intervals); + Node_const_iterator curr = bitstream_tree.begin(); + std::advance(curr,i); + return bitstream_tree.lower(curr); + } + + //! The upper boundary of the \c i th root + virtual Boundary right_boundary(int i) const { + CGAL_assertion(i>=0 && i < number_of_intervals); + Node_const_iterator curr = bitstream_tree.begin(); + std::advance(curr,i); + return bitstream_tree.upper(curr); + } + + //! Returns the polynomial which is isolated + Polynomial polynomial() const { + return f_; + } + + /*! + * \brief When does the isolation algorithm terminate? + * + * This method must be specialised by derived classes + */ + virtual bool termination_condition() { + throw Virtual_method_exception(); + return false; + } + + /*! + * \brief Gives an opportunity to process the nodes after + * the subdivision steps are finished + * + * This method must be specialised by derived classes, but can + * remain empty in many cases. + */ + virtual void process_nodes() { + throw Virtual_method_exception(); + return; + } + + /*! \brief Returns whether the \c i th root is definitely a simple root + * of the isolated polynomial + * + * Must be specialised by derived class + */ + virtual bool is_certainly_simple_root(int i) const { + throw Virtual_method_exception(); + return false; + } + + /*! \brief Returns whether the \c i th root is definitely a multiple root + * of the isolated polynomial + * + * Must be specialised by derived class + */ + virtual bool is_certainly_multiple_root(int i) const { + throw Virtual_method_exception(); + return false; + } + + + virtual int multiplicity_of_root(int i) const { + CGAL_assertion(i>=0 && i < number_of_intervals); + return -1; + } + + virtual int get_upper_bound_for_multiplicity(int i) const { + CGAL_assertion(i>=0 && i < number_of_intervals); + Node_const_iterator curr = bitstream_tree.begin(); + std::advance(curr,i); + return bitstream_tree.min_var(curr); + } + + //! Must be specialized by the derived class + virtual int degree_of_gcd() const { + throw Virtual_method_exception(); + return -1; + } + + //! Must be specialized by the derived class + virtual Polynomial square_free_part() const { + throw Virtual_method_exception(); + return Polynomial(); + } + + //! Must be specialized by the derived class + virtual Handle inverse_transform_isolator() const { + throw Virtual_method_exception(); + return Handle(); + } + + bool is_isolated() const { + return is_isolated_; + } + + Bitstream_descartes_rndl_tree_traits traits() const { + return traits_; + } + + Bitstream_tree get_tree() const { + + return bitstream_tree; + + } + + //! type to distinguish used constructor + Bitstream_descartes_bfs_type type_; + + protected: + + //! Polynomial which is isolated + Polynomial f_; + + //! The traits class + Bitstream_descartes_rndl_tree_traits traits_; + + //! The tree of the Bitstream Descartes method + mutable Bitstream_tree bitstream_tree; + + //! The number of detected isolating intervals + int number_of_intervals; + + //! Has isolation already taken place + mutable bool is_isolated_; + + }; + +/* + * \brief Representation for square free polynomials + */ + template + class Square_free_descartes_bfs_rep + : public Generic_descartes_bfs_rep { + + + public: + + //! Traits type + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + //! The generic representation + typedef Generic_descartes_bfs_rep Base; + + //! Polynomial type + typedef typename Base::Polynomial Polynomial; + + //! Iterator for the leaves in the bitstream tree + typedef typename Base::Node_iterator Node_iterator; + + //! The type of the tree that controls the Bitstream instance + typedef typename Base::Bitstream_tree Bitstream_tree; + + /*! + * \brief Constructor with the square free polynomial f. + */ + Square_free_descartes_bfs_rep( + Polynomial f, + Bitstream_descartes_rndl_tree_traits traits) : + Base(SQUARE_FREE_DESCARTES, f,traits) { + } + + /*! + * \brief Constructor with the square free polynomial f. + */ + Square_free_descartes_bfs_rep( + Polynomial f, + Bitstream_tree tree, + Bitstream_descartes_rndl_tree_traits traits) : + Base(SQUARE_FREE_DESCARTES, f, tree, traits) { + } + + + + //! Needed for reference counting + virtual CGAL::Reference_counted_hierarchy<>* clone() { + return new Square_free_descartes_bfs_rep(*this); + } + + /*! + * \brief Terminates when all detected roots are simple + */ + virtual bool termination_condition() { + for(Node_iterator curr=Base::bitstream_tree.begin(); + curr != Base::bitstream_tree.end();curr++) { + if(Base::bitstream_tree.max_var(curr)!=1) { + return false; + } + } + return true; + } + + //! nothing to do here + virtual void process_nodes() { + return; + } + + //! Polynomial is square free, so gcd is 1 + virtual int degree_of_gcd() const { + return 0; + } + + //! Polynomial is square free + virtual Polynomial square_free_part() const { + return this->f_; + } + + //! Always true + virtual bool is_certainly_simple_root(int i) const { + return true; + } + + //! Always false + virtual bool is_certainly_multiple_root(int i) const { + return false; + } + + }; + +/* + * \brief Representation for polynomials with at most one multiple root + */ + template + class M_k_descartes_bfs_rep + : public Generic_descartes_bfs_rep { + + + public: + + //! Traits class + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + //! Generic representation + typedef Generic_descartes_bfs_rep Base; + + //! Polynomial type + typedef typename Base::Polynomial Polynomial; + + //! Iterator for the leaves of the Bitstream Descartes tree + typedef typename Base::Node_iterator Node_iterator; + + //! Constant iterator for the leaves + typedef typename Base::Node_const_iterator Node_const_iterator; + + //! The interval boundaries are represented in this type + typedef typename Bitstream_descartes_rndl_tree_traits::Boundary + Boundary; + + //! The type of the tree that controls the Bitstream instance + typedef typename Base::Bitstream_tree Bitstream_tree; + + /*! + * \brief Constructor for a polynomial f, not necessarily square + * free + * + * The values m + * and k need to be the exact number of real roots of f + * counted without multiplicity and the degree of the greatest common + * divisor of f with its partial derivative, respectively. + */ + M_k_descartes_bfs_rep(Polynomial f,int m, int k, + Bitstream_descartes_rndl_tree_traits traits) : + Base(M_K_DESCARTES, f,traits), + number_of_roots(m), + gcd_degree(k), + index_of_multiple(-1) { + } + + M_k_descartes_bfs_rep(Polynomial f,int m, int k, + Bitstream_tree tree, + Bitstream_descartes_rndl_tree_traits traits) : + Base(M_K_DESCARTES, f, tree, traits), + number_of_roots(m), + gcd_degree(k), + index_of_multiple(-1) { + } + + //! Default constructor + M_k_descartes_bfs_rep() { + } + + //! Needed for reference counting + virtual CGAL::Reference_counted_hierarchy<>* clone() { + return new M_k_descartes_bfs_rep(*this); + } + + /*! + * \brief Termination condition + * + * If m-1 simple and one more leaf is detected, the Bitstream + * Descartes method is stopped. If the minimal sign + * variation drops under k in each leaf, a + * \c Non_generic_position_exception is thrown. + */ + virtual bool termination_condition() { + int counted_simple_roots=0; + int max_max_var = 0; + for(Node_iterator curr=Base::bitstream_tree.begin(); + curr != Base::bitstream_tree.end();curr++) { + int max_var=Base::bitstream_tree.max_var(curr); + if(max_var > max_max_var) { + max_max_var=max_var; + } + if(max_var==1) { // && Base::bitstream_tree.max_var(curr)==1) { + ++counted_simple_roots; + } + } + //AcX_DSTREAM("Situation: " << this->number_of_intervals << " intervals " << this->number_of_roots << " are expected" << std::endl); + if(this->number_of_intervals==this->number_of_roots + && counted_simple_roots>=number_of_roots-1) { + return true; + } + if(max_max_var<=gcd_degree) { + throw CGAL::CGALi::Non_generic_position_exception(); + } + + return false; + + } + + //! The index of the (possibly) multiple root is computed here. + virtual void process_nodes() { + int i=0; + for(Node_iterator curr=Base::bitstream_tree.begin(); + curr != Base::bitstream_tree.end();curr++) { + if(Base::bitstream_tree.max_var(curr) > 1 ) { + index_of_multiple=i; + return; + } + else { + ++i; + } + } + return; + } + + //! Returns k + virtual int degree_of_gcd() const { + return gcd_degree; + } + + //! True for all roots except for the candidate + virtual bool is_certainly_simple_root(int i) const { + return (i!=index_of_multiple); + } + + //! Always false + virtual bool is_certainly_multiple_root(int i) const { + return false; + } + + + protected: + + //! The "m" + int number_of_roots; + + //! The "k" + int gcd_degree; + + //! The candidate's index + int index_of_multiple; + + }; + + +/* + * \brief Representation for polynomials + * with square free polynomial known as well + */ + template + class Exchange_descartes_bfs_rep + : public Generic_descartes_bfs_rep { + + + public: + + //! Traits class + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + //! Generic representation + typedef Generic_descartes_bfs_rep Base; + + //! Polynomial type + typedef typename Base::Polynomial Polynomial; + + //! Iterator for the leaves of the Bitstream Descartes tree + typedef typename Base::Node_iterator Node_iterator; + + //! Constant iterator for the leaves + typedef typename Base::Node_const_iterator Node_const_iterator; + + //! The interval boundaries are represented in this type + typedef typename Bitstream_descartes_rndl_tree_traits::Boundary + Boundary; + + //! The type of the tree that controls the Bitstream instance + typedef typename Base::Bitstream_tree Bitstream_tree; + + /*! + * \brief Constructor for a polynomial f, not necessarily square + * free + * + * \e sq_free_f is the square free part of the polynomial f + */ + Exchange_descartes_bfs_rep(Polynomial f,Polynomial sq_free_f, + Bitstream_descartes_rndl_tree_traits traits) + : Base(EXCHANGE_DESCARTES,f,traits),sq_free_f(sq_free_f), + sq_free_descartes(Square_free_descartes_tag(),sq_free_f,traits) + { + } + + Exchange_descartes_bfs_rep(Polynomial f,Polynomial sq_free_f, + Bitstream_tree tree, + Bitstream_descartes_rndl_tree_traits traits) + : Base(EXCHANGE_DESCARTES,f,tree,traits),sq_free_f(sq_free_f), + sq_free_descartes(Square_free_descartes_tag(),sq_free_f,traits) + { + } + + //! Default constructor + Exchange_descartes_bfs_rep() + {} + + //! Needed for reference counting + virtual CGAL::Reference_counted_hierarchy<>* clone() { + return new Exchange_descartes_bfs_rep(*this); + } + + virtual void isolate() { + + // TODO: More clever isolation + + //AcX_DSTREAM("Starting isolation" << std::endl); + + Node_iterator curr = this->bitstream_tree.begin(),dummy,new_curr; + + if(curr==this->bitstream_tree.end()) { + this->is_isolated_=true; + return; + } + + int newly_created; + + while(! this->termination_condition()) { + + + if(curr==this->bitstream_tree.end()) { + curr=this->bitstream_tree.begin(); + } + + newly_created = this->bitstream_tree.subdivide(curr,dummy,new_curr); + this->number_of_intervals+=newly_created-1; + curr=new_curr; + //AcX_DSTREAM("done" << std::endl); + } + + this->process_nodes(); + this->is_isolated_ = true; + } + + + /*! + * \brief Termination condition + * + */ + virtual bool termination_condition() { + + int n = this->number_of_intervals; + + if( n != sq_free_descartes.number_of_real_roots() ) { + return false; + } + + Node_iterator curr = this->bitstream_tree.begin(); + + for( int i = 0; i < n-1; i++ ) { + if(this->bitstream_tree.upper(curr) > + sq_free_descartes.left_boundary(i+1)) { + return false; + } + curr++; + } + + curr = this->bitstream_tree.begin(); + curr++; + + for( int i = 1; i < n; i++ ) { + if(this->bitstream_tree.lower(curr) < + sq_free_descartes.right_boundary(i-1)) { + return false; + } + curr++; + } + + return true; + + } + + //! The index of the (possibly) multiple root is computed here. + virtual void process_nodes() { + + } + + //! Polynomial is square free, so gcd is 1 + virtual int degree_of_gcd() const { + return this->f_.degree() - sq_free_f.degree(); + } + + //! Polynomial is square free + virtual Polynomial square_free_part() const { + return sq_free_f; + } + + //! True for all with sign variation one + virtual bool is_certainly_simple_root(int i) const { + CGAL_assertion(i>=0 && i < this->number_of_intervals); + Node_const_iterator curr = this->bitstream_tree.begin(); + std::advance(curr,i); + return this->bitstream_tree.max_var(curr) == 1; + } + + //! True at least for roots with even sign variation + virtual bool is_certainly_multiple_root(int i) const { + CGAL_assertion(i>=0 && i < this->number_of_intervals); + Node_const_iterator curr = this->bitstream_tree.begin(); + std::advance(curr,i); + return this->bitstream_tree.max_var(curr) % 2 == 0; + } + + + protected: + + //! The square free part of the polynomial + Polynomial sq_free_f; + + //! The Isolator for the square free part + Bitstream_descartes_bfs + sq_free_descartes; + + + }; + + +/* + * \brief Representation for polynomials with several multiple roots + */ + template + class M_ki_descartes_bfs_rep : public + Generic_descartes_bfs_rep { + + public: + + //! Traits class + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + //! Generic representation + typedef Generic_descartes_bfs_rep Base; + + //! Polynomial type + typedef typename Base::Polynomial Polynomial; + + //! Iterator for the leaves of the Bitstream Descartes tree + typedef typename Base::Node_iterator Node_iterator; + + //! Constant iterator for the leaves + typedef typename Base::Node_const_iterator Node_const_iterator; + + //! The Bitstream trees are represented in this type + typedef typename Base::Bitstream_tree Bitstream_tree; + + //! The integers are represented in this type + typedef typename Bitstream_descartes_rndl_tree_traits::Integer + Integer; + + //! The interval boundaries are represented in this type + typedef typename Bitstream_descartes_rndl_tree_traits::Boundary + Boundary; + + /*! + * \brief Constructor for a polynomial f, not necessarily square + * free + */ + M_ki_descartes_bfs_rep(Polynomial f, Polynomial sqf, + Bitstream_descartes_rndl_tree_traits traits) : + Base(M_KI_DESCARTES,sqf,traits), + _m_traits(traits), + _m_f_square_free(sqf) { + } + + //! Default constructor + M_ki_descartes_bfs_rep() { + } + + + //! Needed for reference counting + virtual CGAL::Reference_counted_hierarchy<>* clone() { + return new M_ki_descartes_bfs_rep(*this); + } + + /*! + * \brief Termination condition + * + * Each node of the tree either represents a single root + * of the polynomial, or + * \f[ \deg(f/\gcd(f,f')) \geq \deg(f) + 1 - min_var(node) \f] and + * the boundaries of the node form a single root in \f[ f/\gcd(f,f') \f] + */ + virtual bool termination_condition() { + + //std::cout << "check termination" << std::endl; + + const int fd = this->f_.degree(); + const int sqfd = this->_m_f_square_free.degree(); + + if (sqfd == 1) { + return true; + } + + for (Node_iterator curr=Base::bitstream_tree.begin(); + curr != Base::bitstream_tree.end(); curr++) { + + if (Base::bitstream_tree.min_var(curr) == 0) { + continue; + } + if (Base::bitstream_tree.max_var(curr) <= 1) { + continue; + } + + // else for higher "multiplicities" we check +#if 0 + std::cout << "fd: " << fd << std::endl; + std::cout << "sfd: " << sqfd << std::endl; + std::cout << "min: " << Base::bitstream_tree.min_var(curr) + << std::endl; + std::cout << "max: " << Base::bitstream_tree.max_var(curr) + << std::endl; +#endif + // TODO test max or min? + if (sqfd < fd && + Base::bitstream_tree.min_var(curr) > sqfd) { + //std::cout << "false 1" << std::endl; + return false; + } + + if (sqfd <= fd + 1 - Base::bitstream_tree.min_var(curr)) { + + typename Base::Integer clow, cupp; + long log_bdry_den; + + Base::bitstream_tree.boundaries( + curr, clow, cupp, log_bdry_den + ); + + typename Base::Bitstream_tree::Monomial_basis_tag tag; + + typedef typename Base::Bitstream_tree Tree; + + Tree square_free_tree = + Tree(clow, cupp, + log_bdry_den, + this->_m_f_square_free.begin(), + this->_m_f_square_free.end(), + tag, + this->_m_traits); + + Node_iterator begin = square_free_tree.begin(); + + if (begin == square_free_tree.end()) { + //std::cout << "false 2" << std::endl; + return false; + } else { + if (square_free_tree.min_var(begin) < 1 || + square_free_tree.max_var(begin) > 1) { + //std::cout << "false 3" << std::endl; + return false; + } + } + } else { + //std::cout << "false 4" << std::endl; + return false; + } + } + +#if !NDEBUG + std::cerr << "MKI-Descartes: Termination condition not sufficient" + << std::endl; +#endif + return true; + } + + //! empty + virtual void process_nodes() { + // empty + return; + } + + //! Polynomial is square free, so gcd is 1 + virtual int degree_of_gcd() const { + return _m_f_degree - _m_f_square_free_degree; + } + + //! Polynomial is square free + virtual Polynomial square_free_part() const { + return this->_m_f_square_free; + } + + virtual bool is_certainly_simple_root(int i) const { + Node_iterator curr = Base::bitstream_tree.begin(); + std::advance(curr, i); + return (Base::bitstream_tree.min_var(curr) == 1); + } + + virtual bool is_certainly_multiple_root(int i) const { + Node_iterator curr = Base::bitstream_tree.begin(); + std::advance(curr, i); + return (Base::bitstream_tree.min_var(curr) > 1); + } + + virtual int multiplicity_of_root(int i) const { + Node_iterator curr = Base::bitstream_tree.begin(); + std::advance(curr, i); + return (Base::bitstream_tree.min_var(curr)); + } + + protected: + + Bitstream_descartes_rndl_tree_traits _m_traits; + + int _m_f_degree; + + Polynomial _m_f_square_free; + + int _m_f_square_free_degree; + }; + + template + class Backshear_descartes_bfs_rep + : public Generic_descartes_bfs_rep { + + + public: + + typedef EventRefinement Event_refinement; + + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + typedef Generic_descartes_bfs_rep Base; + + typedef typename Base::Polynomial Polynomial; + + typedef typename Base::Node_iterator Node_iterator; + + typedef std::list::iterator Marking_iterator; + + typedef std::list::const_iterator Marking_const_iterator; + + typedef typename Base::Node_const_iterator Node_const_iterator; + + typedef typename Base::Boundary Boundary; + + Backshear_descartes_bfs_rep( + Polynomial f, + int number_of_non_event_points, + int number_of_events, + Event_refinement event_refinement, + Bitstream_descartes_rndl_tree_traits traits) : + Base(BACKSHEAR_DESCARTES,f,traits), + number_of_non_event_points(number_of_non_event_points), + number_of_events(number_of_events), + event_refinement(event_refinement) { + } + + Backshear_descartes_bfs_rep() { + } + + virtual CGAL::Reference_counted_hierarchy<>* clone() { + return new Backshear_descartes_bfs_rep(*this); + } + + virtual void isolate() { + + Node_iterator curr = Base::bitstream_tree.begin(),sub_begin,new_curr; + + if(curr==Base::bitstream_tree.end()) { + this->is_isolated_ = true; + return; + } + markings.clear(); + markings.push_back(this->check_marking(curr)); + + Marking_iterator curr_mark=markings.begin(),mark_helper; + + int newly_created; + + while(! this->termination_condition()) { + //AcX_DSTREAM("Subdivision..." << number_of_intervals << std::endl); + if(curr==Base::bitstream_tree.end()) { + curr=Base::bitstream_tree.begin(); + CGAL_assertion(curr_mark==markings.end()); + curr_mark=markings.begin(); + } + if(Base::bitstream_tree.max_var(curr)==1) { + ++curr; + ++curr_mark; + //AcX_DSTREAM("nothing happend" << std::endl); + } + else { + newly_created = + Base::bitstream_tree.subdivide(curr,sub_begin,new_curr); + mark_helper=markings.erase(curr_mark); + curr_mark=mark_helper; + for(Node_iterator tmp_curr=sub_begin;tmp_curr!=new_curr;tmp_curr++) { + markings.insert(curr_mark,check_marking(tmp_curr)); + } + Base::number_of_intervals+=newly_created-1; + curr=new_curr; + //AcX_DSTREAM(newly_created << " new intervals, marking size: " << markings.size() << std::endl); + + } + } + this->process_nodes(); + this->is_isolated_ = true; + } + + + + virtual bool termination_condition() { + int marked_intervals=0; + int unmarked_odd_intervals=0; + Node_iterator curr=Base::bitstream_tree.begin(); + Marking_iterator curr_mark=markings.begin(); + for(;curr!=Base::bitstream_tree.end();curr++) { + if((*curr_mark) >=0) { + ++marked_intervals; + } + else { + if(Base::bitstream_tree.min_var(curr)%2==1) { // odd + ++unmarked_odd_intervals; + } + } + ++curr_mark; + } + CGAL_assertion(curr_mark==markings.end()); + return ((marked_intervals==number_of_events) + && (unmarked_odd_intervals==number_of_non_event_points)); + } + + virtual void process_nodes() { + Node_iterator curr=Base::bitstream_tree.begin(),curr_helper; + Marking_iterator curr_mark=markings.begin(); + while(curr!=Base::bitstream_tree.end()) { + if(((*curr_mark)==-1) && (Base::bitstream_tree.min_var(curr)%2==0)) { + ++curr; + curr_helper=curr; + curr_helper--; + Base::bitstream_tree.erase(curr_helper); + curr_mark=markings.erase(curr_mark); + Base::number_of_intervals--; + } else { + ++curr_mark; + ++curr; + } + } + CGAL_assertion(curr_mark==markings.end()); + + //AcX_DSTREAM(markings.size() << " " << number_of_non_event_points << " " << number_of_events << std::endl); + CGAL_assertion(static_cast(markings.size()) + ==number_of_non_event_points+number_of_events); + return; + } + + virtual bool is_certainly_simple_root(int i) const { + CGAL_assertion(i>=0 && i < Base::number_of_intervals); + Node_const_iterator curr=Base::bitstream_tree.begin(); + std::advance(curr,i); + return (Base::bitstream_tree.max_var(curr)==1); + } + + virtual bool is_certainly_multiple_root(int i) const { + CGAL_assertion(i>=0 && i < Base::number_of_intervals); + Marking_const_iterator curr=markings.begin(); + std::advance(curr,i); + return (*curr>=0); + } + + + protected: + + int number_of_non_event_points; + + int number_of_events; + + Event_refinement event_refinement; + + std::list markings; + + protected: + + int check_marking(Node_iterator node) { + Boundary lower = Base::bitstream_tree.lower(node), + upper=Base::bitstream_tree.upper(node); + for(int i=0;i + class Inverse_transform_descartes_bfs_rep + : public Generic_descartes_bfs_rep { + + public: + + //! The traits class for approximations + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + //! The Coeeficient type of the input polynomial + typedef typename Bitstream_descartes_rndl_tree_traits::Coefficient + Coefficient; + + //! The polynomial type + typedef CGAL::Polynomial Polynomial; + + typedef Inverse_transform_descartes_bfs_rep + Self; + + //! The used integer type + typedef typename Bitstream_descartes_rndl_tree_traits::Integer Integer; + + //! How the boundaries of the isolating intervals are represented + typedef typename Bitstream_descartes_rndl_tree_traits::Boundary + Boundary; + + typedef Generic_descartes_bfs_rep + Base; + + //! The type for the inverse isolator + typedef typename Base::Handle Handle; + + + + + /*! + * \brief Constructor + */ + Inverse_transform_descartes_bfs_rep(const Polynomial& f, + Handle inv_descartes, + Boundary translation) : + Base(INVERSE_TRANSFORM_DESCARTES), + q(translation), + inv_descartes(inv_descartes) + { + CGAL_assertion(inv_descartes.is_isolated()); + this->is_isolated_ = true; + this->f_ = f; + this->traits_ = inv_descartes.traits(); + this->number_of_intervals + = inv_descartes.number_of_real_roots() - 1; + number_of_roots_smaller_zero = 0; + while(inv_descartes.right_boundary + (number_of_roots_smaller_zero) + < 0) { + number_of_roots_smaller_zero++; + } + + } + + //! Destructor (does nothing) + virtual ~Inverse_transform_descartes_bfs_rep() { + } + + //! Needed for the referencing counting mechanism + virtual CGAL::Reference_counted_hierarchy<>* clone() { + return new Inverse_transform_descartes_bfs_rep(*this); + } + + virtual void refine_interval(int i) const { + inv_descartes.refine_interval(transform_index(i)); + } + + virtual void isolate() const { + } + + + //! The lower boundary of the \c i th root + virtual Boundary left_boundary(int i) const { + return 1/inv_descartes.right_boundary(transform_index(i)) + q; + } + + //! The upper boundary of the \c i th root + virtual Boundary right_boundary(int i) const { + return 1/inv_descartes.left_boundary(transform_index(i)) + q; + } + + + /*! \brief Returns whether the \c i th root is definitely a simple root + * of the isolated polynomial + * + * Must be specialised by derived class + */ + virtual bool is_certainly_simple_root(int i) const { + return inv_descartes.is_certainly_simple_root(transform_index(i)); + } + + /*! \brief Returns whether the \c i th root is definitely a multiple root + * of the isolated polynomial + * + * Must be specialised by derived class + */ + virtual bool is_certainly_multiple_root(int i) const { + return inv_descartes.is_certainly_multiple_root(transform_index(i)); + } + + virtual Handle inverse_transform_isolator() const { + + return inv_descartes; + + } + + + protected: + + //! Translation factor q + Boundary q; + + //! The isolator for f(1/x + q)) + Handle inv_descartes; + + //! Roots <=q of inv_descartes + int number_of_roots_smaller_zero; + + int transform_index(int i) const { + CGAL_assertion(i>=0 && i < this->number_of_intervals); + int ret_value; + if ( i < number_of_roots_smaller_zero ) { + ret_value = number_of_roots_smaller_zero - i - 1; + } else { + ret_value = + this->number_of_intervals - + (i - number_of_roots_smaller_zero); + } + CGAL_assertion(ret_value >= 0 && i < this->number_of_intervals); + return ret_value; + } + + }; + +/* + * \brief Adaptor for roots of a vert line + * (needed as dummy in surface analysis) + * + */ + template + class Vert_line_adapter_descartes_bfs_rep + : public Generic_descartes_bfs_rep { + + public: + + //! The traits class for approximations + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + //! type of vert line + typedef VertLine Vert_line; + + //! type of Curve_analysis_2 + typedef typename Vert_line::Curve_analysis_2 Curve_analysis_2; + + //! type of Curve_kernel_2; + typedef typename Curve_analysis_2::Algebraic_curve_kernel_2 + Curve_kernel_2; + + //! The Coeeficient type of the input polynomial + typedef typename Bitstream_descartes_rndl_tree_traits::Coefficient + Coefficient; + + //! The polynomial type + typedef CGAL::Polynomial Polynomial; + + typedef Inverse_transform_descartes_bfs_rep + Self; + + //! The used integer type + typedef typename Bitstream_descartes_rndl_tree_traits::Integer Integer; + + //! How the boundaries of the isolating intervals are represented + typedef typename Bitstream_descartes_rndl_tree_traits::Boundary + Boundary; + + typedef Generic_descartes_bfs_rep + Base; + + //! The type for the inverse isolator + typedef typename Base::Handle Handle; + + /*! + * \brief Constructor + */ + template + Vert_line_adapter_descartes_bfs_rep(InputIterator begin, + InputIterator end, + Bitstream_descartes_rndl_tree_traits traits) + : Base(VERT_LINE_ADAPTER_DESCARTES) + { + std::copy( begin, end, std::back_inserter(root_vec) ); + this->is_isolated_ = true; + this->traits_ = traits; + this->f_ = Polynomial(0); + this->number_of_intervals + = static_cast(root_vec.size()); + // Isolate all real roots until intervals are disjoint: + for(int i = 1; i < this->number_of_real_roots(); i++ ){ + while( left_boundary(i) < right_boundary(i-1) ) { + if( right_boundary(i)-left_boundary(i) < + right_boundary(i-1)-left_boundary(i-1) ) { + refine_interval(i-1); + } else { + refine_interval(i); + } + } + + } + } + + + //! Destructor (does nothing) + virtual ~Vert_line_adapter_descartes_bfs_rep() { + + } + + //! Needed for the referencing counting mechanism + virtual CGAL::Reference_counted_hierarchy<>* clone() { + return new Vert_line_adapter_descartes_bfs_rep(*this); + } + + virtual void refine_interval(int i) const { + typename Curve_kernel_2::Refine_y_2 + refine_y; // TODO call _object + refine_y(root_vec[i].first.algebraic_real_2(root_vec[i].second)); + } + + virtual void isolate() const { + } + + + //! The lower boundary of the \c i th root + virtual Boundary left_boundary(int i) const { + typename Curve_kernel_2::Lower_boundary_y_2 + lower_boundary_y; // TODO call _object + return lower_boundary_y( + root_vec[i].first.algebraic_real_2(root_vec[i].second) + ); + } + + //! The upper boundary of the \c i th root + virtual Boundary right_boundary(int i) const { + typename Curve_kernel_2::Upper_boundary_y_2 + upper_boundary_y; // TODO call _object + return upper_boundary_y( + root_vec[i].first.algebraic_real_2(root_vec[i].second) + ); + } + + + /*! \brief Returns whether the \c i th root is definitely a simple root + * of the isolated polynomial + * + */ + virtual bool is_certainly_simple_root(int i) const { + return false; + } + + /*! \brief Returns whether the \c i th root is definitely + * a multiple root + * of the isolated polynomial + * + */ + virtual bool is_certainly_multiple_root(int i) const { + return false; + } + + protected: + + //! Roots stored as pair of a AcX::Vert_line and an integer denoting the + //! index + std::vector > root_vec; + + }; + + /*! + * \brief Class for the Bitstream Descartes method + * + * Class for the real root isolation of polynomials, using the Bitstream + * Descartes method. The polynomials coefficient type is arbitrary, the + * approximations of the coefficient type are obtained with the + * \c BitstreamDescartesRndlTreeTraits parameter. For the requirements + * of this traits class, see the documentation of + * CGAL::Bitstream_descartes_rndl_tree. + * + * Internally, an instance of CGAL::Bitstream_descartes_rndl_tree is explored + * in a specific way. That exploration strategy depends on the constructor + * that is used to create the object. A tag is passed that defines the + * variant of the Bitstream Descartes method: The Square_free_descartes_tag + * starts the usual Bitstream method for square free integer polynomials. + * With the M_k_descartes tag, it is able to handle one multiple root in + * favourable situations, the Backshear_descartes_tag allows to isolate + * even more complicated polynomials, if the multiple roots with even + * multiplicity can be refined from outside. See the corresponding + * constructors for more information. + * + */ + template + class Bitstream_descartes_bfs + : ::CGAL::Handle_with_policy< + CGAL::CGALi::Generic_descartes_bfs_rep > { + + public: + + //! Traits class + typedef BitstreamDescartesRndlTreeTraits + Bitstream_descartes_rndl_tree_traits; + + // The generic representation class + typedef + CGAL::CGALi::Generic_descartes_bfs_rep Rep; + + // The Handle type + typedef ::CGAL::Handle_with_policy Base; + + //! The coefficients of the polynomial + typedef typename Bitstream_descartes_rndl_tree_traits::Coefficient + Coefficient; + + //! The polynomial's type + typedef CGAL::Polynomial Polynomial; + + typedef Bitstream_descartes_bfs + Self; + + // Type for the Bitstream Descartes tree +#if BITSTREAM_USES_E08_TREE + typedef NiX::Bitstream_descartes_E08_tree + + Bitstream_tree; +#else + typedef CGAL::CGALi::Bitstream_descartes_rndl_tree + + Bitstream_tree; +#endif + + //! Type for Integers + typedef typename Bitstream_descartes_rndl_tree_traits::Integer Integer; + + //! Iterator type for the leaves of the Descartes tree + typedef typename Bitstream_tree::Node_iterator + Node_iterator; + + //! Const iterator for the leaves + typedef typename Bitstream_tree::Node_const_iterator + Node_const_iterator; + + //! Type for the interval boundaries of the isolating intervals + typedef typename Bitstream_descartes_rndl_tree_traits::Boundary + Boundary; + + //! Default constructor + Bitstream_descartes_bfs() : Base(new Rep()) {} + + /*! + * \brief Constructor for the square free Descartes method + * + * The polynomial \c f must not have multiple real roots. The + * Bitstream Descartes tree is traversed in a bfs manner until + * all leaves have sign variation zero or one. + */ + Bitstream_descartes_bfs(Square_free_descartes_tag t, + Polynomial f, + Bitstream_descartes_rndl_tree_traits traits + = Bitstream_descartes_rndl_tree_traits(), + bool isolate=true) + : Base(new CGAL::CGALi::Square_free_descartes_bfs_rep + (f,traits)) + { + if(isolate) { + this->isolate(); + } + } + + /*! + * \brief Constructor for the square free Descartes method, + * using a precomputed tree + * + * The polynomial \c f must not have multiple real roots. The + * Bitstream Descartes tree is traversed in a bfs manner until + * all leaves have sign variation zero or one. + * The tree must be adequate for the polynomial. + * Use that constructor only if you know what you're doing! + */ + Bitstream_descartes_bfs(Square_free_descartes_tag t, + Polynomial f, + Bitstream_tree tree, + Bitstream_descartes_rndl_tree_traits traits + = Bitstream_descartes_rndl_tree_traits(), + bool isolate=true) + : Base(new CGAL::CGALi::Square_free_descartes_bfs_rep + (f, tree, traits)) + { + if(isolate) { + this->isolate(); + } + } + + /*! + * \brief Constructor for the m-k-Descartes method + * + * The polynomial \c f must have exactly \c m real roots, counted without + * multiplicity, and the degree of gcd(f,f') must be \c k. In this + * case, the constructor either isolates the real roots of \c f sucessfully + * or a Non_generic_position_exception is thrown. Such an exception + * certainly occurs if \c f has more than one multiple real root. If \c f + * has at most one multiple root over the complex numbers, the roots are + * certainly isolated with success. + */ + Bitstream_descartes_bfs(M_k_descartes_tag t, + Polynomial f,int m,int k, + Bitstream_descartes_rndl_tree_traits traits + = Bitstream_descartes_rndl_tree_traits(), + bool isolate=true) + : Base(new CGAL::CGALi::M_k_descartes_bfs_rep + (f,m,k,traits)) + { + if(isolate) { + this->isolate(); + } + } + + + Bitstream_descartes_bfs(M_k_descartes_tag t, + Polynomial f,int m,int k, + Bitstream_tree tree, + Bitstream_descartes_rndl_tree_traits traits + = Bitstream_descartes_rndl_tree_traits(), + bool isolate=true) + : Base(new CGAL::CGALi::M_k_descartes_bfs_rep + (f,m,k,tree,traits)) + { + if(isolate) { + this->isolate(); + } + } + + + + + /*! + * \brief Constructor for the Exchange-Descartes method + * + * The polynomial \c f can have several multiple roots. + * The value of \c must represent the its square-free counterpart. + */ + Bitstream_descartes_bfs(Exchange_descartes_tag t, + Polynomial f, + Polynomial sqf, + Bitstream_descartes_rndl_tree_traits traits = + Bitstream_descartes_rndl_tree_traits(), + bool isolate = true) + : Base(new CGAL::CGALi::Exchange_descartes_bfs_rep + (f,sqf,traits)) + { + if(isolate) { + this->isolate(); + } + } + + Bitstream_descartes_bfs(Exchange_descartes_tag t, + Polynomial f, + Polynomial sqf, + Bitstream_tree tree, + Bitstream_descartes_rndl_tree_traits traits = + Bitstream_descartes_rndl_tree_traits(), + bool isolate = true) + : Base(new CGAL::CGALi::Exchange_descartes_bfs_rep + (f,sqf,tree,traits)) + { + if(isolate) { + this->isolate(); + } + } + + + + /*! + * \brief Constructor for the m-ki-Descartes method + * + * The polynomial \c f can have several multiple roots. + * The value of \c must represent the its square-free counterpart. + */ + Bitstream_descartes_bfs(M_ki_descartes_tag t, + Polynomial f, + Polynomial sqf, + Bitstream_descartes_rndl_tree_traits traits = + Bitstream_descartes_rndl_tree_traits(), + bool isolate = true) + : Base(new CGAL::CGALi::M_ki_descartes_bfs_rep + (f,sqf,traits)) + { + if(isolate) { + this->isolate(); + } + } + + + /*! + * \brief Constructor for the Backshear-Decartes method + * + * The polynomial \c f must have exactly \c number_of_real_roots + * many real roots, counted without multiplicity. Additionally, a set of + * \c number_of_events root can be refined to arbitrary precision with the + * \c event_refinement object. This must support three operations + * for each 0<=i: + *
  • lower_boundary(i), upper_boundary(i) gives an interval (not + * necessarily isolating) of some root of \c f
    • + *
  • refine(i) refines the corresponding interval
+ * Note that the roots in \c event_refinement need not be sorted. All roots + * which are not covered by \c event_refinement must have odd multiplicity. + */ + template + Bitstream_descartes_bfs(Backshear_descartes_tag t, + Polynomial f, + int number_of_real_roots, + int number_of_events, + EventRefinement event_refinement, + Bitstream_descartes_rndl_tree_traits traits + = Bitstream_descartes_rndl_tree_traits(), + bool isolate=true) + : Base(new + CGAL::CGALi::Backshear_descartes_bfs_rep + + (f,number_of_real_roots-number_of_events, + number_of_events,event_refinement,traits)) + { + if(isolate) { + this->isolate(); + } + } + + + /*! + * \brief Constructor for the Inverse-Transformation-Descartes method + * + */ + Bitstream_descartes_bfs(Inverse_transform_descartes_tag t, + Polynomial f, + Self inv_descartes, + Boundary translation = Boundary(0)) + : Base(new CGAL::CGALi::Inverse_transform_descartes_bfs_rep + (f,inv_descartes, + translation)) + { + // No isolation necessary + } + + /*! + * \brief Constructor for the Vert-line-adapter-Descartes method + * + */ + template + Bitstream_descartes_bfs(Vert_line_adapter_descartes_tag t, + InputIterator begin, + InputIterator end, + Bitstream_descartes_rndl_tree_traits traits) + : Base(new CGAL::CGALi::Vert_line_adapter_descartes_bfs_rep + + (begin, end, traits) ) + { + // No isolation necessary + } + + //! return the type of the used descartes method + Bitstream_descartes_bfs_type type() const { + return this->ptr()->type_; + } + + //! Return the polynomial + Polynomial polynomial() const { + CGAL_assertion(is_isolated()); + return this->ptr()->polynomial(); + } + + //! Returns the traits class + Bitstream_descartes_rndl_tree_traits traits() const { + return this->ptr()->traits(); + } + + //! Number of real roots of the polynomial + int number_of_real_roots() const { + CGAL_assertion(is_isolated()); + return this->ptr()->number_of_real_roots(); + } + + //! Refine the ith isolating interval + void refine_interval(int i) const { + CGAL_assertion(is_isolated()); + this->ptr()->refine_interval(i); + } + + //! The left boundary of the ith isolating interval + Boundary left_boundary(int i) const { + CGAL_assertion(is_isolated()); + return this->ptr()->left_boundary(i); + } + + //! The right boundary of the ith isolating interval + Boundary right_boundary(int i) const { + CGAL_assertion(is_isolated()); + return this->ptr()->right_boundary(i); + } + + //! The length of the ith isolating interval + Boundary length(int i) const { + CGAL_assertion(is_isolated()); + return (this->ptr()->right_boundary(i) - + this->ptr()->left_boundary(i)); + } + + /*! + * \brief Returns true if the ith root is known to be a simple + * root of the curve. + */ + bool is_certainly_simple_root(int i) const { + CGAL_assertion(is_isolated()); + return this->ptr()->is_certainly_simple_root(i); + } + + /*! + * \brief Returns true if the ith root is known to be a multiple + * root of the curve. + */ + bool is_certainly_multiple_root(int i) const { + CGAL_assertion(is_isolated()); + return this->ptr()->is_certainly_multiple_root(i); + } + + + /*! + * \brief Returns the multiplicity of the root if know, otherwise -1 + */ + int multiplicity_of_root(int i) const { + CGAL_assertion(is_isolated()); + return this->ptr()->multiplicity_of_root(i); + } + + /*! + * Returns an upper bound for the multiplicity of the ith root + */ + int get_upper_bound_for_multiplicity(int i) const { + CGAL_assertion(is_isolated()); + return this->ptr()->get_upper_bound_for_multiplicity(i); + } + + /*! + * \brief Returns the isolator of the polynomial f(1/x + q), if known + */ + Self inverse_transform_isolator() const { + return this->ptr()->inverse_transform_isolator(); + } + + + public: + + //! Starts the isolation of the real roots. + void isolate() { + CGAL_assertion(! is_isolated()); + this->ptr()->isolate(); + } + + bool is_isolated() const { + return this->ptr()->is_isolated(); + } + + Bitstream_tree get_tree() const { + + return this->ptr()->get_tree(); + + } + + //! returns the degree of the gcd of f and its derivative, if known + int degree_of_gcd() const { + return this->ptr()->degree_of_gcd(); + } + + //! returns the square free part of f, if known + Polynomial square_free_part() const { + return this->ptr()->square_free_part(); + } + + + }; + + + + } // namespace CGALi + +CGAL_END_NAMESPACE + +#endif diff --git a/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_traits_on_vert_line.h b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_traits_on_vert_line.h new file mode 100644 index 00000000000..2b073a58445 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_traits_on_vert_line.h @@ -0,0 +1,687 @@ +// TODO: Add licence +// +// 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) : Michael Kerber +// +// ============================================================================ + +#ifndef CGAL_BITSTREAM_DESCARTES_TRAITS_ON_VERT_LINE +#define CGAL_BITSTREAM_DESCARTES_TRAITS_ON_VERT_LINE 1 + +#include +#include +#include +#include +#include + +CGAL_BEGIN_NAMESPACE + + // Workaround for min/max. This is just temporary! + // Necessary to prevent compiler confusion about std::min/max + // and CGAL::min/max. The better solution is to put the CGAL functions into + // some seperate namespace and include it into CGAL-namepsace with the + // using command. The workaround works only for Sqrt_extensions +template< class NT, class Root > +inline Sqrt_extension min BOOST_PREVENT_MACRO_SUBSTITUTION +( const Sqrt_extension& x , + const Sqrt_extension& y) { + return Min >() (x,y); +} + +template< class NT, class Root > +inline Sqrt_extension max BOOST_PREVENT_MACRO_SUBSTITUTION +( const Sqrt_extension& x , + const Sqrt_extension& y) { + return Min >() (x,y); +} + +namespace CGALi { + + +// Representation of the Bitstream Descartes traits class +template +class Bitstream_descartes_traits_on_vert_line_rep { + +public: + + typedef AlgebraicReal Algebraic_real; + typedef Integer_ Integer; + typedef Coefficient_ Coefficient; + + // Coefficient is a polynomial type + typedef typename Coefficient::NT Coeff_NT; + + typedef typename Algebraic_real::Rational Boundary; + + typedef typename Algebraic_real::Rational Rational; + + typedef typename CGAL::Coercion_traits::Type + Coercion; + + typedef typename CGAL::Coercion_traits::Cast + Rational_to_coercion_cast; + + typedef typename + CGAL::Get_arithmetic_kernel::Arithmetic_kernel:: + Bigfloat_interval BFI; + + typedef typename + CGAL::Get_arithmetic_kernel::Arithmetic_kernel:: + Bigfloat BF; + + typedef CGAL::CGALi::Best_approximation_cache + Best_approximation_cache; + + typedef Bitstream_descartes_traits_on_vert_line_rep + Self; + + Bitstream_descartes_traits_on_vert_line_rep + (Algebraic_real alpha,bool use_approx_memory) + : alpha_(alpha), use_approx_memory(use_approx_memory) + { + } + + + /* + Bitstream_descartes_traits_on_vert_line_rep + (const Self& s) + : alpha_(s.alpha) + { + } + */ + + bool uses_approx() const { + return use_approx_memory; + } + + Best_approximation_cache approx_mem() const { + return approximations_; + } + + + Algebraic_real alpha() const { + return alpha_; + } + + class Approximator { + + private: + + Algebraic_real alpha_; + + Best_approximation_cache approximations_; + + bool use_approx_mem; + + long curr_prec; + + public: + Approximator(Algebraic_real alpha) : alpha_(alpha),use_approx_mem(true), curr_prec(4) {}; + + Approximator(Algebraic_real alpha, + Best_approximation_cache approximations) + : alpha_(alpha), approximations_(approximations), use_approx_mem(true), curr_prec(4) {}; + + Approximator() {}; + + template + class Coeff_to_bfi_functor { + + public: + + typedef NT argument_type; + typedef BFI result_type; + + BFI operator() (NT x) const { + + return CGAL::convert_to_bfi(x); + + } + + }; + + + typedef boost::numeric::interval Interval; + + long log_width(Algebraic_real alpha) const { + CGAL_assertion(! alpha.is_rational()); + double d = CGAL::to_double(alpha.high()- alpha.low()); + if(d != 0.0) { + return CGAL::CGALi::ceil_log2_abs(d); + } + return 0; + } + +#if AcX_USE_NO_BFI_APPROX_IN_BITSTREAM_TRAITS +#warning uses no bfi-approx! + +Integer operator() (Coefficient f, long p) { + + //AcX_DSTREAM("Called approximator with " << f << " and " << p << std::endl); + //AcX_DSTREAM("Alpha: " << (CGAL::to_double(alpha_.high()-alpha_.low())) << " " << rational_to_double(alpha_.high()) << std::endl); + + + + if(use_approx_mem) { + if(approximations_.approximation_exists(f)) { + long prec; + Integer approx; + approximations_.get_best_approximation(f,prec,approx); + if(prec>=p) { + //AcX_DSTREAM("use older approx.." << std::flush); + Integer divisor = CGAL::ipower((Integer)2,prec-p); + Integer return_value = CGAL::div(approx,divisor); + //AcX_DSTREAM("done" << std::endl); + return return_value; + } + } + } + + //AcX_DSTREAM("Have to approximate..." << std::endl); + + Coercion bound; + + Rational_to_coercion_cast rational_to_coercion_cast; + + if(p<0) { + bound = rational_to_coercion_cast + (Rational(CGAL::ipower((Integer)2,-p))); + } + else { + typename CGAL::Fraction_traits::Compose compose; + bound = + rational_to_coercion_cast + (compose(1,CGAL::ipower(Integer(2),p))); + } + Interval f_alpha_iv; + + + + if(alpha_.is_rational()) { + Coercion f_alpha=f.evaluate(alpha_.low()); + f_alpha_iv=Interval(f_alpha,f_alpha); + + } + else { + typedef Approximate_arithmetic_controller + + Approximate_controller; + + Approximate_controller approx_controller(f,alpha_); + + while(true) { + + f_alpha_iv = approx_controller.interval_approximation(); + + //AcX_DSTREAM("done" << std::endl); + if(CGAL::compare + (f_alpha_iv.upper()-f_alpha_iv.lower(),bound) == + CGAL::SMALLER) { + break; + } + //AcX_DSTREAM("refine..." << std::flush); + + approx_controller.refine_value(); + + //AcX_DSTREAM("done" << std::endl); + } + } + + + Integer return_value; + typename CGAL::Fraction_traits::Decompose decompose; + typedef typename + CGAL::Fraction_traits::Numerator_type Numerator; + typedef typename + CGAL::Fraction_traits::Denominator_type Denominator; + if(sign(f_alpha_iv.lower())==CGAL::ZERO) { + return_value = Integer(0); + } + else if(CGAL::sign(f_alpha_iv.lower()) != + CGAL::sign(f_alpha_iv.upper())) { + return_value = Integer(0); + } + else { + Numerator num; + Denominator denom; + + decompose(f_alpha_iv.upper(),num,denom); + + + if(p<0) { + denom *= CGAL::ipower(Integer(2),-p); + } + else { + num *= CGAL::ipower(Integer(2),p); + } + return_value = CGAL::div(num,denom); + + } + + if(use_approx_mem) { + approximations_.update_approximation(f,p,return_value); + } + //AcX_DSTREAM(" returning " << " " << return_value << std::endl); + + return return_value; + } + + +#else + + CGAL::Polynomial _convert_polynomial_to_bfi(Coefficient f) const { + std::vector coeffs; + for(int i = 0; i <= f.degree(); i++) { + coeffs.push_back(CGAL::convert_to_bfi(f[i])); + } + return CGAL::Polynomial(coeffs.begin(), coeffs.end()); + } + + + Integer operator() (Coefficient f, long p) { + + typename CGAL::CGALi::Float_traits::Get_exponent get_exp; + typename CGAL::CGALi::Float_traits::Get_mantissa get_m; + + long old_prec = CGAL::get_precision(BFI()); + + CGAL::Polynomial f_bfi; + BFI alpha_bfi, f_alpha_bfi; + + long prec = 16; + + Integer return_value; + + long wbit = 0; + + while(true) { + CGAL::set_precision(BFI(),prec); + + f_bfi = _convert_polynomial_to_bfi(f); + alpha_bfi = CGAL::convert_to_bfi(alpha_); + + f_alpha_bfi = f_bfi.evaluate(alpha_bfi); + + if(!CGAL::singleton(f_alpha_bfi)) { + long ceil = CGAL::CGALi::ceil_log2_abs(f_alpha_bfi); + long signi = CGAL::get_significant_bits(f_alpha_bfi); + wbit = ceil - signi + p; + + } + + if(wbit<-5 || CGAL::singleton(f_alpha_bfi)) { + break; + } else { + prec*=2; + } + } + BF lower = CGAL::lower(f_alpha_bfi); + + long shift = - (p + get_exp(lower)); + Integer bfi_m(get_m(lower)); + if( shift > 0 ){ + while(shift>63) { + bfi_m = (bfi_m >> 63); + shift-=63; + } + bfi_m = (bfi_m >> shift); + }else{ + // add 0 bits + bfi_m = (bfi_m << -shift); + } + return_value = bfi_m; + CGAL::set_precision(BFI(),old_prec); + + return return_value; + } +#endif + + }; + + + class Lower_bound_log2_abs { + + private: + Algebraic_real alpha_; + + public: + Lower_bound_log2_abs(Algebraic_real alpha) : alpha_(alpha) {} + Lower_bound_log2_abs() {}; + + typedef boost::numeric::interval Interval; + + long operator() (Coefficient f) { + //std::cout << "Called lower_bound_log2_abs with " + // << f << std::flush; + + CGAL_assertion(! alpha_.is_root_of(f)); + + Interval f_alpha_iv; + Coercion abs_lower,abs_upper; + + typedef CGAL::CGALi::Approximate_arithmetic_controller + Approximate_controller; + + Approximate_controller approx_controller(f,alpha_); + + while(true) { + f_alpha_iv = approx_controller.interval_approximation(); + CGAL::Sign lower_sign = CGAL::sign(f_alpha_iv.lower()); + if(CGAL::sign(f_alpha_iv.upper())==lower_sign) { + if(lower_sign==CGAL::POSITIVE) { + abs_lower=f_alpha_iv.lower(); + abs_upper=f_alpha_iv.upper(); + } + else { + abs_lower=CGAL::abs(f_alpha_iv.upper()); + abs_upper=CGAL::abs(f_alpha_iv.lower()); + } + + BFI bfi_low = CGAL::convert_to_bfi(abs_lower), + bfi_high = CGAL::convert_to_bfi(abs_upper); + long lower_bound = CGAL::CGALi::floor_log2_abs(bfi_low), + upper_bound = CGAL::CGALi::ceil_log2_abs(bfi_high); + + CGAL_assertion(upper_bound>=lower_bound); + if(upper_bound-lower_bound <=2) { + //std::cout << "returning " << lower_bound << std::endl; + return lower_bound; + } + } + approx_controller.refine_value(); + } + } + + }; + + class Upper_bound_log2_abs_approximator { + + private: + Algebraic_real alpha_; + + std::vector coefficients_for_alpha; + + int refinements_of_alpha; + + bool zero_test_enabled; + + int refinement_limit; + + // Stores id of polynomials which are known to vanish (or not to + // vanish) at alpha + std::vector zeroes,non_zeroes; + + public: + Upper_bound_log2_abs_approximator(Algebraic_real alpha) + : alpha_(alpha),refinements_of_alpha(0),zero_test_enabled(false) + {} + + Upper_bound_log2_abs_approximator() {}; + + typedef boost::numeric::interval Boundary_interval; + typedef boost::numeric::interval Interval; + + bool initial_upper_bound + (Coefficient f, long& ub_log2_abs,bool& is_certainly_zero) { + return improve_upper_bound(f,ub_log2_abs,is_certainly_zero); + } + + bool improve_upper_bound + (const Coefficient f, long& ub_log2_abs,bool& is_certainly_zero) { + //AcX_DSTREAM("improve upper bound.." << f << std::flush); + if(std::find(coefficients_for_alpha.begin(), + coefficients_for_alpha.end(), + f.id())!=coefficients_for_alpha.end()) { + coefficients_for_alpha.clear(); + //AcX_DSTREAM("refine.." << std::flush); + alpha_.bisect(); + ++refinements_of_alpha; + if(refinements_of_alpha>=refinement_limit) { + zero_test_enabled=true; + } + //AcX_DSTREAM("done.." << std::flush); + } + coefficients_for_alpha.push_back(f.id()); + if(zero_test_enabled) { + if(std::find(zeroes.begin(), + zeroes.end(), + f.id())!=zeroes.end()) { + is_certainly_zero=true; + return true; + } + else if(std::find(non_zeroes.begin(), + non_zeroes.end(), + f.id())!=non_zeroes.end()) { + is_certainly_zero=false; + } + else { + bool zero = CGAL::CGALi::is_root_of(alpha_,f); + if(zero) { + zeroes.push_back(f.id()); + is_certainly_zero=true; + return true; + } + else { + non_zeroes.push_back(f.id()); + is_certainly_zero=false; + } + } + } + else { + is_certainly_zero=false; + } + Boundary_interval alpha_iv(alpha_.low(),alpha_.high()); + Interval f_alpha_iv = evaluate_iv(f,alpha_iv); + CGAL::Sign lower_sign = CGAL::sign(f_alpha_iv.lower()); + bool iv_contains_zero = lower_sign + != CGAL::sign(f_alpha_iv.upper()); + Coercion abs_upper + =(CGAL::abs(f_alpha_iv.lower())= lb_log2_abs); + //AcX_DSTREAM("Upper: " << ub_log2_abs << " Lower: " << lb_log2_abs << std::endl); + //AcX_DSTREAM(((ub_log2_abs - lb_log2_abs) <= 2) << std::endl); + + return ((ub_log2_abs - lb_log2_abs) <= 2); + } + else { + //AcX_DSTREAM("Upper: " << ub_log2_abs << std::endl); + return false; + } + } + }; + +private: + + Algebraic_real alpha_; + + Best_approximation_cache approximations_; + + bool use_approx_memory; + +}; + +/*! + * \brief Traits class for the Bitstream Descartes method over algebraic + * extensions. + * + * Approximates coefficients of the polynomial + * g=f|x=a + * using interval arithmetic. The \c AlgebraicReal object \c a is fixed + * for the traits class. Each coefficient is represented by an univariate + * integer polynomial. For a more detailled specification of its member + * functions, see the documentation of + * \c CGAL::CGALi::Bitstream_descartes_rndl_tree + * + * Using this traits class, bivariate polynomials can be isolated with + * AcX::Bitstream_descartes_bfs, the inner variable \c x is interpreted as + * the algebraic number \c a. + */ +template +class Bitstream_descartes_traits_on_vert_line + : public ::CGAL::Handle_with_policy + > { + +public: + + //! Coefficient type + typedef Coefficient_ Coefficient; + + //! Algebraic real type + typedef AlgebraicReal Algebraic_real; + + //! Integer type + typedef Integer_ Integer; + + typedef Bitstream_descartes_traits_on_vert_line + Self; + + typedef CGAL::CGALi::Bitstream_descartes_traits_on_vert_line_rep + Rep; + + typedef ::CGAL::Handle_with_policy Handle; + + typedef typename Rep::Boundary Boundary; + + //! Functor type for the approximation of coefficients + typedef typename Rep::Approximator Approximator; + + //! Functor type for lower bounds on coefficients + typedef typename Rep::Lower_bound_log2_abs Lower_bound_log2_abs; + + //! Functor type for upper bounds on coefficients + typedef typename Rep::Upper_bound_log2_abs_approximator + Upper_bound_log2_abs_approximator; + + typedef typename Rep::Rational Rational; + + // We can just take these from NT_traits + typedef typename CGAL::Real_embeddable_traits::Sign + Sign; + typedef typename CGAL::CGALi::Real_embeddable_extension + ::Ceil_log2_abs Ceil_log2_abs_Integer; + typedef typename CGAL::CGALi::Real_embeddable_extension + ::Ceil_log2_abs Ceil_log2_abs_long; + + + /*! + * \brief Constructor for the polynomial f|x=a + * + * Initialises an instance such that all bivariate polynomials are + * interpreted as univariate polynomials in the algebraic extension + * Z[a]. + * If the \c use_approx_mem flag is set, the best known approximation + * of each coefficient is stored in a + * \c CGAL::CGALi::Intern::Best_approximation_cache instance. + */ + Bitstream_descartes_traits_on_vert_line + (AlgebraicReal alpha=AlgebraicReal(),bool use_approx_mem = true) + : Handle(alpha,use_approx_mem) + { + } + + /* + Bitstream_descartes_traits_on_vert_line + (const Self& s) + : Handle(s) + { + } + */ + + Algebraic_real alpha() const { + return this->ptr()->alpha(); + } + + //! See the documentation of CGAL::CGALi::Bitstream_descartes_rndl_tree. + Approximator approximator_object() const { + if(this->ptr()->uses_approx()) { + return Approximator(this->ptr()->alpha(),this->ptr()->approx_mem()); + } + else { + return Approximator(this->ptr()->alpha()); + } + } + + //! See the documentation of CGAL::CGALi::Bitstream_descartes_rndl_tree. + Upper_bound_log2_abs_approximator + upper_bound_log2_abs_approximator_object() const { + return Upper_bound_log2_abs_approximator(this->ptr()->alpha()); + } + + //! See the documentation of CGAL::CGALi::Bitstream_descartes_rndl_tree. + Lower_bound_log2_abs lower_bound_log2_abs_object() const { + return Lower_bound_log2_abs(this->ptr()->alpha()); + } + + + //! See the documentation of CGAL::CGALi::Bitstream_descartes_rndl_tree. + class Boundary_creator { + + public: + + Boundary_creator() {} + + Boundary operator() (Integer x,long p) { + Integer num=x, denom,two(2),q,r; + if(p <0) { + CGAL::div_mod(num,two,q,r); + while(r==Integer(0) && p<0) { + num=q; + p++; + CGAL::div_mod(num,two,q,r); + } + denom = CGAL::ipower(Integer(2),-p); + } + else { + num*=CGAL::ipower(Integer(2),p); + denom=1; + } + Rational b(num,denom); + CGAL::simplify(b); + return b; + } + }; + +}; + +} // namespace CGALi + +CGAL_END_NAMESPACE + +#endif diff --git a/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Non_generic_position_exception.h b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Non_generic_position_exception.h new file mode 100644 index 00000000000..77d17d82c81 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Non_generic_position_exception.h @@ -0,0 +1,44 @@ +// TODO: Add licence +// +// 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) : Michael Kerber +// +// ============================================================================ + + + +#ifndef CGAL_NON_GENERIC_POSITION_EXCEPTION +#define CGAL_NON_GENERIC_POSITION_EXCEPTION + + +CGAL_BEGIN_NAMESPACE + + namespace CGALi { + + /*! + * \brief Exception class for not sufficiently generic positions. + * + * Must be thrown whenever a curve cannot be analysed because its position + * is not "good enough". + */ + class Non_generic_position_exception { + + public: + + //! Default constructible + Non_generic_position_exception() {} + + }; + + } // namespace CGALi + +CGAL_END_NAMESPACE + + +#endif diff --git a/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/alg_real_utils.h b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/alg_real_utils.h new file mode 100644 index 00000000000..36fa6bf0d8e --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/alg_real_utils.h @@ -0,0 +1,691 @@ +// TODO: Add licence +// +// 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) : Michael Kerber +// +// ============================================================================ + + +#ifndef CGAL_ALG_REAL_UTILS +#define CGAL_ALG_REAL_UTILS 1 + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + +/* + * \brief Function for merging two sets + * + * This function is similar to the \c std::union_set operation. + * Additionally, it provides a sequence of AcX::Three_valued + * providing information to which input set the corresponding root + * in the merged sequence belonged + * + * The BinaryFunction must have the Result type CGAL::Comparison_result. + */ +template +std::pair +set_union_with_source(InputIterator1 first_begin, + InputIterator1 first_end, + InputIterator2 second_begin, + InputIterator2 second_end, + OutputIterator1 merged_values, + OutputIterator2 merged_values_info, + BinaryFunction compare) { + + InputIterator1 first_it=first_begin; + InputIterator2 second_it=second_begin; + + while((first_it != first_end) || (second_it!=second_end)) { + if(first_it == first_end) { + *merged_values=*second_it++; + ++merged_values; + *merged_values_info++ = CGAL::CGALi::ROOT_OF_SECOND_SET; + continue; + } + if(second_it == second_end) { + *merged_values++=*first_it++; + *merged_values_info++ = CGAL::CGALi::ROOT_OF_FIRST_SET; + continue; + } + + CGAL::Comparison_result c = compare(*first_it,*second_it); + + if(c==CGAL::EQUAL) { + *merged_values++=*first_it++; + ++second_it; + *merged_values_info++ = CGAL::CGALi::ROOT_OF_BOTH_SETS; + continue; + } + if(c==CGAL::SMALLER) { + *merged_values++=*first_it++; + *merged_values_info++ = CGAL::CGALi::ROOT_OF_FIRST_SET; + continue; + } + if(c==CGAL::LARGER) { + *merged_values++=*second_it++; + *merged_values_info++ = CGAL::CGALi::ROOT_OF_SECOND_SET; + continue; + } + } + return std::make_pair(merged_values,merged_values_info); +} + + + +/*! \brief Tries to find a SIMPLE rational q with aq + * is a power of two, and is not too big. There is no guarantee to find + * the rational value between a and b of minimal + * bit size. + */ +template +typename Algebraic_real::Rational +simple_rational_between(const Algebraic_real& a, + const Algebraic_real&b) { + + //srb.start(); + typedef typename Algebraic_real::Rational Rational; + typename CGAL::Fraction_traits::Compose compose; + typedef typename + CGAL::Get_arithmetic_kernel::Arithmetic_kernel AK; + typedef typename AK::Bigfloat Bigfloat; + typedef typename AK::Bigfloat_interval Bigfloat_interval; + typedef typename AK::Integer Integer; + + long old_prec = CGAL::get_precision(Bigfloat_interval()); + + CGAL_assertion(a= b.low()) { + Rational size_a=a.high()-a.low(), + size_b=b.high() - b.low(); + while(a.high() >= b.low()) { + if(size_a < size_b) { + b.refine(); + size_b=b.high() - b.low(); + } else { + a.refine(); + size_a=a.high()-a.low(); + } + } + } + //srb_a.stop(); + + //srb_b.start(); + Bigfloat x=CGAL::upper(CGAL::convert_to_bfi(a.high())), + y=CGAL::lower(CGAL::convert_to_bfi(b.low())); + if(x>=y) { + Rational size_a=a.high() - a.low(), + size_b=b.high() - b.low(), + size_max = size_a>size_b ? size_a : size_b, + size_int = b.low()-a.high(); + while(x>=y) { + //std::cout << "x and y: " << x << " and " << y << std::endl; + //std::cout << "sizes: " << CGAL::to_double(size_int) << " " << CGAL::to_double(size_max) << std::endl; + if(size_int>size_max) { + CGAL::set_precision(Bigfloat_interval(), + 2*CGAL::get_precision(Bigfloat_interval())); + x=CGAL::upper(CGAL::convert_to_bfi(a.high())); + y=CGAL::lower(CGAL::convert_to_bfi(b.low())); + } else { + if(size_a < size_b) { + b.refine(); + size_b=b.high() - b.low(); + y=CGAL::lower(CGAL::convert_to_bfi(b.low())); + } else { + a.refine(); + size_a=a.high()-a.low(); + x=CGAL::upper(CGAL::convert_to_bfi(a.high())); + } + size_max = size_a>size_b ? size_a : size_b; + size_int = b.low()-a.high(); + } + } + } + //srb_b.stop(); + //std::cout << "Intermediate2: " << x << " " << y << std::endl; + typename CGAL::CGALi::Float_traits::Get_mantissa mantissa; + typename CGAL::CGALi::Float_traits::Get_exponent exponent; + + Integer x_m = mantissa(x), + y_m=mantissa(y); + long x_e = exponent(x), y_e = exponent(y); + //std::cout << "Floats1: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl; + if((x_m > 0 && y_m < 0) || x_m < 0 && y_m > 0) { + //srb.stop(); + return Rational(0); + } + bool negative=false; + if(x_m<=0 && y_m <=0) { + x_m=-x_m; + y_m=-y_m; + std::swap(x_m,y_m); + std::swap(x_e,y_e); + negative=true; + } + // Now, we have that (x_m,x_e) represents a number smaller than (y_m,y_e) + //srb_c.start(); + //std::cout << "Floats2: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl; + + // As long as the mantissa is even, simplify + while(x_m != 0 && (x_m & 1)==0 ) { + x_m=x_m >> 1; + x_e++; + } + while(y_m != 0 && (y_m & 1)==0 ) { + y_m=y_m >> 1; + y_e++; + } + //srb_c.stop(); + //std::cout << "Floats3: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl; + + // Bring both numbers to a common exponent + //srb_d.start(); + long min_e = x_e < y_e ? x_e : y_e; + while(x_e > min_e) { + x_m=x_m << 1; + x_e--; + } + while(y_e > min_e) { + y_m=y_m << 1; + y_e--; + } + //srb_d.stop(); + CGAL_assertion(y_e==x_e && x_e==min_e); + + //std::cout << "Floats4: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl; + + // Avoid mantissas to have difference one + if(y_m-x_m==Integer(1)) { + x_m=x_m << 1; + y_m=y_m << 1; + x_e--; + y_e--; + min_e--; + } + //std::cout << "Floats5: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl; + Integer final_mantissa(0); + //srb_e.start(); + long x_log = x_m==Integer(0) ? -1 : CGAL::CGALi::floor_log2_abs(x_m), + y_log = y_m==Integer(0) ? -1 : CGAL::CGALi::floor_log2_abs(y_m), + old_log = y_log; + while(x_log==y_log) { + //std::cout << "here" << std::endl; + while(old_log > y_log) { + final_mantissa = final_mantissa << 1; + old_log--; + } + CGAL_assertion((x_m & ((Integer(1) << x_log) - 1)) == x_m - CGAL::ipower(Integer(2),x_log)); + x_m = x_m & ((Integer(1) << x_log) - 1); // x_m - CGAL::ipower(Integer(2),x_log); + y_m = y_m & ((Integer(1) << y_log) - 1); // y_m - CGAL::ipower(Integer(2),y_log); + + final_mantissa++; + old_log=y_log; + x_log = x_m==0 ? -1 : CGAL::CGALi::floor_log2_abs(x_m); + y_log = y_m==0 ? -1 : CGAL::CGALi::floor_log2_abs(y_m); + } + //srb_e.stop(); + // Now, x_log != y_log, in fact, y_log is greater + CGAL_assertion(x_log y_log) { + final_mantissa = final_mantissa << 1; + old_log--; + } + if((y_m & ((Integer(1) << y_log) - 1 ))==0) { // y_m - CGAL::ipower(Integer(2),y_log)==0) { + // Now, the constructed value would be equal to + while(y_log!=0 && x_log==y_log-1) { + final_mantissa = final_mantissa << 1; + final_mantissa++; + y_log--; + x_m = x_m==0 ? 0 : x_m & ((Integer(1) << x_log) - 1); //x_m - CGAL::ipower(Integer(2),x_log); + x_log = x_m==0 ? -1 : CGAL::CGALi::floor_log2_abs(x_m); + } + final_mantissa = final_mantissa << 1; + final_mantissa++; + y_log--; + } else { + final_mantissa++; + } + //srb_f.stop(); + min_e += y_log; + Rational rat_between; + //std::cout << "Min_e: " << min_e << std::endl; + if(min_e > 0) { + rat_between = compose(final_mantissa << min_e, + Integer(1)); + } else { + rat_between = compose(final_mantissa, Integer(1) << -min_e); + } + if(negative) { + rat_between = -rat_between; + } + //std::cout << "Result: " << a.high() << " " << rat_between << " " << b.low() << std::endl; + CGAL_assertion(a.high() < rat_between); + CGAL_assertion(b.low() > rat_between); + CGAL::set_precision(Bigfloat_interval(),old_prec); + //srb.stop(); + return rat_between; +} + + +/*! + * \brief Produces intermediate rational values for a list of + * algebraic reals. + * + * For a list of Algebraic real values with \c n elements, a list with + * n+1 elements of rational values is given such that the + * ith element is + * between the ith and the (i+1)th element of the input list + * + * The input list must be in increasing order + */ +template +OutputIterator find_intermediate_values(InputIterator start, + InputIterator end, + OutputIterator output) { + typedef typename std::iterator_traits::value_type + Alg_real; + typedef typename Alg_real::Rational Rational; + if(start==end) { + // Empty vector, create one element + *output=Rational(0); + ++output; + return output; + } + switch( CGAL::sign( start->low() ) ) { + case(CGAL::ZERO): { + *output = Rational(-1); + break; + } + case(CGAL::POSITIVE): { + *output = Rational(0); + break; + } + case(CGAL::NEGATIVE): { + Alg_real small_value = Alg_real(Rational(2)*start->low()); + *output + = simple_rational_between(small_value,*start); + } + } + output++; + + ++output; + InputIterator it_1(start),it_2(start); + ++it_2; + while(it_2 != end) { + CGAL_assertion(it_1->compare(*it_2)==CGAL::SMALLER); + Rational beta + = simple_rational_between(*it_1,*it_2); + *output=beta; + ++output; + ++it_1; + ++it_2; + } + switch( CGAL::sign( it_1->high() ) ) { + case(CGAL::ZERO): { + *output = Rational(1); + break; + } + case(CGAL::NEGATIVE): { + *output = Rational(0); + break; + } + case(CGAL::POSITIVE): { + Alg_real big_value = Alg_real(Rational(2)*it_1->high()); + *output + = simple_rational_between(*it_1,big_value); + } + } + output++; + return output; +} + +/*! + * \brief finds a Rational value left of an Algebraic real alpha + */ +template typename AlgebraicReal::Rational +simple_rational_left_of(AlgebraicReal ar) { + + typedef AlgebraicReal Algebraic_real; + typedef typename Algebraic_real::Rational Rational; + + switch( CGAL::sign( ar ) ) { + case(CGAL::ZERO): { + return Rational(-1); + break; + } + case(CGAL::POSITIVE): { + return Rational(0); + break; + } + case(CGAL::NEGATIVE): { + Algebraic_real small_value + = Algebraic_real(Rational(2)*ar.low()); + + return simple_rational_between(small_value,ar); + // = small_value.rational_between(ar); + //= ar.low()-1; + } + } + // never reached + return Rational(0); +} + +/*! + * \brief finds a Rational value rightt of an Algebraic real alpha + */ +template typename AlgebraicReal::Rational +simple_rational_right_of(AlgebraicReal ar) { + + typedef AlgebraicReal Algebraic_real; + typedef typename Algebraic_real::Rational Rational; + + return -simple_rational_left_of(-ar); + +} + +/*! + * \brief Tries to get the (non-zero) sign of f(alpha) with interval + * arithmetic. + * + * Performs refinements of alpha and evaluates \c f at \c alpha. If + * zero is still in the interval when the isolating interval of \c alpha + * is smaller then \c eps, CGAL::ZERO + * is returned, denoting that the sign could not be computed. Otherwise, + * the correct sign of \c f(alpha) is returned + * + * If \c n is set to zero (default), \c alpha is refined + * until the sign can be determined correctly. This may lead to an + * infinite loop, if f(alpha)=0 + * + */ +template +CGAL::Sign estimate_sign_at(AlgebraicReal alpha, + const Poly_1_& f, + const Boundary_ eps + =Boundary_(0)) { + typedef Boundary_ Boundary; + typedef Poly_1_ Poly_1; + typedef AlgebraicReal Algebraic_real; + typedef typename CGAL::Coercion_traits< + typename Poly_1::NT, + Boundary>::Type Interval_boundary; + typename CGAL::Coercion_traits::Cast boundary_cast; + typedef boost::numeric::interval Interval; + Interval ev; + CGAL::Sign sign=CGAL::ZERO; + + typedef CGAL::CGALi::Approximate_arithmetic_controller + Approximate_controller; + + Approximate_controller approx_controller(f,alpha); + + while(sign==CGAL::ZERO) { + ev = approx_controller.interval_approximation(); + + if(! boost::numeric::in_zero(ev)) { + sign=CGAL::compare(ev.lower(),0); + CGAL_assertion(sign!=CGAL::ZERO); + } + else { + if(CGAL::compare(alpha.high()-alpha.low(), + boundary_cast(eps))==CGAL::SMALLER) { + break; + } + approx_controller.refine_value(); + + } + } + return sign; +} + + +/*! + * \brief Evaluates poynomial at intervals with the Horner scheme + */ +template +boost::numeric::interval +::Type > +evaluate_iv(Poly1_ f,boost::numeric::interval iv) { + + typedef Boundary_ Boundary; + typedef Poly1_ Poly1; + typedef typename + CGAL::Coercion_traits::Type Coercion; + //typename CGAL::Coercion_traits::Cast cast; + CGAL_assertion(f.degree()>=0); +#if !AcX_USE_AFFINE_ARITHMETIC + typename CGAL::Coercion_traits::Cast cast; + typedef boost::numeric::interval Coercion_interval; + Coercion_interval iv_cast(cast(iv.lower()),cast(iv.upper())); + int n=f.degree(); + Coercion_interval ret(cast(f[n]),cast(f[n])); + for(int i=n-1;i>=0;i--) { + ret *= iv_cast; + ret += Coercion_interval(cast(f[i]),cast(f[i])); + } + return ret; +#else + typedef boost::numeric::interval Coercion_interval; + int n=f.degree(); + Boundary x0 = (iv.lower() + iv.upper()) / 2; + Boundary x1 = iv.upper() - x0; + + //! Now here four test-implementations of Affine Arithmetic + // AF1 is the first form, which should avoid the fact, that affine forms + // grow by doing non-affine operations + if(AcX_USE_AFFINE_ARITHMETIC == 1) + { + //! AF 1 + //std::cout << "Using Affine Arithmetic -> AF1" << std::endl; + Boundary y0 = f[n]; + Boundary y1 = 0, y2 = 0, spread = 0; + for(int i = n-1; i >= 0; i--) + { + y2 = abs(x0)*y2 + abs(x1)*(abs(y1)+abs(y2)); + y1 = x0*y1 + x1*y0; + y0 = x0*y0 + f[i]; + } + spread = abs(y1) + abs(y2); + return Coercion_interval(y0-spread, y0+spread); + } + // AF2 is the second form, which should avoid the fact described in AF1 + // and to decrease the difficulty of even powers + else if(AcX_USE_AFFINE_ARITHMETIC == 2) + { + //! AF 2 + //std::cout << "Using Affine Arithmetic -> AF2" << std::endl; + Boundary y0 = f[n]; + Boundary y1 = 0, k2 = 0, k3 = 0, e1 = 0, e2 = 0, e3 = 0, low = 0, high = 0, tmp = 0; + for(int i = n-1; i >= 0; i--) + { + if(x0 > 0) + { + k2 = x0 * e2; + k3 = x0 * e3; + } + else + { + k2 = -x0 * e3; + k3 = -x0 * e2; + } + tmp = x1 * y1; + if(tmp > 0) + { + k2 = k2 + tmp; + } + else + { + k3 = k3 + abs(tmp); + } + e1 = abs(x0)*e1 + abs(x1)*(e1 + e2 + e3); + e2 = k2; + e3 = k3; + y1 = x0*y1 + x1*y0; + y0 = x0*y0 + f[i]; + } + low = abs(y1) + e1 + e3; + high = abs(y1) + e1 + e2; + return Coercion_interval(y0-low, y0+high); + } + // QF (quadratic form) is the third form, which should avoid the facts of + // AF1 and decreases better the problem of even powers than AF2 + else if(AcX_USE_AFFINE_ARITHMETIC == 3) + { + //! QF + //std::cout << "Using Affine Arithmetic -> QF" << std::endl; + Boundary y0 = f[n]; + Boundary y1 = 0, z1 = 0, e1 = 0, spread = 0, low = 0, high = 0; + for(int i = n-1; i >= 0; i--) + { + e1 = abs(x0)*e1 + abs(x1)*e1 + abs(x1*z1); + z1 = x0*z1 + x1*y1; + y1 = x0*y1 + x1*y0; + y0 = x0*y0 + f[i]; + } + spread = abs(y1) + e1; + if(z1 > 0) + { + high = spread + z1; + } + else + { + low = spread - z1; + } + return Coercion_interval(y0-low, y0+high); + } + // Here we create a own version of the implementation for Affine Arithmetic + // We sort the powers to null, even and odd and compute the intervals + else if(AcX_USE_AFFINE_ARITHMETIC == 4) + { + //! AF own creation + //std::cout << "Using Affine Arithmetic -> own AF" << std::endl; + Boundary mid_point = (iv.lower() + iv.upper()) / 2; + Boundary width = iv.upper() - mid_point; + typedef CGAL::Polynomial Rat_poly_1; + typedef CGAL::Polynomial Coer_poly_1; + Rat_poly_1 affine_poly(mid_point, width); + Coer_poly_1 poly_aa = f.evaluate(affine_poly); + int n = poly_aa.degree(); + CGAL_assertion(n>=0); + Coercion low = poly_aa[0]; + Coercion high = poly_aa[0]; + for(int i = 1; i <= n ; i++) + { + CGAL::Sign sign = CGAL::sign(poly_aa[i]); + if(sign == CGAL::ZERO) + { + continue; + } + if(i%2 == 0) + { + if(sign == CGAL::POSITIVE) + { + high = high + poly_aa[i]; + } + else + { + low = low + poly_aa[i]; + } + } + else + { + if(sign == CGAL::POSITIVE) + { + low = low - poly_aa[i]; + high = high + poly_aa[i]; + } + else + { + low = low + poly_aa[i]; + high = high - poly_aa[i]; + } + } + } + return boost::numeric::interval(low, high); + } +#endif +} + +//! \brief replaces the \c NiX::is_root_of function, using \c AcX::ntl_gcd +template +bool is_root_of(const AlgebraicReal & x,Polynomial p) { + if(x.is_rational()) { + typedef typename AlgebraicReal::Rational Rational; + Rational exact_value = x.rational(); + typedef typename CGAL::Coercion_traits::Type + Type; + return CGAL::sign(p.evaluate(exact_value))==CGAL::ZERO; + } + else { + CGAL_precondition(p.degree()>=0); + + if(p.degree()==0) return p.is_zero(); + Polynomial g=CGAL::CGALi::gcd_utcf(p,x.polynomial()); + return g.sign_at(x.low())!=g.sign_at(x.high()); + + } + return false; +} + +/* + * \brief Removes the leading term of the polynomial \c f as long as it + * vanishes at \c alpha + * + */ +template +Poly_2 poly_non_vanish_leading_term(const Poly_2& pol,Algebraic_real alpha) { + Poly_2 f(pol); + while(true) { + if(CGAL::CGALi::is_root_of(alpha,f.lcoeff())) { + typename Poly_2::const_iterator poly_end = f.end(); + if(f.begin()==poly_end) { + break; + } + poly_end--; + f=Poly_2(f.begin(),poly_end); + } + else { + break; + } + } + return f; +} + +} // namespace CGALi + + +CGAL_END_NAMESPACE + +#endif // CGAL_ALG_REAL_UTILS diff --git a/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/enums.h b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/enums.h new file mode 100644 index 00000000000..da52a44cfd5 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/enums.h @@ -0,0 +1,41 @@ +// TODO: Add licence +// +// 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) : Michael Kerber +// +// ============================================================================ + + +#ifndef AcX_ENUMS_H +#define AcX_ENUMS_H 1 + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + + enum Three_valued { + + ROOT_OF_FIRST_SET = 1, + ROOT_OF_BOTH_SETS = 0, + ROOT_OF_SECOND_SET=-1, + + }; + + enum Degeneracy_strategy { + + SHEAR_STRATEGY = 0, + EXCEPTION_STRATEGY = 1, + + }; + +} // namespace CGALi + +CGAL_END_NAMESPACE + +#endif