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_canonicalized_pair.h $ 7 // $Id: Rational_function_canonicalized_pair.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_CANONICALIZED_PAIR_H 14 #define CGAL_RATIONAL_FUNCTION_CANONICALIZED_PAIR_H 15 16 #include <CGAL/license/Arrangement_on_surface_2.h> 17 18 19 #include <CGAL/Arr_rat_arc/Base_rational_arc_ds_1.h> 20 #include <CGAL/Arr_rat_arc/Rational_function.h> 21 #include <CGAL/Handle_with_policy.h> 22 #include <CGAL/assertions.h> 23 24 namespace CGAL { 25 namespace Arr_rational_arc { 26 27 template <typename AlgebraicKernel_d_1> 28 class Rational_function_canonicalized_pair_rep: 29 public Base_rational_arc_ds_1<AlgebraicKernel_d_1> 30 { 31 public: 32 typedef AlgebraicKernel_d_1 Algebraic_kernel_d_1; 33 typedef Base_rational_arc_ds_1<Algebraic_kernel_d_1> Base; 34 35 typedef CGAL::Arr_rational_arc::Rational_function<Algebraic_kernel_d_1> 36 Rational_function; 37 typedef typename Base::Polynomial_1 Polynomial_1; 38 typedef typename Base::Algebraic_real_1 Algebraic_real_1; 39 typedef typename Base::Algebraic_vector Algebraic_vector; 40 typedef typename Base::Multiplicity Multiplicity; 41 typedef typename Base::Multiplicity_vector Multiplicity_vector; 42 typedef typename Base::Root_multiplicity_vector Root_multiplicity_vector; 43 typedef typename Base::Solve_1 Solve_1; 44 typedef typename Base::Bound Bound; 45 typedef typename Base::Coefficient Coefficient; 46 47 public: Rational_function_canonicalized_pair_rep(const Rational_function & f,const Rational_function & g,Algebraic_kernel_d_1 * ak_ptr)48 Rational_function_canonicalized_pair_rep(const Rational_function& f, 49 const Rational_function& g, 50 Algebraic_kernel_d_1* ak_ptr) 51 :_f(f),_g(g),_ak_ptr(ak_ptr) 52 { 53 CGAL_precondition(_ak_ptr != nullptr); 54 //canonicalized representation 55 if ( !(f.id() < g.id()) ) 56 std::swap(_f,_g); 57 _resultant = (_f.numer() * _g.denom() - _f.denom() * _g.numer()); 58 59 //f and g are not the same... 60 CGAL_precondition(CGAL::is_zero(_resultant) == false); 61 62 Solve_1 solve_1(_ak_ptr->solve_1_object()); 63 Root_multiplicity_vector rm_vec; 64 65 #if 1 66 solve_1(_resultant,std::back_inserter(rm_vec)); 67 #else 68 CGAL_assertion(false); 69 // don't use this code yet since g and resultant/g may still have a 70 // common root 71 if( CGAL::internal::may_have_common_factor(_f.denom(),_g.denom())){ 72 Polynomial_1 g = CGAL::gcd(_f.denom(), _g.denom()); 73 solve_1(CGAL::integral_division(_resultant, g), 74 std::back_inserter(rm_vec)); 75 solve_1(g, std::back_inserter(rm_vec)); 76 std::sort(rm_vec.begin(), rm_vec.end()); 77 }else{ 78 solve_1(_resultant, std::back_inserter(rm_vec)); 79 } 80 #endif 81 82 _roots.reserve(rm_vec.size()); 83 _multiplicities.reserve(rm_vec.size()); 84 85 for (typename Root_multiplicity_vector::iterator it = rm_vec.begin(); 86 it != rm_vec.end() ; 87 ++it) 88 { 89 _roots.push_back(it->first); 90 _multiplicities.push_back(it->second); 91 } 92 93 Algebraic_vector tmp; 94 tmp.reserve(_f.poles().size() + _g.poles().size()); 95 _event_roots.reserve(rm_vec.size() + _f.poles().size() + _g.poles().size()); 96 //merge the roots of f,g & resultant 97 std::merge( _f.poles().begin(), _f.poles().end(), 98 _g.poles().begin(), _g.poles().end(), 99 std::back_inserter(tmp)); 100 std::merge(tmp.begin(), tmp.end(), 101 _roots.begin(), _roots.end(), 102 std::back_inserter(_event_roots)); 103 104 //remove duplicate entries 105 typename Algebraic_vector::iterator new_end,old_end; 106 new_end = std::unique(_event_roots.begin(), old_end=_event_roots.end()); 107 _event_roots.erase(new_end,old_end); 108 109 //compute the values of _is_above 110 bool curr_is_above = get_is_above_near_minus_infinity(); 111 _is_above.push_back(curr_is_above); 112 113 // TBD: is_sorted is provided both by stl and boost. The interface of 114 // boost version changed, and as of version ? it accepts a single 115 // argument, that is a model of the SinglePassRange concept. The stl 116 // version, on the other hand, is on it's way out. 117 // (The name is wrong-should be is_ordered). 118 // 119 // CGAL_precondition(std::is_sorted(_f.poles().begin(),_f.poles().end())); 120 // CGAL_precondition(std::is_sorted(_g.poles().begin(),_g.poles().end())); 121 // CGAL_precondition(std::is_sorted(_roots.begin(),_roots.end())); 122 123 typename Algebraic_vector ::const_iterator it_f_pole = _f.poles().begin(); 124 typename Multiplicity_vector::const_iterator it_f_mult = 125 _f.pole_multiplicities().begin(); 126 typename Algebraic_vector ::const_iterator it_g_pole = _g.poles().begin(); 127 typename Multiplicity_vector::const_iterator it_g_mult = 128 _g.pole_multiplicities().begin(); 129 typename Algebraic_vector ::const_iterator it_r_root = _roots.begin(); 130 typename Multiplicity_vector::const_iterator it_r_mult = 131 _multiplicities.begin() ; 132 while ( (it_f_pole != _f.poles().end()) || 133 (it_g_pole != _g.poles().end()) || 134 (it_r_root != _roots.end())) 135 { 136 //std::cout << "_f._poles.size() " << _f._poles.size() <<std::endl; 137 //std::cout << "_g._poles.size() " << _g._poles.size()<< std::endl; 138 //std::cout << "_roots.size() " << _roots.size() << std::endl; 139 //if current event is only a pole of f 140 if ((it_f_pole != _f.poles().end()) && 141 ((it_g_pole == _g.poles().end()) || (*it_f_pole < *it_g_pole)) && 142 ((it_r_root == _roots.end()) || (*it_f_pole < *it_r_root))) 143 { 144 if (*it_f_mult % 2 == 1) 145 curr_is_above = (curr_is_above == true ) ? false : true; 146 _is_above.push_back(curr_is_above); 147 ++it_f_pole; 148 ++it_f_mult; 149 } 150 //if current event is only a pole of g 151 else if ((it_g_pole != _g.poles().end()) && 152 ((it_f_pole == _f.poles().end()) || (*it_g_pole < *it_f_pole)) && 153 ((it_r_root == _roots.end()) || (*it_g_pole < *it_r_root))) 154 { 155 if (*it_g_mult % 2 == 1) 156 curr_is_above = (curr_is_above == true ) ? false : true; 157 _is_above.push_back(curr_is_above); 158 ++it_g_pole; 159 ++it_g_mult; 160 } 161 //if current event is only a pole of r 162 else if ((it_r_root != _roots.end()) && 163 ((it_f_pole == _f.poles().end()) || (*it_r_root < *it_f_pole)) && 164 ((it_g_pole == _g.poles().end()) || (*it_r_root < *it_g_pole))) 165 { 166 if (*it_r_mult % 2 == 1) 167 curr_is_above = (curr_is_above == true ) ? false : true; 168 _is_above.push_back(curr_is_above); 169 ++it_r_root; 170 ++it_r_mult; 171 } 172 //if current event is a pole of f g and r 173 else if ( (it_f_pole != _f.poles().end()) && 174 (it_g_pole != _g.poles().end()) && 175 (it_r_root != _roots.end()) && 176 (*it_r_root == *it_f_pole) && 177 (*it_r_root == *it_g_pole)) 178 { 179 //both functions switch signs 180 if (_f.sign_at(*it_r_root, CGAL::NEGATIVE) == 181 _g.sign_at(*it_r_root, CGAL::NEGATIVE)) 182 curr_is_above = !curr_is_above; 183 if (*it_r_mult % 2 == 1) 184 curr_is_above = !curr_is_above; 185 if (_f.sign_at(*it_r_root, CGAL::POSITIVE) == 186 _g.sign_at(*it_r_root, CGAL::POSITIVE)) 187 curr_is_above = !curr_is_above; 188 _is_above.push_back(curr_is_above); 189 ++it_f_pole; 190 ++it_f_mult; 191 ++it_g_pole; 192 ++it_g_mult; 193 ++it_r_root; 194 ++it_r_mult; 195 } 196 else 197 { 198 //should not be reached 199 CGAL_postcondition_msg(false,"invalid case in computing _is_above"); 200 } 201 } 202 //std::cout << "_is_above.size() " << _is_above.size() <<std::endl; 203 //std::cout << " _event_roots.size() + 1 " << _event_roots.size() + 1 << std::endl; 204 CGAL_postcondition(_is_above.size() == _event_roots.size() + 1); 205 //check for validity using explicit computation 206 CGAL_postcondition_code(std::vector<bool> tmp_is_above = 207 compute_is_above_explicitly(); 208 ); 209 CGAL_postcondition(_is_above == tmp_is_above); 210 } 211 212 Comparison_result compare_f_g_at(const Algebraic_real_1& x, 213 CGAL::Sign epsilon = CGAL::ZERO) const 214 { 215 //f and g must be different 216 CGAL_precondition(!CGAL::is_zero(_resultant)); 217 218 //find interval 219 typename Algebraic_vector::const_iterator iter = 220 std::lower_bound(_event_roots.begin(), _event_roots.end(),x); 221 222 //case of a value larger than largest root 223 if (iter == _event_roots.end()) 224 return (_is_above.back() ? CGAL:: LARGER : CGAL::SMALLER); 225 226 typename Algebraic_vector::iterator::difference_type dist = 227 iter - _event_roots.begin(); 228 229 //if x is not a root, ignore epsilons 230 if (*iter != x){ 231 return (_is_above[dist] ? CGAL:: LARGER : CGAL::SMALLER); 232 } 233 234 //x is a root 235 if (epsilon == CGAL::ZERO){ 236 CGAL_precondition(_f.poles().end() == 237 std::find(_f.poles().begin(), _f.poles().end(),x)); 238 CGAL_precondition(_g.poles().end() == 239 std::find(_g.poles().begin(), _g.poles().end(),x)); 240 return (CGAL::EQUAL); 241 } 242 243 if (epsilon == CGAL::NEGATIVE) 244 return (_is_above[dist] ? CGAL:: LARGER : CGAL::SMALLER); 245 else // CGAL::POSITIVE 246 return (_is_above[dist+1] ? CGAL:: LARGER : CGAL::SMALLER); 247 } compare_f_g_at(Arr_parameter_space boundary)248 Comparison_result compare_f_g_at(Arr_parameter_space boundary) const 249 { 250 CGAL_precondition((boundary == ARR_LEFT_BOUNDARY) || 251 (boundary == ARR_RIGHT_BOUNDARY) ); 252 253 //f and g are the same... 254 if (CGAL::is_zero(_resultant)) 255 return CGAL::EQUAL; 256 257 if (boundary == ARR_LEFT_BOUNDARY) 258 return _is_above.front() ? CGAL::LARGER : CGAL::SMALLER ; 259 else // boundary = ARR_RIGHT_BOUNDARY 260 return _is_above.back() ? CGAL::LARGER : CGAL::SMALLER ; 261 } 262 is_intersecting_in_range(const Arr_parameter_space left_parameter_space,const Algebraic_real_1 left,const Arr_parameter_space right_parameter_space,const Algebraic_real_1 right)263 bool is_intersecting_in_range(const Arr_parameter_space left_parameter_space, 264 const Algebraic_real_1 left, 265 const Arr_parameter_space right_parameter_space, 266 const Algebraic_real_1 right) const 267 { 268 //the two function intersect iff the left index and the right index 269 //of the array _is_above are different 270 271 //get left index 272 typename Algebraic_vector::const_iterator left_index = 273 (left_parameter_space == ARR_LEFT_BOUNDARY) ? 274 _event_roots.begin(): 275 std::lower_bound(_event_roots.begin(), _event_roots.end(),left); 276 //check if intersect at left index 277 if (*left_index == left) 278 return true; 279 280 //get right index 281 typename Algebraic_vector::const_iterator right_index = 282 (right_parameter_space == ARR_RIGHT_BOUNDARY) ? 283 _event_roots.begin(): 284 std::lower_bound(_event_roots.begin(), _event_roots.end(),right); 285 //check if intersect at right index 286 if (*right_index == right) 287 return true; 288 289 //check if indices are the same 290 return (left_index == right_index); 291 } 292 f()293 const Rational_function& f() const 294 { 295 return _f; 296 } 297 g()298 const Rational_function& g() const 299 { 300 return _g; 301 } 302 roots()303 const Algebraic_vector & roots() const 304 { 305 return _roots; 306 } 307 multiplicities()308 const Multiplicity_vector & multiplicities() const 309 { 310 return _multiplicities; 311 } 312 313 private: get_is_above_near_minus_infinity()314 bool get_is_above_near_minus_infinity() 315 { 316 bool r = _get_is_above_near_minus_infinity(); 317 CGAL_postcondition(r == __get_is_above_near_minus_infinity()); 318 return r; 319 } 320 _get_is_above_near_minus_infinity()321 bool _get_is_above_near_minus_infinity() 322 { 323 int f_deg = CGAL::degree(_f.numer()) - CGAL::degree(_f.denom()); 324 int g_deg = CGAL::degree(_g.numer()) - CGAL::degree(_g.denom()); 325 326 if (f_deg > g_deg) //f is stronger than g 327 { 328 CGAL::Sign sign = _f.sign_near_minus_infinity(); 329 return (sign == CGAL::NEGATIVE) ? false : 330 (sign == CGAL::POSITIVE) ? true : 331 (_g.sign_near_minus_infinity() == CGAL::NEGATIVE); // _f == zero; 332 } 333 if (f_deg < g_deg) //g is stronger than f 334 { 335 CGAL::Sign sign = _g.sign_near_minus_infinity(); 336 return (sign == CGAL::NEGATIVE) ? true : 337 (sign == CGAL::POSITIVE) ? false : 338 (_f.sign_near_minus_infinity() == CGAL::POSITIVE); // _g == zero; 339 } 340 341 //both have the same degree difference, 342 //check who's leading coeeficient ratio is larger 343 Coefficient lead_coeff_ratio = 344 CGAL::abs(CGAL::leading_coefficient(_f.numer())) * 345 CGAL::leading_coefficient(_g.denom()) - 346 CGAL::abs(CGAL::leading_coefficient(_g.numer())) * 347 CGAL::leading_coefficient(_f.denom()); 348 if (lead_coeff_ratio > 0)//f is stronger than g 349 return (_f.sign_near_minus_infinity() == CGAL::NEGATIVE ) ? false : true; 350 if (lead_coeff_ratio < 0) //g is stronger than f 351 return (_g.sign_near_minus_infinity() == CGAL::NEGATIVE ) ? true : false; 352 353 //ratio is the same, instead of continuing with taylor expansion, 354 //compute explicitly 355 356 return __get_is_above_near_minus_infinity(); 357 } 358 __get_is_above_near_minus_infinity()359 bool __get_is_above_near_minus_infinity() 360 { 361 Bound b; 362 if (_event_roots.empty()) 363 b = Bound(0); 364 else 365 { 366 b = (_ak_ptr->approximate_relative_1_object()(_event_roots.front(), 0)).first - 1; //lower bound of first root 367 } 368 return is_above_at(b); 369 } 370 is_above_at(const Bound & b)371 bool is_above_at(const Bound& b) 372 { 373 //return true if f is above g at b which means return (sign == positive) of : 374 // 375 // numer_f * denom _g - numer_f * denom _g 376 //f-g= -------------------------------------- at b 377 // denom _f * denom _g 378 379 //TODO: unnescecary construction of real 380 Algebraic_real_1 x(_ak_ptr->construct_algebraic_real_1_object()(b)); 381 382 CGAL::Sign numer = _ak_ptr->sign_at_1_object()(_resultant, x); 383 CGAL::Sign denom_f = _ak_ptr->sign_at_1_object()(_f.denom(),x); 384 CGAL::Sign denom_g = _ak_ptr->sign_at_1_object()(_g.denom(),x); 385 CGAL::Sign denom = denom_f * denom_g ; 386 CGAL::Sign s = numer*denom; 387 388 CGAL_precondition(s != CGAL::ZERO); 389 390 return (s == CGAL::POSITIVE); 391 } compute_is_above_explicitly()392 std::vector<bool> compute_is_above_explicitly() 393 { 394 std::vector<bool> tmp_is_above; 395 if (_event_roots.size()== 0) 396 { 397 Bound b = 1; //all bound are legal, choose 1 for simplicity 398 tmp_is_above.push_back(is_above_at(b)); 399 return tmp_is_above; 400 } 401 402 tmp_is_above.reserve(_event_roots.size()+1); 403 //left boundary 404 Bound b = (_ak_ptr->approximate_relative_1_object() 405 (_event_roots.front(),0)).first - 1; //lower bound of first root 406 tmp_is_above.push_back(is_above_at(b)); 407 408 //mid intervals 409 typename Algebraic_vector::size_type i; 410 for (i = 0; i < _event_roots.size()-1; ++i) 411 { 412 b = _ak_ptr->bound_between_1_object()(_event_roots[i],_event_roots[i+1]); 413 tmp_is_above.push_back(is_above_at(b)); 414 } 415 416 //right boundary 417 b = (_ak_ptr->approximate_relative_1_object()(_event_roots.back(), 0)).second + 1; //lower bound of last root 418 tmp_is_above.push_back(is_above_at(b)); 419 return tmp_is_above; 420 } 421 422 private: 423 Rational_function _f,_g; 424 Polynomial_1 _resultant; 425 Algebraic_vector _roots; 426 //roots of resultant merged with roots of f & g's denomenators 427 Algebraic_vector _event_roots; 428 Multiplicity_vector _multiplicities; 429 //is f above g in the interval induced by index i of _roots 430 std::vector<bool> _is_above; 431 Algebraic_kernel_d_1* _ak_ptr; 432 }; // Rational_function_canonicalized_pair_rep 433 434 template <typename Algebraic_kernel_> 435 class Rational_function_canonicalized_pair: 436 public Handle_with_policy<Rational_function_canonicalized_pair_rep 437 <Algebraic_kernel_> > 438 { 439 public: 440 typedef Algebraic_kernel_ Algebraic_kernel_d_1; 441 typedef Rational_function_canonicalized_pair_rep<Algebraic_kernel_d_1> 442 Rep; 443 typedef Handle_with_policy<Rep> Base; 444 typedef Rational_function_canonicalized_pair<Algebraic_kernel_d_1> 445 Self; 446 typedef typename Rep::Rational_function Rational_function; 447 typedef typename Rep::Algebraic_real_1 Algebraic_real_1; 448 typedef typename Rep::Polynomial_1 Polynomial_1; 449 typedef typename Rep::Algebraic_vector Algebraic_vector; 450 typedef typename Rep::Multiplicity_vector Multiplicity_vector; 451 typedef typename Rep::Root_multiplicity_vector Root_multiplicity_vector; 452 Rational_function_canonicalized_pair(const Rational_function & f,const Rational_function & g,Algebraic_kernel_d_1 * ak_ptr)453 Rational_function_canonicalized_pair(const Rational_function& f, 454 const Rational_function& g, 455 Algebraic_kernel_d_1* ak_ptr) : 456 Base(f, g, ak_ptr) {} 457 //used to solve VS bug... Rational_function_canonicalized_pair(const Self & p)458 Rational_function_canonicalized_pair(const Self & p) : 459 Base(static_cast<const Base &> (p)) {} 460 461 Comparison_result compare_f_g_at(const Algebraic_real_1& x, 462 CGAL::Sign epsilon = CGAL::ZERO) const 463 { 464 return this->ptr()->compare_f_g_at(x,epsilon); 465 } 466 compare_f_g_at(Arr_parameter_space boundary)467 Comparison_result compare_f_g_at(Arr_parameter_space boundary) const 468 { 469 return this->ptr()->compare_f_g_at(boundary); 470 } 471 is_intersecting_in_range(const Arr_parameter_space left_parameter_space,const Algebraic_real_1 left,const Arr_parameter_space right_parameter_space,const Algebraic_real_1 right)472 bool is_intersecting_in_range(const Arr_parameter_space left_parameter_space, 473 const Algebraic_real_1 left, 474 const Arr_parameter_space right_parameter_space, 475 const Algebraic_real_1 right) const 476 { 477 return this->ptr()->is_intersecting_in_range(left_parameter_space, left, 478 right_parameter_space, right); 479 } 480 f()481 const Rational_function& f() const 482 { 483 return this->ptr()->f(); 484 } 485 g()486 const Rational_function& g() const 487 { 488 return this->ptr()->g(); 489 } 490 roots()491 const Algebraic_vector & roots() const 492 { 493 return this->ptr()->roots(); 494 } 495 multiplicities()496 const Multiplicity_vector & multiplicities() const 497 { 498 return this->ptr()->multiplicities(); 499 } 500 }; //Rational_function_canonicalized_pair 501 502 } //namespace Arr_rational_arc 503 } //namespace CGAL { 504 505 #endif //CGAL_RATIONAL_FUNCTION_CANONICALIZED_PAIR_H 506