1 // Copyright (c) 2006-2009 Max-Planck-Institute Saarbruecken (Germany). 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/Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/Algebraic_real_rep_bfi.h $ 7 // $Id: Algebraic_real_rep_bfi.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Michael Hemmer <hemmer@mpi-inf.mpg.de> 12 // 13 // ============================================================================ 14 15 // TODO: The comments are all original EXACUS comments and aren't adapted. So 16 // they may be wrong now. 17 18 // This is code is expreimental ! 19 20 /*! \file NiX/Algebraic_real_rep_bfi.h 21 \brief Algebraic_real_rep with refinement via interval rep 22 */ 23 24 #ifndef CGAL_ALGEBRAIC_KERNEL_D_ALGEBRAIC_REAL_REP_BFI_H 25 #define CGAL_ALGEBRAIC_KERNEL_D_ALGEBRAIC_REAL_REP_BFI_H 26 27 #include <CGAL/basic.h> 28 #include <CGAL/Arithmetic_kernel.h> 29 #include <CGAL/Polynomial_type_generator.h> 30 #include <CGAL/Algebraic_kernel_d/Algebraic_real_rep.h> 31 #include <CGAL/Interval_traits.h> 32 #include <CGAL/Bigfloat_interval_traits.h> 33 #include <CGAL/convert_to_bfi.h> 34 #include <CGAL/Sqrt_extension_fwd.h> 35 #include <CGAL/Get_arithmetic_kernel.h> 36 37 namespace CGAL { 38 39 // it would be nice to remove the explicit use of Sqrt_extension 40 // in this file. However, it is much more efficient to convert the root 41 // only once. see convert_to_bfi 42 //template <class COEFF, class ROOT> class Sqrt_extension; //use forward declaration instead 43 44 namespace internal { 45 46 // definition of the Algebraic_real_rep_bfi x: 47 48 //IS_GENERAL: 49 // low_ lower bound of x 50 // high_ upper bound of x 51 // polynomial_ a square free polynomial 52 // sign_at_low_ = polynomial_.evaluate(low_) 53 // x is the only root of polynomial_ in the open interval ]low_,high_[ 54 // low_ != x != high 55 // ******************* EXEPTION ******************* 56 // x is rational: in this case low=high=x 57 58 59 template< class Coefficient_, class Field_> 60 class Algebraic_real_rep_bfi 61 : public Algebraic_real_rep<Coefficient_,Field_> { 62 63 typedef Algebraic_real_rep<Coefficient_,Field_> Base; 64 65 typedef Coefficient_ Coefficient; 66 typedef Field_ Field; 67 68 typedef typename CGAL::Get_arithmetic_kernel<Coefficient>::Arithmetic_kernel AT; 69 typedef typename AT::Bigfloat BF; 70 typedef typename AT::Bigfloat_interval BFI; 71 typedef typename AT::Field_with_sqrt FWS; 72 73 typedef typename CGAL::Polynomial_type_generator<Coefficient,1>::Type Poly; 74 typedef typename Poly::const_iterator PIterator; 75 76 typedef Algebraic_real_rep_bfi <Coefficient,Field> Self; 77 78 mutable std::vector<BFI> polynomial_approx; 79 mutable long current_prec; 80 81 private: 82 template <class Polynomial_1, class OI> 83 inline convert_coeffs(const Polynomial_1 & poly,OI it)84 void convert_coeffs(const Polynomial_1& poly, OI it) const { 85 typename CGAL::Polynomial_traits_d<Polynomial_1>::Get_coefficient 86 coeff; 87 for(int i = 0; i <= CGAL::degree(poly); i++){ 88 *it++ = convert_to_bfi(coeff(poly,i)); 89 } 90 } 91 92 template <class COEFF, class ROOT, class ACDE_TAG, class FP_TAG, class OI > 93 inline 94 void convert_coeffs(const CGAL::Polynomial<CGAL::Sqrt_extension<COEFF,ROOT,ACDE_TAG,FP_TAG>> & poly,OI it)95 convert_coeffs( 96 const CGAL::Polynomial< CGAL::Sqrt_extension<COEFF,ROOT,ACDE_TAG,FP_TAG> >& poly, 97 OI it ) const { 98 99 BFI root(0); 100 for(int i = 0; i <= CGAL::degree(poly); i++){ 101 if(poly[i].is_extended()){ 102 // typename Coercion_traits<FWS,ROOT >::Cast cast_root; 103 // root = convert_to_bfi(NiX::sqrt(cast_root(poly[i].root()))); 104 root = CGAL::sqrt(convert_to_bfi(poly[i].root())); 105 break; 106 } 107 } 108 109 for(int i = 0; i <= CGAL::degree(poly); i++){ 110 if(poly[i].is_extended()){ 111 *it++ = 112 convert_to_bfi(poly[i].a0()) + 113 convert_to_bfi(poly[i].a1()) * root ; 114 }else{ 115 *it++ = convert_to_bfi(poly[i].a0()); 116 } 117 } 118 } 119 120 typedef CGAL::Sign TRI_BOOL; 121 122 public: low()123 const Field& low() const {return this->low_;} high()124 const Field& high() const {return this->high_;} polynomial()125 const Poly& polynomial() const {return this->polynomial_;} sign_at_low()126 const TRI_BOOL& sign_at_low() const {return this->sign_at_low_;} 127 128 template <class NTX> strong_refine(const NTX & m)129 void strong_refine(const NTX& m) const{ 130 if(is_rational()) return; 131 132 if( NTX(low()) <= m && m <= NTX(high()) ){ 133 CGAL_precondition(polynomial().sign_at(m)!=CGAL::ZERO); 134 refine(); 135 while( NTX(low()) <= m && m <= NTX(high())) refine(); 136 } 137 } 138 is_rational()139 bool is_rational() const {return this->is_rational_;} rational()140 const Field& rational() const { 141 CGAL_precondition(is_rational()); 142 CGAL_precondition(low()==high()); 143 CGAL_precondition(polynomial().sign_at(low())==CGAL::ZERO); 144 return this->low_; 145 } 146 147 CGAL::Comparison_result 148 compare(const Field& y, bool are_distinct = false) const { 149 if(is_rational()) return CGAL::compare(rational(),y); 150 if(y <= low()) { 151 strong_refine(y); 152 return CGAL::LARGER; 153 } 154 if(high() <= y) { 155 strong_refine(y); 156 return CGAL::SMALLER; 157 } 158 159 // now: low < y < high 160 if(!are_distinct){ 161 if(sign_of_polynomial_at(y)==CGAL::ZERO){ 162 this->learn_from(y); 163 return CGAL::EQUAL ; 164 } 165 }else{ 166 CGAL_precondition(polynomial().sign_at(y)!=CGAL::ZERO); 167 } 168 strong_refine(y); 169 CGAL_postcondition(y < low() || high() < y ); 170 if(y < low()) return CGAL::LARGER; 171 else return CGAL::SMALLER; 172 } 173 174 CGAL::Comparison_result 175 compare (const Algebraic_real_rep_bfi& y, bool are_distinct = false) const{ 176 if( is_rational()) return -y.compare( rational()); 177 if(y.is_rational()) return compare(y.rational()); 178 179 // type of both x and y IS_GENERAL 180 if ( high() <= y.low() ) return CGAL::SMALLER; 181 if ( low() >= y.high()) return CGAL::LARGER; 182 183 // intersection isolating intervals is ]L,R[ 184 Field L = (low() > y.low() ) ? low() : y.low() ; 185 Field R = (high() < y.high()) ? high() : y.high() ; 186 187 // refine to smaller intervals at intersection interval boundaries 188 // this can change type() only to IS_RATIONAL 189 this->refine_at(L); 190 this->refine_at(R); 191 y.refine_at(L); 192 y.refine_at(R); 193 194 if ( is_rational()) return -y.compare( rational()); 195 if (y.is_rational()) return compare(y.rational()); 196 197 // type of both x and y still IS_GENERAL 198 if ( high() <= y.low() ) return CGAL::SMALLER; 199 if ( y.high() <= low() ) return CGAL::LARGER; 200 201 // filter 1 (optional): determine distinctness by refining intervals 202 #if NiX_REFINEMENTS_BEFORE_GCD > 0 203 if (!are_distinct) { 204 // we may want to refine a bit and hope for the best 205 // because computing the gcd is expensive 206 for (int ntries=0; ntries < NiX_REFINEMENTS_BEFORE_GCD; ntries++) { 207 y.refine(); 208 refine(); 209 if (y.is_rational()) return compare(y.rational()); 210 if ( is_rational()) return -y.compare( rational()); 211 if ( high() <= y.low()) return CGAL::SMALLER; 212 if (y.high() <= low()) return CGAL::LARGER; 213 } 214 } 215 #endif 216 // filter 2: probabilistically check coprimality 217 if (!are_distinct) { 218 are_distinct = !(may_have_common_factor(polynomial(), 219 y.polynomial())); 220 } 221 if (!are_distinct) { 222 // OK, filters failed. So we have to do the actual work 223 // and compute the gcd of the defining polynomials. 224 225 // we have ]low(), high()[ == ]y.low(),y.high()[ == ]L,R[ 226 // and let both numbers decide for the gcd or its complement 227 Poly F1,F2,G; 228 G = CGAL::gcd_up_to_constant_factor(polynomial(),y.polynomial()); 229 F1 = CGAL::integral_division_up_to_constant_factor(polynomial(),G); 230 CGAL_postcondition(CGAL::degree(F1)== 231 CGAL::degree(polynomial())-CGAL::degree(G)); 232 F2 = CGAL::integral_division_up_to_constant_factor(y.polynomial(),G); 233 CGAL_postcondition(CGAL::degree(F2)== 234 CGAL::degree(y.polynomial())-CGAL::degree(G)); 235 236 this->learn_from(G,F1); 237 y.learn_from(G,F2); 238 239 // this may simplify them due to degree loss 240 if (y.is_rational()) return compare(y.rational()); 241 if ( is_rational()) return -y.compare( rational()); 242 243 // type of x and y is still IS_GENERAL 244 // check for equality 245 if (G.sign_at(L)!=G.sign_at(R)){ 246 this->introduce(y); 247 return CGAL::EQUAL; 248 } 249 } 250 251 // if we are here, we know the numbers to be distinct 252 // refine to disjointness 253 for (;;) { 254 y.refine(); 255 refine(); 256 if (y.is_rational()) return compare(y.rational()); 257 if ( is_rational()) return -y.compare( rational()); 258 if ( high() <= y.low()) return CGAL::SMALLER; 259 if (y.high() <= low()) return CGAL::LARGER; 260 } 261 } 262 263 264 //! creates the algebraic real from int \a i. 265 explicit Algebraic_real_rep_bfi(int i = 0) Base(i)266 :Base(i){} 267 268 //! creates the algebraic real from Field \a m. Algebraic_real_rep_bfi(const Field & m)269 explicit Algebraic_real_rep_bfi(const Field& m) 270 :Base(m), current_prec(53) {} 271 272 /*! \brief creates the algebraic real as the unique root of \a P 273 in the open interval <var>]low,high[</var>. 274 \pre the polynomial \a P is square free 275 \pre P(low)!=0 276 \pre P(high)<0 277 \pre x is the one and only root in the open interval of \a P. 278 */ Algebraic_real_rep_bfi(const Poly & P,Field LOW,Field HIGH)279 Algebraic_real_rep_bfi(const Poly& P, Field LOW, Field HIGH): 280 Base(P,LOW,HIGH), current_prec(53){}; 281 282 283 //! copy constructor Algebraic_real_rep_bfi(const Self & y)284 Algebraic_real_rep_bfi(const Self& y) 285 : Base(y), current_prec(53){} 286 287 // assignment 288 Algebraic_real_rep_bfi& operator=(const Self& y) { 289 if ( this != & y) { 290 this->erase_from_list(); 291 this->copy_all_members(y); 292 this->introduce(y); 293 current_prec = y->current_prec; 294 } 295 NiX_expensive_postcond(this->self_test()); 296 return *this; 297 } 298 299 evaluate_polynomial_approx(const BFI & x)300 BFI evaluate_polynomial_approx(const BFI& x) const { 301 // std::cout << "eval approx begin"<< std::endl; 302 typedef std::vector<BFI> BFI_VEC; 303 typedef typename BFI_VEC::reverse_iterator RIT; 304 305 BFI result(0); 306 for(RIT rit = polynomial_approx.rbegin(); 307 rit != polynomial_approx.rend(); 308 rit++){ 309 result = result * x + (*rit); 310 } 311 // std::cout << "eval approx end"<< std::endl; 312 return result; 313 } 314 refine_poly_approximation()315 void refine_poly_approximation() const { 316 CGAL_precondition(current_prec > 1); 317 current_prec *= 2; 318 // std::cout <<"ALGREAL: refine approx: "<< current_prec<<std::endl; 319 set_precision(BFI(),current_prec); 320 polynomial_approx.clear(); 321 convert_coeffs( 322 this->polynomial(), 323 std::back_inserter(polynomial_approx)); 324 325 Self const *next = static_cast<const Self *>(this->next); 326 while(this != next){ 327 // std::cout << this << " " << next << std::endl; 328 next->polynomial_approx = polynomial_approx; 329 next->current_prec = current_prec; 330 next = static_cast<const Self *>(next->next); 331 } 332 }; 333 update_poly_approximation()334 void update_poly_approximation() const { 335 long old_prec = set_precision( BFI(), current_prec ); 336 337 polynomial_approx.clear(); 338 convert_coeffs( this->polynomial(), std::back_inserter( polynomial_approx ) ); 339 340 set_precision( BFI(), old_prec ); 341 342 // TODO: Problems if the next block gets executed 343 // Self const *next = static_cast<const Self *>(this->next); 344 // while(this != next){ 345 // // std::cout << this << " " << next << std::endl; 346 // next->polynomial_approx = polynomial_approx; 347 // next->current_prec = current_prec; 348 // next = static_cast<const Self *>(next->next); 349 // } 350 } 351 refine()352 void refine() const{ 353 if(this->is_rational()) return; 354 355 // std::cout << "refine begin ------- "<< std::endl; 356 357 Field m = (this->low()+this->high())/Field(2); 358 359 // Currently, sign_of_polynomial_at performs exactly the needed 360 // refinement. 361 // TODO: But what if this changes? 362 sign_of_polynomial_at( m ); 363 //std::cout << "refine end ----------------- "<<current_prec<< std::endl; 364 } 365 366 protected: sign_of_polynomial_at(const Field & f)367 virtual CGAL::Sign sign_of_polynomial_at( const Field& f ) const { 368 //return polynomial().sign_at( f ); 369 long old_prec = set_precision(BFI(),current_prec); 370 371 Field m = f; 372 CGAL::simplify(m); 373 374 // if polynomial has changed, the approximated polynomial gets 375 // updated 376 if ( (int) polynomial_approx.size() != 377 CGAL::degree(this->polynomial_)+1) { 378 update_poly_approximation(); 379 } 380 381 CGAL_postcondition(polynomial_approx.size() > 0); 382 383 BFI eval = evaluate_polynomial_approx(convert_to_bfi(m)); 384 385 CGAL::Sign s = CGAL::sign(CGAL::lower(eval)); 386 387 // correct sign if needed 388 if( s*CGAL::sign(CGAL::upper(eval) ) != CGAL::POSITIVE ){ 389 390 //std::cout << "APPROX FAILED-------------------------------"<<std::endl; 391 s = this->polynomial().sign_at(m); 392 if ( s != CGAL::ZERO ) { 393 refine_poly_approximation(); 394 } 395 } 396 397 CGAL_postcondition(s == this->polynomial_.sign_at(m)); 398 399 if ( s == CGAL::ZERO ) { 400 this->learn_from(m); 401 }else{ 402 if ( s == this->sign_at_low() ) this->low_ = m; 403 else this->high_ = m; 404 } 405 406 set_precision(BFI(),old_prec); 407 408 return s; 409 } 410 411 }; 412 413 }//namespace internal 414 415 } //namespace CGAL 416 417 #endif //CGAL_ALGEBRAIC_KERNEL_D_ALGEBRAIC_REAL_REP_BFI_H 418