songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Code refactoring. 

The obsolete and unused code was removed. The interface was cleaned. The
memory leaks dissapeared (according to valgrind), because the pointers
to RS memory were removed.

For test purposes, the rational interface is not tested.
This commit is contained in:
Luis Peñaranda 2013-11-19 15:43:29 -02:00
parent da1791759c
commit cbdca2c35d
11 changed files with 1997 additions and 547 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

179
Algebraic_kernel_d/include/CGAL/RS/ak_1.h Normal file
View File

@ -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 <luis.penaranda@gmx.com>
#ifndef CGAL_RS_AK_1_H
#define CGAL_RS_AK_1_H
#include <cstddef> // included only to define size_t
#include <CGAL/Polynomial_traits_d.h>
#include "algebraic_1.h"
#include "comparator_1.h"
#include "signat_1.h"
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,
class Bound_,
class Isolator_,
class Refiner_,
class Ptraits_=CGAL::Polynomial_traits_d<Polynomial_> >
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<Polynomial_1,Bound>
Signat;
typedef CGAL::RS_AK1::Simple_comparator_1<Polynomial_1,
Bound,
Refiner,
Signat,
Ptraits>
Comparator;
public:
typedef CGAL::RS_AK1::Algebraic_1<Polynomial_1,
Bound,
Refiner,
Comparator,
Ptraits>
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,
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,
Bound,
Algebraic_real_1,
Isolator,
Comparator,
Signat,
Ptraits>
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,
Bound,
Algebraic_real_1,
Isolator,
Signat,
Ptraits>
Solve_1;
typedef CGAL::RS_AK1::Number_of_solutions_1<Polynomial_1,Isolator>
Number_of_solutions_1;
typedef CGAL::RS_AK1::Sign_at_1<Polynomial_1,
Bound,
Algebraic_real_1,
Refiner,
Signat,
Ptraits>
Sign_at_1;
typedef CGAL::RS_AK1::Is_zero_at_1<Polynomial_1,
Bound,
Algebraic_real_1,
Refiner,
Signat,
Ptraits>
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_1_H

418
Algebraic_kernel_d/include/CGAL/RS/algebraic_1.h
View File

