1 // Copyright (c) 2008 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/Polynomial/include/CGAL/Polynomial_traits_d.h $ 7 // $Id: Polynomial_traits_d.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@informatik.uni-mainz.de> 12 // Sebastian Limbach <slimbach@mpi-inf.mpg.de> 13 // 14 // ============================================================================ 15 #ifndef CGAL_POLYNOMIAL_TRAITS_D_H 16 #define CGAL_POLYNOMIAL_TRAITS_D_H 17 18 #include <CGAL/disable_warnings.h> 19 20 #include <CGAL/basic.h> 21 #include <functional> 22 #include <list> 23 #include <vector> 24 #include <utility> 25 26 #include <CGAL/Polynomial/fwd.h> 27 #include <CGAL/Polynomial/misc.h> 28 #include <CGAL/Polynomial/Polynomial_type.h> 29 #include <CGAL/Polynomial/Monomial_representation.h> 30 #include <CGAL/Polynomial/Degree.h> 31 #include <CGAL/polynomial_utils.h> 32 #include <CGAL/Polynomial/square_free_factorize.h> 33 #include <CGAL/Polynomial/modular_filter.h> 34 #include <CGAL/extended_euclidean_algorithm.h> 35 #include <CGAL/Polynomial/resultant.h> 36 #include <CGAL/Polynomial/subresultants.h> 37 #include <CGAL/Polynomial/sturm_habicht_sequence.h> 38 39 #include <CGAL/boost/iterator/transform_iterator.hpp> 40 #include <CGAL/tss.h> 41 42 43 #define CGAL_POLYNOMIAL_TRAITS_D_BASE_TYPEDEFS \ 44 private: \ 45 typedef Polynomial_traits_d< Polynomial< Coefficient_type_ > > PT; \ 46 typedef Polynomial_traits_d< Coefficient_type_ > PTC; \ 47 \ 48 typedef Polynomial<Coefficient_type_> Polynomial_d; \ 49 typedef Coefficient_type_ Coefficient_type; \ 50 \ 51 typedef typename Innermost_coefficient_type<Polynomial_d>::Type \ 52 Innermost_coefficient_type; \ 53 static const int d = Dimension<Polynomial_d>::value; \ 54 \ 55 \ 56 typedef std::pair< Exponent_vector, Innermost_coefficient_type > \ 57 Exponents_coeff_pair; \ 58 typedef std::vector< Exponents_coeff_pair > Monom_rep; \ 59 \ 60 typedef CGAL::Recursive_const_flattening< d-1, \ 61 typename CGAL::Polynomial<Coefficient_type>::const_iterator > \ 62 Coefficient_const_flattening; \ 63 \ 64 typedef typename \ 65 Coefficient_const_flattening::Recursive_flattening_iterator \ 66 Innermost_coefficient_const_iterator; \ 67 \ 68 typedef typename Polynomial_d::const_iterator \ 69 Coefficient_const_iterator; \ 70 \ 71 typedef std::pair<Innermost_coefficient_const_iterator, \ 72 Innermost_coefficient_const_iterator> \ 73 Innermost_coefficient_const_iterator_range; \ 74 \ 75 typedef std::pair<Coefficient_const_iterator, \ 76 Coefficient_const_iterator> \ 77 Coefficient_const_iterator_range; \ 78 79 80 81 namespace CGAL { 82 83 namespace internal { 84 85 // Base class for functors depending on the algebraic category of the 86 // innermost coefficient 87 template< class Coefficient_type_, class ICoeffAlgebraicCategory > 88 class Polynomial_traits_d_base_icoeff_algebraic_category { 89 public: 90 typedef Null_functor Multivariate_content; 91 }; 92 93 // Specializations 94 template< class Coefficient_type_ > 95 class Polynomial_traits_d_base_icoeff_algebraic_category< 96 Polynomial< Coefficient_type_ >, Integral_domain_without_division_tag > 97 : public Polynomial_traits_d_base_icoeff_algebraic_category< 98 Polynomial< Coefficient_type_ >, Null_tag > {}; 99 100 template< class Coefficient_type_ > 101 class Polynomial_traits_d_base_icoeff_algebraic_category< 102 Polynomial< Coefficient_type_ >, Integral_domain_tag > 103 : public Polynomial_traits_d_base_icoeff_algebraic_category< 104 Polynomial< Coefficient_type_ >, Integral_domain_without_division_tag > {}; 105 106 template< class Coefficient_type_ > 107 class Polynomial_traits_d_base_icoeff_algebraic_category< 108 Polynomial< Coefficient_type_ >, Unique_factorization_domain_tag > 109 : public Polynomial_traits_d_base_icoeff_algebraic_category< 110 Polynomial< Coefficient_type_ >, Integral_domain_tag > { 111 CGAL_POLYNOMIAL_TRAITS_D_BASE_TYPEDEFS 112 113 public: 114 115 116 struct Multivariate_content 117 : public CGAL::cpp98::unary_function< Polynomial_d , Innermost_coefficient_type >{ 118 Innermost_coefficient_type operatorMultivariate_content119 operator()(const Polynomial_d& p) const { 120 typedef Innermost_coefficient_const_iterator IT; 121 Innermost_coefficient_type content(0); 122 typename PT::Construct_innermost_coefficient_const_iterator_range range; 123 for (IT it = range(p).first; it != range(p).second; it++){ 124 content = CGAL::gcd(content, *it); 125 if(CGAL::is_one(content)) break; 126 } 127 return content; 128 } 129 }; 130 }; 131 132 template< class Coefficient_type_ > 133 class Polynomial_traits_d_base_icoeff_algebraic_category< 134 Polynomial< Coefficient_type_ >, Euclidean_ring_tag > 135 : public Polynomial_traits_d_base_icoeff_algebraic_category< 136 Polynomial< Coefficient_type_ >, Unique_factorization_domain_tag > 137 {}; 138 139 template< class Coefficient_type_ > 140 class Polynomial_traits_d_base_icoeff_algebraic_category< 141 Polynomial< Coefficient_type_ >, Field_tag > 142 : public Polynomial_traits_d_base_icoeff_algebraic_category< 143 Polynomial< Coefficient_type_ >, Integral_domain_tag > { 144 CGAL_POLYNOMIAL_TRAITS_D_BASE_TYPEDEFS 145 146 public: 147 148 // Multivariate_content; 149 struct Multivariate_content 150 : public CGAL::cpp98::unary_function< Polynomial_d , Innermost_coefficient_type >{ operatorMultivariate_content151 Innermost_coefficient_type operator()(const Polynomial_d& p) const { 152 if( CGAL::is_zero(p) ) 153 return Innermost_coefficient_type(0); 154 else 155 return Innermost_coefficient_type(1); 156 } 157 }; 158 }; 159 160 template< class Coefficient_type_ > 161 class Polynomial_traits_d_base_icoeff_algebraic_category< 162 Polynomial< Coefficient_type_ >, Field_with_sqrt_tag > 163 : public Polynomial_traits_d_base_icoeff_algebraic_category< 164 Polynomial< Coefficient_type_ >, Field_tag > {}; 165 166 template< class Coefficient_type_ > 167 class Polynomial_traits_d_base_icoeff_algebraic_category< 168 Polynomial< Coefficient_type_ >, Field_with_kth_root_tag > 169 : public Polynomial_traits_d_base_icoeff_algebraic_category< 170 Polynomial< Coefficient_type_ >, Field_with_sqrt_tag > {}; 171 172 template< class Coefficient_type_ > 173 class Polynomial_traits_d_base_icoeff_algebraic_category< 174 Polynomial< Coefficient_type_ >, Field_with_root_of_tag > 175 : public Polynomial_traits_d_base_icoeff_algebraic_category< 176 Polynomial< Coefficient_type_ >, Field_with_kth_root_tag > {}; 177 178 // Base class for functors depending on the algebraic category of the 179 // Polynomial type 180 template< class Coefficient_type_, class PolynomialAlgebraicCategory > 181 class Polynomial_traits_d_base_polynomial_algebraic_category { 182 public: 183 typedef Null_functor Univariate_content; 184 typedef Null_functor Square_free_factorize; 185 }; 186 187 // Specializations 188 template< class Coefficient_type_ > 189 class Polynomial_traits_d_base_polynomial_algebraic_category< 190 Polynomial< Coefficient_type_ >, Integral_domain_without_division_tag > 191 : public Polynomial_traits_d_base_polynomial_algebraic_category< 192 Polynomial< Coefficient_type_ >, Null_tag > {}; 193 194 template< class Coefficient_type_ > 195 class Polynomial_traits_d_base_polynomial_algebraic_category< 196 Polynomial< Coefficient_type_ >, Integral_domain_tag > 197 : public Polynomial_traits_d_base_polynomial_algebraic_category< 198 Polynomial< Coefficient_type_ >, Integral_domain_without_division_tag > {}; 199 200 template< class Coefficient_type_ > 201 class Polynomial_traits_d_base_polynomial_algebraic_category< 202 Polynomial< Coefficient_type_ >, Unique_factorization_domain_tag > 203 : public Polynomial_traits_d_base_polynomial_algebraic_category< 204 Polynomial< Coefficient_type_ >, Integral_domain_tag > { 205 CGAL_POLYNOMIAL_TRAITS_D_BASE_TYPEDEFS 206 207 public: 208 209 // Univariate_content 210 struct Univariate_content 211 : public CGAL::cpp98::unary_function< Polynomial_d , Coefficient_type>{ operatorUnivariate_content212 Coefficient_type operator()(const Polynomial_d& p) const { 213 return p.content(); 214 } 215 }; 216 217 // Square_free_factorize; 218 struct Square_free_factorize{ 219 220 template < class OutputIterator > operatorSquare_free_factorize221 OutputIterator operator()( const Polynomial_d& p, OutputIterator oi) const { 222 std::vector<Polynomial_d> factors; 223 std::vector<int> mults; 224 225 square_free_factorize 226 ( p, std::back_inserter(factors), std::back_inserter(mults) ); 227 228 CGAL_postcondition( factors.size() == mults.size() ); 229 for(unsigned int i = 0; i < factors.size(); i++){ 230 *oi++=std::make_pair(factors[i],mults[i]); 231 } 232 233 return oi; 234 } 235 236 template< class OutputIterator > operatorSquare_free_factorize237 OutputIterator operator()( 238 const Polynomial_d& p , 239 OutputIterator oi, 240 Innermost_coefficient_type& a ) const { 241 242 if( CGAL::is_zero(p) ) { 243 a = Innermost_coefficient_type(0); 244 return oi; 245 } 246 247 typedef Polynomial_traits_d< Polynomial_d > PT; 248 typename PT::Innermost_leading_coefficient ilcoeff; 249 typename PT::Multivariate_content mcontent; 250 a = CGAL::unit_part( ilcoeff( p ) ) * mcontent( p ); 251 252 return (*this)( p/Polynomial_d(a), oi); 253 } 254 }; 255 }; 256 257 template< class Coefficient_type_ > 258 class Polynomial_traits_d_base_polynomial_algebraic_category< 259 Polynomial< Coefficient_type_ >, Euclidean_ring_tag > 260 : public Polynomial_traits_d_base_polynomial_algebraic_category< 261 Polynomial< Coefficient_type_ >, Unique_factorization_domain_tag > {}; 262 263 264 // Polynomial_traits_d_base class connecting the two base classes which depend 265 // on the algebraic category of the innermost coefficient type and the poly- 266 // nomial type. 267 268 // First the general base class for the innermost coefficient 269 template< class InnermostCoefficient_type, 270 class ICoeffAlgebraicCategory, class PolynomialAlgebraicCategory > 271 class Polynomial_traits_d_base { 272 typedef InnermostCoefficient_type ICoeff; 273 public: 274 static const int d = 0; 275 276 typedef ICoeff Polynomial_d; 277 typedef ICoeff Coefficient_type; 278 typedef ICoeff Innermost_coefficient_type; 279 280 struct Degree 281 : public CGAL::cpp98::unary_function< ICoeff , int > { operatorDegree282 int operator()(const ICoeff&) const { return 0; } 283 }; 284 struct Total_degree 285 : public CGAL::cpp98::unary_function< ICoeff , int > { operatorTotal_degree286 int operator()(const ICoeff&) const { return 0; } 287 }; 288 289 typedef Null_functor Construct_polynomial; 290 typedef Null_functor Get_coefficient; 291 typedef Null_functor Leading_coefficient; 292 typedef Null_functor Univariate_content; 293 typedef Null_functor Multivariate_content; 294 typedef Null_functor Shift; 295 typedef Null_functor Negate; 296 typedef Null_functor Invert; 297 typedef Null_functor Translate; 298 typedef Null_functor Translate_homogeneous; 299 typedef Null_functor Scale_homogeneous; 300 typedef Null_functor Differentiate; 301 302 struct Is_square_free 303 : public CGAL::cpp98::unary_function< ICoeff, bool > { operatorIs_square_free304 bool operator()( const ICoeff& ) const { 305 return true; 306 } 307 }; 308 309 struct Make_square_free 310 : public CGAL::cpp98::unary_function< ICoeff, ICoeff>{ operatorMake_square_free311 ICoeff operator()( const ICoeff& x ) const { 312 if (CGAL::is_zero(x)) return x ; 313 else return ICoeff(1); 314 } 315 }; 316 317 typedef Null_functor Square_free_factorize; 318 typedef Null_functor Pseudo_division; 319 typedef Null_functor Pseudo_division_remainder; 320 typedef Null_functor Pseudo_division_quotient; 321 322 struct Gcd_up_to_constant_factor 323 : public CGAL::cpp98::binary_function< ICoeff, ICoeff, ICoeff >{ operatorGcd_up_to_constant_factor324 ICoeff operator()(const ICoeff& x, const ICoeff& y) const { 325 if (CGAL::is_zero(x) && CGAL::is_zero(y)) 326 return ICoeff(0); 327 else 328 return ICoeff(1); 329 } 330 }; 331 332 typedef Null_functor Integral_division_up_to_constant_factor; 333 334 struct Univariate_content_up_to_constant_factor 335 : public CGAL::cpp98::unary_function< ICoeff, ICoeff >{ operatorUnivariate_content_up_to_constant_factor336 ICoeff operator()(const ICoeff& ) const { 337 // TODO: Why not return 0 if argument is 0 ? 338 return ICoeff(1); 339 } 340 }; 341 342 typedef Null_functor Square_free_factorize_up_to_constant_factor; 343 typedef Null_functor Resultant; 344 typedef Null_functor Canonicalize; 345 typedef Null_functor Evaluate_homogeneous; 346 347 struct Innermost_leading_coefficient 348 :public CGAL::cpp98::unary_function <ICoeff, ICoeff>{ operatorInnermost_leading_coefficient349 const ICoeff& operator()(const ICoeff& x){return x;} 350 }; 351 352 struct Degree_vector{ 353 typedef Exponent_vector result_type; 354 typedef Coefficient_type argument_type; 355 // returns the exponent vector of inner_most_lcoeff. operatorDegree_vector356 result_type operator()(const Coefficient_type&) const{ 357 return Exponent_vector(); 358 } 359 }; 360 361 struct Get_innermost_coefficient 362 : public CGAL::cpp98::binary_function< ICoeff, Polynomial_d, Exponent_vector > { operatorGet_innermost_coefficient363 const ICoeff& operator()( const Polynomial_d& p, Exponent_vector ) { 364 return p; 365 } 366 }; 367 368 typedef Null_functor Evaluate ; 369 370 struct Substitute{ 371 public: 372 template <class Input_iterator> 373 typename 374 CGAL::Coercion_traits< 375 typename std::iterator_traits<Input_iterator>::value_type, 376 Innermost_coefficient_type>::Type operatorSubstitute377 operator()( 378 const Innermost_coefficient_type& p, 379 Input_iterator CGAL_precondition_code(begin), 380 Input_iterator CGAL_precondition_code(end) ) const { 381 CGAL_precondition(end == begin); 382 typedef typename std::iterator_traits<Input_iterator>::value_type 383 value_type; 384 typedef CGAL::Coercion_traits<Innermost_coefficient_type,value_type> CT; 385 return typename CT::Cast()(p); 386 } 387 }; 388 389 struct Substitute_homogeneous{ 390 public: 391 // this is the end of the recursion 392 // begin contains the homogeneous variabel 393 // hdegree is the remaining degree 394 template <class Input_iterator> 395 typename 396 CGAL::Coercion_traits< 397 typename std::iterator_traits<Input_iterator>::value_type, 398 Innermost_coefficient_type>::Type operatorSubstitute_homogeneous399 operator()( 400 const Innermost_coefficient_type& p, 401 Input_iterator begin, 402 Input_iterator CGAL_precondition_code(end), 403 int hdegree) const { 404 405 typedef typename std::iterator_traits<Input_iterator>::value_type 406 value_type; 407 typedef CGAL::Coercion_traits<Innermost_coefficient_type,value_type> CT; 408 typename CT::Type result = 409 typename CT::Cast()(CGAL::ipower(*begin++,hdegree)) 410 * typename CT::Cast()(p); 411 412 CGAL_precondition(end == begin); 413 CGAL_precondition(hdegree >= 0); 414 return result; 415 } 416 }; 417 }; 418 419 // Now the version for the polynomials with all functors provided by all 420 // polynomials 421 template< class Coefficient_type_, 422 class ICoeffAlgebraicCategory, class PolynomialAlgebraicCategory > 423 class Polynomial_traits_d_base< Polynomial< Coefficient_type_ >, 424 ICoeffAlgebraicCategory, PolynomialAlgebraicCategory > 425 : public Polynomial_traits_d_base_icoeff_algebraic_category< 426 Polynomial< Coefficient_type_ >, ICoeffAlgebraicCategory >, 427 public Polynomial_traits_d_base_polynomial_algebraic_category< 428 Polynomial< Coefficient_type_ >, PolynomialAlgebraicCategory > { 429 430 typedef Polynomial_traits_d< Polynomial< Coefficient_type_ > > PT; 431 typedef Polynomial_traits_d< Coefficient_type_ > PTC; 432 433 public: 434 typedef Polynomial<Coefficient_type_> Polynomial_d; 435 typedef Coefficient_type_ Coefficient_type; 436 437 typedef typename internal::Innermost_coefficient_type<Polynomial_d>::Type 438 Innermost_coefficient_type; 439 static const int d = Dimension<Polynomial_d>::value; 440 441 private: 442 typedef std::pair< Exponent_vector, Innermost_coefficient_type > 443 Exponents_coeff_pair; 444 typedef std::vector< Exponents_coeff_pair > Monom_rep; 445 446 typedef CGAL::Recursive_const_flattening< d-1, 447 typename CGAL::Polynomial<Coefficient_type>::const_iterator > 448 Coefficient_const_flattening; 449 450 public: 451 typedef typename Coefficient_const_flattening::Recursive_flattening_iterator 452 Innermost_coefficient_const_iterator; 453 typedef typename Polynomial_d::const_iterator Coefficient_const_iterator; 454 455 typedef std::pair<Innermost_coefficient_const_iterator, 456 Innermost_coefficient_const_iterator> 457 Innermost_coefficient_const_iterator_range; 458 459 typedef std::pair<Coefficient_const_iterator, 460 Coefficient_const_iterator> 461 Coefficient_const_iterator_range; 462 463 464 // We use our own Strict Weak Ordering predicate in order to avoid 465 // problems when calling sort for a Exponents_coeff_pair where the 466 // coeff type has no comparison operators available. 467 private: 468 struct Compare_exponents_coeff_pair 469 : public CGAL::cpp98::binary_function< 470 std::pair< Exponent_vector, Innermost_coefficient_type >, 471 std::pair< Exponent_vector, Innermost_coefficient_type >, 472 bool > 473 { operatorCompare_exponents_coeff_pair474 bool operator()( 475 const std::pair< Exponent_vector, Innermost_coefficient_type >& p1, 476 const std::pair< Exponent_vector, Innermost_coefficient_type >& p2 ) const { 477 // TODO: Precondition leads to an error within test_translate in 478 // Polynomial_traits_d test 479 // CGAL_precondition( p1.first != p2.first ); 480 return p1.first < p2.first; 481 } 482 }; 483 484 public: 485 486 // 487 // Functors as defined in the reference manual 488 // (with sometimes slightly extended functionality) 489 490 // Construct_polynomial; 491 struct Construct_polynomial { 492 493 typedef Polynomial_d result_type; 494 operatorConstruct_polynomial495 Polynomial_d operator()() const { 496 return Polynomial_d(0); 497 } 498 499 template <class T> operatorConstruct_polynomial500 Polynomial_d operator()( T a ) const { 501 return Polynomial_d(a); 502 } 503 504 //! construct the constant polynomial a0 operatorConstruct_polynomial505 Polynomial_d operator() (const Coefficient_type& a0) const 506 {return Polynomial_d(a0);} 507 508 //! construct the polynomial a0 + a1*x operatorConstruct_polynomial509 Polynomial_d operator() ( 510 const Coefficient_type& a0, const Coefficient_type& a1) const 511 {return Polynomial_d(a0,a1);} 512 513 //! construct the polynomial a0 + a1*x + a2*x^2 operatorConstruct_polynomial514 Polynomial_d operator() ( 515 const Coefficient_type& a0, const Coefficient_type& a1, 516 const Coefficient_type& a2) const 517 {return Polynomial_d(a0,a1,a2);} 518 519 //! construct the polynomial a0 + a1*x + ... + a3*x^3 operatorConstruct_polynomial520 Polynomial_d operator() ( 521 const Coefficient_type& a0, const Coefficient_type& a1, 522 const Coefficient_type& a2, const Coefficient_type& a3) const 523 {return Polynomial_d(a0,a1,a2,a3);} 524 525 //! construct the polynomial a0 + a1*x + ... + a4*x^4 operatorConstruct_polynomial526 Polynomial_d operator() ( 527 const Coefficient_type& a0, const Coefficient_type& a1, 528 const Coefficient_type& a2, const Coefficient_type& a3, 529 const Coefficient_type& a4) const 530 {return Polynomial_d(a0,a1,a2,a3,a4);} 531 532 //! construct the polynomial a0 + a1*x + ... + a5*x^5 operatorConstruct_polynomial533 Polynomial_d operator() ( 534 const Coefficient_type& a0, const Coefficient_type& a1, 535 const Coefficient_type& a2, const Coefficient_type& a3, 536 const Coefficient_type& a4, const Coefficient_type& a5) const 537 {return Polynomial_d(a0,a1,a2,a3,a4,a5);} 538 539 //! construct the polynomial a0 + a1*x + ... + a6*x^6 operatorConstruct_polynomial540 Polynomial_d operator() ( 541 const Coefficient_type& a0, const Coefficient_type& a1, 542 const Coefficient_type& a2, const Coefficient_type& a3, 543 const Coefficient_type& a4, const Coefficient_type& a5, 544 const Coefficient_type& a6) const 545 {return Polynomial_d(a0,a1,a2,a3,a4,a5,a6);} 546 547 //! construct the polynomial a0 + a1*x + ... + a7*x^7 operatorConstruct_polynomial548 Polynomial_d operator() ( 549 const Coefficient_type& a0, const Coefficient_type& a1, 550 const Coefficient_type& a2, const Coefficient_type& a3, 551 const Coefficient_type& a4, const Coefficient_type& a5, 552 const Coefficient_type& a6, const Coefficient_type& a7) const 553 {return Polynomial_d(a0,a1,a2,a3,a4,a5,a6,a7);} 554 555 //! construct the polynomial a0 + a1*x + ... + a8*x^8 operatorConstruct_polynomial556 Polynomial_d operator() ( 557 const Coefficient_type& a0, const Coefficient_type& a1, 558 const Coefficient_type& a2, const Coefficient_type& a3, 559 const Coefficient_type& a4, const Coefficient_type& a5, 560 const Coefficient_type& a6, const Coefficient_type& a7, 561 const Coefficient_type& a8) const 562 {return Polynomial_d(a0,a1,a2,a3,a4,a5,a6,a7,a8);} 563 564 #if 1 565 private: 566 template <class Input_iterator, class NT> Polynomial_d construct_value_typeConstruct_polynomial567 construct_value_type(Input_iterator begin, Input_iterator end, NT) const { 568 typedef CGAL::Coercion_traits<NT,Coefficient_type> CT; 569 CGAL_static_assertion((boost::is_same<typename CT::Type,Coefficient_type>::value)); 570 typename CT::Cast cast; 571 return Polynomial_d( 572 boost::make_transform_iterator(begin,cast), 573 boost::make_transform_iterator(end,cast)); 574 } 575 576 template <class Input_iterator, class NT> Polynomial_d construct_value_typeConstruct_polynomial577 construct_value_type(Input_iterator begin, Input_iterator end, std::pair<Exponent_vector,NT>) const { 578 return (*this)(begin,end,false);// construct from non sorted monom rep 579 } 580 581 public: 582 template< class Input_iterator > operatorConstruct_polynomial583 Polynomial_d operator()( Input_iterator begin, Input_iterator end) const { 584 if(begin == end ) return Polynomial_d(0); 585 typedef typename std::iterator_traits<Input_iterator>::value_type value_type; 586 return construct_value_type(begin,end,value_type()); 587 } 588 589 template< class Input_iterator > operatorConstruct_polynomial590 Polynomial_d operator()( Input_iterator begin, Input_iterator end, bool is_sorted) const { 591 // Avoid compiler warning 592 (void)is_sorted; 593 if(begin == end ) return Polynomial_d(0); 594 Monom_rep monom_rep(begin,end); 595 // if(!is_sorted) 596 std::sort(monom_rep.begin(),monom_rep.end(),Compare_exponents_coeff_pair()); 597 return Create_polynomial_from_monom_rep<Coefficient_type>()(monom_rep.begin(),monom_rep.end()); 598 } 599 #else 600 601 // Construct from Coefficient type 602 template< class Input_iterator > 603 inline Polynomial_d constructConstruct_polynomial604 construct( Input_iterator begin, Input_iterator end, Tag_true) const { 605 if(begin == end ) return Polynomial_d(0); 606 return Polynomial_d(begin,end); 607 } 608 // Construct from momom rep 609 template< class Input_iterator > 610 inline Polynomial_d constructConstruct_polynomial611 construct( Input_iterator begin, Input_iterator end, Tag_false) const { 612 // construct from non sorted monom rep 613 return (*this)(begin,end,false); 614 } 615 616 template< class Input_iterator > 617 Polynomial_d operatorConstruct_polynomial618 operator()( Input_iterator begin, Input_iterator end ) const { 619 if(begin == end ) return Polynomial_d(0); 620 typedef typename std::iterator_traits<Input_iterator>::value_type value_type; 621 typedef Boolean_tag<boost::is_same<value_type,Coefficient_type>::value> 622 Is_coeff; 623 std::vector<value_type> vec(begin,end); 624 return construct(vec.begin(),vec.end(),Is_coeff()); 625 } 626 627 628 template< class Input_iterator > 629 Polynomial_d operatorConstruct_polynomial630 operator()(Input_iterator begin, Input_iterator end , bool is_sorted) const{ 631 if(!is_sorted) 632 std::sort(begin,end,Compare_exponents_coeff_pair()); 633 return Create_polynomial_from_monom_rep< Coefficient_type >()(begin,end); 634 } 635 #endif 636 637 public: 638 639 template< class T > 640 class Create_polynomial_from_monom_rep { 641 public: 642 template <class Monom_rep_iterator> operatorConstruct_polynomial643 Polynomial_d operator()( 644 Monom_rep_iterator begin, 645 Monom_rep_iterator end) const { 646 647 Innermost_coefficient_type zero(0); 648 std::vector< Innermost_coefficient_type > coefficients; 649 for(Monom_rep_iterator it = begin; it != end; it++){ 650 int current_exp = it->first[0]; 651 if((int) coefficients.size() < current_exp) 652 coefficients.resize(current_exp,zero); 653 coefficients.push_back(it->second); 654 } 655 return Polynomial_d(coefficients.begin(),coefficients.end()); 656 } 657 }; 658 659 template< class T > 660 class Create_polynomial_from_monom_rep< Polynomial < T > > { 661 public: 662 template <class Monom_rep_iterator> operatorConstruct_polynomial663 Polynomial_d operator()( 664 Monom_rep_iterator begin, 665 Monom_rep_iterator end) const { 666 667 typedef Polynomial_traits_d<Coefficient_type> PT; 668 typename PT::Construct_polynomial construct; 669 670 CGAL_static_assertion(PT::d != 0); // Coefficient_type is a Polynomial 671 std::vector<Coefficient_type> coefficients; 672 673 Coefficient_type zero(0); 674 while(begin != end){ 675 int current_exp = begin->first[PT::d]; 676 // fill up with zeros until current exp is reached 677 if((int) coefficients.size() < current_exp) 678 coefficients.resize(current_exp,zero); 679 680 // select range for coefficient of current exponent 681 Monom_rep_iterator coeff_end = begin; 682 while( coeff_end != end && coeff_end->first[PT::d] == current_exp ){ 683 ++coeff_end; 684 } 685 coefficients.push_back(construct(begin, coeff_end)); 686 begin = coeff_end; 687 } 688 return Polynomial_d(coefficients.begin(),coefficients.end()); 689 } 690 }; 691 }; 692 693 // Get_coefficient; 694 struct Get_coefficient 695 : public CGAL::cpp98::binary_function<Polynomial_d, int, Coefficient_type > { 696 operatorGet_coefficient697 const Coefficient_type& operator()( const Polynomial_d& p, int i) const { 698 CGAL_STATIC_THREAD_LOCAL_VARIABLE(Coefficient_type, zero, 0); 699 CGAL_precondition( i >= 0 ); 700 typename PT::Degree degree; 701 if( i > degree(p) ) 702 return zero; 703 return p[i]; 704 } 705 }; 706 707 // Get_innermost_coefficient; 708 struct Get_innermost_coefficient 709 : public CGAL::cpp98::binary_function< Polynomial_d, 710 Exponent_vector, 711 Innermost_coefficient_type > 712 { 713 714 const Innermost_coefficient_type& operatorGet_innermost_coefficient715 operator()( const Polynomial_d& p, Exponent_vector ev ) const { 716 CGAL_precondition( !ev.empty() ); 717 typename PTC::Get_innermost_coefficient gic; 718 typename PT::Get_coefficient gc; 719 int exponent = ev.back(); 720 ev.pop_back(); 721 return gic( gc( p, exponent ), ev ); 722 }; 723 }; 724 725 typedef CGAL::internal::Monomial_representation<Polynomial_d> Monomial_representation; 726 727 // Swap variable x_i with x_j 728 struct Swap { 729 typedef Polynomial_d result_type; 730 typedef Polynomial_d first_argument_type; 731 typedef int second_argument_type; 732 typedef int third_argument_type; 733 734 public: 735 operatorSwap736 Polynomial_d operator()(const Polynomial_d& p, int i, int j ) const { 737 CGAL_precondition(0 <= i && i < d); 738 CGAL_precondition(0 <= j && j < d); 739 typedef std::pair< Exponent_vector, Innermost_coefficient_type > 740 Exponents_coeff_pair; 741 Monomial_representation gmr; 742 Construct_polynomial construct; 743 typedef std::vector< Exponents_coeff_pair > Monom_vector; 744 typedef typename Monom_vector::iterator MVIterator; 745 Monom_vector monoms; 746 gmr( p, std::back_inserter( monoms ) ); 747 for( MVIterator it = monoms.begin(); it != monoms.end(); ++it ) { 748 std::swap(it->first[i],it->first[j]); 749 } 750 // sort only once ! 751 std::sort(monoms.begin(), monoms.end(),Compare_exponents_coeff_pair()); 752 return construct(monoms.begin(), monoms.end(),true); 753 } 754 }; 755 756 757 // Move; 758 // move variable x_i to position of x_j 759 // order of other variables remains 760 // default j = d makes x_i the othermost variable 761 struct Move { 762 typedef Polynomial_d result_type; 763 typedef Polynomial_d first_argument_type; 764 typedef int second_argument_type; 765 typedef int third_argument_type; 766 767 Polynomial_d operatorMove768 operator()(const Polynomial_d& p, int i, int j = (d-1) ) const { 769 CGAL_precondition(0 <= i && i < d); 770 CGAL_precondition(0 <= j && j < d); 771 typedef std::pair< Exponent_vector, Innermost_coefficient_type > 772 Exponents_coeff_pair; 773 typedef std::vector< Exponents_coeff_pair > Monom_rep; 774 Monomial_representation gmr; 775 Construct_polynomial construct; 776 Monom_rep mon_rep; 777 gmr( p, std::back_inserter( mon_rep ) ); 778 for( typename Monom_rep::iterator it = mon_rep.begin(); 779 it != mon_rep.end(); 780 ++it ) { 781 // this is as good as std::rotate since it uses swap also 782 if (i < j) 783 for( int k = i; k < j; k++ ) 784 std::swap(it->first[k],it->first[k+1]); 785 else 786 for( int k = i; k > j; k-- ) 787 std::swap(it->first[k],it->first[k-1]); 788 789 } 790 return construct( mon_rep.begin(), mon_rep.end() ); 791 } 792 }; 793 794 795 struct Permute { 796 typedef Polynomial_d result_type; operatorPermute797 template <typename Input_iterator> Polynomial_d operator() 798 (const Polynomial_d& p, Input_iterator first, Input_iterator last) const { 799 Construct_polynomial construct; 800 Monomial_representation gmr; 801 Monom_rep mon_rep; 802 gmr( p, std::back_inserter( mon_rep )); 803 std::vector<int> on_place, number_is; 804 int i= 0; 805 for (Input_iterator iter = first ; iter != last ; ++iter) 806 number_is.push_back (i++); 807 on_place = number_is; 808 int rem_place = 0, rem_number = i= 0; 809 for(Input_iterator iter = first ; iter != last ; ++iter){ 810 for( typename Monom_rep::iterator it = mon_rep.begin(); it != 811 mon_rep.end(); ++it ) 812 std::swap(it->first[number_is[i]],it->first[(*iter)]); 813 814 815 rem_place= number_is[i]; 816 rem_number= on_place[(*iter)]; 817 818 on_place[(*iter)] = i; 819 on_place[rem_place]=rem_number; 820 number_is[rem_number]=rem_place; 821 number_is[i++]= (*iter); 822 } 823 return construct( mon_rep.begin(), mon_rep.end() ); 824 } 825 }; 826 827 //Degree; 828 typedef CGAL::internal::Degree<Polynomial_d> Degree; 829 830 // Total_degree; 831 struct Total_degree : public CGAL::cpp98::unary_function< Polynomial_d , int >{ operatorTotal_degree832 int operator()(const Polynomial_d& p) const { 833 typedef Polynomial_traits_d<Coefficient_type> COEFF_POLY_TRAITS; 834 typename COEFF_POLY_TRAITS::Total_degree total_degree; 835 Degree degree; 836 CGAL_precondition( degree(p) >= 0); 837 838 int result = 0; 839 for(int i = 0; i <= degree(p) ; i++){ 840 if( ! CGAL::is_zero( p[i]) ) 841 result = (std::max)(result , total_degree(p[i]) + i ); 842 } 843 return result; 844 } 845 }; 846 847 // Leading_coefficient; 848 struct Leading_coefficient 849 : public CGAL::cpp98::unary_function< Polynomial_d , Coefficient_type>{ operatorLeading_coefficient850 const Coefficient_type& operator()(const Polynomial_d& p) const { 851 return p.lcoeff(); 852 } 853 }; 854 855 // Innermost_leading_coefficient; 856 struct Innermost_leading_coefficient 857 : public CGAL::cpp98::unary_function< Polynomial_d , Innermost_coefficient_type>{ 858 const Innermost_coefficient_type& operatorInnermost_leading_coefficient859 operator()(const Polynomial_d& p) const { 860 typename PTC::Innermost_leading_coefficient ilcoeff; 861 typename PT::Leading_coefficient lcoeff; 862 return ilcoeff(lcoeff(p)); 863 } 864 }; 865 866 867 //return a canonical representative of all constant multiples. 868 struct Canonicalize 869 : public CGAL::cpp98::unary_function<Polynomial_d, Polynomial_d>{ 870 871 private: canonicalize_Canonicalize872 inline Polynomial_d canonicalize_(Polynomial_d p, CGAL::Tag_true) const 873 { 874 typedef typename Polynomial_traits_d<Polynomial_d>::Innermost_coefficient_type IC; 875 typename Polynomial_traits_d<Polynomial_d>::Innermost_leading_coefficient ilcoeff; 876 typename Algebraic_extension_traits<IC>::Normalization_factor nfac; 877 878 IC tmp = nfac(ilcoeff(p)); 879 if(tmp != IC(1)){ 880 p *= Polynomial_d(tmp); 881 } 882 remove_scalar_factor(p); 883 p /= p.unit_part(); 884 p.simplify_coefficients(); 885 886 CGAL_postcondition(nfac(ilcoeff(p)) == IC(1)); 887 return p; 888 }; 889 canonicalize_Canonicalize890 inline Polynomial_d canonicalize_(Polynomial_d p, CGAL::Tag_false) const 891 { 892 remove_scalar_factor(p); 893 p /= p.unit_part(); 894 p.simplify_coefficients(); 895 return p; 896 }; 897 898 public: 899 Polynomial_d operatorCanonicalize900 operator()( const Polynomial_d& p ) const { 901 if (CGAL::is_zero(p)) return p; 902 903 typedef Innermost_coefficient_type IC; 904 typedef typename Algebraic_extension_traits<IC>::Is_extended Is_extended; 905 return canonicalize_(p, Is_extended()); 906 } 907 }; 908 909 // Differentiate; 910 struct Differentiate 911 : public CGAL::cpp98::unary_function<Polynomial_d, Polynomial_d>{ 912 Polynomial_d operatorDifferentiate913 operator()(Polynomial_d p, int i = (d-1)) const { 914 if (i == (d-1) ){ 915 p.diff(); 916 }else{ 917 Swap swap; 918 p = swap(p,i,d-1); 919 p.diff(); 920 p = swap(p,i,d-1); 921 } 922 return p; 923 } 924 }; 925 926 // Evaluate; 927 struct Evaluate 928 :public CGAL::cpp98::binary_function<Polynomial_d,Coefficient_type,Coefficient_type>{ 929 // Evaluate with respect to one variable 930 Coefficient_type operatorEvaluate931 operator()(const Polynomial_d& p, const Coefficient_type& x) const { 932 return p.evaluate(x); 933 } 934 #define ICOEFF typename First_if_different<Innermost_coefficient_type, Coefficient_type>::Type operatorEvaluate935 Coefficient_type operator() 936 ( const Polynomial_d& p, const ICOEFF& x) const 937 { 938 return p.evaluate(x); 939 } 940 #undef ICOEFF 941 }; 942 943 // Evaluate_homogeneous; 944 struct Evaluate_homogeneous{ 945 typedef Coefficient_type result_type; 946 typedef Polynomial_d first_argument_type; 947 typedef Coefficient_type second_argument_type; 948 typedef Coefficient_type third_argument_type; 949 operatorEvaluate_homogeneous950 Coefficient_type operator()( 951 const Polynomial_d& p, const Coefficient_type& a, const Coefficient_type& b) const 952 { 953 return p.evaluate_homogeneous(a,b); 954 } 955 #define ICOEFF typename First_if_different<Innermost_coefficient_type, Coefficient_type>::Type operatorEvaluate_homogeneous956 Coefficient_type operator() 957 ( const Polynomial_d& p, const ICOEFF& a, const ICOEFF& b) const 958 { 959 return p.evaluate_homogeneous(a,b); 960 } 961 #undef ICOEFF 962 963 }; 964 965 // Is_zero_at; 966 struct Is_zero_at { 967 private: 968 typedef Algebraic_structure_traits<Innermost_coefficient_type> AST; 969 typedef typename AST::Is_zero::result_type BOOL; 970 public: 971 typedef BOOL result_type; 972 973 template< class Input_iterator > operatorIs_zero_at974 BOOL operator()( 975 const Polynomial_d& p, 976 Input_iterator begin, 977 Input_iterator end ) const { 978 typename PT::Substitute substitute; 979 return( CGAL::is_zero( substitute( p, begin, end ) ) ); 980 } 981 }; 982 983 // Is_zero_at_homogeneous; 984 struct Is_zero_at_homogeneous { 985 private: 986 typedef Algebraic_structure_traits<Innermost_coefficient_type> AST; 987 typedef typename AST::Is_zero::result_type BOOL; 988 public: 989 typedef BOOL result_type; 990 991 template< class Input_iterator > operatorIs_zero_at_homogeneous992 BOOL operator() 993 ( const Polynomial_d& p, Input_iterator begin, Input_iterator end ) const 994 { 995 typename PT::Substitute_homogeneous substitute_homogeneous; 996 return( CGAL::is_zero( substitute_homogeneous( p, begin, end ) ) ); 997 } 998 }; 999 1000 // Sign_at, Sign_at_homogeneous, Compare 1001 // define XXX_ even though ICoeff may not be Real_embeddable 1002 // select propoer XXX among XXX_ or Null_functor using ::boost::mpl::if_ 1003 private: 1004 struct Sign_at_ { 1005 private: 1006 typedef Real_embeddable_traits<Innermost_coefficient_type> RT; 1007 public: 1008 typedef typename RT::Sign result_type; 1009 1010 template< class Input_iterator > operatorSign_at_1011 result_type operator()( 1012 const Polynomial_d& p, 1013 Input_iterator begin, 1014 Input_iterator end ) const 1015 { 1016 typename PT::Substitute substitute; 1017 return CGAL::sign( substitute( p, begin, end ) ); 1018 } 1019 }; 1020 1021 struct Sign_at_homogeneous_ { 1022 typedef Real_embeddable_traits<Innermost_coefficient_type> RT; 1023 public: 1024 typedef typename RT::Sign result_type; 1025 1026 template< class Input_iterator > operatorSign_at_homogeneous_1027 result_type operator()( 1028 const Polynomial_d& p, 1029 Input_iterator begin, 1030 Input_iterator end) const { 1031 typename PT::Substitute_homogeneous substitute_homogeneous; 1032 return CGAL::sign( substitute_homogeneous( p, begin, end ) ); 1033 } 1034 }; 1035 1036 typedef Real_embeddable_traits<Innermost_coefficient_type> RET_IC; 1037 typedef typename RET_IC::Is_real_embeddable IC_is_real_embeddable; 1038 public: 1039 typedef typename ::boost::mpl::if_<IC_is_real_embeddable,Sign_at_,Null_functor>::type Sign_at; 1040 typedef typename ::boost::mpl::if_<IC_is_real_embeddable,Sign_at_homogeneous_,Null_functor>::type Sign_at_homogeneous; 1041 typedef typename Real_embeddable_traits<Polynomial_d>::Compare Compare; 1042 1043 1044 struct Construct_coefficient_const_iterator_range 1045 : public CGAL::cpp98::unary_function< Polynomial_d, 1046 Coefficient_const_iterator_range> { 1047 Coefficient_const_iterator_range operatorConstruct_coefficient_const_iterator_range1048 operator () (const Polynomial_d& p) const { 1049 return std::make_pair( p.begin(), p.end() ); 1050 } 1051 }; 1052 1053 struct Construct_innermost_coefficient_const_iterator_range 1054 : public CGAL::cpp98::unary_function< Polynomial_d, 1055 Innermost_coefficient_const_iterator_range> { 1056 Innermost_coefficient_const_iterator_range operatorConstruct_innermost_coefficient_const_iterator_range1057 operator () (const Polynomial_d& p) const { 1058 return std::make_pair( 1059 typename Coefficient_const_flattening::Flatten()(p.end(),p.begin()), 1060 typename Coefficient_const_flattening::Flatten()(p.end(),p.end())); 1061 } 1062 }; 1063 1064 struct Is_square_free 1065 : public CGAL::cpp98::unary_function< Polynomial_d, bool >{ operatorIs_square_free1066 bool operator()( const Polynomial_d& p ) const { 1067 if( !internal::may_have_multiple_factor( p ) ) 1068 return true; 1069 1070 Gcd_up_to_constant_factor gcd_utcf; 1071 Univariate_content_up_to_constant_factor ucontent_utcf; 1072 Integral_division_up_to_constant_factor idiv_utcf; 1073 Differentiate diff; 1074 1075 Coefficient_type content = ucontent_utcf( p ); 1076 typename PTC::Is_square_free isf; 1077 1078 if( !isf( content ) ) 1079 return false; 1080 1081 Polynomial_d regular_part = idiv_utcf( p, Polynomial_d( content ) ); 1082 1083 Polynomial_d g = gcd_utcf(regular_part,diff(regular_part)); 1084 return ( g.degree() == 0 ); 1085 } 1086 }; 1087 1088 1089 struct Make_square_free 1090 : public CGAL::cpp98::unary_function< Polynomial_d, Polynomial_d >{ 1091 Polynomial_d operatorMake_square_free1092 operator()(const Polynomial_d& p) const { 1093 if (CGAL::is_zero(p)) return p; 1094 Gcd_up_to_constant_factor gcd_utcf; 1095 Univariate_content_up_to_constant_factor ucontent_utcf; 1096 Integral_division_up_to_constant_factor idiv_utcf; 1097 Differentiate diff; 1098 typename PTC::Make_square_free msf; 1099 1100 Coefficient_type content = ucontent_utcf(p); 1101 Polynomial_d result = Polynomial_d(msf(content)); 1102 1103 Polynomial_d regular_part = idiv_utcf(p,Polynomial_d(content)); 1104 Polynomial_d g = gcd_utcf(regular_part,diff(regular_part)); 1105 1106 1107 result *= idiv_utcf(regular_part,g); 1108 return Canonicalize()(result); 1109 1110 } 1111 }; 1112 1113 // Pseudo_division; 1114 struct Pseudo_division { 1115 typedef Polynomial_d result_type; 1116 void operatorPseudo_division1117 operator()( 1118 const Polynomial_d& f, const Polynomial_d& g, 1119 Polynomial_d& q, Polynomial_d& r, Coefficient_type& D) const { 1120 Polynomial_d::pseudo_division(f,g,q,r,D); 1121 } 1122 }; 1123 1124 struct Pseudo_division_quotient 1125 :public CGAL::cpp98::binary_function<Polynomial_d, Polynomial_d, Polynomial_d> { 1126 1127 Polynomial_d operatorPseudo_division_quotient1128 operator()(const Polynomial_d& f, const Polynomial_d& g) const { 1129 Polynomial_d q,r; 1130 Coefficient_type D; 1131 Polynomial_d::pseudo_division(f,g,q,r,D); 1132 return q; 1133 } 1134 }; 1135 1136 struct Pseudo_division_remainder 1137 :public CGAL::cpp98::binary_function<Polynomial_d, Polynomial_d, Polynomial_d> { 1138 1139 Polynomial_d operatorPseudo_division_remainder1140 operator()(const Polynomial_d& f, const Polynomial_d& g) const { 1141 Polynomial_d q,r; 1142 Coefficient_type D; 1143 Polynomial_d::pseudo_division(f,g,q,r,D); 1144 return r; 1145 } 1146 }; 1147 1148 struct Gcd_up_to_constant_factor 1149 :public CGAL::cpp98::binary_function<Polynomial_d, Polynomial_d, Polynomial_d> { 1150 Polynomial_d operatorGcd_up_to_constant_factor1151 operator()(const Polynomial_d& p, const Polynomial_d& q) const { 1152 if(p==q) return CGAL::canonicalize(p); 1153 if (CGAL::is_zero(p) && CGAL::is_zero(q)){ 1154 return Polynomial_d(0); 1155 } 1156 // apply modular filter first 1157 if (internal::may_have_common_factor(p,q)){ 1158 return internal::gcd_utcf_(p,q); 1159 }else{ 1160 return Polynomial_d(1); 1161 } 1162 } 1163 }; 1164 1165 struct Integral_division_up_to_constant_factor 1166 :public CGAL::cpp98::binary_function<Polynomial_d, Polynomial_d, Polynomial_d> { 1167 1168 1169 1170 Polynomial_d operatorIntegral_division_up_to_constant_factor1171 operator()(const Polynomial_d& p, const Polynomial_d& q) const { 1172 typedef Innermost_coefficient_type IC; 1173 1174 typename PT::Construct_polynomial construct; 1175 typename PT::Innermost_leading_coefficient ilcoeff; 1176 typename PT::Construct_innermost_coefficient_const_iterator_range range; 1177 typedef Algebraic_extension_traits<Innermost_coefficient_type> AET; 1178 typename AET::Denominator_for_algebraic_integers dfai; 1179 typename AET::Normalization_factor nfac; 1180 1181 1182 IC ilcoeff_q = ilcoeff(q); 1183 // this factor is needed in case IC is an Algebraic extension 1184 IC dfai_q = dfai(range(q).first, range(q).second); 1185 // make dfai_q a 'scalar' 1186 ilcoeff_q *= dfai_q * nfac(dfai_q); 1187 1188 Polynomial_d result = (p * construct(ilcoeff_q)) / q; 1189 1190 return Canonicalize()(result); 1191 } 1192 }; 1193 1194 struct Univariate_content_up_to_constant_factor 1195 :public CGAL::cpp98::unary_function<Polynomial_d, Coefficient_type> { 1196 Coefficient_type operatorUnivariate_content_up_to_constant_factor1197 operator()(const Polynomial_d& p) const { 1198 typename PTC::Gcd_up_to_constant_factor gcd_utcf; 1199 1200 if(CGAL::is_zero(p)) return Coefficient_type(0); 1201 if(PT::d == 1) return Coefficient_type(1); 1202 1203 Coefficient_type result(0); 1204 for(typename Polynomial_d::const_iterator it = p.begin(); 1205 it != p.end(); 1206 it++){ 1207 result = gcd_utcf(*it,result); 1208 } 1209 return result; 1210 1211 } 1212 }; 1213 1214 struct Square_free_factorize_up_to_constant_factor { 1215 private: 1216 typedef Coefficient_type Coeff; 1217 typedef Innermost_coefficient_type ICoeff; 1218 1219 // rsqff_utcf computes the sqff recursively for Coeff 1220 // end of recursion: ICoeff 1221 1222 template < class OutputIterator > rsqff_utcfSquare_free_factorize_up_to_constant_factor1223 OutputIterator rsqff_utcf ( ICoeff , OutputIterator oi) const{ 1224 return oi; 1225 } 1226 1227 template < class OutputIterator > rsqff_utcfSquare_free_factorize_up_to_constant_factor1228 OutputIterator rsqff_utcf ( 1229 typename First_if_different<Coeff,ICoeff>::Type c, 1230 OutputIterator oi) const { 1231 1232 typename PTC::Square_free_factorize_up_to_constant_factor sqff; 1233 std::vector<std::pair<Coefficient_type,int> > fac_mul_pairs; 1234 sqff(c,std::back_inserter(fac_mul_pairs)); 1235 1236 for(unsigned int i = 0; i < fac_mul_pairs.size(); i++){ 1237 Polynomial_d factor(fac_mul_pairs[i].first); 1238 int mult = fac_mul_pairs[i].second; 1239 *oi++=std::make_pair(factor,mult); 1240 } 1241 return oi; 1242 } 1243 1244 public: 1245 template < class OutputIterator> 1246 OutputIterator operatorSquare_free_factorize_up_to_constant_factor1247 operator()(Polynomial_d p, OutputIterator oi) const { 1248 if (CGAL::is_zero(p)) return oi; 1249 1250 Univariate_content_up_to_constant_factor ucontent_utcf; 1251 Integral_division_up_to_constant_factor idiv_utcf; 1252 Coefficient_type c = ucontent_utcf(p); 1253 1254 p = idiv_utcf( p , Polynomial_d(c)); 1255 std::vector<Polynomial_d> factors; 1256 std::vector<int> mults; 1257 square_free_factorize_utcf( 1258 p, std::back_inserter(factors), std::back_inserter(mults)); 1259 for(unsigned int i = 0; i < factors.size() ; i++){ 1260 *oi++=std::make_pair(factors[i],mults[i]); 1261 } 1262 if (CGAL::total_degree(c) == 0) 1263 return oi; 1264 else 1265 return rsqff_utcf(c,oi); 1266 } 1267 }; 1268 1269 struct Shift 1270 : public CGAL::cpp98::binary_function< Polynomial_d,int,Polynomial_d >{ 1271 1272 Polynomial_d operatorShift1273 operator()(const Polynomial_d& p, int e, int i = (d-1)) const { 1274 Construct_polynomial construct; 1275 Monomial_representation gmr; 1276 Monom_rep monom_rep; 1277 gmr(p,std::back_inserter(monom_rep)); 1278 for(typename Monom_rep::iterator it = monom_rep.begin(); 1279 it != monom_rep.end(); 1280 it++){ 1281 it->first[i]+=e; 1282 } 1283 return construct(monom_rep.begin(), monom_rep.end()); 1284 } 1285 }; 1286 1287 struct Negate 1288 : public CGAL::cpp98::unary_function< Polynomial_d, Polynomial_d >{ 1289 operatorNegate1290 Polynomial_d operator()(const Polynomial_d& p, int i = (d-1)) const { 1291 Construct_polynomial construct; 1292 Monomial_representation gmr; 1293 Monom_rep monom_rep; 1294 gmr(p,std::back_inserter(monom_rep)); 1295 for(typename Monom_rep::iterator it = monom_rep.begin(); 1296 it != monom_rep.end(); 1297 it++){ 1298 if (it->first[i] % 2 != 0) 1299 it->second = - it->second; 1300 } 1301 return construct(monom_rep.begin(), monom_rep.end()); 1302 } 1303 }; 1304 1305 struct Invert 1306 : public CGAL::cpp98::unary_function< Polynomial_d , Polynomial_d >{ operatorInvert1307 Polynomial_d operator()(Polynomial_d p, int i = (PT::d-1)) const { 1308 if (i == (d-1)){ 1309 p.reversal(); 1310 }else{ 1311 p = Swap()(p,i,PT::d-1); 1312 p.reversal(); 1313 p = Swap()(p,i,PT::d-1); 1314 } 1315 return p ; 1316 } 1317 }; 1318 1319 struct Translate 1320 : public CGAL::cpp98::binary_function< Polynomial_d , Innermost_coefficient_type, 1321 Polynomial_d >{ 1322 Polynomial_d operatorTranslate1323 operator()( 1324 Polynomial_d p, 1325 const Innermost_coefficient_type& c, 1326 int i = (d-1)) 1327 const { 1328 if (i == (d-1) ){ 1329 p.translate(Coefficient_type(c)); 1330 }else{ 1331 Swap swap; 1332 p = swap(p,i,d-1); 1333 p.translate(Coefficient_type(c)); 1334 p = swap(p,i,d-1); 1335 } 1336 return p; 1337 } 1338 }; 1339 1340 struct Translate_homogeneous{ 1341 typedef Polynomial_d result_type; 1342 typedef Polynomial_d first_argument_type; 1343 typedef Innermost_coefficient_type second_argument_type; 1344 typedef Innermost_coefficient_type third_argument_type; 1345 1346 Polynomial_d operatorTranslate_homogeneous1347 operator()(Polynomial_d p, 1348 const Innermost_coefficient_type& a, 1349 const Innermost_coefficient_type& b, 1350 int i = (d-1) ) const { 1351 if (i == (d-1) ){ 1352 p.translate(Coefficient_type(a),Coefficient_type(b)); 1353 }else{ 1354 Swap swap; 1355 p = swap(p,i,d-1); 1356 p.translate(Coefficient_type(a),Coefficient_type(b)); 1357 p = swap(p,i,d-1); 1358 } 1359 return p; 1360 } 1361 }; 1362 1363 struct Scale 1364 : public CGAL::cpp98::binary_function< Polynomial_d, 1365 Innermost_coefficient_type, 1366 Polynomial_d > 1367 { operatorScale1368 Polynomial_d operator()( Polynomial_d p, const Innermost_coefficient_type& c, 1369 int i = (PT::d-1) ) const { 1370 CGAL_precondition( i <= d-1 ); 1371 CGAL_precondition( i >= 0 ); 1372 typename PT::Scale_homogeneous scale_homogeneous; 1373 return scale_homogeneous( p, c, Innermost_coefficient_type(1), i ); 1374 } 1375 1376 }; 1377 1378 struct Scale_homogeneous{ 1379 typedef Polynomial_d result_type; 1380 typedef Polynomial_d first_argument_type; 1381 typedef Innermost_coefficient_type second_argument_type; 1382 typedef Innermost_coefficient_type third_argument_type; 1383 1384 Polynomial_d operatorScale_homogeneous1385 operator()( 1386 Polynomial_d p, 1387 const Innermost_coefficient_type& a, 1388 const Innermost_coefficient_type& b, 1389 int i = (d-1) ) const { 1390 1391 CGAL_precondition( ! CGAL::is_zero(b) ); 1392 CGAL_precondition( i <= d-1 ); 1393 CGAL_precondition( i >= 0 ); 1394 1395 if (i != (d-1) ) p = Swap()(p,i,d-1); 1396 1397 if(CGAL::is_one(b)) 1398 p.scale_up(Coefficient_type(a)); 1399 else 1400 if(CGAL::is_one(a)) 1401 p.scale_down(Coefficient_type(b)); 1402 else 1403 p.scale(Coefficient_type(a),Coefficient_type(b)); 1404 1405 if (i != (d-1) ) p = Swap()(p,i,d-1); 1406 1407 return p; 1408 } 1409 }; 1410 1411 struct Resultant 1412 : public CGAL::cpp98::binary_function<Polynomial_d, Polynomial_d, Coefficient_type>{ 1413 1414 Coefficient_type operatorResultant1415 operator()( 1416 const Polynomial_d& p, 1417 const Polynomial_d& q) const { 1418 return internal::resultant(p,q); 1419 } 1420 }; 1421 1422 // Polynomial subresultants (aka subresultant polynomials) 1423 struct Polynomial_subresultants { 1424 1425 template<typename OutputIterator> operatorPolynomial_subresultants1426 OutputIterator operator()( 1427 const Polynomial_d& p, 1428 const Polynomial_d& q, 1429 OutputIterator out, 1430 int i = (d-1) ) const { 1431 if(i == (d-1) ) 1432 return CGAL::internal::polynomial_subresultants<PT>(p,q,out); 1433 else 1434 return CGAL::internal::polynomial_subresultants<PT>(Move()(p,i), 1435 Move()(q,i), 1436 out); 1437 } 1438 }; 1439 1440 // principal subresultants (aka scalar subresultants) 1441 struct Principal_subresultants { 1442 1443 template<typename OutputIterator> operatorPrincipal_subresultants1444 OutputIterator operator()( 1445 const Polynomial_d& p, 1446 const Polynomial_d& q, 1447 OutputIterator out, 1448 int i = (d-1) ) const { 1449 if(i == (d-1) ) 1450 return CGAL::internal::principal_subresultants<PT>(p,q,out); 1451 else 1452 return CGAL::internal::principal_subresultants<PT>(Move()(p,i), 1453 Move()(q,i), 1454 out); 1455 } 1456 }; 1457 1458 // Subresultants with cofactors 1459 struct Polynomial_subresultants_with_cofactors { 1460 1461 template<typename OutputIterator1, 1462 typename OutputIterator2, 1463 typename OutputIterator3> operatorPolynomial_subresultants_with_cofactors1464 OutputIterator1 operator()( 1465 const Polynomial_d& p, 1466 const Polynomial_d& q, 1467 OutputIterator1 out_sres, 1468 OutputIterator2 out_co_p, 1469 OutputIterator3 out_co_q, 1470 int i = (d-1) ) const { 1471 if(i == (d-1) ) 1472 return CGAL::internal::polynomial_subresultants_with_cofactors<PT> 1473 (p,q,out_sres,out_co_p,out_co_q); 1474 else 1475 return CGAL::internal::polynomial_subresultants_with_cofactors<PT> 1476 (Move()(p,i),Move()(q,i),out_sres,out_co_p,out_co_q); 1477 } 1478 }; 1479 1480 // Sturm-Habicht sequence (aka signed subresultant sequence) 1481 struct Sturm_habicht_sequence { 1482 1483 template<typename OutputIterator> operatorSturm_habicht_sequence1484 OutputIterator operator()( 1485 const Polynomial_d& p, 1486 OutputIterator out, 1487 int i = (d-1) ) const { 1488 if(i == (d-1) ) 1489 return CGAL::internal::sturm_habicht_sequence<PT>(p,out); 1490 else 1491 return CGAL::internal::sturm_habicht_sequence<PT>(Move()(p,i), 1492 out); 1493 } 1494 }; 1495 1496 // Sturm-Habicht sequence with cofactors 1497 struct Sturm_habicht_sequence_with_cofactors { 1498 1499 template<typename OutputIterator1, 1500 typename OutputIterator2, 1501 typename OutputIterator3> operatorSturm_habicht_sequence_with_cofactors1502 OutputIterator1 operator()( 1503 const Polynomial_d& p, 1504 OutputIterator1 out_stha, 1505 OutputIterator2 out_f, 1506 OutputIterator3 out_fx, 1507 int i = (d-1) ) const { 1508 if(i == (d-1) ) 1509 return CGAL::internal::sturm_habicht_sequence_with_cofactors<PT> 1510 (p,out_stha,out_f,out_fx); 1511 else 1512 return CGAL::internal::sturm_habicht_sequence_with_cofactors<PT> 1513 (Move()(p,i),out_stha,out_f,out_fx); 1514 } 1515 }; 1516 1517 // Principal Sturm-Habicht sequence (formal leading coefficients 1518 // of Sturm-Habicht sequence) 1519 struct Principal_sturm_habicht_sequence { 1520 1521 template<typename OutputIterator> operatorPrincipal_sturm_habicht_sequence1522 OutputIterator operator()( 1523 const Polynomial_d& p, 1524 OutputIterator out, 1525 int i = (d-1) ) const { 1526 if(i == (d-1) ) 1527 return CGAL::internal::principal_sturm_habicht_sequence<PT>(p,out); 1528 else 1529 return CGAL::internal::principal_sturm_habicht_sequence<PT> 1530 (Move()(p,i),out); 1531 } 1532 }; 1533 1534 1535 // returns the Exponten_vector of the innermost leading coefficient 1536 struct Degree_vector{ 1537 typedef Exponent_vector result_type; 1538 typedef Polynomial_d argument_type; 1539 1540 // returns the exponent vector of inner_most_lcoeff. operatorDegree_vector1541 result_type operator()(const Polynomial_d& polynomial) const{ 1542 typename PTC::Degree_vector degree_vector; 1543 Exponent_vector result = degree_vector(polynomial.lcoeff()); 1544 result.push_back(polynomial.degree()); 1545 return result; 1546 } 1547 }; 1548 1549 // substitute every variable by its new value in the iterator range 1550 // begin refers to the innermost/first variable 1551 struct Substitute{ 1552 public: 1553 template <class Input_iterator> 1554 typename CGAL::Coercion_traits< 1555 typename std::iterator_traits<Input_iterator>::value_type, 1556 Innermost_coefficient_type 1557 >::Type operatorSubstitute1558 operator()( 1559 const Polynomial_d& p, 1560 Input_iterator begin, 1561 Input_iterator end) const { 1562 typedef typename std::iterator_traits<Input_iterator> ITT; 1563 typedef typename ITT::iterator_category Category; 1564 return (*this)(p,begin,end,Category()); 1565 } 1566 1567 template <class Input_iterator> 1568 typename CGAL::Coercion_traits< 1569 typename std::iterator_traits<Input_iterator>::value_type, 1570 Innermost_coefficient_type 1571 >::Type operatorSubstitute1572 operator()( 1573 const Polynomial_d& p, 1574 Input_iterator begin, 1575 Input_iterator end, 1576 std::forward_iterator_tag) const { 1577 typedef typename std::iterator_traits<Input_iterator> ITT; 1578 std::list<typename ITT::value_type> list(begin,end); 1579 return (*this)(p,list.begin(),list.end()); 1580 } 1581 1582 template <class Input_iterator> 1583 typename 1584 CGAL::Coercion_traits 1585 <typename std::iterator_traits<Input_iterator>::value_type, 1586 Innermost_coefficient_type>::Type operatorSubstitute1587 operator()( 1588 const Polynomial_d& p, 1589 Input_iterator begin, 1590 Input_iterator end, 1591 std::bidirectional_iterator_tag) const { 1592 1593 typedef typename std::iterator_traits<Input_iterator>::value_type 1594 value_type; 1595 typedef CGAL::Coercion_traits<Innermost_coefficient_type,value_type> CT; 1596 typename PTC::Substitute subs; 1597 1598 typename CT::Type x = typename CT::Cast()(*(--end)); 1599 1600 int i = Degree()(p); 1601 typename CT::Type y = 1602 subs(Get_coefficient()(p,i),begin,end); 1603 1604 while (--i >= 0){ 1605 y *= x; 1606 y += subs(Get_coefficient()(p,i),begin,end); 1607 } 1608 return y; 1609 } 1610 }; // substitute every variable by its new value in the iterator range 1611 1612 1613 1614 // begin refers to the innermost/first variable 1615 struct Substitute_homogeneous{ 1616 1617 template<typename Input_iterator> 1618 struct Result_type{ 1619 typedef std::iterator_traits<Input_iterator> ITT; 1620 typedef typename ITT::value_type value_type; 1621 typedef Coercion_traits<value_type, Innermost_coefficient_type> CT; 1622 typedef typename CT::Type Type; 1623 }; 1624 1625 public: 1626 1627 template <class Input_iterator> 1628 typename Result_type<Input_iterator>::Type operatorSubstitute_homogeneous1629 operator()( const Polynomial_d& p, Input_iterator begin, Input_iterator end) const{ 1630 int hdegree = Total_degree()(p); 1631 1632 typedef std::iterator_traits<Input_iterator> ITT; 1633 std::list<typename ITT::value_type> list(begin,end); 1634 1635 // make the homogeneous variable the first in the list 1636 list.push_front(list.back()); 1637 list.pop_back(); 1638 1639 // reverse and begin with the outermost variable 1640 return (*this)(p, list.rbegin(), list.rend(), hdegree); 1641 } 1642 1643 // this operator is undcoumented and for internal use: 1644 // the iterator range starts with the outermost variable 1645 // and ends with the homogeneous variable 1646 template <class Input_iterator> 1647 typename Result_type<Input_iterator>::Type operatorSubstitute_homogeneous1648 operator()( 1649 const Polynomial_d& p, 1650 Input_iterator begin, 1651 Input_iterator end, 1652 int hdegree) const{ 1653 1654 1655 typedef std::iterator_traits<Input_iterator> ITT; 1656 typedef typename ITT::value_type value_type; 1657 typedef Coercion_traits<value_type, Innermost_coefficient_type> CT; 1658 1659 typename PTC::Substitute_homogeneous subsh; 1660 1661 typename CT::Type x = typename CT::Cast()(*begin++); 1662 1663 1664 int i = Degree()(p); 1665 typename CT::Type y = subsh(Get_coefficient()(p,i),begin,end, hdegree-i); 1666 1667 while (--i >= 0){ 1668 y *= x; 1669 y += subsh(Get_coefficient()(p,i),begin,end,hdegree-i); 1670 } 1671 return y; 1672 } 1673 }; 1674 1675 }; 1676 1677 } // namespace internal 1678 1679 // Definition of Polynomial_traits_d 1680 // 1681 // In order to determine the algebraic category of the innermost coefficient, 1682 // the Polynomial_traits_d_base class with "Null_tag" is used. 1683 1684 template< class Polynomial > 1685 class Polynomial_traits_d 1686 : public internal::Polynomial_traits_d_base< Polynomial, 1687 typename Algebraic_structure_traits< 1688 typename internal::Innermost_coefficient_type<Polynomial>::Type >::Algebraic_category, 1689 typename Algebraic_structure_traits< Polynomial >::Algebraic_category > , 1690 public Algebraic_structure_traits<Polynomial>{ 1691 1692 //------------ Rebind ----------- 1693 private: 1694 template <class T, int d> 1695 struct Gen_polynomial_type{ 1696 typedef CGAL::Polynomial<typename Gen_polynomial_type<T,d-1>::Type> Type; 1697 }; 1698 template <class T> 1699 struct Gen_polynomial_type<T,0>{ typedef T Type; }; 1700 1701 public: 1702 template <class T, int d> 1703 struct Rebind{ 1704 typedef Polynomial_traits_d<typename Gen_polynomial_type<T,d>::Type> Other; 1705 }; 1706 //------------ Rebind ----------- 1707 }; 1708 1709 } //namespace CGAL 1710 1711 #include <CGAL/enable_warnings.h> 1712 1713 #endif // CGAL_POLYNOMIAL_TRAITS_D_H 1714