1 // Copyright (c) 2011 Tel-Aviv University (Israel), INRIA Sophia-Antipolis (France). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Arrangement_on_surface_2/include/CGAL/Arr_rat_arc/Rational_function.h $ 7 // $Id: Rational_function.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // Author(s) : Oren Salzman <orenzalz@post.tau.ac.il > 11 // Michael Hemmer <Michael.Hemmer@sophia.inria.fr> 12 13 #ifndef CGAL_RATIONAL_FUNCTION_H 14 #define CGAL_RATIONAL_FUNCTION_H 15 16 #include <CGAL/license/Arrangement_on_surface_2.h> 17 18 19 #include <CGAL/tss.h> 20 #include <CGAL/Arr_rat_arc/Base_rational_arc_ds_1.h> 21 #include <CGAL/Handle_with_policy.h> 22 namespace CGAL { 23 namespace Arr_rational_arc { 24 25 template <class AlgebraicKernel_d_1 > 26 class Rational_function_rep : public Base_rational_arc_ds_1<AlgebraicKernel_d_1> 27 { 28 public: 29 typedef AlgebraicKernel_d_1 Algebraic_kernel_d_1; 30 typedef Base_rational_arc_ds_1<Algebraic_kernel_d_1> Base; 31 32 typedef typename Base::Polynomial_1 Polynomial_1; 33 typedef typename Base::Algebraic_real_1 Algebraic_real_1; 34 typedef typename Base::Algebraic_vector Algebraic_vector; 35 typedef typename Base::Multiplicity Multiplicity; 36 typedef typename Base::Multiplicity_vector Multiplicity_vector; 37 typedef typename Base::Root_multiplicity_vector Root_multiplicity_vector; 38 typedef typename Base::Solve_1 Solve_1; 39 40 public: Rational_function_rep()41 Rational_function_rep() : _ak_ptr(nullptr){} Rational_function_rep(const Polynomial_1 & numer,const Polynomial_1 & denom,Algebraic_kernel_d_1 * ak_ptr)42 Rational_function_rep(const Polynomial_1& numer, 43 const Polynomial_1& denom, 44 Algebraic_kernel_d_1* ak_ptr): 45 _numer(numer), _denom(denom),_ak_ptr(ak_ptr) 46 { 47 initialize(); 48 } 49 50 CGAL::Sign sign_at (const Algebraic_real_1& x, 51 CGAL::Sign epsilon = CGAL::ZERO) const 52 { 53 //find interval 54 typename Algebraic_vector::const_iterator iter = 55 std::lower_bound(_event_roots.begin(), _event_roots.end(),x); 56 57 //case of a value larger than largest root 58 if (iter == _event_roots.end()) 59 return (_sign.back()); 60 61 typename Algebraic_vector::iterator::difference_type dist = 62 iter - _event_roots.begin(); 63 64 //if x is not a root, ignore epsilons 65 if (*iter != x) 66 return (_sign[dist]); 67 68 //x is a root 69 if (epsilon == CGAL::ZERO) 70 return (CGAL::EQUAL); 71 else if (epsilon == CGAL::NEGATIVE) 72 return (_sign[dist] ); 73 else // CGAL::POSITIVE 74 return (_sign[dist+1]); 75 } 76 sign_near_minus_infinity()77 CGAL::Sign sign_near_minus_infinity() const 78 { 79 return _sign.front(); 80 } 81 82 bool operator==(const Rational_function_rep& other) const 83 { 84 return ((this->_numer == other.numer()) && 85 (this->_denom == other.denom()) ); 86 } 87 numer()88 const Polynomial_1& numer() const 89 { 90 return _numer; 91 } 92 denom()93 const Polynomial_1& denom() const 94 { 95 return _denom; 96 } 97 poles()98 const Algebraic_vector& poles() const 99 { 100 return _poles; 101 } 102 pole_multiplicities()103 const Multiplicity_vector& pole_multiplicities() const 104 { 105 return _pole_multiplicities; 106 } 107 108 private: initialize()109 void initialize() 110 { 111 CGAL_precondition(_ak_ptr != nullptr); 112 CGAL_precondition(CGAL::is_zero(_denom) == false); 113 if (CGAL::is_zero(_numer)) 114 { 115 //function does not change sign 116 _sign.push_back(CGAL::ZERO); 117 return; 118 } 119 120 Solve_1 solve_1 (_ak_ptr->solve_1_object()); 121 Root_multiplicity_vector rm_poles_vec,rm_intersctions_vec; 122 solve_1(_denom, std::back_inserter(rm_poles_vec)); //poles 123 solve_1(_numer, std::back_inserter(rm_intersctions_vec)); //intersections with zero 124 125 //reserve memory 126 typename Root_multiplicity_vector::size_type num_of_poles = 127 rm_poles_vec.size(); 128 typename Root_multiplicity_vector::size_type num_of_intersections = 129 rm_intersctions_vec.size(); 130 131 _poles.reserve(num_of_poles); 132 _pole_multiplicities.reserve(num_of_poles); 133 134 _event_roots.reserve(num_of_poles + num_of_intersections); 135 _pole_multiplicities.reserve(num_of_poles + num_of_intersections); 136 137 //initialize poles 138 for ( typename Root_multiplicity_vector::iterator it = rm_poles_vec.begin(); 139 it != rm_poles_vec.end() ; 140 ++it) 141 { 142 _poles.push_back(it->first); 143 _pole_multiplicities.push_back(it->second); 144 } 145 146 //initialize events 147 Root_multiplicity_vector events_mult_vec; 148 std::merge( rm_poles_vec.begin() , rm_poles_vec.end() , 149 rm_intersctions_vec.begin() , rm_intersctions_vec.end() , 150 std::back_inserter(events_mult_vec)); 151 152 typename Root_multiplicity_vector::iterator it; 153 for (it = events_mult_vec.begin(); it != events_mult_vec.end(); ++it) 154 { 155 _event_roots.push_back(it->first); 156 _event_multiplicities.push_back(it->second); 157 } 158 159 //initialize left most interval (at -oo) 160 //for (ax^n+.../x^m+...) the sign is: sign(a) * (-1)^(n +m) 161 CGAL::Sign curr_sign = CGAL::sign(CGAL::leading_coefficient(_numer)); 162 if ((CGAL::degree(_numer) + CGAL::degree(_denom)) % 2 == 1) 163 curr_sign = curr_sign * CGAL::NEGATIVE; 164 _sign.push_back(curr_sign); 165 166 typename Multiplicity_vector::iterator it2; 167 for (it2 = _event_multiplicities.begin(); 168 it2 != _event_multiplicities.end(); 169 ++it2) 170 { 171 if (*it2 % 2 == 1) 172 curr_sign = curr_sign * CGAL::NEGATIVE; 173 _sign.push_back(curr_sign); 174 } 175 } 176 177 private: 178 Polynomial_1 _numer; 179 Polynomial_1 _denom; 180 Algebraic_vector _poles; //roots of the denominator 181 Multiplicity_vector _pole_multiplicities; //multiplicities of the poles 182 Algebraic_vector _event_roots; //poles and intersection points with y=0, function can change signs in these events 183 Multiplicity_vector _event_multiplicities; //multiplicities of the events 184 std::vector<CGAL::Sign> _sign; //function's sign in the corresponding interval induced by _event_roots (if no roots then only one value) 185 mutable Algebraic_kernel_d_1* _ak_ptr; 186 187 };//Rational_function_rep 188 189 template < class Algebraic_kernel_ > 190 class Rational_function: 191 public Handle_with_policy<Rational_function_rep<Algebraic_kernel_> > 192 { 193 public: 194 typedef Algebraic_kernel_ Algebraic_kernel_d_1; 195 typedef Handle_with_policy<Rational_function_rep<Algebraic_kernel_> > 196 Base; 197 typedef Rational_function<Algebraic_kernel_d_1> Self; 198 typedef Rational_function_rep<Algebraic_kernel_> Rep; 199 typedef typename Rep::Algebraic_real_1 Algebraic_real_1; 200 typedef typename Rep::Polynomial_1 Polynomial_1; 201 typedef typename Rep::Algebraic_vector Algebraic_vector; 202 typedef typename Rep::Multiplicity_vector Multiplicity_vector; 203 204 typedef typename Base::Id_type Id_type; 205 private: get_default_instance()206 static Self& get_default_instance() 207 { 208 CGAL_STATIC_THREAD_LOCAL_VARIABLE_0(Algebraic_kernel_d_1, kernel); 209 CGAL_STATIC_THREAD_LOCAL_VARIABLE_3(Self, x, Polynomial_1(0), Polynomial_1(1), &kernel); 210 return x; 211 } 212 public: Rational_function(const Polynomial_1 & numer,const Polynomial_1 & denom,Algebraic_kernel_d_1 * ak_ptr)213 Rational_function(const Polynomial_1& numer, 214 const Polynomial_1& denom, 215 Algebraic_kernel_d_1* ak_ptr) : 216 Base(numer,denom,ak_ptr) {} 217 218 //used to solve VS bug... Rational_function()219 Rational_function () : 220 Base(static_cast<const Base &> (get_default_instance())) {} 221 222 // explicit copy-constructor, required by VC9 Rational_function(const Self & r)223 Rational_function (const Self & r) 224 : Base(static_cast<const Base &> (r)) {} 225 226 Self& operator=(const Self&)=default; 227 228 CGAL::Sign sign_at(const Algebraic_real_1& x, 229 CGAL::Sign epsilon = CGAL::ZERO) const 230 { 231 return this->ptr()->sign_at(x,epsilon); 232 } 233 sign_near_minus_infinity()234 CGAL::Sign sign_near_minus_infinity() const 235 { 236 return this->ptr()->sign_near_minus_infinity (); 237 } 238 bool operator== (const Self & other) const 239 { 240 if (this->is_identical (other)) 241 return true; 242 return (*(this->ptr()) == *(other.ptr())); 243 } 244 numer()245 const Polynomial_1& numer () const 246 { 247 return this->ptr()->numer(); 248 } 249 denom()250 const Polynomial_1& denom () const 251 { 252 return this->ptr()->denom(); 253 } 254 poles()255 const Algebraic_vector& poles () const 256 { 257 return this->ptr()->poles(); 258 } 259 pole_multiplicities()260 const Multiplicity_vector& pole_multiplicities () const 261 { 262 return this->ptr()->pole_multiplicities(); 263 } 264 }; //Rational_function 265 266 } //namespace Arr_rational_arc 267 } //namespace CGAL { 268 269 #endif //CGAL_RATIONAL_FUNCTION_H 270