@ -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 // 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 // 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, // published by the Free Software Foundation; version 2.1 of the License.
// or (at your option) any later version. // See the file LICENSE.LGPL distributed with CGAL.
// //
// Licensees holding a valid commercial license may use this file in // Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software. // accordance with the commercial license agreement provided with the software.
@ -19,123 +19,337 @@
#ifndef CGAL_RS_ALGEBRAIC_1_H #ifndef CGAL_RS_ALGEBRAIC_1_H
#define CGAL_RS_ALGEBRAIC_1_H #define CGAL_RS_ALGEBRAIC_1_H
#include <CGAL/Gmpz.h>
#include <CGAL/Gmpq.h>
#include <CGAL/Gmpfr.h>
#include <CGAL/Gmpfi.h>
#include <boost/operators.hpp> #include <boost/operators.hpp>
#include <CGAL/Real_embeddable_traits.h>
#include <CGAL/Gmpq.h>
#include <iostream>
namespace CGAL{ namespace CGAL{
namespace RS_AK1{
class RS_polynomial_1; // This class represents the simplest algebraic number one can think about.
class Algebraic_1; // 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
bool operator<(const Algebraic_1&,const Algebraic_1&); // root. Polynomial type and bound type are the first two template
bool operator==(const Algebraic_1&,const Algebraic_1&); // parameters.
//
// representation of algebraic numbers // The third template parameter is a refiner, a function object that
class Algebraic_1_rep{ // receives the polynomial and the bounds defining an algebraic number and
public: // an integer p, and modifies the two bounds until the difference between
mutable mpfi_t interval; // the two bounds is less than x*2^(-p), where x is the value of the
RS_polynomial_1 *polynomial; // represented algebraic number. The signature of a refiner must be:
int numroot; // void
int multiplicity; // Refiner_()(const Polynomial_&,Bound_&,Bound_&,int p);
mutable Sign lefteval; //
// The fourth template argument is a comparator, a function object that
Algebraic_1_rep(): // receives the polynomials and bounds defining two algebraic numbres and
polynomial(NULL),numroot(-1),multiplicity(-1),lefteval(ZERO){} // just compares them, returning a CGAL::Comparison_result. The signature
~Algebraic_1_rep(){} // 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).
template <class Polynomial_,
class Bound_,
class Refiner_,
class Comparator_,
class Ptraits_>
class Algebraic_1:
boost::totally_ordered<Algebraic_1<Polynomial_,
Bound_,
Refiner_,
Comparator_,
Ptraits_>,
double>{
private: private:
Algebraic_1_rep(const Algebraic_1_rep&); typedef Polynomial_ Polynomial;
Algebraic_1_rep& operator=(const Algebraic_1_rep&); 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<Polynomial,
Bound,
Refiner,
Comparator,
Ptraits>
Algebraic;
Polynomial pol;
mutable Bound left,right;
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)||(left<d&&right>d));
}
// 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 <class T>
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(){}
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);
}; };
// The class of the algebraic numbers. #define CGAL_RS_COMPARE_ALGEBRAIC_TYPE(_t) \
class Algebraic_1 : Handle_for<Algebraic_1_rep>, bool operator<(_t t)const \
boost::totally_ordered1<Algebraic_1>{ {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;}
typedef Handle_for<Algebraic_1_rep> Base; 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<Bound> 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<double,double> to_interval()const{
typedef Real_embeddable_traits<Bound> 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 Polynomial_,
class Bound_,
class Refiner_,
class Comparator_,
class Ptraits_>
class Real_embeddable_traits<RS_AK1::Algebraic_1<Polynomial_,
Bound_,
Refiner_,
Comparator_,
Ptraits_> >:
public INTERN_RET::Real_embeddable_traits_base<
RS_AK1::Algebraic_1<Polynomial_,
Bound_,
Refiner_,
Comparator_,
Ptraits_>,
CGAL::Tag_true>{
typedef Polynomial_ P;
typedef Bound_ B;
typedef Refiner_ R;
typedef Comparator_ C;
typedef Ptraits_ T;
public: public:
Algebraic_1(const Algebraic_1&,mpfr_prec_t,mpfr_prec_t); typedef RS_AK1::Algebraic_1<P,B,R,C,T> Type;
Algebraic_1(); typedef CGAL::Tag_true Is_real_embeddable;
Algebraic_1(int); typedef bool Boolean;
Algebraic_1(unsigned); typedef CGAL::Sign Sign;
Algebraic_1(long); typedef CGAL::Comparison_result Comparison_result;
Algebraic_1(unsigned long);
Algebraic_1(double);
Algebraic_1(mpz_srcptr);
Algebraic_1(mpq_srcptr);
Algebraic_1(mpfr_srcptr);
Algebraic_1(mpfi_srcptr);
Algebraic_1(const Gmpz&);
Algebraic_1(const Gmpq&);
Algebraic_1(const Gmpfr&);
// the only interesting constructor typedef INTERN_RET::Real_embeddable_traits_base<Type,CGAL::Tag_true>
Algebraic_1(mpfi_srcptr, Base;
const RS_polynomial_1&, typedef typename Base::Compare Compare;
int,
int
//,mpfi_ptr,mpfi_ptr
);
// the another interesting variant class Sgn:public std::unary_function<Type,CGAL::Sign>{
Algebraic_1(mpfi_srcptr, public:
const RS_polynomial_1&, CGAL::Sign operator()(const Type &a)const{
int, return Compare()(a,Type(0));
int, }
//mpfi_ptr,mpfi_ptr, };
CGAL::Sign);
// functions related to the member data class To_double:public std::unary_function<Type,double>{
mpfi_srcptr mpfi()const; public:
mpfi_ptr mpfi(); double operator()(Type a)const{return a.to_double();}
Gmpfi interval()const; };
Gmpfr inf()const;
Gmpfr sup()const;
const RS_polynomial_1& pol()const;
int nr()const;
int mult()const;
void set_mpfi(mpfi_srcptr);
void set_mpfi_ptr(mpfi_srcptr);
void clear_pol();
void set_pol(const RS_polynomial_1 &);
void set_nr(int);
void set_mult(int);
void set_prec(mp_prec_t);
void set_lefteval(Sign)const;
mp_prec_t get_prec()const;
Sign lefteval()const;
mpfr_srcptr left()const;
mpfr_srcptr right()const;
bool is_consistent()const;
bool is_point()const; // are endpoints equal?
bool overlaps(const Algebraic_1&)const;
bool contains(int n)const;
bool contains(mpfr_srcptr)const;
bool contains(const Gmpz&)const;
Algebraic_1 operator+()const; class To_interval:
Algebraic_1 operator-()const; public std::unary_function<Type,std::pair<double,double> >{
public:
std::pair<double,double> operator()(const Type &a)const{
return a.to_interval();
}
};
bool is_valid()const; class Is_zero:public std::unary_function<Type,Boolean>{
bool is_finite()const; public:
//template<class> bool operator()(const Type &a)const{
double to_double()const; return Sgn()(a)==CGAL::ZERO;
std::pair<double,double> to_interval() const; }
Algebraic_1 sqrt()const; };
}; // class Algebraic_1
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 B,class R,class C,class T>
inline
RS_AK1::Algebraic_1<P,B,R,C,T> min
BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_1<P,B,R,C,T> a,
RS_AK1::Algebraic_1<P,B,R,C,T> b){
return(a<b?a:b);
}
template <class P,class B,class R,class C,class T>
inline
RS_AK1::Algebraic_1<P,B,R,C,T> max
BOOST_PREVENT_MACRO_SUBSTITUTION(RS_AK1::Algebraic_1<P,B,R,C,T> a,
RS_AK1::Algebraic_1<P,B,R,C,T> b){
return(a>b?a:b);
}
template <class P,class B,class R,class C,class T>
inline
std::ostream& operator<<(std::ostream &o,
const RS_AK1::Algebraic_1<P,B,R,C,T> &a){
return(o<<'['<<a.get_pol()<<','<<
a.get_left()<<','<<
a.get_right()<<']');
}
template <class P,class B,class R,class C,class T>
inline
std::istream& operator>>(std::istream &i,
RS_AK1::Algebraic_1<P,B,R,C,T> &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<P,B,R,C,T>(pol,lb,rb);
return i;
}
} // namespace CGAL } // namespace CGAL
#include <CGAL/RS/algebraic_1_constructors.h>
#include <CGAL/RS/algebraic_1_operators.h>
#include <CGAL/RS/algebraic_1_member.h> // other member functions
#include <CGAL/RS/algebraic_1_comparisons.h>
#include <CGAL/RS/algebraic_1_other.h> // related non-member functions
#include <CGAL/RS/algebraic_1_real_embeddable.h>
#endif // CGAL_RS_ALGEBRAIC_1_H #endif // CGAL_RS_ALGEBRAIC_1_H

