cleanup & passed kernel as argument

This commit is contained in:
Efi Fogel 2011-06-20 09:58:35 +00:00
parent 3b19f90449
commit e9eda73bcc
3 changed files with 482 additions and 342 deletions

View File

@ -21,61 +21,72 @@ namespace Arr_rational_arc {
//-------------------
//Algebraic_point_2_rep
//-------------------
template < class Algebraic_kernel_ >
class Algebraic_point_2_rep: public Base_rational_arc_ds_1<Algebraic_kernel_>
template <typename Algebraic_kernel_>
class Algebraic_point_2_rep : public Base_rational_arc_ds_1<Algebraic_kernel_>
{
public:
typedef Algebraic_kernel_ Algebraic_kernel;
typedef Base_rational_arc_ds_1<Algebraic_kernel> Base;
typedef Algebraic_kernel_ Algebraic_kernel;
typedef Base_rational_arc_ds_1<Algebraic_kernel> Base;
typedef CGAL::Arr_rational_arc::Rational_function<Algebraic_kernel> Rational_function;
typedef CGAL::Arr_rational_arc::Rational_function_pair<Algebraic_kernel> Rational_function_pair;
typedef CGAL::Arr_rational_arc::Rational_function<Algebraic_kernel>
Rational_function;
typedef CGAL::Arr_rational_arc::Rational_function_pair<Algebraic_kernel>
Rational_function_pair;
typedef typename Base::Algebraic_real_1 Algebraic_real_1;
typedef typename Algebraic_kernel::Bound Bound;
typedef typename Base::Integer Integer ;
typedef typename Base::Rational Rational ;
typedef typename Base::Polynomial_1 Polynomial_1;
typedef typename Base::Algebraic_real_1 Algebraic_real_1;
typedef typename Algebraic_kernel::Bound Bound;
typedef typename Base::Integer Integer ;
typedef typename Base::Rational Rational ;
typedef typename Base::Polynomial_1 Polynomial_1;
typedef typename Base::Root_multiplicity_vector Root_multiplicity_vector;
typedef typename Base::Root_multiplicity_vector Root_multiplicity_vector;
typedef typename Get_arithmetic_kernel<Algebraic_real_1>::Arithmetic_kernel AK;
typedef typename AK::Bigfloat_interval BFI;
typedef Bigfloat_interval_traits<BFI> BFI_traits;
typedef CGAL::Polynomial<BFI> BFI_polynomial;
typedef CGAL::Polynomial_traits_d<BFI_polynomial> BFI_polynomial_traits;
typedef typename Get_arithmetic_kernel<Algebraic_real_1>::Arithmetic_kernel
AK;
typedef typename AK::Bigfloat_interval BFI;
typedef Bigfloat_interval_traits<BFI> BFI_traits;
typedef CGAL::Polynomial<BFI> BFI_polynomial;
typedef CGAL::Polynomial_traits_d<BFI_polynomial> BFI_polynomial_traits;
typedef typename Base::FT_rat_1 FT_rat_1;
typedef typename Base::Polynomial_traits_1 Polynomial_traits_1;
typedef typename Base::FT_rat_1 FT_rat_1;
typedef typename Base::Polynomial_traits_1 Polynomial_traits_1;
typedef CGAL::Arr_rational_arc::Cache<Algebraic_kernel_> Cache;
typedef CGAL::Arr_rational_arc::Cache<Algebraic_kernel_>
Cache;
public:
Algebraic_point_2_rep(){}
Algebraic_point_2_rep(const Rational_function& rational_function,const Algebraic_real_1& x_coordinate):
_rational_function(rational_function),_x_coordinate(x_coordinate) {}
Algebraic_point_2_rep(const Rational_function& rational_function,
const Algebraic_real_1& x_coordinate) :
_rational_function(rational_function),
_x_coordinate(x_coordinate) {}
//assignment oparator
Algebraic_point_2_rep & operator = (const Algebraic_point_2_rep & other)
Algebraic_point_2_rep& operator=(const Algebraic_point_2_rep& other)
{
if (this != &other) // protect against invalid self-assignment
{
_rational_function = other._rational_function;
_x_coordinate = other._x_coordinate;
}
{
_rational_function = other._rational_function;
_x_coordinate = other._x_coordinate;
}
return *this;
}
Comparison_result compare_xy_2 (const Algebraic_point_2_rep & other,Cache& cache) const
Comparison_result compare_xy_2(const Algebraic_point_2_rep& other,
Cache& cache, Algebraic_kernel& kernel) const
{
Comparison_result comp = CGAL::compare (_x_coordinate, other.x());
Comparison_result comp = CGAL::compare(_x_coordinate, other.x());
if (comp != EQUAL)
return comp;
if (this->_rational_function == other.rational_function())
if (_rational_function == other.rational_function())
return EQUAL;
Rational_function_pair rat_func_pair = cache.get_rational_pair( this->_rational_function,
other.rational_function());
Rational_function_pair rat_func_pair =
cache.get_rational_pair(_rational_function, other.rational_function(),
kernel);
return rat_func_pair.compare_f_g_at(_x_coordinate);
}
Algebraic_real_1& x()
{
return _x_coordinate;
@ -85,6 +96,7 @@ public:
{
return _x_coordinate;
}
const Rational_function& rational_function() const
{
return _rational_function;
@ -95,33 +107,36 @@ public:
Algebraic_real_1 y() const
{
typedef CGAL::Polynomial<Polynomial_1> Polynomial_2;
//converting the defining polynomial of x and the rational function to bivariate polynomials
//converting the defining polynomial of x and the rational function to
//bivariate polynomials
Polynomial_2 f(_algebraic_kernel.compute_polynomial_1_object()(_x_coordinate));
Polynomial_2 y(CGAL::shift(Polynomial_2(1),1));
Polynomial_2 g(_rational_function.numer() - y * _rational_function.denom());
f=CGAL::swap(f,0,1); //swap x and y in the polynomial f
g=CGAL::swap(g,0,1); //swap x and y in the polynomial g
Polynomial_1 r(CGAL::resultant(f,g)); //compute the resultant in x (polynomial in y)
f=CGAL::swap(f, 0, 1); //swap x and y in the polynomial f
g=CGAL::swap(g, 0, 1); //swap x and y in the polynomial g
//compute the resultant in x (polynomial in y)
Polynomial_1 r(CGAL::resultant(f,g));
//solve for all roots of resultant
std::list<Algebraic_real_1> roots;
_algebraic_kernel.solve_1_object()(r,false,std::back_inserter(roots));
_algebraic_kernel.solve_1_object()(r, false, std::back_inserter(roots));
//isolate the right root
unsigned int initial_precision = 16;
int error_bound = 2;
while (roots.size() > 1)
{
std::pair<Bound,Bound> y_bounds (this->approximate_absolute_y(error_bound,initial_precision));
std::pair<Bound, Bound>
y_bounds(this->approximate_absolute_y(error_bound,initial_precision));
while (CGAL::compare(roots.front(),y_bounds.first) == SMALLER)
roots.pop_front();
while (CGAL::compare(y_bounds.second,roots.back()) == SMALLER)
roots.pop_back();
error_bound*=2;
error_bound *= 2;
}
CGAL_postcondition (roots.size() == 1);
@ -158,8 +173,10 @@ public:
set_precision(precision);
BFI x_bfi(convert_to_bfi(_x_coordinate));
BFI_polynomial numer_bfi(convert_to_bfi_extended(_rational_function.numer()));
BFI_polynomial denom_bfi(convert_to_bfi_extended(_rational_function.denom()));
BFI_polynomial
numer_bfi(convert_to_bfi_extended(_rational_function.numer()));
BFI_polynomial
denom_bfi(convert_to_bfi_extended(_rational_function.denom()));
BFI y_numer_bfi(evaluate(numer_bfi,x_bfi));
BFI y_denom_bfi(evaluate(denom_bfi,x_bfi));
@ -167,9 +184,11 @@ public:
if (CGAL::zero_in(y_denom_bfi) == false)
{
BFI y_bfi(y_numer_bfi/y_denom_bfi);
if (CGAL::compare (CGAL::width(y_bfi),Rational(CGAL::lower(CGAL::abs(y_bfi)))*error_bound ) == SMALLER)
if (CGAL::compare(CGAL::width(y_bfi),
Rational(CGAL::lower(CGAL::abs(y_bfi))) * error_bound )
== SMALLER)
return std::make_pair(Bound(CGAL::lower(y_bfi)),
Bound(CGAL::upper(y_bfi)) );
Bound(CGAL::upper(y_bfi)));
}
else precision*=2;
}
@ -232,16 +251,21 @@ private:
else precision*=2;
}
}
template <class NTX>
static typename CGAL::Coercion_traits<NTX,typename Get_arithmetic_kernel<NTX>::Arithmetic_kernel::Bigfloat_interval>::Type
template <typename NTX>
static typename
CGAL::Coercion_traits<NTX,
typename Get_arithmetic_kernel<NTX>::
Arithmetic_kernel::Bigfloat_interval>::Type
convert_to_bfi_extended(const NTX& x)
{
typedef typename Get_arithmetic_kernel<NTX>::Arithmetic_kernel AK;
typedef typename AK::Bigfloat_interval BFI;
typedef CGAL::Coercion_traits<NTX,BFI> CT;
typedef CGAL::Coercion_traits<NTX, BFI> CT;
return typename CT::Cast()(x);
}
double evaluate_at(const Polynomial_1& poly,const double x) const
double evaluate_at(const Polynomial_1& poly, const double x) const
{
double x_val = 1;
double ret_val(0);
@ -252,6 +276,7 @@ private:
}
return ret_val;
}
private:
Rational_function _rational_function; //supporting rational function
Algebraic_real_1 _x_coordinate;
@ -261,13 +286,15 @@ private:
template <class Algebraic_kernel_ >
class Algebraic_point_2: public Handle_with_policy<Algebraic_point_2_rep<Algebraic_kernel_> >
template <typename Algebraic_kernel_ >
class Algebraic_point_2 :
public Handle_with_policy<Algebraic_point_2_rep<Algebraic_kernel_> >
{
public:
typedef Algebraic_kernel_ Algebraic_kernel;
typedef Handle_with_policy<Algebraic_point_2_rep<Algebraic_kernel_> > Base;
typedef Algebraic_kernel_ Algebraic_kernel;
typedef Handle_with_policy<Algebraic_point_2_rep<Algebraic_kernel_> >
Base;
typedef Algebraic_point_2<Algebraic_kernel> Self;
typedef Algebraic_point_2_rep<Algebraic_kernel_> Rep;
typedef typename Rep::Rational Rational;
@ -285,7 +312,8 @@ private:
static Rational_function rational_function(numer,denom);
static Algebraic_kernel kernel;
static Algebraic_real_1 x_coordinate = kernel.construct_algebraic_real_1_object()(Rational(0));
static Algebraic_real_1 x_coordinate =
kernel.construct_algebraic_real_1_object()(Rational(0));
static Self default_instance(rational_function,x_coordinate);
@ -294,20 +322,25 @@ private:
/*static Self x = Self(Rational(0),Rational(0),_cache);
return x; */
}
public:
explicit Algebraic_point_2( const Rational_function& rational_function,
const Algebraic_real_1& x_coordinate)
: Base(rational_function,x_coordinate) {}
explicit Algebraic_point_2(const Rational_function& rational_function,
const Algebraic_real_1& x_coordinate) :
Base(rational_function,x_coordinate) {}
//used to solve VS bug...
Algebraic_point_2 (const Self & p = get_default_instance()) : Base(static_cast<const Base &> (p)) {}
Algebraic_point_2(const Self & p = get_default_instance()) :
Base(static_cast<const Base &> (p)) {}
Comparison_result compare_xy_2 (const Self & other,Cache& cache) const
Comparison_result compare_xy_2(const Algebraic_point_2& other,
Cache& cache, Algebraic_kernel& kernel) const
{
if (this->is_identical (other))
return CGAL::EQUAL;
return this->ptr()->compare_xy_2(*other.ptr(),cache);
return this->ptr()->compare_xy_2(*other.ptr(), cache, kernel);
}
Algebraic_real_1& x()
{
if (this->is_shared())
@ -324,45 +357,52 @@ public:
{
return this->ptr()->rational_function();
}
Algebraic_real_1 y() const
{
return this->ptr()->y();
}
std::pair<double,double> to_double() const
{
return this->ptr()->to_double();
}
std::pair<Bound,Bound> approximate_absolute_x( int a) const
{
return this->ptr()->approximate_absolute_x(a);
}
std::pair<Bound,Bound> approximate_absolute_y( int a) const
{
return this->ptr()->approximate_absolute_y(a);
}
std::pair<Bound,Bound> approximate_relative_x( int r) const
{
return this->ptr()->approximate_relative_x(r);
}
std::pair<Bound,Bound> approximate_relative_y( int r) const
{
return this->ptr()->approximate_relative_y(r);
}
std::ostream& print (std::ostream& os) const
std::ostream& print(std::ostream& os) const
{
return this->ptr()->print(os);
}
}; //Algebraic_point_2
template < class Algebraic_kernel_ >
std::ostream&
operator<< (std::ostream& os,
const Algebraic_point_2<Algebraic_kernel_> & p)
template < typename Algebraic_kernel_>
std::ostream& operator<<(std::ostream& os,
const Algebraic_point_2<Algebraic_kernel_> & p)
{
return (p.print (os));
return (p.print(os));
}
} //namespace Arr_rational_arc
} //namespace CGAL {
#endif //CGAL_ALBERAIC_POINT_D_1_H

