diff --git a/Algebraic_kernel_d/include/CGAL/RS/ak_z_1.h b/Algebraic_kernel_d/include/CGAL/RS/ak_z_1.h new file mode 100644 index 00000000000..fad00697678 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/ak_z_1.h @@ -0,0 +1,211 @@ +// 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_Z_1_H +#define CGAL_RS_AK_Z_1_H + +#include // included only to define size_t +#include +#include "algebraic_z_1.h" +#include "comparator_1.h" +#include "signat_1.h" +#include "functors_z_1.h" + +// This file defines the "Z-algebraic kernel". This kind of kernel makes +// all the internal operations with an integer polynomial type (the name +// "Z" comes from there). For this, a converter functor (passed as a +// template parameter) is used, which converts the input polynomial to the +// integer representation. + +namespace CGAL{ +namespace RS_AK1{ + +template , + class ZPtraits_=CGAL::Polynomial_traits_d > +class Algebraic_kernel_1{ + public: + typedef ExtPolynomial_ Polynomial_1; + typedef IntPolynomial_ ZPolynomial_1; + typedef PolConverter_ PolConverter; + typedef typename Polynomial_1::NT Coefficient; + typedef Bound_ Bound; + private: + typedef Isolator_ Isolator; + typedef Refiner_ Refiner; + typedef Ptraits_ Ptraits; + typedef ZPtraits_ ZPtraits; + typedef CGAL::RS_AK1::Signat_1 + Signat; + typedef CGAL::RS_AK1::Simple_comparator_1 + Comparator; + public: + typedef CGAL::RS_AK1::Algebraic_z_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_Z_1_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/algebraic_z_1.h b/Algebraic_kernel_d/include/CGAL/RS/algebraic_z_1.h new file mode 100644 index 00000000000..f75bb84b8f4 --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/algebraic_z_1.h @@ -0,0 +1,353 @@ +// 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_ALGEBRAIC_Z_1_H +#define CGAL_RS_ALGEBRAIC_Z_1_H + +#include "algebraic_1.h" + +namespace CGAL{ +namespace RS_AK1{ + +template +class Algebraic_z_1: +public Algebraic_1, +boost::totally_ordered, + double>{ + protected: + typedef Polynomial_ Polynomial; + typedef ZPolynomial_ ZPolynomial; + typedef Bound_ Bound; + typedef ZRefiner_ ZRefiner; + typedef ZComparator_ ZComparator; + typedef ZPtraits_ ZPtraits; + typedef typename ZPtraits::Coefficient_type ZCoefficient; + typedef typename ZPtraits::Scale ZScale; + typedef Ptraits_ Ptraits; + typedef typename Ptraits::Coefficient_type Coefficient; + typedef typename Ptraits::Scale Scale; + typedef Algebraic_1 + Algebraic_base; + typedef Algebraic_z_1 + Algebraic; + ZPolynomial zpol; + + public: + Algebraic_z_1(){}; + Algebraic_z_1(const Polynomial &p, + const ZPolynomial &zp, + const Bound &l, + const Bound &r):Algebraic_base(p,l,r),zpol(zp){ + CGAL_assertion(l<=r); + } + Algebraic_z_1(const Algebraic &a): + Algebraic_base(a.pol,a.left,a.right),zpol(a.zpol){} + Algebraic_z_1(const Bound &b){ + typedef typename Ptraits::Shift Shift; + typedef typename ZPtraits::Shift ZShift; + this->left=b; + this->right=b; + Gmpq q(b); + this->pol=Coefficient(mpq_denref(q.mpq()))* + Shift()(Polynomial(1),1,0)- + Coefficient(mpq_numref(q.mpq())); + zpol=ZCoefficient(mpq_denref(q.mpq()))* + ZShift()(ZPolynomial(1),1,0)- + ZCoefficient(mpq_numref(q.mpq())); + CGAL_assertion(this->left==this->right); + CGAL_assertion(this->left==b); + } + template + Algebraic_z_1(const T &t){ + typedef typename Ptraits::Shift Shift; + typedef typename ZPtraits::Shift ZShift; + CGAL::Gmpq q(t); + this->pol=Coefficient(mpq_denref(q.mpq()))* + Shift()(Polynomial(1),1,0)- + Coefficient(mpq_numref(q.mpq())); + zpol=ZCoefficient(mpq_denref(q.mpq()))* + ZShift()(ZPolynomial(1),1,0)- + ZCoefficient(mpq_numref(q.mpq())); + this->left=Bound(t,std::round_toward_neg_infinity); + this->right=Bound(t,std::round_toward_infinity); + CGAL_assertion(this->left<=t&&this->right>=t); + } + Algebraic_z_1(const CGAL::Gmpq &q){ + typedef typename Ptraits::Shift Shift; + typedef typename ZPtraits::Shift ZShift; + this->pol=Coefficient(mpq_denref(q.mpq()))* + Shift()(Polynomial(1),1,0)- + Coefficient(mpq_numref(q.mpq())); + zpol=ZCoefficient(mpq_denref(q.mpq()))* + ZShift()(ZPolynomial(1),1,0)- + ZCoefficient(mpq_numref(q.mpq())); + this->left=Bound(); + this->right=Bound(); + mpfr_t b; + mpfr_init(b); + mpfr_set_q(b,q.mpq(),GMP_RNDD); + mpfr_swap(b,this->left.fr()); + mpfr_set_q(b,q.mpq(),GMP_RNDU); + mpfr_swap(b,this->right.fr()); + mpfr_clear(b); + CGAL_assertion(this->left<=q&&this->right>=q); + } + ~Algebraic_z_1(){} + + ZPolynomial get_zpol()const{return zpol;} + void set_zpol(const ZPolynomial &p){zpol=p;} + + Algebraic operator-()const{ + return Algebraic(Scale()(this->get_pol(),Coefficient(-1)), + ZScale()(get_zpol(),ZCoefficient(-1)), + -this->right, + -this->left); + } + +#define CGAL_RS_COMPARE_ALGEBRAIC_Z(_a) \ + (ZComparator()(get_zpol(),this->get_left(),this->get_right(), \ + (_a).get_zpol(),(_a).get_left(),(_a).get_right())) + + Comparison_result compare(Algebraic a)const{ + return CGAL_RS_COMPARE_ALGEBRAIC_Z(a); + }; + +#define CGAL_RS_COMPARE_ALGEBRAIC_Z_TYPE(_t) \ + bool operator<(_t t)const \ + {Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::SMALLER;} \ + bool operator>(_t t)const \ + {Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::LARGER;} \ + bool operator==(_t t)const \ + {Algebraic a(t);return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::EQUAL;} + + bool operator==(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::EQUAL;} + bool operator!=(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)!=CGAL::EQUAL;} + bool operator<(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::SMALLER;} + bool operator<=(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)!=CGAL::LARGER;} + bool operator>(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)==CGAL::LARGER;} + bool operator>=(Algebraic a)const + {return CGAL_RS_COMPARE_ALGEBRAIC_Z(a)!=CGAL::SMALLER;} + + CGAL_RS_COMPARE_ALGEBRAIC_Z_TYPE(double) + +#undef CGAL_RS_COMPARE_ALGEBRAIC_Z_TYPE +#undef CGAL_RS_COMPARE_ALGEBRAIC_Z + +#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; + ZRefiner()(get_zpol(), + this->get_left(), + this->get_right(), + CGAL_RS_DBL_PREC); + CGAL_assertion(TD()(this->get_left())==TD()(this->get_right())); + return TD()(this->get_left()); + } + std::pair to_interval()const{ + typedef Real_embeddable_traits RT; + typedef typename RT::To_interval TI; + return std::make_pair(TI()(this->get_left()).first, + TI()(this->get_right()).second); + } +#undef CGAL_RS_DBL_PREC +}; // class Algebraic_z_1 + +} // namespace RS_AK1 + +// We define Algebraic_z_1 as real-embeddable +template +class Real_embeddable_traits >: +public INTERN_RET::Real_embeddable_traits_base< + RS_AK1::Algebraic_z_1, + CGAL::Tag_true>{ + typedef Polynomial_ P; + typedef ZPolynomial_ ZP; + typedef Bound_ B; + typedef ZRefiner_ R; + typedef ZComparator_ C; + typedef Ptraits_ T; + typedef ZPtraits_ ZT; + + public: + + typedef RS_AK1::Algebraic_z_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; + } + }; +}; + +template +inline +RS_AK1::Algebraic_z_1 min +BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_z_1 a, + RS_AK1::Algebraic_z_1 b){ + return(a +inline +RS_AK1::Algebraic_z_1 max +BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_z_1 a, + RS_AK1::Algebraic_z_1 b){ + return(a>b?a:b); +} + +template +inline +std::ostream& operator<<(std::ostream &o, + const RS_AK1::Algebraic_z_1 &a){ + return(o<<'['< +inline +std::istream& operator>>(std::istream &i, + RS_AK1::Algebraic_z_1 &a){ + std::istream::int_type c; + P pol; + ZP zpol; + 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>>zpol; + 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_z_1(pol,zpol,lb,rb); + return i; +} + +} // namespace CGAL + +#endif // CGAL_RS_ALGEBRAIC_Z_1_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/functors_z_1.h b/Algebraic_kernel_d/include/CGAL/RS/functors_z_1.h new file mode 100644 index 00000000000..6534a38f76e --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/RS/functors_z_1.h @@ -0,0 +1,709 @@ +// 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_FUNCTORS_Z_1_H +#define CGAL_RS_FUNCTORS_Z_1_H + +#include +#include + +namespace CGAL{ +namespace RS_AK1{ + +template +struct Construct_algebraic_real_1{ + typedef Polynomial_ Polynomial; + typedef ZPolynomial_ ZPolynomial; + typedef PolConverter_ PolConverter; + typedef Algebraic_ Algebraic; + typedef Bound_ Bound; + typedef Coefficient_ Coefficient; + typedef Isolator_ Isolator; + typedef Algebraic result_type; + + template + Algebraic operator()(const T &a)const{ + return Algebraic(a); + } + + Algebraic operator()(const Polynomial &p,size_t i)const{ + ZPolynomial zp=PolConverter()(p); + Isolator isol(zp); + return Algebraic(p, + zp, + isol.left_bound(i), + isol.right_bound(i)); + } + + Algebraic operator()(const Polynomial &p, + const Bound &l, + const Bound &r)const{ + return Algebraic(p,PolConverter()(p),l,r); + } + +}; // struct Construct_algebraic_1 + +template +struct Compute_polynomial_1: +public std::unary_function{ + typedef Polynomial_ Polynomial; + typedef Algebraic_ Algebraic; + Polynomial operator()(const Algebraic &x)const{ + return x.get_pol(); + } +}; // struct Compute_polynomial_1 + +template +struct Is_coprime_1: +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; + } +}; // struct Is_coprime_1 + +template +struct Make_coprime_1{ + 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; + } +}; // struct Make_coprime_1 + +template +struct Solve_1{ + typedef Polynomial_ Polynomial_1; + typedef ZPolynomial_ ZPolynomial_1; + typedef PolConverter_ PolConverter; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Isolator_ Isolator; + typedef Signat_ ZSignat; + 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; + typedef ZPtraits_ ZPtraits; + typedef typename ZPtraits::Gcd_up_to_constant_factor ZGcd; + typedef typename ZPtraits::Square_free_factorize_up_to_constant_factor + ZSqfr; + typedef typename ZPtraits::Degree ZDegree; + typedef typename ZPtraits::Make_square_free ZSfpart; + + template + OutputIterator operator()(const Polynomial_1 &p, + OutputIterator res)const{ + typedef std::pair polmult; + typedef std::vector sqvec; + typedef std::pair zpolmult; + typedef std::vector zsqvec; + + ZPolynomial_1 zp=PolConverter()(p); + Polynomial_1 sfp=Sfpart()(p); + ZPolynomial_1 zsfp=PolConverter()(sfp); + zsqvec zsfv; + ZSqfr()(zp,std::back_inserter(zsfv)); + Isolator isol(zsfp); + int *m=(int*)calloc(isol.number_of_real_roots(),sizeof(int)); + for(typename zsqvec::iterator i=zsfv.begin();i!=zsfv.end();++i){ + int k=ZDegree()(i->first); + ZSignat 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{ + ZPolynomial_1 zp=PolConverter()(p); + Isolator isol(zp); + 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 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: + size_t 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){ + ZPolynomial_1 zifirst=PolConverter()(i->first); + int k=ZDegree()(zifirst); + ZSignat signof(zifirst); + for(size_t j=index_l;k&&jsecond; + --k; + } + } + } + } + for(size_t 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 ZPolynomial_ ZPolynomial_1; + typedef PolConverter_ PolConverter; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Refiner_ Refiner; + typedef Signat_ ZSignat; + typedef Ptraits_ Ptraits; + typedef ZPtraits_ ZPtraits; + + 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_zpol(), + 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::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_zpol(), + 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; + ZSignat sign_at_gcd(PolConverter()(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 ZPolynomial_ ZPolynomial_1; + typedef PolConverter_ PolConverter; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Refiner_ Refiner; + typedef Signat_ ZSignat; + typedef Ptraits_ Ptraits; + typedef ZPtraits_ ZPtraits; + + 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::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_zpol(), + 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; + ZSignat sign_at_gcd(PolConverter()(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 ZPolynomial_ ZPolynomial_1; + typedef PolConverter_ PolConverter; + typedef Isolator_ Isolator; + size_t operator()(const Polynomial_1 &p)const{ + // TODO: make sure that p is square free (precondition) + Isolator isol(PolConverter()(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_zpol(),al,ar, + b.get_zpol(),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_zpol(),al,ar, + balg.get_zpol(),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_zpol(), + al, + ar, + balg.get_zpol(), + 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 ZPolynomial_ ZPolynomial_1; + typedef PolConverter_ PolConverter; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Isolator_ Isolator; + typedef Comparator_ Comparator; + typedef Signat_ Signat; + typedef Ptraits_ Ptraits; + typedef ZPtraits_ ZPtraits; + + 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 >{ + typedef Polynomial_ Polynomial_1; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Refiner_ Refiner; + + // 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_zpol(),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 + +template +struct Approximate_relative_1: +public std::binary_function >{ + typedef Polynomial_ Polynomial_1; + typedef Bound_ Bound; + typedef Algebraic_ Algebraic; + typedef Refiner_ Refiner; + + 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_zpol(), + 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 + +} // namespace RS_AK1 +} // namespace CGAL + +#endif // CGAL_RS_FUNCTORS_Z_1_H