127
Algebraic_kernel_d/include/CGAL/RS/bisection_refiner_1.h Normal file
View File

@ -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 <luis.penaranda@gmx.com>
// 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 <CGAL/Polynomial_traits_d.h>
#include "signat_1.h"
namespace CGAL{
template <class Polynomial_,class Bound_>
struct Bisection_refiner_1{
typedef CGAL::RS_AK1::Signat_1<Polynomial_,Bound_> Signat;
void operator()(const Polynomial_&,Bound_&,Bound_&,int);
}; // class Bisection_refiner_1
// TODO: Write in a generic way, if possible.
template <class Polynomial_,class Bound_>
void
Bisection_refiner_1<Polynomial_,Bound_>::
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<Polynomial<Gmpz>,Gmpfr>::
operator()(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec){
typedef Polynomial<Gmpz> Polynomial;
typedef Polynomial_traits_d<Polynomial> Ptraits;
typedef Ptraits::Degree Degree;
typedef Ptraits::Make_square_free Sfpart;
typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
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 ["<<left<<","<<right<<"]"<<std::endl;
#ifndef CGAL_GMPFR_NO_REFCOUNT
// Make sure the endpoints do not share references. If some of them
// does, copy it.
if(!left.is_unique()){
Gmpfr new_left(0,left.get_precision());
mpfr_set(new_left.fr(),left.fr(),GMP_RNDN);
left=new_left;
CGAL_assertion_code(new_left=Gmpfr();)
CGAL_assertion(left.is_unique());
}
if(!right.is_unique()){
Gmpfr new_right(0,right.get_precision());
mpfr_set(new_right.fr(),right.fr(),GMP_RNDN);
right=new_right;
CGAL_assertion_code(new_right=Gmpfr();)
CGAL_assertion(right.is_unique());
}
#endif // CGAL_GMPFR_NO_REFCOUNT
Polynomial sfpp=Sfpart()(pol);
Signat signof(sfpp);
CGAL::Sign sl,sc;
mp_prec_t pl,pc;
mpfr_t center;
sl=signof(left);
if(sl==ZERO)
return;
pl=left.get_precision();
pc=right.get_precision();
pc=(pl>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<prec;++i){
CGAL_assertion_code(round=)
mpfr_add(center,left.fr(),right.fr(),GMP_RNDN);
CGAL_assertion(!round);
CGAL_assertion_code(round=)
mpfr_div_2ui(center,center,1,GMP_RNDN);
CGAL_assertion(!round);
sc=signof(Gmpfr(center));
if(sc==ZERO){ // we have a root
CGAL_assertion_code(round=)
mpfr_set(left.fr(),center,GMP_RNDN);
CGAL_assertion(!round);
mpfr_swap(right.fr(),center);
break;
}
if(sc==sl)
mpfr_swap(left.fr(),center);
else
mpfr_swap(right.fr(),center);
}
mpfr_clear(center);
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
} // namespace CGAL
#endif // CGAL_RS_BISECTION_REFINER_1_H

89
Algebraic_kernel_d/include/CGAL/RS/comparator_1.h Normal file
View File

@ -0,0 +1,89 @@
// 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_COMPARATOR_1_H
#define CGAL_RS_COMPARATOR_1_H
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,
class Bound_,
class Refiner_,
class Signat_,
class Ptraits_>
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(r1<l2)
return SMALLER;
}else{
if(r2<l1)
return LARGER;
}
Polynomial G=Gcd()(p1,p2);
if(Degree()(G)==0)
return compare_unequal(p1,l1,r1,p2,l2,r2);
Signat sg(G);
CGAL::Sign sleft=sg(l1>l2?l1:l2);
if(sleft==ZERO)
return EQUAL;
CGAL::Sign sright=sg(r1<r2?r1:r2);
if(sleft!=sright)
return EQUAL;
else
return compare_unequal(p1,l1,r1,p2,l2,r2);
}
// This function compares two algebraic numbers, assuming that they
// are not equal.
CGAL::Comparison_result
compare_unequal(const Polynomial &p1,Bound &l1,Bound &r1,
const Polynomial &p2,Bound &l2,Bound &r2)const{
CGAL_precondition(l1<=r1&&l2<=r2);
int prec=CGAL::max(
CGAL::max(l1.get_precision(),
r1.get_precision()),
CGAL::max(l2.get_precision(),
r2.get_precision()));
do{
prec*=2;
Refiner()(p1,l1,r1,prec);
Refiner()(p2,l2,r2,prec);
CGAL_assertion(l1<=r1&&l2<=r2);
}while(l1<=l2?r1>=l2:r2>=l1);
return (r1<l2?SMALLER:LARGER);
}
}; // struct Simple_comparator_1
} // namespace RS_AK1
} // namespace CGAL
#endif // CGAL_RS_COMPARATOR_1_H