View File

@ -8,10 +8,11 @@
namespace CGAL {
namespace Arr_rational_arc {
//-------------------
//Cache
//-------------------
template < class Algebraic_kernel_ >
template <typename Algebraic_kernel_>
class Cache : public Base_rational_arc_ds_1<Algebraic_kernel_>
{
public:
@ -26,24 +27,31 @@ public:
typedef typename Base::FT_rat_1 FT_rat_1;
typedef typename Base::Polynomial_traits_1 Polynomial_traits_1;
typedef CGAL::Arr_rational_arc::Rational_function<Algebraic_kernel> Rational_function;
typedef CGAL::Arr_rational_arc::Rational_function_canonicalized_pair<Algebraic_kernel> Rational_function_canonicalized_pair;
typedef CGAL::Arr_rational_arc::Rational_function_pair<Algebraic_kernel> Rational_function_pair;
typedef CGAL::Arr_rational_arc::Rational_function<Algebraic_kernel>
Rational_function;
typedef CGAL::Arr_rational_arc::Rational_function_canonicalized_pair<Algebraic_kernel>
Rational_function_canonicalized_pair;
typedef CGAL::Arr_rational_arc::Rational_function_pair<Algebraic_kernel>
Rational_function_pair;
typedef std::pair<Polynomial_1,Polynomial_1> Rational_function_key;
class Less_compare_rational_function_key
{
public:
Comparison_result compare(const Rational_function_key& k1, const Rational_function_key& k2) const
Comparison_result compare(const Rational_function_key& k1,
const Rational_function_key& k2) const
{
Comparison_result cr = typename Polynomial_traits_1::Compare()(k1.first,k2.first);
Comparison_result cr =
typename Polynomial_traits_1::Compare()(k1.first, k2.first);
if (cr != CGAL::EQUAL)
return (cr);
cr = typename Polynomial_traits_1::Compare()(k1.second,k2.second);
cr = typename Polynomial_traits_1::Compare()(k1.second, k2.second);
return cr;
}
bool operator()(const Rational_function_key& k1, const Rational_function_key& k2) const
bool operator()(const Rational_function_key& k1,
const Rational_function_key& k2) const
{
Comparison_result cr = compare(k1,k2);
return ((cr == CGAL::LARGER) ? true : false);
@ -56,7 +64,7 @@ public:
Rational_function_map;
typedef typename Rational_function::Id_type Rational_function_id_type;
typedef std::pair<Rational_function_id_type,Rational_function_id_type>
typedef std::pair<Rational_function_id_type, Rational_function_id_type>
Rational_function_canonicalized_pair_key;
class Less_compare_rational_function_pair_key
@ -70,11 +78,14 @@ public:
};
typedef typename std::map< Rational_function_canonicalized_pair_key,
Rational_function_canonicalized_pair,
Less_compare_rational_function_pair_key> Rational_function_canonicalized_pair_map;
Rational_function_canonicalized_pair,
Less_compare_rational_function_pair_key>
Rational_function_canonicalized_pair_map;
public:
Cache() : _rat_func_map_watermark(128),_rat_pair_map_watermark(128){};
void initialize(const Self& other, Algebraic_kernel& kernel = Algebraic_kernel())
void initialize(const Self& other,
Algebraic_kernel& kernel = Algebraic_kernel())
{
//copy rational function map
typename Rational_function_map::const_iterator iter1;
@ -85,8 +96,9 @@ public:
if (iter1->second.is_shared())
{
Rational_function_key key = iter1->first;
Rational_function f(iter1->second.numer(),
iter1->second.denom(),kernel); //construct new instance
//construct new instance
Rational_function f(iter1->second.numer(), iter1->second.denom(),
kernel);
_rat_func_map.insert(std::make_pair(key,f));
}
}
@ -102,8 +114,9 @@ public:
if (iter2->second.is_shared())
{
Rational_function_canonicalized_pair_key key = iter2->first;
//construct new instance
Rational_function_canonicalized_pair p(iter2->second.f(),
iter2->second.g(),kernel); //construct new instance
iter2->second.g(), kernel);
_rat_pair_map.insert(std::make_pair(key,p));
}
}
@ -120,9 +133,10 @@ public:
return _rat_pair_map;
}
const Rational_function& get_rational_function(const Polynomial_1& numer,
const Polynomial_1& denom,
Algebraic_kernel& kernel = Algebraic_kernel())
const Rational_function& get_rational_function(const Polynomial_1& numer,
const Polynomial_1& denom,
Algebraic_kernel& kernel =
Algebraic_kernel())
{
Rational_function_key key = get_key(numer,denom);
@ -141,7 +155,8 @@ public:
//then insert the new element
Rational_function f(numer,denom,kernel);
typename Rational_function_map::iterator it2 = _rat_func_map.insert(it,std::make_pair(key,f));
typename Rational_function_map::iterator it2 =
_rat_func_map.insert(it,std::make_pair(key,f));
return it2->second;
}
}
@ -151,22 +166,26 @@ public:
Integer numer,denom;
typename FT_rat_1::Decompose()(rat,numer,denom);
Polynomial_1 numer_poly = typename Polynomial_traits_1::Construct_polynomial()(numer);
Polynomial_1 denom_poly = typename Polynomial_traits_1::Construct_polynomial()(denom);
Polynomial_1 numer_poly =
typename Polynomial_traits_1::Construct_polynomial()(numer);
Polynomial_1 denom_poly =
typename Polynomial_traits_1::Construct_polynomial()(denom);
return get_rational_function (numer_poly,denom_poly,kernel);
}
const Rational_function_pair get_rational_pair ( const Rational_function& f,
const Rational_function& g,
Algebraic_kernel& kernel = Algebraic_kernel())
const Rational_function_pair get_rational_pair(const Rational_function& f,
const Rational_function& g,
Algebraic_kernel& kernel =
Algebraic_kernel())
{
CGAL_precondition(!(f==g));
Rational_function_canonicalized_pair_key key = get_key(f,g);
bool is_opposite = (f.id() < g.id()) ? false : true ;
//look if element exists in cache already
typename Rational_function_canonicalized_pair_map::iterator it = _rat_pair_map.lower_bound(key);
typename Rational_function_canonicalized_pair_map::iterator it =
_rat_pair_map.lower_bound(key);
if(it != _rat_pair_map.end() && !(_rat_pair_map.key_comp()(key, it->first)))
{
@ -192,15 +211,17 @@ public:
rat_pair_map_clean_up();
}
private:
Rational_function_key get_key(const Polynomial_1& numer, const Polynomial_1& denom) const
Rational_function_key get_key(const Polynomial_1& numer,
const Polynomial_1& denom) const
{
return Rational_function_key(numer, denom);
}
Rational_function_key get_key(const Rational_function& f) const
{
return get_key( f.numer(),f.denom());
return get_key(f.numer(), f.denom());
}
Rational_function_canonicalized_pair_key get_key(const Rational_function& f, const Rational_function& g) const
Rational_function_canonicalized_pair_key
get_key(const Rational_function& f, const Rational_function& g) const
{
return (f.id() < g.id()) ?
Rational_function_canonicalized_pair_key( f.id(),g.id() ):