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