957
Algebraic_kernel_d/include/CGAL/RS/functors_1.h
View File

File diff suppressed because it is too large Load Diff

123
Algebraic_kernel_d/include/CGAL/RS/rs23_k_isolator_1.h Normal file
View File

@ -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 <luis.penaranda@gmx.com>
#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 <CGAL/Gmpfi.h>
#include <vector>
namespace CGAL{
template <class Polynomial_,class Bound_>
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<Interval> _real_roots;
};
template <class Polynomial_,class Bound_>
RS23_k_isolator_1<Polynomial_,Bound_>::
RS23_k_isolator_1(const Polynomial_ &p){
CGAL_error_msg("not implemented for these polynomial/bound types");
}
template <>
RS23_k_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,Gmpfr>::
RS23_k_isolator_1(const CGAL::Polynomial<CGAL::Gmpz> &p):_polynomial(p){
typedef CGAL::Polynomial<CGAL::Gmpz> Pol;
typedef CGAL::Gmpfr Bound;
typedef typename CGAL::RS3::RS3_k_refiner_1<Pol,Bound> KRefiner;
int numsols;
unsigned int degree=p.degree();
mpz_t *coeffs=(mpz_t*)malloc((degree+1)*sizeof(mpz_t));
std::vector<Gmpfi> 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<numsols;++j){
Gmpfr left(intervals[j].inf());
Gmpfr right(intervals[j].sup());
KRefiner()(p,left,right,53);
_real_roots.push_back(intervals[j]);
}
free(coeffs);
}
template <class Polynomial_,class Bound_>
Polynomial_
RS23_k_isolator_1<Polynomial_,Bound_>::polynomial()const{
return _polynomial;
}
template <class Polynomial_,class Bound_>
int
RS23_k_isolator_1<Polynomial_,Bound_>::number_of_real_roots()const{
return _real_roots.size();
}
template <class Polynomial_,class Bound_>
bool
RS23_k_isolator_1<Polynomial_,Bound_>::is_exact_root(int i)const{
return _real_roots[i].inf()==_real_roots[i].sup();
}
template <class Polynomial_,class Bound_>
Bound_
RS23_k_isolator_1<Polynomial_,Bound_>::left_bound(int i)const{
return _real_roots[i].inf();
}
template <class Polynomial_,class Bound_>
Bound_
RS23_k_isolator_1<Polynomial_,Bound_>::right_bound(int i)const{
return _real_roots[i].sup();
}
} // namespace CGAL
#endif // CGAL_RS_RS23_K_ISOLATOR_1_H

