From cbdca2c35d689ccbf705d7cb04036b255a53a9df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Luis=20Pe=C3=B1aranda?= Date: Tue, 19 Nov 2013 15:43:29 -0200 Subject: [PATCH] Code refactoring. The obsolete and unused code was removed. The interface was cleaned. The memory leaks dissapeared (according to valgrind), because the pointers to RS memory were removed. For test purposes, the rational interface is not tested. --- Algebraic_kernel_d/include/CGAL/RS/ak_1.h | 179 +++ .../include/CGAL/RS/algebraic_1.h | 428 +++++-- .../include/CGAL/RS/bisection_refiner_1.h | 127 ++ .../include/CGAL/RS/comparator_1.h | 89 ++ .../include/CGAL/RS/functors_1.h | 1059 ++++++++++------- .../include/CGAL/RS/rs23_k_isolator_1.h | 123 ++ .../include/CGAL/RS/rs2_calls.h | 134 +++ .../include/CGAL/RS/rs2_isolator_1.h | 105 ++ .../include/CGAL/RS/rs3_k_refiner_1.h | 108 ++ .../include/CGAL/RS/rs3_refiner_1.h | 108 ++ Algebraic_kernel_d/include/CGAL/RS/signat_1.h | 84 ++ 11 files changed, 1997 insertions(+), 547 deletions(-) create mode 100644 Algebraic_kernel_d/include/CGAL/RS/ak_1.h create mode 100644 Algebraic_kernel_d/include/CGAL/RS/bisection_refiner_1.h create mode 100644 Algebraic_kernel_d/include/CGAL/RS/comparator_1.h create mode 100644 Algebraic_kernel_d/include/CGAL/RS/rs23_k_isolator_1.h create mode 100644 Algebraic_kernel_d/include/CGAL/RS/rs2_calls.h create mode 100644 Algebraic_kernel_d/include/CGAL/RS/rs2_isolator_1.h create mode 100644 Algebraic_kernel_d/include/CGAL/RS/rs3_k_refiner_1.h create mode 100644 Algebraic_kernel_d/include/CGAL/RS/rs3_refiner_1.h create mode 100644 Algebraic_kernel_d/include/CGAL/RS/signat_1.h diff --git a/Algebraic_kernel_d/include/CGAL/RS/ak_1.h b/Algebraic_kernel_d/include/CGAL/RS/ak_1.h new file mode 100644 index 00000000000..e75651879ff --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/ak_1.h @@ -0,0 +1,179 @@ +// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author: Luis Peñaranda + +#ifndef CGAL_RS_AK_1_H +#define CGAL_RS_AK_1_H + +#include // included only to define size_t +#include +#include "algebraic_1.h" +#include "comparator_1.h" +#include "signat_1.h" + +namespace CGAL{ +namespace RS_AK1{ + +template > +class Algebraic_kernel_1{ + public: + typedef Polynomial_ Polynomial_1; + typedef typename Polynomial_1::NT Coefficient; + typedef Bound_ Bound; + private: + typedef Isolator_ Isolator; + typedef Refiner_ Refiner; + typedef Ptraits_ Ptraits; + typedef CGAL::RS_AK1::Signat_1 + Signat; + typedef CGAL::RS_AK1::Simple_comparator_1 + Comparator; + public: + typedef CGAL::RS_AK1::Algebraic_1 + Algebraic_real_1; + typedef size_t size_type; + typedef unsigned Multiplicity_type; + + // default constructor and destructor + public: + Algebraic_kernel_1(){}; + ~Algebraic_kernel_1(){}; + + // functors from the CGAL concept + public: + typedef CGAL::RS_AK1::Construct_algebraic_real_1 + Construct_algebraic_real_1; + typedef CGAL::RS_AK1::Compute_polynomial_1 + Compute_polynomial_1; + typedef CGAL::RS_AK1::Isolate_1 + Isolate_1; + typedef typename Ptraits::Is_square_free Is_square_free_1; + typedef typename Ptraits::Make_square_free Make_square_free_1; + typedef typename Ptraits::Square_free_factorize + Square_free_factorize_1; + typedef CGAL::RS_AK1::Is_coprime_1 + Is_coprime_1; + typedef CGAL::RS_AK1::Make_coprime_1 + Make_coprime_1; + typedef CGAL::RS_AK1::Solve_1 + Solve_1; + typedef CGAL::RS_AK1::Number_of_solutions_1 + Number_of_solutions_1; + + typedef CGAL::RS_AK1::Sign_at_1 + Sign_at_1; + typedef CGAL::RS_AK1::Is_zero_at_1 + Is_zero_at_1; + typedef CGAL::RS_AK1::Compare_1 + Compare_1; + typedef CGAL::RS_AK1::Bound_between_1 + Bound_between_1; + typedef CGAL::RS_AK1::Approximate_absolute_1 + Approximate_absolute_1; + typedef CGAL::RS_AK1::Approximate_relative_1 + Approximate_relative_1; + +#define CREATE_FUNCTION_OBJECT(T,N) \ + T N##_object()const{return T();} + CREATE_FUNCTION_OBJECT(Construct_algebraic_real_1, + construct_algebraic_real_1) + CREATE_FUNCTION_OBJECT(Compute_polynomial_1, + compute_polynomial_1) + CREATE_FUNCTION_OBJECT(Isolate_1, + isolate_1) + CREATE_FUNCTION_OBJECT(Is_square_free_1, + is_square_free_1) + CREATE_FUNCTION_OBJECT(Make_square_free_1, + make_square_free_1) + CREATE_FUNCTION_OBJECT(Square_free_factorize_1, + square_free_factorize_1) + CREATE_FUNCTION_OBJECT(Is_coprime_1, + is_coprime_1) + CREATE_FUNCTION_OBJECT(Make_coprime_1, + make_coprime_1) + CREATE_FUNCTION_OBJECT(Solve_1, + solve_1) + CREATE_FUNCTION_OBJECT(Number_of_solutions_1, + number_of_solutions_1) + CREATE_FUNCTION_OBJECT(Sign_at_1, + sign_at_1) + CREATE_FUNCTION_OBJECT(Is_zero_at_1, + is_zero_at_1) + CREATE_FUNCTION_OBJECT(Compare_1, + compare_1) + CREATE_FUNCTION_OBJECT(Bound_between_1, + bound_between_1) + CREATE_FUNCTION_OBJECT(Approximate_absolute_1, + approximate_absolute_1) + CREATE_FUNCTION_OBJECT(Approximate_relative_1, + approximate_relative_1) +#undef CREATE_FUNCTION_OBJECT + +}; // class Algebraic_kernel_1 + +} // namespace RS_AK1 +} // namespace CGAL + +#endif // CGAL_RS_AK_1_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1.h b/Algebraic_kernel_d/include/CGAL/RS/algebraic_1.h index 6fa85138678..cb0aa8ab565 100644 --- a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1.h +++ b/Algebraic_kernel_d/include/CGAL/RS/algebraic_1.h @@ -1,9 +1,9 @@ -// Copyright (c) 2006-2008 Inria Lorraine (France). All rights reserved. +// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved. // // This file is part of CGAL (www.cgal.org); you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public License as -// published by the Free Software Foundation; either version 3 of the License, -// or (at your option) any later version. +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. // // Licensees holding a valid commercial license may use this file in // accordance with the commercial license agreement provided with the software. @@ -19,123 +19,337 @@ #ifndef CGAL_RS_ALGEBRAIC_1_H #define CGAL_RS_ALGEBRAIC_1_H -#include -#include -#include -#include #include +#include +#include +#include namespace CGAL{ +namespace RS_AK1{ -class RS_polynomial_1; -class Algebraic_1; +// This class represents the simplest algebraic number one can think about. +// One algebraic number is represented by the polynomial of which it is +// root and the two endpoints of an interval that contains it, and no other +// root. Polynomial type and bound type are the first two template +// parameters. +// +// The third template parameter is a refiner, a function object that +// receives the polynomial and the bounds defining an algebraic number and +// an integer p, and modifies the two bounds until the difference between +// the two bounds is less than x*2^(-p), where x is the value of the +// represented algebraic number. The signature of a refiner must be: +// void +// Refiner_()(const Polynomial_&,Bound_&,Bound_&,int p); +// +// The fourth template argument is a comparator, a function object that +// receives the polynomials and bounds defining two algebraic numbres and +// just compares them, returning a CGAL::Comparison_result. The signature +// of a comparator must be: +// CGAL::Comparison_result +// Comparator_()(const Polynomial_&,Bound_&,Bound_&, +// const Polynomial_&,Bound_&,Bound_&); +// The comparator can modify the bounds, with the condition that the +// algebraic numbers remain consistent (one and only one root on each +// interval). -bool operator<(const Algebraic_1&,const Algebraic_1&); -bool operator==(const Algebraic_1&,const Algebraic_1&); +template +class Algebraic_1: +boost::totally_ordered, + double>{ + private: + typedef Polynomial_ Polynomial; + typedef Bound_ Bound; + typedef Refiner_ Refiner; + typedef Comparator_ Comparator; + typedef Ptraits_ Ptraits; + typedef typename Ptraits::Coefficient_type Coefficient; + typedef typename Ptraits::Scale Scale; + typedef Algebraic_1 + Algebraic; -// representation of algebraic numbers -class Algebraic_1_rep{ -public: - mutable mpfi_t interval; - RS_polynomial_1 *polynomial; - int numroot; - int multiplicity; - mutable Sign lefteval; + Polynomial pol; + mutable Bound left,right; - Algebraic_1_rep(): - polynomial(NULL),numroot(-1),multiplicity(-1),lefteval(ZERO){} - ~Algebraic_1_rep(){} + public: + Algebraic_1(){}; + Algebraic_1(const Polynomial &p, + const Bound &l, + const Bound &r):pol(p),left(l),right(r){ + CGAL_assertion(l<=r); + } + Algebraic_1(const Algebraic &a): + pol(a.pol),left(a.left),right(a.right){} + // This assumes that Gmpq is constructible from bound type and + // that polynomial coefficient type is constructible from mpz_t. + Algebraic_1(const Bound &b):left(b),right(b){ + typedef typename Ptraits::Shift shift; + Gmpq q(b); + pol=Coefficient(mpq_denref(q.mpq()))* + shift()(Polynomial(1),1,0)- + Coefficient(mpq_numref(q.mpq())); + CGAL_assertion(left==right&&left==b); + } + // TODO: make this constructor generic, the coefficient type is + // assumed to be constructible from mpz_t (rewrite in terms of + // Gmpq/Gmpz) + Algebraic_1(double d){ + typedef typename Ptraits::Shift shift; + Gmpq q(d); + pol=Coefficient(mpq_denref(q.mpq()))* + shift()(Polynomial(1),1,0)- + Coefficient(mpq_numref(q.mpq())); + left=Bound(d/*,std::round_toward_neg_infinity*/); + right=Bound(d/*,std::round_toward_infinity*/); + CGAL_assertion((left==right&&left==d)||(leftd)); + } + // TODO: Constructors from types such as int, unsigned and long. This + // implementation assumes that the bound type is Gmpfr and that T can + // be exactly converted to Gmpq. + template + Algebraic_1(const T &t){ + typedef typename Ptraits::Shift shift; + CGAL::Gmpq q(t); + pol=Coefficient(mpq_denref(q.mpq()))* + shift()(Polynomial(1),1,0)- + Coefficient(mpq_numref(q.mpq())); + left=Bound(t,std::round_toward_neg_infinity); + right=Bound(t,std::round_toward_infinity); + CGAL_assertion(left<=t&&right>=t); + } + ~Algebraic_1(){} -private: - Algebraic_1_rep(const Algebraic_1_rep&); - Algebraic_1_rep& operator=(const Algebraic_1_rep&); + Algebraic_1& operator=(const Algebraic &a){ + pol=a.get_pol(); + left=a.get_left(); + right=a.get_right(); + } + + Polynomial get_pol()const{return pol;} + Bound& get_left()const{return left;} + Bound& get_right()const{return right;} + + Algebraic operator-()const{ + return Algebraic(Scale()(get_pol(),Coefficient(-1)), + -right, + -left); + } + +#define CGAL_RS_COMPARE_ALGEBRAIC(_a) \ + (Comparator()(get_pol(),get_left(),get_right(), \ + (_a).get_pol(),(_a).get_left(),(_a).get_right())) + + Comparison_result compare(Algebraic a)const{ + return CGAL_RS_COMPARE_ALGEBRAIC(a); + }; + +#define CGAL_RS_COMPARE_ALGEBRAIC_TYPE(_t) \ + bool operator<(_t t)const \ + {Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::SMALLER;} \ + bool operator>(_t t)const \ + {Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::LARGER;} \ + bool operator==(_t t)const \ + {Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::EQUAL;} + + bool operator==(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::EQUAL;} + bool operator!=(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC(a)!=CGAL::EQUAL;} + bool operator<(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::SMALLER;} + bool operator<=(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC(a)!=CGAL::LARGER;} + bool operator>(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC(a)==CGAL::LARGER;} + bool operator>=(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC(a)!=CGAL::SMALLER;} + + CGAL_RS_COMPARE_ALGEBRAIC_TYPE(double) + +#undef CGAL_RS_COMPARE_ALGEBRAIC_TYPE +#undef CGAL_RS_COMPARE_ALGEBRAIC + +#ifdef IEEE_DBL_MANT_DIG +#define CGAL_RS_DBL_PREC IEEE_DBL_MANT_DIG +#else +#define CGAL_RS_DBL_PREC 53 +#endif + double to_double(){ + typedef Real_embeddable_traits RT; + typedef typename RT::To_double TD; + Refiner()(get_pol(),get_left(),get_right(),CGAL_RS_DBL_PREC); + CGAL_assertion(TD()(get_left())==TD()(get_right())); + return TD()(get_left()); + } + std::pair to_interval()const{ + typedef Real_embeddable_traits RT; + typedef typename RT::To_interval TI; + return std::make_pair(TI()(get_left()).first, + TI()(get_right()).second); + } +#undef CGAL_RS_DBL_PREC + + void set_left(const Bound &l)const{ + left=l; + } + void set_right(const Bound &r)const{ + right=r; + } + void set_pol(const Polynomial &p){ + pol=p; + } + +}; // class Algebraic_1 + +} // namespace RS_AK1 + +// We define Algebraic_1 as real-embeddable +template +class Real_embeddable_traits >: +public INTERN_RET::Real_embeddable_traits_base< + RS_AK1::Algebraic_1, + CGAL::Tag_true>{ + typedef Polynomial_ P; + typedef Bound_ B; + typedef Refiner_ R; + typedef Comparator_ C; + typedef Ptraits_ T; + + public: + + typedef RS_AK1::Algebraic_1 Type; + typedef CGAL::Tag_true Is_real_embeddable; + typedef bool Boolean; + typedef CGAL::Sign Sign; + typedef CGAL::Comparison_result Comparison_result; + + typedef INTERN_RET::Real_embeddable_traits_base + Base; + typedef typename Base::Compare Compare; + + class Sgn:public std::unary_function{ + public: + CGAL::Sign operator()(const Type &a)const{ + return Compare()(a,Type(0)); + } + }; + + class To_double:public std::unary_function{ + public: + double operator()(Type a)const{return a.to_double();} + }; + + class To_interval: + public std::unary_function >{ + public: + std::pair operator()(const Type &a)const{ + return a.to_interval(); + } + }; + + class Is_zero:public std::unary_function{ + public: + bool operator()(const Type &a)const{ + return Sgn()(a)==CGAL::ZERO; + } + }; + + class Is_finite:public std::unary_function{ + public: + bool operator()(const Type&)const{return true;} + }; + + class Abs:public std::unary_function{ + public: + Type operator()(const Type &a)const{ + return Sgn()(a)==CGAL::NEGATIVE?-a:a; + } + }; }; -// The class of the algebraic numbers. -class Algebraic_1 : Handle_for, - boost::totally_ordered1{ +template +inline +RS_AK1::Algebraic_1 min +BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_1 a, + RS_AK1::Algebraic_1 b){ + return(a Base; +template +inline +RS_AK1::Algebraic_1 max +BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_1 a, + RS_AK1::Algebraic_1 b){ + return(a>b?a:b); +} -public: +template +inline +std::ostream& operator<<(std::ostream &o, + const RS_AK1::Algebraic_1 &a){ + return(o<<'['< - double to_double()const; - std::pair to_interval() const; - Algebraic_1 sqrt()const; -}; // class Algebraic_1 +template +inline +std::istream& operator>>(std::istream &i, + RS_AK1::Algebraic_1 &a){ + // TODO: cleanly write this function + std::istream::int_type c; + P pol; + B lb,rb; + c=i.get(); + if(c!='['){ + CGAL_error_msg("error reading istream, \'[\' expected"); + return i; + } + i>>pol; + c=i.get(); + if(c!=','){ + CGAL_error_msg("error reading istream, \',\' expected"); + return i; + } + i>>lb; + c=i.get(); + if(c!=','){ + CGAL_error_msg("error reading istream, \',\' expected"); + return i; + } + i>>rb; + c=i.get(); + if(c!=']'){ + CGAL_error_msg("error reading istream, \']\' expected"); + return i; + } + a=RS_AK1::Algebraic_1(pol,lb,rb); + return i; +} } // namespace CGAL -#include -#include -#include // other member functions -#include -#include // related non-member functions -#include - -#endif // CGAL_RS_ALGEBRAIC_1_H +#endif // CGAL_RS_ALGEBRAIC_1_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/bisection_refiner_1.h b/Algebraic_kernel_d/include/CGAL/RS/bisection_refiner_1.h new file mode 100644 index 00000000000..05f5a91bf34 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/bisection_refiner_1.h @@ -0,0 +1,127 @@ +// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author: Luis Peñaranda + +// This file contains the simplest refiner, that bisects the interval a given +// number of times. + +#ifndef CGAL_RS_BISECTION_REFINER_1_H +#define CGAL_RS_BISECTION_REFINER_1_H + +#include +#include "signat_1.h" + +namespace CGAL{ + +template +struct Bisection_refiner_1{ + typedef CGAL::RS_AK1::Signat_1 Signat; + void operator()(const Polynomial_&,Bound_&,Bound_&,int); +}; // class Bisection_refiner_1 + +// TODO: Write in a generic way, if possible. +template +void +Bisection_refiner_1:: +operator()(const Polynomial_&,Bound_&,Bound_&,int){ + CGAL_error_msg("bisection refiner not implemented for these types"); + return; +} + +// TODO: Optimize this function. +template<> +void +Bisection_refiner_1,Gmpfr>:: +operator()(const Polynomial &pol,Gmpfr &left,Gmpfr &right,int prec){ + typedef Polynomial Polynomial; + typedef Polynomial_traits_d Ptraits; + typedef Ptraits::Degree Degree; + typedef Ptraits::Make_square_free Sfpart; + typedef CGAL::RS_AK1::Signat_1 + Signat; + CGAL_precondition(left<=right); + // TODO: add precondition to check whether the interval is a point + // or the evaluations on its endpoints have different signs + //std::cout<<"refining ["<pc?pl:pc)+(mp_prec_t)prec; + mpfr_init2(center,pc); + CGAL_assertion_code(int round=) + mpfr_prec_round(left.fr(),pc,GMP_RNDN); + CGAL_assertion(!round); + CGAL_assertion_code(round=) + mpfr_prec_round(right.fr(),pc,GMP_RNDN); + CGAL_assertion(!round); + for(int i=0;i + +#ifndef CGAL_RS_COMPARATOR_1_H +#define CGAL_RS_COMPARATOR_1_H + +namespace CGAL{ +namespace RS_AK1{ + +template +struct Simple_comparator_1{ + typedef Polynomial_ Polynomial; + typedef Bound_ Bound; + typedef Refiner_ Refiner; + typedef Signat_ Signat; + typedef Ptraits_ Ptraits; + typedef typename Ptraits::Gcd_up_to_constant_factor Gcd; + typedef typename Ptraits::Degree Degree; + + CGAL::Comparison_result + operator()(const Polynomial &p1,Bound &l1,Bound &r1, + const Polynomial &p2,Bound &l2,Bound &r2)const{ + CGAL_precondition(l1<=r1&&l2<=r2); + if(l1<=l2){ + if(r1l2?l1:l2); + if(sleft==ZERO) + return EQUAL; + CGAL::Sign sright=sg(r1=l2:r2>=l1); + return (r1 -#ifndef CGAL_RS_FUNCTORS_H -#define CGAL_RS_FUNCTORS_H +#ifndef CGAL_RS_FUNCTORS_1_H +#define CGAL_RS_FUNCTORS_1_H -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include -#ifdef IEEE_DBL_MANT_DIG -# define CGAL_RS_FUNCTORS_DBL_PREC IEEE_DBL_MANT_DIG -#else -# define CGAL_RS_FUNCTORS_DBL_PREC 53 -#endif +namespace CGAL{ +namespace RS_AK1{ -namespace RSFunctors{ +template +struct Construct_algebraic_real_1{ + typedef Polynomial_ Polynomial; + typedef Algebraic_ Algebraic; + typedef Bound_ Bound; + typedef Coefficient_ Coefficient; + typedef Isolator_ Isolator; + typedef Algebraic result_type; -typedef CGAL::RS_polynomial_1 Polynomial; -typedef CGAL::Algebraic_1 Algebraic; -typedef CGAL::Gmpfr Bound; -typedef int Multiplicity; + template + Algebraic operator()(const T &a)const{ + return Algebraic(a); + } -template + Algebraic operator()(const Polynomial &p,size_t i)const{ + Isolator isol(p); + return Algebraic(p,isol.left_bound(i),isol.right_bound(i)); + } + + Algebraic operator()(const Polynomial &p, + const Bound &l, + const Bound &r)const{ + return Algebraic(p,l,r); + } + +}; // struct Construct_algebraic_1 + +template struct Compute_polynomial_1: -public std::unary_function{ - typedef PolynomialType P; - typedef CGAL::from_rs_poly

