mirror of https://github.com/CGAL/cgal
added the z-kernel, for types different than Gmpz
This commit is contained in:
Luis Peñaranda
parent
89c5170ee0
commit
d33d8bf2cc
|
|
@ -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 <luis.penaranda@gmx.com>
|
||||||
|
|
||||||
|
#ifndef CGAL_RS_AK_Z_1_H
|
||||||
|
#define CGAL_RS_AK_Z_1_H
|
||||||
|
|
||||||
|
#include <cstddef> // included only to define size_t
|
||||||
|
#include <CGAL/Polynomial_traits_d.h>
|
||||||
|
#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 ExtPolynomial_,
|
||||||
|
class IntPolynomial_,
|
||||||
|
class PolConverter_,
|
||||||
|
class Bound_,
|
||||||
|
class Isolator_,
|
||||||
|
class Refiner_,
|
||||||
|
class Ptraits_=CGAL::Polynomial_traits_d<ExtPolynomial_>,
|
||||||
|
class ZPtraits_=CGAL::Polynomial_traits_d<IntPolynomial_> >
|
||||||
|
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<ZPolynomial_1,Bound>
|
||||||
|
Signat;
|
||||||
|
typedef CGAL::RS_AK1::Simple_comparator_1<ZPolynomial_1,
|
||||||
|
Bound,
|
||||||
|
Refiner,
|
||||||
|
Signat,
|
||||||
|
ZPtraits>
|
||||||
|
Comparator;
|
||||||
|
public:
|
||||||
|
typedef CGAL::RS_AK1::Algebraic_z_1<Polynomial_1,
|
||||||
|
ZPolynomial_1,
|
||||||
|
Bound,
|
||||||
|
Refiner,
|
||||||
|
Comparator,
|
||||||
|
Ptraits,
|
||||||
|
ZPtraits>
|
||||||
|
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<Polynomial_1,
|
||||||
|
ZPolynomial_1,
|
||||||
|
PolConverter,
|
||||||
|
Algebraic_real_1,
|
||||||
|
Bound,
|
||||||
|
Coefficient,
|
||||||
|
Isolator>
|
||||||
|
Construct_algebraic_real_1;
|
||||||
|
typedef CGAL::RS_AK1::Compute_polynomial_1<Polynomial_1,
|
||||||
|
Algebraic_real_1>
|
||||||
|
Compute_polynomial_1;
|
||||||
|
typedef CGAL::RS_AK1::Isolate_1<Polynomial_1,
|
||||||
|
ZPolynomial_1,
|
||||||
|
PolConverter,
|
||||||
|
Bound,
|
||||||
|
Algebraic_real_1,
|
||||||
|
Isolator,
|
||||||
|
Comparator,
|
||||||
|
Signat,
|
||||||
|
Ptraits,
|
||||||
|
ZPtraits>
|
||||||
|
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<Polynomial_1,Ptraits>
|
||||||
|
Is_coprime_1;
|
||||||
|
typedef CGAL::RS_AK1::Make_coprime_1<Polynomial_1,Ptraits>
|
||||||
|
Make_coprime_1;
|
||||||
|
typedef CGAL::RS_AK1::Solve_1<Polynomial_1,
|
||||||
|
ZPolynomial_1,
|
||||||
|
PolConverter,
|
||||||
|
Bound,
|
||||||
|
Algebraic_real_1,
|
||||||
|
Isolator,
|
||||||
|
Signat,
|
||||||
|
Ptraits,
|
||||||
|
ZPtraits>
|
||||||
|
Solve_1;
|
||||||
|
typedef CGAL::RS_AK1::Number_of_solutions_1<Polynomial_1,
|
||||||
|
ZPolynomial_1,
|
||||||
|
PolConverter,
|
||||||
|
Isolator>
|
||||||
|
Number_of_solutions_1;
|
||||||
|
|
||||||
|
typedef CGAL::RS_AK1::Sign_at_1<Polynomial_1,
|
||||||
|
ZPolynomial_1,
|
||||||
|
PolConverter,
|
||||||
|
Bound,
|
||||||
|
Algebraic_real_1,
|
||||||
|
Refiner,
|
||||||
|
Signat,
|
||||||
|
Ptraits,
|
||||||
|
ZPtraits>
|
||||||
|
Sign_at_1;
|
||||||
|
typedef CGAL::RS_AK1::Is_zero_at_1<Polynomial_1,
|
||||||
|
ZPolynomial_1,
|
||||||
|
PolConverter,
|
||||||
|
Bound,
|
||||||
|
Algebraic_real_1,
|
||||||
|
Refiner,
|
||||||
|
Signat,
|
||||||
|
Ptraits,
|
||||||
|
ZPtraits>
|
||||||
|
Is_zero_at_1;
|
||||||
|
typedef CGAL::RS_AK1::Compare_1<Algebraic_real_1,
|
||||||
|
Bound,
|
||||||
|
Comparator>
|
||||||
|
Compare_1;
|
||||||
|
typedef CGAL::RS_AK1::Bound_between_1<Algebraic_real_1,
|
||||||
|
Bound,
|
||||||
|
Comparator>
|
||||||
|
Bound_between_1;
|
||||||
|
typedef CGAL::RS_AK1::Approximate_absolute_1<Polynomial_1,
|
||||||
|
Bound,
|
||||||
|
Algebraic_real_1,
|
||||||
|
Refiner>
|
||||||
|
Approximate_absolute_1;
|
||||||
|
typedef CGAL::RS_AK1::Approximate_relative_1<Polynomial_1,
|
||||||
|
Bound,
|
||||||
|
Algebraic_real_1,
|
||||||
|
Refiner>
|
||||||
|
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
|
||||||
|
|
@ -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 <luis.penaranda@gmx.com>
|
||||||
|
|
||||||
|
#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 Polynomial_,
|
||||||
|
class ZPolynomial_,
|
||||||
|
class Bound_,
|
||||||
|
class ZRefiner_,
|
||||||
|
class ZComparator_,
|
||||||
|
class Ptraits_,
|
||||||
|
class ZPtraits_>
|
||||||
|
class Algebraic_z_1:
|
||||||
|
public Algebraic_1<Polynomial_,Bound_,ZRefiner_,ZComparator_,ZPtraits_>,
|
||||||
|
boost::totally_ordered<Algebraic_z_1<Polynomial_,
|
||||||
|
ZPolynomial_,
|
||||||
|
Bound_,
|
||||||
|
ZRefiner_,
|
||||||
|
ZComparator_,
|
||||||
|
Ptraits_,
|
||||||
|
ZPtraits_>,
|
||||||
|
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<Polynomial,Bound,ZRefiner,ZComparator,ZPtraits>
|
||||||
|
Algebraic_base;
|
||||||
|
typedef Algebraic_z_1<Polynomial,
|
||||||
|
ZPolynomial,
|
||||||
|
Bound,
|
||||||
|
ZRefiner,
|
||||||
|
ZComparator,
|
||||||
|
Ptraits,
|
||||||
|
ZPtraits>
|
||||||
|
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 <class T>
|
||||||
|
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<Bound> 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<double,double> to_interval()const{
|
||||||
|
typedef Real_embeddable_traits<Bound> 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 Polynomial_,
|
||||||
|
class ZPolynomial_,
|
||||||
|
class Bound_,
|
||||||
|
class ZRefiner_,
|
||||||
|
class ZComparator_,
|
||||||
|
class Ptraits_,
|
||||||
|
class ZPtraits_>
|
||||||
|
class Real_embeddable_traits<RS_AK1::Algebraic_z_1<Polynomial_,
|
||||||
|
ZPolynomial_,
|
||||||
|
Bound_,
|
||||||
|
ZRefiner_,
|
||||||
|
ZComparator_,
|
||||||
|
Ptraits_,
|
||||||
|
ZPtraits_> >:
|
||||||
|
public INTERN_RET::Real_embeddable_traits_base<
|
||||||
|
RS_AK1::Algebraic_z_1<Polynomial_,
|
||||||
|
ZPolynomial_,
|
||||||
|
Bound_,
|
||||||
|
ZRefiner_,
|
||||||
|
ZComparator_,
|
||||||
|
Ptraits_,
|
||||||
|
ZPtraits_>,
|
||||||
|
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<P,ZP,B,R,C,T,ZT> 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<Type,CGAL::Tag_true>
|
||||||
|
Base;
|
||||||
|
typedef typename Base::Compare Compare;
|
||||||
|
|
||||||
|
class Sgn:public std::unary_function<Type,CGAL::Sign>{
|
||||||
|
public:
|
||||||
|
CGAL::Sign operator()(const Type &a)const{
|
||||||
|
return Compare()(a,Type(0));
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
class To_double:public std::unary_function<Type,double>{
|
||||||
|
public:
|
||||||
|
double operator()(Type a)const{return a.to_double();}
|
||||||
|
};
|
||||||
|
|
||||||
|
class To_interval:
|
||||||
|
public std::unary_function<Type,std::pair<double,double> >{
|
||||||
|
public:
|
||||||
|
std::pair<double,double> operator()(const Type &a)const{
|
||||||
|
return a.to_interval();
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
class Is_zero:public std::unary_function<Type,Boolean>{
|
||||||
|
public:
|
||||||
|
bool operator()(const Type &a)const{
|
||||||
|
return Sgn()(a)==CGAL::ZERO;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
class Is_finite:public std::unary_function<Type,Boolean>{
|
||||||
|
public:
|
||||||
|
bool operator()(const Type&)const{return true;}
|
||||||
|
};
|
||||||
|
|
||||||
|
class Abs:public std::unary_function<Type,Type>{
|
||||||
|
public:
|
||||||
|
Type operator()(const Type &a)const{
|
||||||
|
return Sgn()(a)==CGAL::NEGATIVE?-a:a;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
};
|
||||||
|
|
||||||
|
template <class P,class ZP,class B,class R,class C,class T,class ZT>
|
||||||
|
inline
|
||||||
|
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> min
|
||||||
|
BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> a,
|
||||||
|
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> b){
|
||||||
|
return(a<b?a:b);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class P,class ZP,class B,class R,class C,class T,class ZT>
|
||||||
|
inline
|
||||||
|
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> max
|
||||||
|
BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> a,
|
||||||
|
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> b){
|
||||||
|
return(a>b?a:b);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class P,class ZP,class B,class R,class C,class T,class ZT>
|
||||||
|
inline
|
||||||
|
std::ostream& operator<<(std::ostream &o,
|
||||||
|
const RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> &a){
|
||||||
|
return(o<<'['<<a.get_pol()<<','<<
|
||||||
|
a.get_zpol()<<','<<
|
||||||
|
a.get_left()<<','<<
|
||||||
|
a.get_right()<<']');
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class P,class ZP,class B,class R,class C,class T,class ZT>
|
||||||
|
inline
|
||||||
|
std::istream& operator>>(std::istream &i,
|
||||||
|
RS_AK1::Algebraic_z_1<P,ZP,B,R,C,T,ZT> &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<P,ZP,B,R,C,T,ZT>(pol,zpol,lb,rb);
|
||||||
|
return i;
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace CGAL
|
||||||
|
|
||||||
|
#endif // CGAL_RS_ALGEBRAIC_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 <luis.penaranda@gmx.com>
|
||||||
|
|
||||||
|
#ifndef CGAL_RS_FUNCTORS_Z_1_H
|
||||||
|
#define CGAL_RS_FUNCTORS_Z_1_H
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
#include <CGAL/Gmpfi.h>
|
||||||
|
|
||||||
|
namespace CGAL{
|
||||||
|
namespace RS_AK1{
|
||||||
|
|
||||||
|
template <class Polynomial_,
|
||||||
|
class ZPolynomial_,
|
||||||
|
class PolConverter_,
|
||||||
|
class Algebraic_,
|
||||||
|
class Bound_,
|
||||||
|
class Coefficient_,
|
||||||
|
class Isolator_>
|
||||||
|
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 <class T>
|
||||||
|
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 <class Polynomial_,class Algebraic_>
|
||||||
|
struct Compute_polynomial_1:
|
||||||
|
public std::unary_function<Algebraic_,Polynomial_>{
|
||||||
|
typedef Polynomial_ Polynomial;
|
||||||
|
typedef Algebraic_ Algebraic;
|
||||||
|
Polynomial operator()(const Algebraic &x)const{
|
||||||
|
return x.get_pol();
|
||||||
|
}
|
||||||
|
}; // struct Compute_polynomial_1
|
||||||
|
|
||||||
|
template <class Polynomial_,class Ptraits_>
|
||||||
|
struct Is_coprime_1:
|
||||||
|
public std::binary_function<Polynomial_,Polynomial_,bool>{
|
||||||
|
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 <class Polynomial_,class Ptraits_>
|
||||||
|
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 <class Polynomial_,
|
||||||
|
class ZPolynomial_,
|
||||||
|
class PolConverter_,
|
||||||
|
class Bound_,
|
||||||
|
class Algebraic_,
|
||||||
|
class Isolator_,
|
||||||
|
class Signat_,
|
||||||
|
class Ptraits_,
|
||||||
|
class ZPtraits_>
|
||||||
|
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 <class OutputIterator>
|
||||||
|
OutputIterator operator()(const Polynomial_1 &p,
|
||||||
|
OutputIterator res)const{
|
||||||
|
typedef std::pair<Polynomial_1,int> polmult;
|
||||||
|
typedef std::vector<polmult> sqvec;
|
||||||
|
typedef std::pair<ZPolynomial_1,int> zpolmult;
|
||||||
|
typedef std::vector<zpolmult> 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&&j<isol.number_of_real_roots();++j){
|
||||||
|
if(!m[j]){
|
||||||
|
CGAL::Sign sg_l=
|
||||||
|
signof(isol.left_bound(j));
|
||||||
|
CGAL::Sign sg_r=
|
||||||
|
signof(isol.right_bound(j));
|
||||||
|
if((sg_l!=sg_r)||
|
||||||
|
((sg_l==CGAL::ZERO)&&
|
||||||
|
(sg_r==CGAL::ZERO))){
|
||||||
|
m[j]=i->second;
|
||||||
|
--k;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
for(int l=0;l<isol.number_of_real_roots();++l)
|
||||||
|
*res++=std::make_pair(Algebraic(sfp,
|
||||||
|
zsfp,
|
||||||
|
isol.left_bound(l),
|
||||||
|
isol.right_bound(l)),
|
||||||
|
m[l]);
|
||||||
|
free(m);
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class OutputIterator>
|
||||||
|
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<isol.number_of_real_roots();++l)
|
||||||
|
*res++=Algebraic(p,
|
||||||
|
zp,
|
||||||
|
isol.left_bound(l),
|
||||||
|
isol.right_bound(l));
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class OutputIterator>
|
||||||
|
OutputIterator operator()(const Polynomial_1 &p,
|
||||||
|
const Bound &l,
|
||||||
|
const Bound &u,
|
||||||
|
OutputIterator res)const{
|
||||||
|
typedef std::vector<Algebraic> RV;
|
||||||
|
typedef std::pair<Polynomial_1,int> PM;
|
||||||
|
typedef std::vector<PM> 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<nb_roots&&roots[index_l]<l)
|
||||||
|
++index_l;
|
||||||
|
CGAL_assertion(index_l<=nb_roots);
|
||||||
|
if(index_l==nb_roots)
|
||||||
|
return res;
|
||||||
|
index_u=index_l;
|
||||||
|
while(index_u<nb_roots&&roots[index_u]<u)
|
||||||
|
++index_u;
|
||||||
|
CGAL_assertion(index_u<=nb_roots);
|
||||||
|
if(index_u==index_l)
|
||||||
|
return res;
|
||||||
|
// now, we have to return roots in [index_l,index_u)
|
||||||
|
PMV sfv;
|
||||||
|
Sqfr()(p,std::back_inserter(sfv)); // square-free fact. of p
|
||||||
|
// array to store the multiplicities
|
||||||
|
int *m=(int*)calloc(nb_roots,sizeof(int));
|
||||||
|
// we iterate over all the pairs <root,mult> 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&&j<index_u;++j){
|
||||||
|
if(!m[j]){
|
||||||
|
CGAL::Sign sg_l=
|
||||||
|
signof(roots[j].get_left());
|
||||||
|
CGAL::Sign sg_r=
|
||||||
|
signof(roots[j].get_right());
|
||||||
|
if((sg_l!=sg_r)||
|
||||||
|
((sg_l==CGAL::ZERO)&&
|
||||||
|
(sg_r==CGAL::ZERO))){
|
||||||
|
m[j]=i->second;
|
||||||
|
--k;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
for(size_t l=index_l;l<index_u;++l)
|
||||||
|
*res++=std::make_pair(roots[l],m[l]);
|
||||||
|
free(m);
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class OutputIterator>
|
||||||
|
OutputIterator operator()(const Polynomial_1 &p,
|
||||||
|
bool known_to_be_square_free,
|
||||||
|
const Bound &l,
|
||||||
|
const Bound &u,
|
||||||
|
OutputIterator res)const{
|
||||||
|
typedef std::vector<Algebraic> 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 Polynomial_,
|
||||||
|
class ZPolynomial_,
|
||||||
|
class PolConverter_,
|
||||||
|
class Bound_,
|
||||||
|
class Algebraic_,
|
||||||
|
class Refiner_,
|
||||||
|
class Signat_,
|
||||||
|
class Ptraits_,
|
||||||
|
class ZPtraits_>
|
||||||
|
class Sign_at_1:
|
||||||
|
public std::binary_function<Polynomial_,Algebraic_,CGAL::Sign>{
|
||||||
|
// 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<CGAL::Sign> eval_interv(const Polynomial_1 &p,
|
||||||
|
const Bound &l,
|
||||||
|
const Bound &r)const{
|
||||||
|
typedef typename Ptraits::Substitute Subst;
|
||||||
|
std::vector<CGAL::Gmpfi> 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<CGAL::Sign> 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<CGAL::Sign>::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<CGAL::Sign> unknown=
|
||||||
|
Uncertain<CGAL::Sign>::indeterminate();
|
||||||
|
CGAL::Uncertain<CGAL::Sign> 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 Polynomial_,
|
||||||
|
class ZPolynomial_,
|
||||||
|
class PolConverter_,
|
||||||
|
class Bound_,
|
||||||
|
class Algebraic_,
|
||||||
|
class Refiner_,
|
||||||
|
class Signat_,
|
||||||
|
class Ptraits_,
|
||||||
|
class ZPtraits_>
|
||||||
|
class Is_zero_at_1:
|
||||||
|
public std::binary_function<Polynomial_,Algebraic_,bool>{
|
||||||
|
// 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<CGAL::Sign> eval_interv(const Polynomial_1 &p,
|
||||||
|
const Bound &l,
|
||||||
|
const Bound &r)const{
|
||||||
|
typedef typename Ptraits::Substitute Subst;
|
||||||
|
std::vector<CGAL::Gmpfi> 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<CGAL::Sign> unknown=
|
||||||
|
Uncertain<CGAL::Sign>::indeterminate();
|
||||||
|
CGAL::Uncertain<CGAL::Sign> 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 <class Polynomial_,
|
||||||
|
class ZPolynomial_,
|
||||||
|
class PolConverter_,
|
||||||
|
class Isolator_>
|
||||||
|
struct Number_of_solutions_1:
|
||||||
|
public std::unary_function<Polynomial_,int>{
|
||||||
|
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 <class Algebraic_,
|
||||||
|
class Bound_,
|
||||||
|
class Comparator_>
|
||||||
|
struct Compare_1:
|
||||||
|
public std::binary_function<Algebraic_,Algebraic_,CGAL::Comparison_result>{
|
||||||
|
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 <class T>
|
||||||
|
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 <class Algebraic_,
|
||||||
|
class Bound_,
|
||||||
|
class Comparator_>
|
||||||
|
struct Bound_between_1:
|
||||||
|
public std::binary_function<Algebraic_,Algebraic_,Bound_>{
|
||||||
|
typedef Algebraic_ Algebraic;
|
||||||
|
typedef Bound_ Bound;
|
||||||
|
typedef Comparator_ Comparator;
|
||||||
|
|
||||||
|
Bound operator()(Algebraic a,Algebraic b)const{
|
||||||
|
typedef Compare_1<Algebraic,Bound,Comparator> Compare;
|
||||||
|
typename Bound::Precision_type prec;
|
||||||
|
switch(Compare()(a,b)){
|
||||||
|
case CGAL::LARGER:
|
||||||
|
CGAL_assertion(b.get_right()<a.get_left());
|
||||||
|
prec=CGAL::max(b.get_right().get_precision(),
|
||||||
|
a.get_left().get_precision());
|
||||||
|
return Bound::add(b.get_right(),
|
||||||
|
a.get_left(),
|
||||||
|
1+prec)/2;
|
||||||
|
break;
|
||||||
|
case CGAL::SMALLER:
|
||||||
|
CGAL_assertion(a.get_right()<b.get_left());
|
||||||
|
prec=CGAL::max(a.get_right().get_precision(),
|
||||||
|
b.get_left().get_precision());
|
||||||
|
return Bound::add(a.get_right(),
|
||||||
|
b.get_left(),
|
||||||
|
1+prec)/2;
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
CGAL_error_msg(
|
||||||
|
"bound between two equal numbers");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}; // struct Bound_between_1
|
||||||
|
|
||||||
|
template <class Polynomial_,
|
||||||
|
class ZPolynomial_,
|
||||||
|
class PolConverter_,
|
||||||
|
class Bound_,
|
||||||
|
class Algebraic_,
|
||||||
|
class Isolator_,
|
||||||
|
class Comparator_,
|
||||||
|
class Signat_,
|
||||||
|
class Ptraits_,
|
||||||
|
class ZPtraits_>
|
||||||
|
struct Isolate_1:
|
||||||
|
public std::binary_function<Algebraic_,Polynomial_,std::pair<Bound_,Bound_> >{
|
||||||
|
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<Bound,Bound>
|
||||||
|
operator()(const Algebraic &a,const Polynomial_1 &p)const{
|
||||||
|
std::vector<Algebraic> roots;
|
||||||
|
std::back_insert_iterator<std::vector<Algebraic> > rit(roots);
|
||||||
|
typedef Solve_1<Polynomial_1,
|
||||||
|
ZPolynomial_1,
|
||||||
|
PolConverter,
|
||||||
|
Bound,
|
||||||
|
Algebraic,
|
||||||
|
Isolator,
|
||||||
|
Signat,
|
||||||
|
Ptraits,
|
||||||
|
ZPtraits> Solve;
|
||||||
|
typedef Compare_1<Algebraic,Bound,Comparator> Compare;
|
||||||
|
Solve()(p,false,rit);
|
||||||
|
for(typename std::vector<Algebraic>::size_type i=0;
|
||||||
|
i<roots.size();
|
||||||
|
++i){
|
||||||
|
// we use the comparison functor, that makes both
|
||||||
|
// intervals disjoint iff the algebraic numbers they
|
||||||
|
// represent are not equal
|
||||||
|
Compare()(a,roots[i]);
|
||||||
|
}
|
||||||
|
return std::make_pair(a.left(),a.right());
|
||||||
|
}
|
||||||
|
}; // Isolate_1
|
||||||
|
|
||||||
|
template <class Polynomial_,
|
||||||
|
class Bound_,
|
||||||
|
class Algebraic_,
|
||||||
|
class Refiner_>
|
||||||
|
struct Approximate_absolute_1:
|
||||||
|
public std::binary_function<Algebraic_,int,std::pair<Bound_,Bound_> >{
|
||||||
|
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<Bound,Bound> 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 <class Polynomial_,
|
||||||
|
class Bound_,
|
||||||
|
class Algebraic_,
|
||||||
|
class Refiner_>
|
||||||
|
struct Approximate_relative_1:
|
||||||
|
public std::binary_function<Algebraic_,int,std::pair<Bound_,Bound_> >{
|
||||||
|
typedef Polynomial_ Polynomial_1;
|
||||||
|
typedef Bound_ Bound;
|
||||||
|
typedef Algebraic_ Algebraic;
|
||||||
|
typedef Refiner_ Refiner;
|
||||||
|
|
||||||
|
std::pair<Bound,Bound> 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<unsigned>(
|
||||||
|
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
|
||||||
Loading…
Reference in New Issue