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/Algebraic_point_2.h $ 7 // $Id: Algebraic_point_2.h 4e519a3 2021-05-05T13:15:37+02: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_ALBERAIC_POINT_D_1_H 14 #define CGAL_ALBERAIC_POINT_D_1_H 15 16 #include <CGAL/license/Arrangement_on_surface_2.h> 17 18 19 #include <ostream> 20 #include <CGAL/tss.h> 21 22 #include <CGAL/Arr_rat_arc/Base_rational_arc_ds_1.h> 23 #include <CGAL/Arr_rat_arc/Cache.h> 24 #include <CGAL/Arr_rat_arc/Rational_function.h> 25 #include <CGAL/Arr_rat_arc/Rational_function_canonicalized_pair.h> 26 27 #include <CGAL/Handle_with_policy.h> 28 29 #include <CGAL/Polynomial.h> 30 31 #include <CGAL/Arithmetic_kernel.h> 32 #include <CGAL/convert_to_bfi.h> 33 #include <CGAL/Bigfloat_interval_traits.h> 34 35 namespace CGAL { 36 namespace Arr_rational_arc { 37 38 //------------------- 39 //Algebraic_point_2_rep 40 //------------------- 41 template <typename AlgebraicKernelD_1> 42 class Algebraic_point_2_rep : public Base_rational_arc_ds_1<AlgebraicKernelD_1> 43 { 44 public: 45 typedef AlgebraicKernelD_1 Algebraic_kernel_d_1; 46 typedef Base_rational_arc_ds_1<Algebraic_kernel_d_1> 47 Base; 48 49 typedef CGAL::Arr_rational_arc::Rational_function<Algebraic_kernel_d_1> 50 Rational_function; 51 typedef CGAL::Arr_rational_arc::Rational_function_pair<Algebraic_kernel_d_1> 52 Rational_function_pair; 53 54 typedef typename Base::Algebraic_real_1 Algebraic_real_1; 55 typedef typename Algebraic_kernel_d_1::Bound Bound; 56 typedef typename Base::Integer Integer ; 57 typedef typename Base::Rational Rational ; 58 typedef typename Base::Polynomial_1 Polynomial_1; 59 60 typedef typename Base::Root_multiplicity_vector Root_multiplicity_vector; 61 62 typedef typename Get_arithmetic_kernel<Algebraic_real_1>::Arithmetic_kernel 63 AK; 64 typedef typename AK::Bigfloat_interval BFI; 65 typedef Bigfloat_interval_traits<BFI> BFI_traits; 66 typedef CGAL::Polynomial<BFI> BFI_polynomial; 67 typedef CGAL::Polynomial_traits_d<BFI_polynomial> BFI_polynomial_traits; 68 69 70 typedef typename Base::FT_rat_1 FT_rat_1; 71 typedef typename Base::Polynomial_traits_1 Polynomial_traits_1; 72 73 typedef CGAL::Arr_rational_arc::Cache<Algebraic_kernel_d_1> 74 Cache; 75 76 public: Algebraic_point_2_rep()77 Algebraic_point_2_rep(){} Algebraic_point_2_rep(const Rational_function & rational_function,const Algebraic_real_1 & x_coordinate)78 Algebraic_point_2_rep(const Rational_function& rational_function, 79 const Algebraic_real_1& x_coordinate) : 80 _rational_function(rational_function), 81 _x_coordinate(x_coordinate) {} 82 Algebraic_point_2_rep(const Algebraic_point_2_rep & other)83 Algebraic_point_2_rep(const Algebraic_point_2_rep& other) 84 { 85 if (this != &other) // protect against invalid self-assignment 86 { 87 _rational_function = other._rational_function; 88 _x_coordinate = other._x_coordinate; 89 } 90 } 91 92 //assignment oparator 93 Algebraic_point_2_rep& operator=(const Algebraic_point_2_rep& other) 94 { 95 if (this != &other) // protect against invalid self-assignment 96 { 97 _rational_function = other._rational_function; 98 _x_coordinate = other._x_coordinate; 99 } 100 return *this; 101 } 102 compare_xy_2(const Algebraic_point_2_rep & other,const Cache & cache)103 Comparison_result compare_xy_2(const Algebraic_point_2_rep& other, 104 const Cache& cache) const 105 { 106 Comparison_result comp = CGAL::compare(_x_coordinate, other.x()); 107 if (comp != EQUAL) 108 return comp; 109 if (_rational_function == other.rational_function()) 110 return EQUAL; 111 Rational_function_pair rat_func_pair = 112 cache.get_rational_pair(_rational_function, other.rational_function()); 113 return rat_func_pair.compare_f_g_at(_x_coordinate); 114 } 115 x()116 Algebraic_real_1& x() 117 { 118 return _x_coordinate; 119 } 120 x()121 const Algebraic_real_1& x() const 122 { 123 return _x_coordinate; 124 } 125 rational_function()126 const Rational_function& rational_function() const 127 { 128 return _rational_function; 129 } 130 numerator()131 const Polynomial_1& numerator() const { return _rational_function.numer(); } denominator()132 const Polynomial_1& denominator() const { return _rational_function.denom(); } 133 134 //new functions... y()135 Algebraic_real_1 y() const 136 { 137 typedef CGAL::Polynomial<Polynomial_1> Polynomial_2; 138 //converting the defining polynomial of x and the rational function to 139 //bivariate polynomials 140 Polynomial_2 f(_algebraic_kernel.compute_polynomial_1_object()(_x_coordinate)); 141 142 Polynomial_2 y(CGAL::shift(Polynomial_2(1),1)); 143 Polynomial_2 g(_rational_function.numer() - y * _rational_function.denom()); 144 145 f=CGAL::swap(f, 0, 1); //swap x and y in the polynomial f 146 g=CGAL::swap(g, 0, 1); //swap x and y in the polynomial g 147 //compute the resultant in x (polynomial in y) 148 Polynomial_1 r(CGAL::resultant(f,g)); 149 150 //solve for all roots of resultant 151 std::list<Algebraic_real_1> roots; 152 _algebraic_kernel.solve_1_object()(r, false, std::back_inserter(roots)); 153 //isolate the right root 154 unsigned int initial_precision = 16; 155 int error_bound = 2; 156 157 while (roots.size() > 1) 158 { 159 std::pair<Bound, Bound> 160 y_bounds(this->approximate_absolute_y(error_bound,initial_precision)); 161 while (CGAL::compare(roots.front(),y_bounds.first) == SMALLER) 162 roots.pop_front(); 163 164 while (CGAL::compare(y_bounds.second,roots.back()) == SMALLER) 165 roots.pop_back(); 166 167 error_bound *= 2; 168 } 169 170 CGAL_postcondition (roots.size() == 1); 171 return roots.front(); 172 } to_double()173 std::pair<double,double> to_double() const 174 { 175 double x = CGAL::to_double(_x_coordinate); 176 double numer_val = evaluate_at(_rational_function.numer(),x); 177 double denom_val = evaluate_at(_rational_function.denom(),x); 178 return std::make_pair(x,numer_val/denom_val); 179 } approximate_absolute_x(int a)180 std::pair<Bound,Bound> approximate_absolute_x( int a) const 181 { 182 return _algebraic_kernel.approximate_absolute_1_object()(_x_coordinate,a); 183 } approximate_absolute_y(int a)184 std::pair<Bound,Bound> approximate_absolute_y( int a) const 185 { 186 unsigned int precision = 16; 187 return approximate_absolute_y(a,precision); 188 } approximate_relative_x(int r)189 std::pair<Bound,Bound> approximate_relative_x( int r) const 190 { 191 return _algebraic_kernel.approximate_relative_1_object()(_x_coordinate,r); 192 } approximate_relative_y(int r)193 std::pair<Bound,Bound> approximate_relative_y( int r) const 194 { 195 // return zero if y is zero (the code below would not terminate) 196 // moreover approx relative is actually not well defined 197 if(_rational_function.sign_at(_x_coordinate)==CGAL::ZERO) 198 return std::make_pair(Bound(0),Bound(0)); 199 200 typename BFI_traits::Set_precision set_precision; 201 typename BFI_polynomial_traits::Evaluate evaluate; 202 203 typedef typename BFI_traits::Bound BF; 204 205 long precision = 16; 206 set_precision(precision); 207 BF eps = CGAL::ipower(BF(1)/2,r); 208 209 while (true){ 210 set_precision(precision); 211 BFI x_bfi(convert_to_bfi(_x_coordinate)); 212 213 BFI_polynomial 214 numer_bfi(convert_to_bfi_extended(_rational_function.numer())); 215 BFI_polynomial 216 denom_bfi(convert_to_bfi_extended(_rational_function.denom())); 217 218 BFI y_numer_bfi(evaluate(numer_bfi,x_bfi)); 219 BFI y_denom_bfi(evaluate(denom_bfi,x_bfi)); 220 221 if (CGAL::zero_in(y_denom_bfi) == false) 222 { 223 BFI y_bfi(y_numer_bfi/y_denom_bfi); 224 225 if (CGAL::compare( 226 CGAL::width(y_bfi), 227 CGAL::lower(CGAL::abs(y_bfi)) * eps) 228 == SMALLER) 229 return std::make_pair( 230 Bound(CGAL::lower(y_bfi)), 231 Bound(CGAL::upper(y_bfi))); 232 } 233 else precision*=2; 234 } 235 } 236 print(std::ostream & os)237 std::ostream& print (std::ostream& os) const 238 { 239 std::pair<double,double> double_p; 240 switch(::CGAL::IO::get_mode(os)) 241 { 242 case ::CGAL::IO::PRETTY: 243 double_p = this->to_double(); 244 os <<"(" ; 245 os << double_p.first; 246 os <<" , " ; 247 os << double_p.second; 248 os << ")"; 249 break; 250 251 case ::CGAL::IO::BINARY: 252 std::cerr << "BINARY format not yet implemented" << std::endl; 253 break; 254 255 default: 256 // ASCII 257 os <<"x = " << _x_coordinate<<" "; 258 os <<"rational function = ( " ; 259 Base::print_polynomial(os, _rational_function.numer(), 'x'); 260 os << ") / ("; 261 Base::print_polynomial(os, _rational_function.denom(), 'x'); 262 os << ")"; 263 } 264 265 return os; 266 } 267 private: approximate_absolute_y(int a,unsigned int & precision)268 std::pair<Bound,Bound> approximate_absolute_y( int a,unsigned int& precision ) const 269 { 270 typename BFI_traits::Set_precision set_precision; 271 typename BFI_polynomial_traits::Evaluate evaluate; 272 273 typedef typename BFI_traits::Bound BF; 274 275 BF eps = CGAL::ipower(BF(1)/2,a); 276 while (true) 277 { 278 set_precision(precision); 279 BFI x_bfi(convert_to_bfi(_x_coordinate)); 280 281 BFI_polynomial numer_bfi(convert_to_bfi_extended(_rational_function.numer())); 282 BFI_polynomial denom_bfi(convert_to_bfi_extended(_rational_function.denom())); 283 284 BFI y_numer_bfi(evaluate(numer_bfi,x_bfi)); 285 BFI y_denom_bfi(evaluate(denom_bfi,x_bfi)); 286 287 if (CGAL::zero_in(y_denom_bfi) == false) 288 { 289 BFI y_bfi(y_numer_bfi/y_denom_bfi); 290 if (CGAL::width(y_bfi) < eps ) 291 return std::make_pair( 292 Bound(CGAL::lower(y_bfi)), 293 Bound(CGAL::upper(y_bfi))); 294 295 } 296 else precision*=2; 297 } 298 } 299 300 template <typename NTX> 301 static typename 302 CGAL::Coercion_traits<NTX, 303 typename Get_arithmetic_kernel<NTX>:: 304 Arithmetic_kernel::Bigfloat_interval>::Type convert_to_bfi_extended(const NTX & x)305 convert_to_bfi_extended(const NTX& x) 306 { 307 typedef typename Get_arithmetic_kernel<NTX>::Arithmetic_kernel AK; 308 typedef typename AK::Bigfloat_interval BFI; 309 typedef CGAL::Coercion_traits<NTX, BFI> CT; 310 return typename CT::Cast()(x); 311 } 312 evaluate_at(const Polynomial_1 & poly,const double x)313 double evaluate_at(const Polynomial_1& poly, const double x) const 314 { 315 double x_val = 1; 316 double ret_val(0); 317 for (int i(0); i <= poly.degree(); ++i) 318 { 319 ret_val = ret_val + x_val*CGAL::to_double(poly[i]); 320 x_val = x_val*x; 321 } 322 return ret_val; 323 } 324 325 private: 326 Rational_function _rational_function; //supporting rational function 327 Algebraic_real_1 _x_coordinate; 328 329 Algebraic_kernel_d_1 _algebraic_kernel; 330 }; 331 332 template <typename Algebraic_kernel_ > 333 class Algebraic_point_2 : 334 public Handle_with_policy<Algebraic_point_2_rep<Algebraic_kernel_> > 335 { 336 337 public: 338 typedef Algebraic_kernel_ Algebraic_kernel_d_1; 339 typedef Handle_with_policy<Algebraic_point_2_rep<Algebraic_kernel_> > 340 Base; 341 typedef Algebraic_point_2<Algebraic_kernel_d_1> Self; 342 typedef Algebraic_point_2_rep<Algebraic_kernel_> Rep; 343 typedef typename Rep::Rational Rational; 344 typedef typename Rep::Algebraic_real_1 Algebraic_real_1; 345 typedef typename Rep::Rational_function Rational_function; 346 typedef typename Rep::Bound Bound; 347 typedef typename Rep::Cache Cache; 348 typedef typename Rep::Polynomial_1 Polynomial_1; 349 350 private: get_default_instance()351 static Self& get_default_instance() 352 { 353 CGAL_STATIC_THREAD_LOCAL_VARIABLE_0(Algebraic_kernel_d_1, kernel); 354 typedef typename Rational_function::Polynomial_1 Poly; 355 CGAL_STATIC_THREAD_LOCAL_VARIABLE(Poly, numer,0); 356 CGAL_STATIC_THREAD_LOCAL_VARIABLE(Poly, denom, 1); 357 CGAL_STATIC_THREAD_LOCAL_VARIABLE_3(Rational_function, rational_function, numer, denom, &kernel); 358 359 CGAL_STATIC_THREAD_LOCAL_VARIABLE(Algebraic_real_1, x_coordinate,kernel.construct_algebraic_real_1_object()(Rational(0))); 360 361 CGAL_STATIC_THREAD_LOCAL_VARIABLE_2(Self, default_instance, rational_function, x_coordinate); 362 363 return default_instance; 364 365 /*static Self x = Self(Rational(0),Rational(0),_cache); 366 return x; */ 367 } 368 369 public: Algebraic_point_2(const Rational_function & rational_function,const Algebraic_real_1 & x_coordinate)370 explicit Algebraic_point_2(const Rational_function& rational_function, 371 const Algebraic_real_1& x_coordinate) : 372 Base(rational_function,x_coordinate) {} 373 Algebraic_point_2()374 Algebraic_point_2() : 375 Base(static_cast<const Base &> (get_default_instance())) {} 376 compare_xy_2(const Algebraic_point_2 & other,const Cache & cache)377 Comparison_result compare_xy_2(const Algebraic_point_2& other, 378 const Cache& cache) const 379 { 380 if (this->is_identical (other)) 381 return CGAL::EQUAL; 382 return this->ptr()->compare_xy_2(*other.ptr(), cache); 383 } 384 numerator()385 const Polynomial_1& numerator() const { return this->ptr()->numerator(); } denominator()386 const Polynomial_1& denominator() const { return this->ptr()->denominator(); } 387 x()388 Algebraic_real_1& x() 389 { 390 if (this->is_shared()) 391 this->copy_on_write(); 392 return this->ptr()->x(); 393 } 394 x()395 const Algebraic_real_1& x() const 396 { 397 return this->ptr()->x(); 398 } 399 rational_function()400 const Rational_function& rational_function() const 401 { 402 return this->ptr()->rational_function(); 403 } 404 y()405 Algebraic_real_1 y() const 406 { 407 return this->ptr()->y(); 408 } 409 to_double()410 std::pair<double,double> to_double() const 411 { 412 return this->ptr()->to_double(); 413 } 414 approximate_absolute_x(int a)415 std::pair<Bound,Bound> approximate_absolute_x( int a) const 416 { 417 return this->ptr()->approximate_absolute_x(a); 418 } 419 approximate_absolute_y(int a)420 std::pair<Bound,Bound> approximate_absolute_y( int a) const 421 { 422 return this->ptr()->approximate_absolute_y(a); 423 } 424 approximate_relative_x(int r)425 std::pair<Bound,Bound> approximate_relative_x( int r) const 426 { 427 return this->ptr()->approximate_relative_x(r); 428 } 429 approximate_relative_y(int r)430 std::pair<Bound,Bound> approximate_relative_y( int r) const 431 { 432 return this->ptr()->approximate_relative_y(r); 433 } 434 print(std::ostream & os)435 std::ostream& print(std::ostream& os) const 436 { 437 return this->ptr()->print(os); 438 } 439 }; //Algebraic_point_2 440 441 442 template < typename Algebraic_kernel_> 443 std::ostream& operator<<(std::ostream& os, 444 const Algebraic_point_2<Algebraic_kernel_> & p) 445 { 446 return (p.print(os)); 447 } 448 449 } //namespace Arr_rational_arc 450 } //namespace CGAL { 451 452 #endif //CGAL_ALBERAIC_POINT_D_1_H 453