back; - P operator()(const Algebraic &a)const{ - return back()(a.pol()); +public std::unary_function{ + typedef Polynomial_ Polynomial; + typedef Algebraic_ Algebraic; + Polynomial operator()(const Algebraic &x)const{ + return x.get_pol(); } -}; // Compute_polynomial_1 +}; // struct Compute_polynomial_1 -template -struct Is_square_free_1: -public std::unary_function{ - typedef PolynomialType P; - typedef CGAL::to_rs_poly

convert; - typedef GcdPolicy Gcd; - bool operator()(const P &p)const{ - Polynomial rsp=convert()(p); - return(!(Gcd()(rsp,rsp.derive()).get_degree_static())); - } -}; // Is_square_free_1 - -template -struct Make_square_free_1: -public std::unary_function{ - typedef PolynomialType P; - typedef CGAL::to_rs_poly

convert; - typedef CGAL::from_rs_poly

back; - typedef GcdPolicy Gcd; - P operator()(const P &p)const{ - return back()(CGAL::sfpart_1()(convert()(p))); - } -}; // Make_square_free_1 - -template -struct Square_free_factorize_1{ - typedef PolynomialType P; - typedef CGAL::to_rs_poly

convert; - typedef CGAL::from_rs_poly

back; - typedef GcdPolicy Gcd; - template - OutputIterator operator()(const P &p,OutputIterator oi)const{ - Polynomial rsp=convert()(p); - CGAL::sqfrvec factorization(CGAL::sqfr_1()(rsp)); - for(CGAL::sqfrvec::iterator i=factorization.begin(); - i!=factorization.end(); - ++i){ - *oi++=std::make_pair(back()((*i).first),(*i).second); - } - return oi; - } -}; // Square_free_factorize_1 - -template +template struct Is_coprime_1: -public std::binary_function{ - typedef PolynomialType P; - typedef CGAL::to_rs_poly

