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/Descartes.h $ 7 // $Id: Descartes.h 26355e2 2020-06-25T12:31:21+02:00 Mael Rouxel-Labbé 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 /*! \file NiX/Descartes.h 19 \brief defines class NiX::Descartes. 20 21 Isolate real roots of polynomials. 22 23 This file provides a class to isolate real roots of polynomials, 24 using the algorithm based on the method of Descartes. 25 26 The polynomial has to be a univariat polynomial over any number 27 type which is contained in the real numbers. 28 29 */ 30 31 #ifndef CGAL_ALGEBRAIC_KERNEL_D_DESCARTES_H 32 #define CGAL_ALGEBRAIC_KERNEL_D_DESCARTES_H 33 34 #include <CGAL/basic.h> 35 #include <CGAL/Polynomial.h> 36 37 #include <CGAL/Algebraic_kernel_d/univariate_polynomial_utils.h> 38 #include <CGAL/Algebraic_kernel_d/construct_binary.h> 39 40 #define POLYNOMIAL_REBIND( coeff ) \ 41 typename CGAL::Polynomial_traits_d<Polynomial>::template \ 42 Rebind<coeff,1>::Other::Type 43 44 45 namespace CGAL { 46 47 namespace internal { 48 49 /*! \ingroup NiX_Algebraic_real 50 * \brief A model of concept RealRootIsolator. 51 */ 52 template <class Polynomial_, class Rational_> 53 class Descartes { 54 typedef CGAL::Fraction_traits<Polynomial_> FT_poly; 55 typedef Fraction_traits<Rational_> FT_rat; 56 public: 57 //! First template parameter 58 typedef Polynomial_ Polynomial; 59 //! Second template parameter 60 typedef Rational_ Rational; 61 //! Bound type of the isolating intervals 62 typedef Rational_ Bound; 63 // Integer or Numerator/Denominator type of bound. 64 typedef typename CGAL::Fraction_traits<Rational>::Numerator_type Integer; 65 private: 66 typedef typename Polynomial::NT Coeff; 67 typedef Integer IT; 68 69 Polynomial poly_; 70 int number_of_real_roots_; 71 IT* numerator; 72 IT* denominator_exponent; 73 bool* is_exact; 74 IT LEFT,SCALE,DENOM; 75 bool is_strong_; 76 int k; 77 bool interval_given; 78 79 public: 80 /*! \brief Constructor from univariate square free polynomial. 81 82 The RealRootIsolator provides isolating intervals for the real 83 roots of the polynomial. 84 \pre the polynomial is square free 85 */ 86 Descartes(const Polynomial& P = Polynomial(Coeff(0)), 87 bool is_strong = false, 88 int kk = 2) poly_(P)89 : poly_(P) , 90 is_strong_(is_strong), 91 k(kk), 92 interval_given(false) { 93 94 numerator = new IT[CGAL::degree(P)]; 95 denominator_exponent = new IT[CGAL::degree(P)]; 96 is_exact = new bool[CGAL::degree(P)]; 97 number_of_real_roots_ = 0; 98 if(CGAL::degree(P) == 0) 99 { 100 if(P.is_zero()) number_of_real_roots_ = -1; 101 return; 102 } 103 104 intern_decompose(poly_,typename FT_poly::Is_fraction()); 105 } 106 107 // constructor for coefficient types \c Coeff with given interval 108 // (experimental) 109 Descartes(const Polynomial& P, 110 const Rational& left, 111 const Rational& right, 112 bool is_strong = false, 113 int kk = 2) poly_(P)114 : poly_(P) , 115 is_strong_(is_strong), 116 k(kk), 117 interval_given(true) { 118 119 numerator = new IT[CGAL::degree(P)]; 120 denominator_exponent = new IT[CGAL::degree(P)]; 121 is_exact = new bool[CGAL::degree(P)]; 122 number_of_real_roots_ = 0; 123 if(CGAL::degree(P) == 0) 124 { 125 if(P.is_zero()) number_of_real_roots_ = -1; 126 return; 127 } 128 typename FT_rat::Decompose decompose; 129 typedef typename FT_rat::Numerator Numerator; 130 typedef typename FT_rat::Denominator Denominator; 131 Numerator numleft, numright; 132 Denominator denleft, denright; 133 134 decompose(left,numleft,denleft); 135 decompose(right,numright,denright); 136 137 LEFT = numleft * denright; 138 SCALE = numright * denleft - LEFT; 139 DENOM = denleft * denright; 140 poly_.scale_down(denleft*denright); 141 142 intern_decompose(poly_,typename FT_poly::Is_decomposable()); 143 } 144 145 //! copy constructor Descartes(const Descartes & D)146 Descartes(const Descartes& D) 147 : poly_(D.poly_), 148 number_of_real_roots_(D.number_of_real_roots_), 149 LEFT(D.LEFT), 150 SCALE(D.SCALE), 151 DENOM(D.DENOM), 152 is_strong_(D.is_strong_), 153 k(D.k), 154 interval_given(D.interval_given) { 155 156 numerator = new IT[CGAL::degree(poly_)]; 157 denominator_exponent = new IT[CGAL::degree(poly_)]; 158 is_exact = new bool[CGAL::degree(poly_)]; 159 for(int i=0; i<number_of_real_roots(); i++) 160 { 161 numerator[i] = D.numerator[i]; 162 denominator_exponent[i] = D.denominator_exponent[i]; 163 is_exact[i] = D.is_exact[i]; 164 } 165 } 166 167 // destructor ~Descartes()168 ~Descartes() { 169 delete[] numerator; 170 delete[] denominator_exponent; 171 delete[] is_exact; 172 } 173 174 public: // functions 175 176 /*! \brief returns the defining polynomial*/ polynomial()177 Polynomial polynomial() const { return poly_; } 178 179 //! returns the number of real roots number_of_real_roots()180 int number_of_real_roots() const { return number_of_real_roots_; } 181 182 /*! \brief returns true if the isolating interval is degenerated to a 183 single point. 184 185 If is_exact_root(i) is true, 186 then left_bound(int i) equals \f$root_i\f$. \n 187 If is_exact_root(i) is true, 188 then right_bound(int i) equals \f$root_i\f$. \n 189 */ is_exact_root(int i)190 bool is_exact_root(int i) const { return is_exact[i]; } 191 192 193 194 public: 195 196 left_bound(int i,IT & numerator_,IT & denominator_)197 void left_bound(int i, IT& numerator_, IT& denominator_) const { 198 CGAL_assertion(i >= 0 && i < number_of_real_roots_); 199 construct_binary(denominator_exponent[i], denominator_); 200 numerator_= SCALE * numerator[i] + LEFT * denominator_; 201 denominator_ = denominator_ * DENOM; 202 } 203 204 right_bound(int i,IT & numerator_,IT & denominator_)205 void right_bound(int i,IT& numerator_, IT& denominator_) const { 206 CGAL_assertion(i >= 0 && i < number_of_real_roots_); 207 if(is_exact[i]){ 208 return left_bound(i,numerator_,denominator_); 209 } 210 else{ 211 construct_binary(denominator_exponent[i],denominator_); 212 numerator_= SCALE * (numerator[i]+1) + LEFT * denominator_; 213 denominator_ = denominator_ * DENOM; 214 } 215 } 216 public: 217 218 /*! \brief returns \f${l_i}\f$ the left bound of the isolating interval 219 for root \f$root_{i}\f$. 220 221 In case is_exact_root(i) is true, \f$l_i = root_{i}\f$,\n 222 otherwise: \f$l_i < root_{i}\f$. 223 224 If \f$i-1>=0\f$, then \f$l_i > root_{i-1}\f$. \n 225 If \f$i-1>=0\f$, then \f$l_i >= r_{i-1}\f$, 226 the right bound of \f$root_{i-1}\f$\n 227 228 \pre 0 <= i < number_of_real_roots() 229 */ left_bound(int i)230 Rational left_bound(int i) const { 231 IT numerator_, denominator_; 232 left_bound(i,numerator_,denominator_); 233 return Rational(numerator_) / Rational(denominator_); 234 } 235 236 /*! \brief returns \f${r_i}\f$ the right bound of the isolating interval 237 for root \f$root_{i}\f$. 238 239 In case is_exact_root(i) is true, \f$r_i = root_{i}\f$,\n 240 otherwise: \f$r_i > root_{i}\f$. 241 242 If \f$i+1< n \f$, then \f$r_i < root_{i+1}\f$, 243 where \f$n\f$ is number of real roots.\n 244 If \f$i+1< n \f$, then \f$r_i <= l_{i+1}\f$, 245 the left bound of \f$root_{i+1}\f$\n 246 247 248 \pre 0 <= i < number_of_real_roots() 249 */ right_bound(int i)250 Rational right_bound(int i) const { 251 IT numerator_, denominator_; 252 right_bound(i,numerator_,denominator_); 253 return Rational(numerator_) / Rational(denominator_); 254 } 255 256 private: intern_decompose(Polynomial P_,::CGAL::Tag_true)257 void intern_decompose( Polynomial P_, ::CGAL::Tag_true){ 258 typename FT_poly::Decompose decompose; 259 typename FT_poly::Numerator_type NumP; 260 typename FT_poly::Denominator_type dummy; 261 262 decompose(P_,NumP,dummy); 263 init_with(NumP); 264 } 265 intern_decompose(Polynomial P,::CGAL::Tag_false)266 void intern_decompose( Polynomial P, ::CGAL::Tag_false){ 267 init_with(P); 268 } 269 270 271 template<class Polynomial__> init_with(const Polynomial__ & P)272 void init_with(const Polynomial__& P){ 273 typedef typename Polynomial__::NT Coeff; 274 if(!interval_given) 275 { 276 LEFT = -weak_upper_root_bound<Coeff>(P); 277 SCALE = - LEFT * IT(2); 278 DENOM = IT(1); 279 } 280 Polynomial__ R = ::CGAL::translate(P,Coeff(LEFT)); 281 Polynomial__ Q = ::CGAL::scale_up(R,Coeff(SCALE)); 282 zero_one_descartes<Coeff>(Q,0,0); 283 } 284 285 286 //! returns the polynomial $(1 + x)^n P(1/(1 + x))$. 287 template <class Coeff__> 288 /* 289 typename 290 CGAL::Polynomial_traits_d<Polynomial> 291 ::template Rebind<Coeff__,1>::Other::Type 292 */ POLYNOMIAL_REBIND(Coeff__)293 POLYNOMIAL_REBIND(Coeff__) 294 variation_transformation(const POLYNOMIAL_REBIND(Coeff__)& P) { 295 POLYNOMIAL_REBIND(Coeff__) R = reversal(P); 296 return translate_by_one(R); 297 } 298 299 //! Returns an upper bound on the absolute value of all roots of $P$. 300 /*! The upper bound is a power of two. Only works for univariate 301 * polynomials. 302 */ 303 template <class Coeff__> weak_upper_root_bound(const POLYNOMIAL_REBIND (Coeff__)& P)304 IT weak_upper_root_bound(const POLYNOMIAL_REBIND(Coeff__)& P) { 305 306 typename Real_embeddable_traits<Coeff__>::Abs abs; 307 const int n = CGAL::degree(P); 308 IT r(1); // return value 309 Coeff__ x(1); // needed to "evaluate" the polynomial 310 Coeff__ val; 311 for (;;) { 312 val = -abs(P[n]); 313 for (int i = n-1; i >= 0; i--) { 314 val = val*x + abs(P[i]); 315 } 316 if (val < Coeff__(0)) return r; 317 r *= IT(2); 318 x = Coeff__(r); 319 } 320 } 321 322 //! tests if the polynomial has no root in the interval. 323 template <class Coeff__> not_zero_in_interval(const POLYNOMIAL_REBIND (Coeff__)& P)324 bool not_zero_in_interval(const POLYNOMIAL_REBIND(Coeff__)& P) 325 { 326 if(CGAL::degree(P) == 0) return true; 327 if(internal::sign_variations(variation_transformation<Coeff__>(P)) != 0) 328 return false; 329 return (P[0] != Coeff__(0) && P.evaluate(Coeff__(1)) != Coeff__(0)); 330 } 331 //! Descartes algoritm to determine isolating intervals for the roots 332 //! lying in the interval (0,1). 333 // The parameters $(i,D)$ describe the interval $(i/2^D, (i+1)/2^D)$. 334 // Here $0\leq i < 2^D$. 335 template <class Coeff__> zero_one_descartes(const POLYNOMIAL_REBIND (Coeff__)& P,IT i,IT D)336 void zero_one_descartes(const POLYNOMIAL_REBIND(Coeff__)& P, 337 IT i, IT D) { 338 // Determine the number of sign variations of the transformed 339 // polynomial $(1+x)^nP(1/(1+x))$. This gives the number of 340 // roots of $P$ in $(0,1)$. 341 342 POLYNOMIAL_REBIND(Coeff__) R = variation_transformation<Coeff__>(P); 343 int descarte = sign_variations(R); 344 345 // no root 346 if ( descarte == 0 ) return; 347 348 // exactly one root 349 // Note the termination criterion $P(0)\neq 0$ and $P(1)\neq 0$. 350 // This ensures that the given interval is an isolating interval. 351 if ( descarte == 1 352 && P[0] != Coeff__(0) 353 && P.evaluate(Coeff__(1)) != Coeff__(0) ) { 354 if(is_strong_) { 355 strong_zero_one_descartes<Coeff__>(P,i,D); 356 return; 357 } 358 else { 359 numerator[number_of_real_roots_] = i; 360 denominator_exponent[number_of_real_roots_] = D; 361 is_exact[number_of_real_roots_] = false; 362 number_of_real_roots_++; 363 return; 364 } 365 } 366 367 // more than one root 368 // Refine the interval. 369 i = 2*i; D = D+1; 370 371 // Transform the polynomial such that the first half of the interval 372 // is mapped to the unit interval. 373 POLYNOMIAL_REBIND(Coeff__) Q = scale_down(P,Coeff__(2)); 374 375 // Consider the first half of the interval. 376 zero_one_descartes<Coeff__>(Q,i,D); 377 378 // Test if the polynomial is zero at the midpoint of the interval 379 POLYNOMIAL_REBIND(Coeff__) S = translate_by_one(Q); 380 if ( S[0] == Coeff__(0) ) { 381 numerator[number_of_real_roots_] = i + 1; 382 denominator_exponent[number_of_real_roots_] = D; 383 is_exact[number_of_real_roots_] = true; 384 number_of_real_roots_++; 385 } 386 387 // Consider the second half of the interval. 388 zero_one_descartes<Coeff__>(S,i+1,D); 389 } 390 391 392 //! Strong Descartes algoritm to determine isolating intervals for the 393 //! roots lying in the interval (0,1), where the first 394 //! derivative have no sign change. \pre $P$ has only one root in the 395 //! interval given by $(i,D)$. 396 // The parameters $(i,D)$ describe the interval $(i/2^D, (i+1)/2^D)$. 397 // Here $0\leq i < D$. 398 template <class Coeff__> strong_zero_one_descartes(const POLYNOMIAL_REBIND (Coeff__)& P,IT i,IT D)399 void strong_zero_one_descartes(const POLYNOMIAL_REBIND(Coeff__)& P, 400 IT i, IT D) { 401 402 // Test if the polynomial P' has no roots in the 403 // interval. For further use in Newton, the interval should be not 404 // too large. 405 406 // test if isolating interval is smaller than epsilon 407 // [l,r] -> r-l < epsilon 408 // l = (r-l) * i/2^D + l 409 // r = (r-l) * (i+1)/2^D + l 410 // r-l = (r-l) * 1/2^D 411 // r-l < epsilon = 2^(-k) 412 // <=> (r-l) * 1/2^D < 2^(-k) 413 // <=> 2^D > (r-l) / 2^(-k) 414 // <=> 2^D > (r-l) * 2^k 415 416 POLYNOMIAL_REBIND(Coeff__) PP = CGAL::differentiate(P); 417 if(not_zero_in_interval<Coeff__>(PP)) { // P' 418 IT tmp; 419 construct_binary(D-k, tmp); // tmp = 2^{D-k} 420 if(tmp * DENOM > SCALE ) { 421 numerator[number_of_real_roots_] = i; 422 denominator_exponent[number_of_real_roots_] = D; 423 is_exact[number_of_real_roots_] = false; 424 number_of_real_roots_++; 425 return; 426 } 427 } 428 429 // either $P'$ fails the test, 430 // or the interval is too large 431 // Refine the interval. 432 i = 2*i; D = D+1; 433 434 // Transform the polynomial such that the first half of the interval 435 // is mapped to the unit interval. 436 POLYNOMIAL_REBIND(Coeff__) Q = scale_down(P,Coeff__(2)); 437 438 // Test if the polynomial is zero at the midpoint of the interval 439 POLYNOMIAL_REBIND(Coeff__) S = translate_by_one(Q); 440 if ( S[0] == Coeff__(0) ) { 441 numerator[number_of_real_roots_] = i + 1; 442 denominator_exponent[number_of_real_roots_] = D; 443 is_exact[number_of_real_roots_] = true; 444 number_of_real_roots_++; 445 return; 446 } 447 448 // Consider the first half of the interval. 449 if(sign_variations(variation_transformation<Coeff__>(Q)) == 1) { 450 strong_zero_one_descartes<Coeff__>(Q,i,D); 451 return; 452 } 453 454 // Consider the second half of the interval. 455 strong_zero_one_descartes<Coeff__>(S,i+1,D); 456 return; 457 } 458 }; 459 460 } // namespace internal 461 462 } //namespace CGAL 463 464 #endif // CGAL_ALGEBRAIC_KERNEL_D_DESCARTES_H 465