134
Algebraic_kernel_d/include/CGAL/RS/rs2_calls.h Normal file
View File

@ -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 <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS2_CALLS_H
#define CGAL_RS_RS2_CALLS_H
#include <gmp.h>
#include <CGAL/Gmpfr.h>
#include <CGAL/Gmpfi.h>
#include <rs_exports.h>
#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<nb_elts;++i){
ident_vect=rs_export_list_vect_ibfr_monnode
(ident_node);
CGAL_assertion_msg(rs_export_dim_vect_ibfr
(ident_vect)==1,
"vector dimension must be 1");
ident_elt=rs_export_elt_vect_ibfr(ident_vect,0);
mpfi_ptr root_pointer=
(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt);
mpfi_init2(roots[i],mpfi_get_prec(root_pointer));
mpfi_set(roots[i],root_pointer);
x[i]=roots[i];
// This doesn't work because RS relocates the
// mpfrs that form the mpfi. Nevertheless, the
// mpfi address is not changed.
//x[i]=(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt);
ident_node=rs_export_list_vect_ibfr_nextnode
(ident_node);
}
return nb_elts;
}
template<class OutputIterator>
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;i<nb_elts;++i){
ident_vect=rs_export_list_vect_ibfr_monnode
(ident_node);
CGAL_assertion_msg(rs_export_dim_vect_ibfr
(ident_vect)==1,
"vector dimension must be 1");
ident_elt=rs_export_elt_vect_ibfr(ident_vect,0);
mpfi_ptr root_pointer=
(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt);
mp_prec_t root_prec=mpfi_get_prec(root_pointer);
// Construct Gmpfr's with pointers to endpoints.
Gmpfr left(&(root_pointer->left),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

105
Algebraic_kernel_d/include/CGAL/RS/rs2_isolator_1.h Normal file
View File

@ -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 <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS2_ISOLATOR_1_H
#define CGAL_RS_RS2_ISOLATOR_1_H
#include "rs2_calls.h"
#include <CGAL/Gmpfi.h>
#include <vector>
namespace CGAL{
namespace RS2{
template <class Polynomial_,class Bound_>
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<Interval> _real_roots;
};
template <class Polynomial_,class Bound_>
RS2_isolator_1<Polynomial_,Bound_>::
RS2_isolator_1(const Polynomial_ &p){
CGAL_error_msg("not implemented for these polynomial/bound types");
}
template <>
RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,Gmpfr>::
RS2_isolator_1(const CGAL::Polynomial<CGAL::Gmpz> &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 <class Polynomial_,class Bound_>
Polynomial_
RS2_isolator_1<Polynomial_,Bound_>::polynomial()const{
return _polynomial;
}
template <class Polynomial_,class Bound_>
int
RS2_isolator_1<Polynomial_,Bound_>::number_of_real_roots()const{
return _real_roots.size();
}
template <class Polynomial_,class Bound_>
bool
RS2_isolator_1<Polynomial_,Bound_>::is_exact_root(int i)const{
return _real_roots[i].inf()==_real_roots[i].sup();
}
template <class Polynomial_,class Bound_>
Bound_
RS2_isolator_1<Polynomial_,Bound_>::left_bound(int i)const{
return _real_roots[i].inf();
}
template <class Polynomial_,class Bound_>
Bound_
RS2_isolator_1<Polynomial_,Bound_>::right_bound(int i)const{
return _real_roots[i].sup();
}
} // namespace RS2
} // namespace CGAL
#endif // CGAL_RS_RS2_ISOLATOR_1_H

108
Algebraic_kernel_d/include/CGAL/RS/rs3_k_refiner_1.h Normal file
View File

@ -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 <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS3_K_REFINER_1_H
#define CGAL_RS_RS3_K_REFINER_1_H
#include <CGAL/Polynomial_traits_d.h>
#include "rs2_calls.h"
#include <rs3_fncts.h>
namespace CGAL{
namespace RS3{
template <class Polynomial_,class Bound_>
struct RS3_k_refiner_1{
void operator()(const Polynomial_&,Bound_&,Bound_&,int);
}; // class RS3_k_refiner_1
template <class Polynomial_,class Bound_>
void
RS3_k_refiner_1<Polynomial_,Bound_>::
operator()(const Polynomial_&,Bound_&,Bound_&,int){
CGAL_error_msg("RS3 k-refiner not implemented for these types");
return;
}
template<>
void
RS3_k_refiner_1<Polynomial<Gmpz>,Gmpfr>::
operator()
(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec){
typedef Polynomial<Gmpz> Polynomial;
typedef Polynomial_traits_d<Polynomial> 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 ["<<left<<","<<right<<"]"<<std::endl;
int deg=Degree()(pol);
mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
__mpfi_struct interval;
#ifndef CGAL_GMPFR_NO_REFCOUNT
// Make sure the endpoints do not share references. If some of them
// does, copy it.
if(!left.is_unique()){
Gmpfr new_left(0,left.get_precision());
mpfr_set(new_left.fr(),left.fr(),GMP_RNDN);
left=new_left;
CGAL_assertion_code(new_left=Gmpfr();)
CGAL_assertion(left.is_unique());
}
if(!right.is_unique()){
Gmpfr new_right(0,right.get_precision());
mpfr_set(new_right.fr(),right.fr(),GMP_RNDN);
right=new_right;
CGAL_assertion_code(new_right=Gmpfr();)
CGAL_assertion(right.is_unique());
}
#endif // CGAL_GMPFR_NO_REFCOUNT
interval.left=*(left.fr());
interval.right=*(right.fr());
for(int i=0;i<=deg;++i)
coefficients[i][0]=*(pol[i].mpz());
RS2::RS2_calls::init_solver();
rs3_refine_u_root(&interval,
coefficients,
deg,
prec+CGAL::max(left.get_precision(),
right.get_precision()),
1,
1);
free(coefficients);
mpfr_clear(left.fr());
mpfr_clear(right.fr());
mpfr_custom_init_set(left.fr(),
mpfr_custom_get_kind(&interval.left),
mpfr_custom_get_exp(&interval.left),
mpfr_get_prec(&interval.left),
mpfr_custom_get_mantissa(&interval.left));
mpfr_custom_init_set(right.fr(),
mpfr_custom_get_kind(&interval.right),
mpfr_custom_get_exp(&interval.right),
mpfr_get_prec(&interval.right),
mpfr_custom_get_mantissa(&interval.right));
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
} // namespace RS3
} // namespace CGAL
#endif // CGAL_RS_RS3_K_REFINER_1_H

108
Algebraic_kernel_d/include/CGAL/RS/rs3_refiner_1.h Normal file
View File

@ -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 <luis.penaranda@gmx.com>
#ifndef CGAL_RS_RS3_REFINER_1_H
#define CGAL_RS_RS3_REFINER_1_H
#include <CGAL/Polynomial_traits_d.h>
#include "rs2_calls.h"
#include <rs3_fncts.h>
namespace CGAL{
namespace RS3{
template <class Polynomial_,class Bound_>
struct RS3_refiner_1{
void operator()(const Polynomial_&,Bound_&,Bound_&,int,int=0);
}; // class RS3_refiner_1
template <class Polynomial_,class Bound_>
void
RS3_refiner_1<Polynomial_,Bound_>::
operator()(const Polynomial_&,Bound_&,Bound_&,int,int){
CGAL_error_msg("RS3 refiner not implemented for these types");
return;
}
template<>
void
RS3_refiner_1<Polynomial<Gmpz>,Gmpfr>::
operator()
(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec,int k){
typedef Polynomial<Gmpz> Polynomial;
typedef Polynomial_traits_d<Polynomial> 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 ["<<left<<","<<right<<"]"<<std::endl;
int deg=Degree()(pol);
mpz_t* coefficients=(mpz_t*)malloc((deg+1)*sizeof(mpz_t));
__mpfi_struct interval;
#ifndef CGAL_GMPFR_NO_REFCOUNT
// Make sure the endpoints do not share references. If some of them
// does, copy it.
if(!left.is_unique()){
Gmpfr new_left(0,left.get_precision());
mpfr_set(new_left.fr(),left.fr(),GMP_RNDN);
left=new_left;
CGAL_assertion_code(new_left=Gmpfr();)
CGAL_assertion(left.is_unique());
}
if(!right.is_unique()){
Gmpfr new_right(0,right.get_precision());
mpfr_set(new_right.fr(),right.fr(),GMP_RNDN);
right=new_right;
CGAL_assertion_code(new_right=Gmpfr();)
CGAL_assertion(right.is_unique());
}
#endif // CGAL_GMPFR_NO_REFCOUNT
interval.left=*(left.fr());
interval.right=*(right.fr());
for(int i=0;i<=deg;++i)
coefficients[i][0]=*(pol[i].mpz());
RS2::RS2_calls::init_solver();
rs3_refine_u_root(&interval,
coefficients,
deg,
prec+CGAL::max(left.get_precision(),
right.get_precision()),
k,
k);
free(coefficients);
mpfr_clear(left.fr());
mpfr_clear(right.fr());
mpfr_custom_init_set(left.fr(),
mpfr_custom_get_kind(&interval.left),
mpfr_custom_get_exp(&interval.left),
mpfr_get_prec(&interval.left),
mpfr_custom_get_mantissa(&interval.left));
mpfr_custom_init_set(right.fr(),
mpfr_custom_get_kind(&interval.right),
mpfr_custom_get_exp(&interval.right),
mpfr_get_prec(&interval.right),
mpfr_custom_get_mantissa(&interval.right));
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
} // namespace RS3
} // namespace CGAL
#endif // CGAL_RS_RS3_REFINER_1_H

84
Algebraic_kernel_d/include/CGAL/RS/signat_1.h Normal file
View File

@ -0,0 +1,84 @@
// 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_SIGNAT_1_H
#define CGAL_RS_SIGNAT_1_H
#include <CGAL/Gmpfi.h>
#include <CGAL/Polynomial_traits_d.h>
//#include <boost/mpl/assert.hpp>
namespace CGAL{
namespace RS_AK1{
template <class Polynomial_,class Bound_>
struct Signat_1{
typedef Polynomial_ Polynomial;
typedef Bound_ Bound;
typedef CGAL::Polynomial_traits_d<Polynomial> 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 <class Polynomial_,class Bound_>
inline CGAL::Sign
Signat_1<Polynomial_,Bound_>::operator()(const Bound_ &x)const{
typedef Polynomial_ Polynomial;
typedef Bound_ Bound;
typedef Real_embeddable_traits<Bound> REtraits;
typedef typename REtraits::Sgn BSign;
typedef Algebraic_structure_traits<Bound> 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<AStraits::Is_exact,Tag_true>));
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<Polynomial<Gmpz>,Gmpfr>::operator()(const Gmpfr &x)const{
typedef Signat_1<Polynomial,Gmpq> 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<CGAL::Sign> indet=Uncertain<CGAL::Sign>::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