convert; - typedef GcdPolicy Gcd; - bool operator()(const P &p1,const P &p2)const{ - CGAL::RS_polynomial_1 rsp1=convert()(p1); - CGAL::RS_polynomial_1 rsp2=convert()(p2); - return(!Gcd()(rsp1,rsp2).get_degree_static()); +public std::binary_function{ + typedef Polynomial_ Polynomial; + typedef Ptraits_ Ptraits; + typedef typename Ptraits::Gcd_up_to_constant_factor Gcd; + typedef typename Ptraits::Degree Degree; + inline bool operator()(const Polynomial &p1,const Polynomial &p2)const{ + return Degree()(Gcd()(p1,p2))==0; } -}; // Is_coprime_1 +}; // struct Is_coprime_1 -template +template struct Make_coprime_1{ - typedef PolynomialType P; - typedef CGAL::to_rs_poly

convert; - typedef CGAL::from_rs_poly

back; - typedef GcdPolicy Gcd; - bool operator()(const P &p1,const P &p2,P &g,P &q1,P &q2)const{ - typedef GcdPolicy Gcd; - CGAL::RS_polynomial_1 rsp1=convert()(p1); - CGAL::RS_polynomial_1 rsp2=convert()(p2); - CGAL::RS_polynomial_1 rsg=convert()(g); - rsg=Gcd()(rsp1,rsp2); - g=back()(rsg); - // even when g==1, we calculate q1 and q2 - q1=back()(*CGAL::Ediv_1()(rsp1,rsg)); - q2=back()(*CGAL::Ediv_1()(rsp2,rsg)); - return rsg.get_degree_static()?false:true; + typedef Polynomial_ Polynomial; + typedef Ptraits_ Ptraits; + typedef typename Ptraits::Gcd_up_to_constant_factor Gcd; + typedef typename Ptraits::Degree Degree; + typedef typename Ptraits::Integral_division_up_to_constant_factor + IDiv; + bool operator()(const Polynomial &p1, + const Polynomial &p2, + Polynomial &g, + Polynomial &q1, + Polynomial &q2)const{ + g=Gcd()(p1,p2); + q1=IDiv()(p1,g); + q2=IDiv()(p2,g); + return Degree()(Gcd()(p1,p2))==0; } -}; // Make_coprime_1 +}; // struct Make_coprime_1 -template -struct Solve_RS_1{ - typedef GcdPolicy Gcd; - typedef CGAL::sfpart_1 Sfpart; - typedef CGAL::sqfr_1 Sqfr; - - template - OutputIterator operator()(const Polynomial &p,OutputIterator res)const{ - int nr,*m; - mpfi_ptr *x; - CGAL::sqfrvec sfv=Sqfr()(p); - x=(mpfi_ptr*)malloc(Sfpart()(p).get_degree()*sizeof(mpfi_ptr)); - m=(int*)calloc(Sfpart()(p).get_degree(),sizeof(int)); - nr=solve_1(x,Sfpart()(p)); - CGAL_assertion_msg(nr>=0,"error in resolution"); - for(CGAL::sqfrvec::size_type i=0;ileft)); - CGAL::Sign sg_r= - CGAL::RSSign::signat(sfv[i].first,&(x[j]->right)); - if((sg_l!=sg_r)||((sg_l==CGAL::ZERO)&&(sg_r==CGAL::ZERO))){ - m[j]=sfv[i].second; - --k; - } - } - } - } - for(int i=0;i - OutputIterator operator()(const Polynomial &p, - bool known_to_be_square_free, - OutputIterator res)const{ - int nr,m; - mpfi_ptr *x; - if(known_to_be_square_free){ - p.set_sf(); - x=(mpfi_ptr*)malloc(p.get_degree()*sizeof(mpfi_ptr)); - nr=solve_1(x,p); - CGAL_assertion_msg(nr>=0,"error in resolution"); - m=1; // we know in this case that multiplicity is 1 - }else{ - x=(mpfi_ptr*)malloc(Sfpart()(p).get_degree()*sizeof(mpfi_ptr)); - nr=solve_1(x,Sfpart()(p)); - CGAL_assertion_msg(nr>=0,"error in resolution"); - m=0; // we won't calculate multiplicities - } - for(int i=0;i +template struct Solve_1{ - typedef PolynomialType P; - typedef CGAL::to_rs_poly

convert; - typedef GcdPolicy Gcd; - typedef Solve_RS_1 Solve_RS; + typedef Polynomial_ Polynomial_1; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Isolator_ Isolator; + typedef Signat_ Signat; + typedef Ptraits_ Ptraits; + typedef typename Ptraits::Gcd_up_to_constant_factor Gcd; + typedef typename Ptraits::Square_free_factorize_up_to_constant_factor + Sqfr; + typedef typename Ptraits::Degree Degree; + typedef typename Ptraits::Make_square_free Sfpart; - template - OutputIterator operator()(const P &p,OutputIterator res)const{ - return Solve_RS()(convert()(p),res); - } + template + OutputIterator operator()(const Polynomial_1 &p, + OutputIterator res)const{ + typedef std::pair polmult; + typedef std::vector sqvec; - template - OutputIterator operator()(const P &p, - bool known_to_be_square_free, - OutputIterator res)const{ - return Solve_RS()(convert()(p),known_to_be_square_free,res); - } - - template - OutputIterator operator()(const P &p, - const Bound& lower, - const Bound& upper, - OutputIterator res)const{ - typedef std::vector > RMV; - RMV roots; - this->operator()(p,std::back_inserter(roots)); - for(typename RMV::iterator it = roots.begin(); it != roots.end();it++){ - if(lower <= it->first && it->first <= upper) - *res++=*it; - } - return res; - } - - template - OutputIterator operator()(const P &p, - bool known_to_be_square_free, - const Bound& lower, - const Bound& upper, - OutputIterator res)const{ - typedef std::vector< Algebraic > RV; - RV roots; - this->operator()(p,known_to_be_square_free,std::back_inserter(roots)); - for(typename RV::iterator it = roots.begin(); it != roots.end();it++){ - if(lower <= *it && *it <= upper) - *res++=*it; - } - return res; - } -}; // Solve_1 - -template -struct Construct_alg_1{ - typedef PolynomialType Poly; - typedef CoeffType Coeff; - typedef GcdFunction Gcd; - typedef CGAL::to_rs_poly convert; - typedef Solve_1 Solve; - - Algebraic operator()(int a) const { - return Algebraic(a); - } - - Algebraic operator()(const Bound a) const { - return Algebraic(a); - } - - Algebraic operator()(const Coeff a) const { - return Algebraic(a); - } - - Algebraic operator()(const Poly &p,int i) const { - CGAL_precondition(CGAL::is_square_free(p)); - std::vector roots; - std::back_insert_iterator > rootsit(roots); - Solve()(p,true,rootsit); - return roots[i]; - } - - Algebraic operator()(const Poly &p,Bound l,Bound u) const { - mpfi_t i; - mpfi_init(i); - mpfr_set(&i->left,l.fr(),GMP_RNDD); - mpfr_set(&i->right,u.fr(),GMP_RNDU); - return Algebraic(i,convert()(p),0,0 - //,NULL,NULL - ); - } -}; // Construct_alg_1 - -template -struct Number_of_solutions_1: -public std::unary_function{ - typedef PolynomialType P; - typedef CGAL::to_rs_poly

convert; - - int operator()(const P &p)const{ - int nr; - mpfi_ptr *x; - CGAL::RS_polynomial_1 rspoly=convert()(p); - x=(mpfi_ptr*)malloc(rspoly.get_degree()*sizeof(mpfi_ptr)); - nr=solve_1(x,rspoly); - CGAL_assertion_msg(nr>=0,"error in resolution"); - free(x); - return nr; - } -}; // Number_of_solutions_1 - -template -struct Sign_at_1: -public std::binary_function{ - typedef PolynomialType P; - typedef GcdPolicy Gcd; - typedef CGAL::to_rs_poly

convert; - - CGAL::Sign operator()(const P &p,const Algebraic &a)const{ - return RS3::sign_1(convert()(p),a); - } -}; // Sign_at_1 - -template -struct Is_zero_at_1: -public std::binary_function{ - typedef PolynomialType P; - typedef GcdPolicy Gcd; - typedef Sign_at_1 Sign_at; - - bool operator()(const P &p,const Algebraic &a)const{ - return (Sign_at()(p,a)==CGAL::ZERO); - } -}; // Is_zero_at_1 - -template -struct Compare_1: - public std::binary_function{ - typedef GcdPolicy Gcd; - typedef CGAL::Comparison_result Comparison_result; - typedef CGAL::Gmpz Gmpz; - typedef CGAL::Gmpq Gmpq; - - Comparison_result operator()(const Algebraic &r1,const Algebraic &r2)const{ - return CGAL::RS_COMPARE::compare_1(r1,r2); - } - - Comparison_result operator()(const int &r1, const Algebraic &r2)const{ - return this->operator()(Algebraic(r1),r2);} - Comparison_result operator()(const Bound &r1,const Algebraic &r2)const{ - return this->operator()(Algebraic(r1),r2);} - Comparison_result operator()(const Gmpz &r1, const Algebraic &r2)const{ - return this->operator()(Algebraic(r1),r2);} - Comparison_result operator()(const Gmpq &r1, const Algebraic &r2)const{ - return this->operator()(Algebraic(r1),r2);} - Comparison_result operator()(const Algebraic &r1,const int &r2)const{ - return this->operator()(r1,Algebraic(r2));} - Comparison_result operator()(const Algebraic &r1,const Bound &r2)const{ - return this->operator()(r1,Algebraic(r2));} - Comparison_result operator()(const Algebraic &r1,const Gmpz &r2)const{ - return this->operator()(r1,Algebraic(r2));} - Comparison_result operator()(const Algebraic &r1,const Gmpq &r2)const{ - return this->operator()(r1,Algebraic(r2));} -}; // Compare_1 - -template -struct Isolate_1: -public std::binary_function >{ - typedef PolynomialType Poly; - typedef GcdFunction Gcd; - typedef CGAL::to_rs_poly convert; - typedef Solve_1 Solve; - typedef Compare_1 Compare; - - std::pair operator()(const Algebraic &a,const Poly &p)const{ - std::vector roots; - std::back_insert_iterator > rootsit(roots); - Solve()(p,true,rootsit); - for(std::vector::size_type i=0;i -struct Bound_between_1: - public std::binary_function{ - typedef GcdPolicy Gcd; - Bound operator()(const Algebraic &x1,const Algebraic &x2)const{ - double l,r,m; - switch(CGAL::RS_COMPARE::compare_1(x1,x2)){ - case CGAL::LARGER: - CGAL_assertion(x2.sup() - x1.inf().get_precision()? - 1+x2.sup().get_precision(): - 1+x1.inf().get_precision()))/2; - break; - case CGAL::SMALLER: - CGAL_assertion(x1.sup() - x2.inf().get_precision()? - 1+x1.sup().get_precision(): - 1+x2.inf().get_precision()))/2; - break; - default: - CGAL_error_msg("bound between two equal numbers"); - } + Polynomial_1 sfp=Sfpart()(p); + sqvec sfv; + Sqfr()(p,std::back_inserter(sfv)); + Isolator isol(sfp); + int *m=(int*)calloc(isol.number_of_real_roots(),sizeof(int)); + for(typename sqvec::iterator i=sfv.begin();i!=sfv.end();++i){ + int k=Degree()(i->first); + Signat signof(i->first); + for(int j=0;k&&jsecond; + --k; + } + } + } + } + for(int l=0;l + OutputIterator operator()(const Polynomial_1 &p, + bool known_to_be_square_free, + OutputIterator res)const{ + Isolator isol(p); + for(int l=0;l + OutputIterator operator()(const Polynomial_1 &p, + const Bound &l, + const Bound &u, + OutputIterator res)const{ + typedef std::vector RV; + typedef typename RV::iterator RVI; + typedef std::pair PM; + typedef std::vector PMV; + typedef typename PMV::iterator PMVI; + CGAL_precondition_msg(l<=u, + "left bound must be <= right bound"); + RV roots; // all roots of the polynomial + this->operator()(p,false,std::back_inserter(roots)); + size_t nb_roots=roots.size(); + // indices of the first and last roots to be reported: + int index_l=0,index_u; + while(index_l and match the + // roots in the interval [index_l,index_u) + for(PMVI i=sfv.begin();i!=sfv.end();++i){ + int k=Degree()(i->first); + Signat signof(i->first); + for(int j=index_l;k&&jsecond; + --k; + } + } + } + } + for(int l=index_l;l + OutputIterator operator()(const Polynomial_1 &p, + bool known_to_be_square_free, + const Bound &l, + const Bound &u, + OutputIterator res)const{ + typedef std::vector RV; + typedef typename RV::iterator RVI; + CGAL_precondition_msg(l<=u, + "left bound must be <= right bound"); + RV roots; + this->operator()(p, + known_to_be_square_free, + std::back_inserter(roots)); + for(RVI it=roots.begin();it!=roots.end();it++) + if(*it>=l&&*it<=u) + *res++=*it; + return res; + } + +}; // struct Solve_1 + +template +class Sign_at_1: +public std::binary_function{ + // This implementation will work with any polynomial type whose + // coefficient type is explicit interoperable with Gmpfi. + // TODO: Make this function generic. + public: + typedef Polynomial_ Polynomial_1; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Refiner_ Refiner; + typedef Signat_ Signat; + typedef Ptraits_ Ptraits; + + private: + CGAL::Uncertain eval_interv(const Polynomial_1 &p, + const Bound &l, + const Bound &r)const{ + typedef typename Ptraits::Substitute Subst; + std::vector substitutions; + substitutions.push_back(CGAL::Gmpfi(l,r)); + CGAL::Gmpfi eval=Subst()(p, + substitutions.begin(), + substitutions.end()); + return eval.sign(); + } + + // This function assumes that the sign of the evaluation is not zero, + // it just refines x until having the correct sign. + CGAL::Sign refine_and_return(const Polynomial_1 &p,Algebraic x)const{ + CGAL::Gmpfr xl(x.get_left()); + CGAL::Gmpfr xr(x.get_right()); + CGAL::Uncertain s; + for(;;){ + Refiner()(x.get_pol(), + xl, + xr, + 2*CGAL::max(xl.get_precision(), + xr.get_precision())); + s=eval_interv(p,xl,xr); + if(!s.is_same(Uncertain::indeterminate())){ + x.set_left(xl); + x.set_right(xr); + return s; + } + } + } + + public: + CGAL::Sign operator()(const Polynomial_1 &p,Algebraic x)const{ + typedef typename Ptraits::Substitute Subst; + typedef typename Ptraits::Gcd_up_to_constant_factor Gcd; + typedef typename Ptraits::Make_square_free Sfpart; + typedef typename Ptraits::Degree Degree; + typedef typename Ptraits::Differentiate Deriv; + CGAL::Uncertain unknown= + Uncertain::indeterminate(); + CGAL::Uncertain s=eval_interv(p, + x.get_left(), + x.get_right()); + if(!s.is_same(unknown)) + return s; + // We are not sure about the sign. We calculate the gcd in + // order to know if both polynomials have common roots. + Polynomial_1 sfpp=Sfpart()(p); + Polynomial_1 gcd=Gcd()(sfpp,Sfpart()(x.get_pol())); + if(Degree()(gcd)==0) + return refine_and_return(p,x); + + // At this point, gcd is not 1; we proceed as follows: + // -we refine x until having p monotonic in x's interval (to be + // sure that p has at most one root on that interval), + // -if the gcd has a root on this interval, both roots are + // equal (we return 0), otherwise, we refine until having a + // result. + + // How to assure that p is monotonic in an interval: when its + // derivative is never zero in that interval. + Polynomial_1 dsfpp=Deriv()(sfpp); + CGAL::Gmpfr xl(x.get_left()); + CGAL::Gmpfr xr(x.get_right()); + while(eval_interv(dsfpp,xl,xr).is_same(unknown)){ + Refiner()(x.get_pol(), + xl, + xr, + 2*CGAL::max(xl.get_precision(), + xr.get_precision())); + } + x.set_left(xl); + x.set_right(xr); + + // How to know that the gcd has a root: evaluate endpoints. + CGAL::Sign sleft,sright; + Signat sign_at_gcd(gcd); + if((sleft=sign_at_gcd(x.get_left()))==CGAL::ZERO|| + (sright=sign_at_gcd(x.get_right()))==CGAL::ZERO|| + (sleft!=sright)) + return CGAL::ZERO; + return refine_and_return(p,x); + } +}; // struct Sign_at_1 + +template +class Is_zero_at_1: +public std::binary_function{ + // This implementation will work with any polynomial type whose + // coefficient type is explicit interoperable with Gmpfi. + // TODO: Make this function generic. + public: + typedef Polynomial_ Polynomial_1; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Refiner_ Refiner; + typedef Signat_ Signat; + typedef Ptraits_ Ptraits; + + private: + CGAL::Uncertain eval_interv(const Polynomial_1 &p, + const Bound &l, + const Bound &r)const{ + typedef typename Ptraits::Substitute Subst; + std::vector substitutions; + substitutions.push_back(CGAL::Gmpfi(l,r)); + CGAL::Gmpfi eval=Subst()(p, + substitutions.begin(), + substitutions.end()); + return eval.sign(); + } + + public: + bool operator()(const Polynomial_1 &p,Algebraic x)const{ + typedef typename Ptraits::Substitute Subst; + typedef typename Ptraits::Gcd_up_to_constant_factor Gcd; + typedef typename Ptraits::Make_square_free Sfpart; + typedef typename Ptraits::Degree Degree; + typedef typename Ptraits::Differentiate Deriv; + CGAL::Uncertain unknown= + Uncertain::indeterminate(); + CGAL::Uncertain s=eval_interv(p, + x.get_left(), + x.get_right()); + if(!s.is_same(unknown)) + return (s==CGAL::ZERO); + // We are not sure about the sign. We calculate the gcd in + // order to know if both polynomials have common roots. + Polynomial_1 sfpp=Sfpart()(p); + Polynomial_1 gcd=Gcd()(sfpp,Sfpart()(x.get_pol())); + if(Degree()(gcd)==0) + return false; + + // At this point, gcd is not 1; we proceed as follows: + // -we refine x until having p monotonic in x's interval (to be + // sure that p has at most one root on that interval), + // -if the gcd has a root on this interval, both roots are + // equal (we return 0), otherwise, we refine until having a + // result. + + // How to assure that p is monotonic in an interval: when its + // derivative is never zero in that interval. + Polynomial_1 dsfpp=Deriv()(sfpp); + CGAL::Gmpfr xl(x.get_left()); + CGAL::Gmpfr xr(x.get_right()); + while(eval_interv(dsfpp,xl,xr).is_same(unknown)){ + Refiner()(x.get_pol(), + xl, + xr, + 2*CGAL::max(xl.get_precision(), + xr.get_precision())); + } + x.set_left(xl); + x.set_right(xr); + + // How to know that the gcd has a root: evaluate endpoints. + CGAL::Sign sleft,sright; + Signat sign_at_gcd(gcd); + return((sleft=sign_at_gcd(x.get_left()))==CGAL::ZERO|| + (sright=sign_at_gcd(x.get_right()))==CGAL::ZERO|| + (sleft!=sright)); + } +}; // class Is_zero_at_1 + +// TODO: it says in the manual that this should return a size_type, but test +// programs assume that this is equal to int +template +struct Number_of_solutions_1: +public std::unary_function{ + typedef Polynomial_ Polynomial_1; + typedef Isolator_ Isolator; + size_t operator()(const Polynomial_1 &p)const{ + // TODO: make sure that p is square free (precondition) + Isolator isol(p); + return isol.number_of_real_roots(); + } +}; // struct Number_of_solutions_1 + +// This functor not only compares two algebraic numbers. In case they are +// different, it refines them until they do not overlap. +template +struct Compare_1: +public std::binary_function{ + typedef Algebraic_ Algebraic; + typedef Bound_ Bound; + typedef Comparator_ Comparator; + + CGAL::Comparison_result operator()(Algebraic a,Algebraic b)const{ + Bound al=a.get_left(); + Bound ar=a.get_right(); + Bound bl=b.get_left(); + Bound br=b.get_right(); + CGAL::Comparison_result c=Comparator()(a.get_pol(),al,ar, + b.get_pol(),bl,br); + a.set_left(al); + a.set_right(ar); + b.set_left(bl); + b.set_right(br); + return c; + } + + CGAL::Comparison_result operator()(Algebraic a,const Bound &b)const{ + Bound al=a.get_left(); + Bound ar=a.get_right(); + Algebraic balg(b); + CGAL::Comparison_result c=Comparator()(a.get_pol(),al,ar, + balg.get_pol(),b,b); + a.set_left(al); + a.set_right(ar); + return c; + } + + template + CGAL::Comparison_result operator()(Algebraic a,const T &b)const{ + Bound al=a.get_left(); + Bound ar=a.get_right(); + Algebraic balg(b); + CGAL::Comparison_result c=Comparator()(a.get_pol(), + al, + ar, + balg.get_pol(), + balg.get_left(), + balg.get_right()); + a.set_left(al); + a.set_right(ar); + return c; + } + +}; // class Compare_1 + +template +struct Bound_between_1: +public std::binary_function{ + typedef Algebraic_ Algebraic; + typedef Bound_ Bound; + typedef Comparator_ Comparator; + + Bound operator()(Algebraic a,Algebraic b)const{ + typedef Compare_1 Compare; + typename Bound::Precision_type prec; + switch(Compare()(a,b)){ + case CGAL::LARGER: + CGAL_assertion(b.get_right() +struct Isolate_1: +public std::binary_function >{ + typedef Polynomial_ Polynomial_1; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Isolator_ Isolator; + typedef Comparator_ Comparator; + typedef Signat_ Signat; + typedef Ptraits_ Ptraits; + + std::pair + operator()(const Algebraic &a,const Polynomial_1 &p)const{ + std::vector roots; + std::back_insert_iterator > rit(roots); + typedef Solve_1 Solve; + typedef Compare_1 Compare; + Solve()(p,false,rit); + for(typename std::vector::size_type i=0; + i struct Approximate_absolute_1: - public std::binary_function >{ +public std::binary_function >{ + typedef Polynomial_ Polynomial_1; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Refiner_ Refiner; - std::pair - operator()(const Algebraic& x, int prec) const { - RS3::refine_1(x,CGAL::abs(prec)); + // This implementation assumes that Bound is Gmpfr. + // TODO: make generic. + std::pair operator()(const Algebraic &x,const int &a)const{ + Bound xl(x.get_left()),xr(x.get_right()); + // refsteps=log2(xl-xr) + mpfr_t temp; + mpfr_init(temp); + mpfr_sub(temp,xr.fr(),xl.fr(),GMP_RNDU); + mpfr_log2(temp,temp,GMP_RNDU); + long refsteps=mpfr_get_si(temp,GMP_RNDU); + mpfr_clear(temp); + Refiner()(x.get_pol(),xl,xr,CGAL::abs(refsteps+a)); + x.set_left(xl); + x.set_right(xr); + CGAL_assertion(a>0? + (xr-xl)*CGAL::ipower(Bound(2),a)<=Bound(1): + (xr-xl)<=CGAL::ipower(Bound(2),-a)); + return std::make_pair(xl,xr); + } +}; // struct Approximate_absolute_1 - CGAL_postcondition_code( - CGAL::Gmpfr::Precision_type - subprec=1+ - std::max(x.inf().get_precision(), - x.sup().get_precision()); - CGAL::Gmpfr width=CGAL::Gmpfr::sub(x.sup(),x.inf(),subprec); - ) - CGAL_postcondition_code(if(prec>0)) - CGAL_postcondition(width*CGAL::ipower(Bound(2),prec)<=Bound(1)); - CGAL_postcondition_code(else) - CGAL_postcondition(width<=CGAL::ipower(Bound(2),-prec)); +template +struct Approximate_relative_1: +public std::binary_function >{ + typedef Polynomial_ Polynomial_1; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Refiner_ Refiner; - return std::make_pair(x.inf(),x.sup()); - } -}; + std::pair operator()(const Algebraic &x,const int &a)const{ + if(CGAL::is_zero(x)) + return std::make_pair(Bound(0),Bound(0)); + Bound error=CGAL::ipower(Bound(2),CGAL::abs(a)); + Bound xl(x.get_left()),xr(x.get_right()); + Bound max_b=(CGAL::max)(CGAL::abs(xr),CGAL::abs(xl)); + while(a>0?(xr-xl)*error>max_b:(xr-xl)>error*max_b){ + Refiner()(x.get_pol(), + xl, + xr, + std::max( + CGAL::abs(a), + CGAL::max(xl.get_precision(), + xr.get_precision()))); + max_b=(CGAL::max)(CGAL::abs(xr),CGAL::abs(xl)); + } + x.set_left(xl); + x.set_right(xr); + CGAL_assertion( + a>0? + (xr-xl)*CGAL::ipower(Bound(2),a)<=max_b: + (xr-xl)<=CGAL::ipower(Bound(2),-a)*max_b); + return std::make_pair(xl,xr); + } +}; // struct Approximate_relative_1 -struct Approximate_relative_1 - :public std::binary_function >{ +} // namespace RS_AK1 +} // namespace CGAL - std::pair operator()(const Algebraic &x, int prec) const { - if(CGAL::is_zero(x)) - return std::make_pair(Bound(0),Bound(0)); - - Bound error=CGAL::ipower(Bound(2),CGAL::abs(prec)); - Bound max_b=(CGAL::max)(CGAL::abs(x.sup()),CGAL::abs(x.inf())); - while(prec>0? - (x.sup()-x.inf())*error>max_b: - (x.sup()-x.inf())>error*max_b){ - RS3::refine_1(x,mpfi_get_prec(x.mpfi())+CGAL::abs(prec)); - max_b=(CGAL::max)(CGAL::abs(x.sup()),CGAL::abs(x.inf())); - } - - CGAL_postcondition_code(if(prec>0)) - CGAL_postcondition((x.sup()-x.inf())*CGAL::ipower(Bound(2),prec)<=max_b); - CGAL_postcondition_code(else) - CGAL_postcondition((x.sup()-x.inf())<=CGAL::ipower(Bound(2),-prec)*max_b); - - return std::make_pair(x.inf(),x.sup()); - } -}; - -} // namespace RSFunctors - -#undef CGAL_RS_FUNCTORS_DBL_PREC - -#endif // CGAL_RS_FUNCTORS_H +#endif // CGAL_RS_FUNCTORS_1_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/rs23_k_isolator_1.h b/Algebraic_kernel_d/include/CGAL/RS/rs23_k_isolator_1.h new file mode 100644 index 00000000000..0328b604eae --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/rs23_k_isolator_1.h @@ -0,0 +1,123 @@ +// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author: Luis Peñaranda + +#ifndef CGAL_RS_RS23_K_ISOLATOR_1_H +#define CGAL_RS_RS23_K_ISOLATOR_1_H + +// This file includes an isolator. Its particularity is that is isolates the +// roots with RS2, and the refines them until reaching Kantorovich criterion. +// This can take long, but the later refinements will be extremely fast with +// RS3. The functor is not in RS2 neither in RS3 namespace, because it uses +// functions from both. + +#include "rs2_calls.h" +#include "rs3_k_refiner_1.h" +#include +#include + +namespace CGAL{ + +template +class RS23_k_isolator_1{ + public: + typedef Polynomial_ Polynomial; + typedef Bound_ Bound; + private: + typedef Gmpfi Interval; + public: + RS23_k_isolator_1(const Polynomial&); + Polynomial polynomial()const; + int number_of_real_roots()const; + bool is_exact_root(int i)const; + Bound left_bound(int i)const; + Bound right_bound(int i)const; + private: + Polynomial _polynomial; + std::vector _real_roots; +}; + +template +RS23_k_isolator_1:: +RS23_k_isolator_1(const Polynomial_ &p){ + CGAL_error_msg("not implemented for these polynomial/bound types"); +} + +template <> +RS23_k_isolator_1,Gmpfr>:: +RS23_k_isolator_1(const CGAL::Polynomial &p):_polynomial(p){ + typedef CGAL::Polynomial Pol; + typedef CGAL::Gmpfr Bound; + typedef typename CGAL::RS3::RS3_k_refiner_1 KRefiner; + int numsols; + unsigned int degree=p.degree(); + mpz_t *coeffs=(mpz_t*)malloc((degree+1)*sizeof(mpz_t)); + std::vector intervals; + for(int i=0;i<=degree;++i) + coeffs[i][0]=*(p[i].mpz()); + RS2::RS2_calls::init_solver(); + RS2::RS2_calls::create_rs_upoly(coeffs,degree,rs_get_default_up()); + set_rs_precisol(0); + set_rs_verbose(0); + rs_run_algo((char*)"UISOLE"); + RS2::RS2_calls::insert_roots(std::back_inserter(intervals)); + // RS2 computed the isolating intervals. Now, we use RS3 to refine each + // root until reaching Kantorovich criterion, before adding it to the + // root vector. + numsols=intervals.size(); + for(int j=0;j +Polynomial_ +RS23_k_isolator_1::polynomial()const{ + return _polynomial; +} + +template +int +RS23_k_isolator_1::number_of_real_roots()const{ + return _real_roots.size(); +} + +template +bool +RS23_k_isolator_1::is_exact_root(int i)const{ + return _real_roots[i].inf()==_real_roots[i].sup(); +} + +template +Bound_ +RS23_k_isolator_1::left_bound(int i)const{ + return _real_roots[i].inf(); +} + +template +Bound_ +RS23_k_isolator_1::right_bound(int i)const{ + return _real_roots[i].sup(); +} + +} // namespace CGAL + +#endif // CGAL_RS_RS23_K_ISOLATOR_1_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/rs2_calls.h b/Algebraic_kernel_d/include/CGAL/RS/rs2_calls.h new file mode 100644 index 00000000000..60966539723 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/rs2_calls.h @@ -0,0 +1,134 @@ +// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author: Luis Peñaranda + +#ifndef CGAL_RS_RS2_CALLS_H +#define CGAL_RS_RS2_CALLS_H + +#include +#include +#include +#include + +#ifdef CGAL_RS_OLD_INCLUDES +#define CGALRS_PTR(a) long int a +#else +#define CGALRS_PTR(a) void *a +#endif + +namespace CGAL{ +namespace RS2{ + +struct RS2_calls{ + + static void init_solver(){ + static bool first=true; + if(first){ + first=false; + rs_init_rs(); + rs_reset_all(); + }else + rs_reset_all(); + } + + static void create_rs_upoly(mpz_t *poly, + const int deg, + CGALRS_PTR(ident_pol)){ + CGALRS_PTR(ident_mon); + CGALRS_PTR(ident_coeff); + rs_import_uppring((char*)"T"); + for(int i=0;i<=deg;++i) + if(mpz_sgn(poly[i])){ // don't add if == 0 + ident_mon=rs_export_new_mon_upp_bz(); + ident_coeff=rs_export_new_gmp(); + rs_import_bz_gmp + (ident_coeff,TO_RSPTR_IN(&(poly[i]))); + rs_dset_mon_upp_bz(ident_mon,ident_coeff,i); + rs_dappend_list_mon_upp_bz(ident_pol, + ident_mon); + } + } + + static int affiche_sols_eqs(mpfi_ptr *x){ + CGALRS_PTR(ident_sols_eqs); + CGALRS_PTR(ident_node); + CGALRS_PTR(ident_vect); + CGALRS_PTR(ident_elt); + int nb_elts; + ident_sols_eqs=rs_get_default_sols_eqs(); + nb_elts=rs_export_list_vect_ibfr_nb(ident_sols_eqs); + ident_node=rs_export_list_vect_ibfr_firstnode(ident_sols_eqs); + mpfi_t *roots=(mpfi_t*)malloc(nb_elts*sizeof(mpfi_t)); + for(int i=0;i + static OutputIterator insert_roots(OutputIterator x){ + CGALRS_PTR(ident_sols_eqs); + CGALRS_PTR(ident_node); + CGALRS_PTR(ident_vect); + CGALRS_PTR(ident_elt); + int nb_elts; + ident_sols_eqs=rs_get_default_sols_eqs(); + nb_elts=rs_export_list_vect_ibfr_nb(ident_sols_eqs); + ident_node=rs_export_list_vect_ibfr_firstnode(ident_sols_eqs); + for(int i=0;ileft),root_prec); + Gmpfr right(&(root_pointer->right),root_prec); + // Copy them, to have the data out of RS memory. + *x++=Gmpfi(left,right,root_prec+1); + ident_node=rs_export_list_vect_ibfr_nextnode + (ident_node); + } + return x; + } + +}; // struct RS2_calls + +} // namespace RS2 +} // namespace CGAL + +#endif // CGAL_RS_RS2_CALLS_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/rs2_isolator_1.h b/Algebraic_kernel_d/include/CGAL/RS/rs2_isolator_1.h new file mode 100644 index 00000000000..a09beefa5c1 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/rs2_isolator_1.h @@ -0,0 +1,105 @@ +// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author: Luis Peñaranda + +#ifndef CGAL_RS_RS2_ISOLATOR_1_H +#define CGAL_RS_RS2_ISOLATOR_1_H + +#include "rs2_calls.h" +#include +#include + +namespace CGAL{ +namespace RS2{ + +template +class RS2_isolator_1{ + public: + typedef Polynomial_ Polynomial; + typedef Bound_ Bound; + private: + typedef Gmpfi Interval; + public: + RS2_isolator_1(const Polynomial&); + Polynomial polynomial()const; + int number_of_real_roots()const; + bool is_exact_root(int i)const; + Bound left_bound(int i)const; + Bound right_bound(int i)const; + private: + Polynomial _polynomial; + std::vector _real_roots; +}; + +template +RS2_isolator_1:: +RS2_isolator_1(const Polynomial_ &p){ + CGAL_error_msg("not implemented for these polynomial/bound types"); +} + +template <> +RS2_isolator_1,Gmpfr>:: +RS2_isolator_1(const CGAL::Polynomial &p):_polynomial(p){ + unsigned int degree=p.degree(); + mpz_t *coeffs=(mpz_t*)malloc((degree+1)*sizeof(mpz_t)); + mpfi_ptr *intervals_mpfi=(mpfi_ptr*)malloc(degree*sizeof(mpfi_ptr)); + for(int i=0;i<=degree;++i) + coeffs[i][0]=*(p[i].mpz()); + RS2::RS2_calls::init_solver(); + RS2::RS2_calls::create_rs_upoly(coeffs,degree,rs_get_default_up()); + free(coeffs); + set_rs_precisol(0); + set_rs_verbose(0); + rs_run_algo((char*)"UISOLE"); + RS2::RS2_calls::insert_roots(std::back_inserter(_real_roots)); + free(intervals_mpfi); +} + +template +Polynomial_ +RS2_isolator_1::polynomial()const{ + return _polynomial; +} + +template +int +RS2_isolator_1::number_of_real_roots()const{ + return _real_roots.size(); +} + +template +bool +RS2_isolator_1::is_exact_root(int i)const{ + return _real_roots[i].inf()==_real_roots[i].sup(); +} + +template +Bound_ +RS2_isolator_1::left_bound(int i)const{ + return _real_roots[i].inf(); +} + +template +Bound_ +RS2_isolator_1::right_bound(int i)const{ + return _real_roots[i].sup(); +} + +} // namespace RS2 +} // namespace CGAL + +#endif // CGAL_RS_RS2_ISOLATOR_1_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/rs3_k_refiner_1.h b/Algebraic_kernel_d/include/CGAL/RS/rs3_k_refiner_1.h new file mode 100644 index 00000000000..11508c37915 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/rs3_k_refiner_1.h @@ -0,0 +1,108 @@ +// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author: Luis Peñaranda + +#ifndef CGAL_RS_RS3_K_REFINER_1_H +#define CGAL_RS_RS3_K_REFINER_1_H + +#include +#include "rs2_calls.h" +#include + +namespace CGAL{ +namespace RS3{ + +template +struct RS3_k_refiner_1{ + void operator()(const Polynomial_&,Bound_&,Bound_&,int); +}; // class RS3_k_refiner_1 + +template +void +RS3_k_refiner_1:: +operator()(const Polynomial_&,Bound_&,Bound_&,int){ + CGAL_error_msg("RS3 k-refiner not implemented for these types"); + return; +} + +template<> +void +RS3_k_refiner_1,Gmpfr>:: +operator() +(const Polynomial &pol,Gmpfr &left,Gmpfr &right,int prec){ + typedef Polynomial Polynomial; + typedef Polynomial_traits_d Ptraits; + typedef Ptraits::Degree Degree; + CGAL_precondition(left<=right); + // TODO: add precondition to check whether the interval is a point + // or the evaluations on its endpoints have different signs + //std::cout<<"refining ["< + +#ifndef CGAL_RS_RS3_REFINER_1_H +#define CGAL_RS_RS3_REFINER_1_H + +#include +#include "rs2_calls.h" +#include + +namespace CGAL{ +namespace RS3{ + +template +struct RS3_refiner_1{ + void operator()(const Polynomial_&,Bound_&,Bound_&,int,int=0); +}; // class RS3_refiner_1 + +template +void +RS3_refiner_1:: +operator()(const Polynomial_&,Bound_&,Bound_&,int,int){ + CGAL_error_msg("RS3 refiner not implemented for these types"); + return; +} + +template<> +void +RS3_refiner_1,Gmpfr>:: +operator() +(const Polynomial &pol,Gmpfr &left,Gmpfr &right,int prec,int k){ + typedef Polynomial Polynomial; + typedef Polynomial_traits_d Ptraits; + typedef Ptraits::Degree Degree; + CGAL_precondition(left<=right); + // TODO: add precondition to check whether the interval is a point + // or the evaluations on its endpoints have different signs + //std::cout<<"refining ["< + +#ifndef CGAL_RS_SIGNAT_1_H +#define CGAL_RS_SIGNAT_1_H + +#include +#include +//#include + +namespace CGAL{ +namespace RS_AK1{ + +template +struct Signat_1{ + typedef Polynomial_ Polynomial; + typedef Bound_ Bound; + typedef CGAL::Polynomial_traits_d PT; + typedef typename PT::Degree Degree; + Polynomial pol; + Signat_1(const Polynomial &p):pol(p){}; + CGAL::Sign operator()(const Bound&)const; +}; // struct Signat_1 + +template +inline CGAL::Sign +Signat_1::operator()(const Bound_ &x)const{ + typedef Polynomial_ Polynomial; + typedef Bound_ Bound; + typedef Real_embeddable_traits REtraits; + typedef typename REtraits::Sgn BSign; + typedef Algebraic_structure_traits AStraits; + // This generic signat works only when Bound_ is an exact type. For + // non-exact types, an implementation must be provided. + //BOOST_MPL_ASSERT((boost::is_same)); + int d=Degree()(pol); + Bound h(pol[d]); + for(int i=1;i<=d;++i) + h=h*x+pol[d-i]; + return BSign()(h); +} + +template <> +inline CGAL::Sign +Signat_1,Gmpfr>::operator()(const Gmpfr &x)const{ + typedef Signat_1 Exact_sign; + // This seems to work faster for small polynomials: + // return Exact_sign(pol)(x); + int d=Degree()(pol); + if(d==0) + return pol[0].sign(); + Gmpfi h(pol[d],x.get_precision()+2*d); + Uncertain indet=Uncertain::indeterminate(); + if(h.sign().is_same(indet)) + return Exact_sign(pol)(x); + for(int i=1;i<=d;++i){ + h*=x; + h+=pol[d-i]; + if(h.sign().is_same(indet)) + return Exact_sign(pol)(x); + } + CGAL_assertion(!h.sign().is_same(indet)); + return h.sign(); +} + +} // namespace RS_AK1 +} // namespace CGAL + +#endif // CGAL_RS_SIGNAT_1_H