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/Curve_analysis_2.h $ 7 // $Id: Curve_analysis_2.h 4e519a3 2021-05-05T13:15:37+02:00 Sébastien Loriot 8 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Michael Kerber <mkerber@mpi-inf.mpg.de> 12 // 13 // ============================================================================ 14 15 #ifndef CGAL_ALGEBRAIC_CURVE_KERNEL_CURVE_ANALYSIS_2_ALCIX_H 16 #define CGAL_ALGEBRAIC_CURVE_KERNEL_CURVE_ANALYSIS_2_ALCIX_H 17 18 #include <CGAL/disable_warnings.h> 19 20 #include <vector> 21 #include <set> 22 #include <map> 23 24 #include <boost/mpl/has_xxx.hpp> 25 #include <boost/type_traits/is_base_of.hpp> 26 #include <boost/mpl/and.hpp> 27 #include <boost/mpl/logical.hpp> 28 #include <boost/type_traits/is_same.hpp> 29 30 31 #include <CGAL/basic.h> 32 #include <CGAL/assertions.h> 33 #include <CGAL/Cache.h> 34 #include <CGAL/function_objects.h> 35 #include <CGAL/Handle_with_policy.h> 36 37 #include <CGAL/Algebraic_kernel_d/Bitstream_descartes.h> 38 #include <CGAL/Algebraic_kernel_d/Interval_evaluate_1.h> 39 #include <CGAL/Algebraic_kernel_d/Bitstream_descartes_rndl_tree_traits.h> 40 #include <CGAL/Algebraic_kernel_d/macros.h> 41 #include <CGAL/Algebraic_kernel_d/exceptions.h> 42 #include <CGAL/Algebraic_kernel_d/enums.h> 43 #include <CGAL/Algebraic_kernel_d/algebraic_curve_kernel_2_tools.h> 44 #include <CGAL/Algebraic_kernel_d/Status_line_CA_1.h> 45 #include <CGAL/Algebraic_kernel_d/Event_line_builder.h> 46 #include <CGAL/Algebraic_kernel_d/Shear_controller.h> 47 #include <CGAL/Algebraic_kernel_d/Shear_transformation.h> 48 #include <CGAL/Algebraic_kernel_d/Bitstream_coefficient_kernel_at_alpha.h> 49 #include <CGAL/Algebraic_kernel_d/shear.h> 50 51 #include <CGAL/Polynomial_traits_d.h> 52 53 54 55 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 56 // put includes here 57 #endif 58 59 60 namespace CGAL { 61 62 template<typename AlgebraicKernelWithAnalysis_2, 63 typename Rep_> 64 class Curve_analysis_2; 65 66 namespace internal { 67 68 template<typename Comparable,bool has_template_typedefs> 69 struct Is_derived_from_Handle_with_policy { 70 typedef boost::false_type Tag; 71 }; 72 73 template<typename Comparable> 74 struct Is_derived_from_Handle_with_policy<Comparable,true> { 75 76 typedef typename 77 boost::is_base_of< CGAL::Handle_with_policy 78 < typename Comparable::T, 79 typename Comparable::Handle_policy, 80 typename Comparable::Allocator >, 81 Comparable 82 >::type Tag; 83 }; 84 85 86 template<typename Comparable,typename Tag> struct Compare_for_vert_line_map_ 87 { 88 bool operator() (const Comparable& a, const Comparable& b) const { 89 return a<b; 90 } 91 }; 92 93 template<typename Comparable> 94 struct Compare_for_vert_line_map_<Comparable,boost::true_type> { 95 96 bool operator() (const Comparable& a, const Comparable& b) const { 97 return CGAL::Handle_id_less_than< Comparable >()(a,b); 98 } 99 }; 100 101 template<typename Comparable> struct Compare_for_vert_line_map 102 : public CGAL::cpp98::binary_function<Comparable,Comparable,bool> { 103 104 BOOST_MPL_HAS_XXX_TRAIT_DEF(T) 105 BOOST_MPL_HAS_XXX_TRAIT_DEF(Handle_policy) 106 BOOST_MPL_HAS_XXX_TRAIT_DEF(Allocator) 107 108 typedef typename CGAL::internal::Is_derived_from_Handle_with_policy 109 < Comparable, 110 has_T<Comparable>::value && 111 has_Handle_policy<Comparable>::value && 112 has_Allocator<Comparable>::value>::Tag Tag; 113 114 public: 115 116 bool operator() (const Comparable& a, const Comparable& b) const { 117 118 return eval(a,b); 119 } 120 121 private: 122 123 Compare_for_vert_line_map_<Comparable,Tag> eval; 124 125 126 }; 127 128 129 // \brief Representation class for algebraic curves. 130 template< typename AlgebraicKernelWithAnalysis_2> 131 class Curve_analysis_2_rep { 132 133 public: 134 //! this instance's template parameter 135 typedef AlgebraicKernelWithAnalysis_2 Algebraic_kernel_with_analysis_2; 136 137 //! the class itself 138 typedef Curve_analysis_2_rep Self; 139 140 //! The handle class 141 typedef CGAL::Curve_analysis_2 142 <Algebraic_kernel_with_analysis_2,Self> Handle; 143 144 //protected: 145 public: 146 147 typedef int size_type; 148 149 CGAL_ACK_SNAP_ALGEBRAIC_CURVE_KERNEL_2_TYPEDEFS(Handle); 150 151 typedef std::map< Bound, Status_line_1 > 152 Vert_line_at_rational_map; 153 154 typedef 155 std::map< Algebraic_real_1, 156 Status_line_1, 157 internal::Compare_for_vert_line_map<Algebraic_real_1> > 158 Vert_line_map; 159 160 //!\name Constructors 161 //!@{ 162 163 //! Default constructor 164 Curve_analysis_2_rep() 165 { 166 } 167 168 //! Constructor with polynomial 169 Curve_analysis_2_rep(Algebraic_kernel_with_analysis_2 *kernel, 170 Polynomial_2 poly, 171 CGAL::Degeneracy_strategy strategy) : 172 _m_kernel(kernel), f(poly), degeneracy_strategy(strategy) 173 { 174 } 175 176 //!@} 177 178 private: 179 180 typedef internal::LRU_hashed_map< 181 Bound, 182 std::vector<Algebraic_real_1>, 183 internal::To_double_hasher > Intermediate_cache; 184 185 Intermediate_cache intermediate_cache; 186 187 typedef internal::Event_line_builder<Algebraic_kernel_with_analysis_2> 188 Event_line_builder; 189 190 191 // Internal information struct about x-coordinates 192 struct Event_coordinate_1 { 193 Event_coordinate_1(){} //added to solve a compilation error of gcc-3.4 (bug?) 194 Algebraic_real_1 val; 195 size_type mult_of_prim_res_root; 196 size_type index_of_prim_res_root; 197 size_type mult_of_content_root; 198 size_type index_of_content_root; 199 size_type mult_of_prim_lcoeff_root; 200 size_type index_of_prim_lcoeff_root; 201 boost::optional<Status_line_1> stack; 202 }; 203 204 // Functor to get the X_coordinate of an Event_coordinate 205 struct Val_functor { 206 typedef Event_coordinate_1 argument_type; 207 typedef Algebraic_real_1 result_type; 208 result_type operator() (argument_type event) const { 209 return event.val; 210 } 211 }; 212 213 214 //! The object holding the information about events, as an optional 215 mutable boost::optional<std::vector<Event_coordinate_1> > 216 event_coordinates; 217 218 //! The algebraic kernel to use 219 Algebraic_kernel_with_analysis_2* _m_kernel; 220 221 //! The polynomial defining the curve 222 boost::optional<Polynomial_2> f; 223 224 //! How degenerate situations are handled 225 CGAL::Degeneracy_strategy degeneracy_strategy; 226 227 /*! 228 * \brief The polynomial without its content (the gcd of the coeffs). 229 * 230 * The content is the greatest common divisor of the coefficients of \c f 231 * considered as polynomial <tt>y</tt>. \c The polynomial f_primitive is 232 * \c f/cont(f). The corresponding curve is equal to the curve of \c f, 233 * only without vertical line components. 234 */ 235 mutable boost::optional<Polynomial_2> f_primitive; 236 237 //! the polynomial containing all roots of the resultant of the primitive 238 //! part of f and its y-derivative 239 mutable boost::optional<Polynomial_1> 240 resultant_of_primitive_and_derivative_y; 241 242 //! the polynomial containing all roots of the resultant of the primitive 243 //! part of f and its x-derivative 244 mutable boost::optional<Polynomial_1> 245 resultant_of_primitive_and_derivative_x; 246 247 //! The Sturm-Habicht polynomials of f 248 mutable boost::optional<std::vector<Polynomial_2> > 249 sturm_habicht_of_primitive; 250 251 //! The content of f 252 mutable boost::optional<Polynomial_1> content; 253 254 //! The non-working shear factors, as far as known 255 mutable std::set<Integer> bad_shears; 256 257 //! The already known shear factors 258 mutable std::map<Integer,Handle> sheared_curves; 259 260 //! Has the curve vertical line components 261 mutable boost::optional<bool> has_vertical_component; 262 263 //! The intermediate values 264 mutable boost::optional<std::vector<boost::optional<Bound> > > 265 intermediate_values; 266 267 //! stores Y_values at rational coordinate 268 mutable Vert_line_at_rational_map vert_line_at_rational_map; 269 270 //! stores vert_lines 271 mutable Vert_line_map vert_line_map; 272 273 /**! \brief Information about whether arcs at +/- infty 274 * are asymptotic to y=beta, 275 * or go to +/- infty also in y-direction 276 */ 277 mutable boost::optional<std::vector<CGAL::Object> > 278 horizontal_asymptotes_left, horizontal_asymptotes_right; 279 280 //! friends 281 friend class ::CGAL::Curve_analysis_2 282 <Algebraic_kernel_with_analysis_2,Self>; 283 284 }; // class Curve_analysis_2_rep 285 } // namespace internal 286 287 288 /*! 289 * \brief Analysis for algebraic curves of arbitrary degree. 290 * 291 * This class constitutes a model for the concept 292 * AlgebraicKernelWithAnalysis_d_2::CurveAnalysis_2. 293 * For a square-free bivariate polynomial \c f, a topologic-geometrical 294 * analysis of the algebraic curve defined by the vanishing set of \c f 295 * is provided. This means, one can ask for the total number, and the position 296 * of the critical x-coordinates of the curve, and for each x-coordinate, 297 * geometric information about the curve can be obtained. This data 298 * is capsuled into an object of type \c Curve_analysis_2::Status_line_1, 299 * which is in fact a \c Status_line_CA_1 object. 300 * 301 * The restriction to square-free curves is a weak one, since the curves 302 * can be made square-free before passed to the analysis. 303 * The \c Construct_curve_2 functor of \c Algebraic_curve_kernel_2 is 304 * doing so, thus it accepts arbitrary bivariate polynomials. 305 * 306 * The analysis is implemented in a "lazy" fashion. This means, when 307 * created, the analysis delays all computations until the information 308 * is queried for the first time. This means, if only parts of the curves 309 * are of interest, only a partial analysis is performed. 310 * We remark that nevertheless, the global \e projection \e step 311 * (i.e., computing the (sub)resultants) must be done once a \c Status_line_1 312 * is queried. Often, this step forms the bottleneck in the whole computation. 313 * 314 * For more details of the algorithm, consult the reference: 315 * A.Eigenwillig, M.Kerber, N.Wolpert: Fast and Exact Geometric Analysis of 316 * Real Algebraic Plane Curves. Proceedings of the International Symposium 317 * on Symbolic and Algebraic Computation (ISSAC 2007), pp. 151-158 318 */ 319 template<typename AlgebraicKernelWithAnalysis_2, 320 typename Rep_ 321 = internal::Curve_analysis_2_rep< AlgebraicKernelWithAnalysis_2> 322 > 323 class Curve_analysis_2 : public ::CGAL::Handle_with_policy< Rep_ > { 324 325 //! \name typedefs 326 //! @{ 327 328 public: 329 //! this instance' first template parameter 330 typedef AlgebraicKernelWithAnalysis_2 Algebraic_kernel_with_analysis_2; 331 332 //! this instance' second template parameter 333 typedef Rep_ Rep; 334 335 private: 336 337 //! The internal type for event coordinates 338 typedef typename Rep::Event_coordinate_1 Event_coordinate_1; 339 340 // Internal class to build lines at events 341 typedef typename Rep::Event_line_builder Event_line_builder; 342 343 // Base class 344 typedef ::CGAL::Handle_with_policy<Rep> Base; 345 346 // This type 347 typedef CGAL::Curve_analysis_2<Algebraic_kernel_with_analysis_2,Rep> Self; 348 349 public: 350 351 //! Indexing type 352 typedef typename Rep::size_type size_type; 353 354 CGAL_ACK_SNAP_ALGEBRAIC_CURVE_KERNEL_2_TYPEDEFS(Self); 355 356 //! Required by the CurveKernel_2 concept 357 typedef Algebraic_real_1 Coordinate_1; 358 359 //! Traits type for Polynomial_2 360 typedef CGAL::Polynomial_traits_d<Polynomial_2> Polynomial_traits_2; 361 362 private: 363 364 /*! 365 * \brief Coercion between the coefficient type of the polynomial 366 * and the bound type of the curve analysis 367 * 368 * Interoperability of both types is required 369 */ 370 typedef CGAL::Coercion_traits<Bound, Coefficient> Coercion; 371 372 /*! 373 * \brief The common supertype that both the coefficient and the bound 374 * type are convertible to 375 */ 376 typedef typename Coercion::Type Coercion_type; 377 378 //! Polynomial over the \c Coercion_type 379 typedef typename CGAL::Polynomial_traits_d<Polynomial_2> 380 ::template Rebind<Coercion_type,1>::Other::Type Poly_coer_1; 381 382 public: 383 384 //! Type to represent points on curves 385 typedef typename Algebraic_kernel_with_analysis_2::Algebraic_real_2 386 Algebraic_real_2; 387 388 //! Required by the CurveKernel_2 concept 389 typedef Algebraic_real_2 Coordinate_2; 390 391 //! type for horizontal asymtote values 392 typedef CGAL::Object Asymptote_y; 393 394 //! @} 395 396 private: 397 398 //! \name Helping structs 399 // @{ 400 401 struct Event_functor { 402 Event_functor(const Self* curve) : curve(curve) {} 403 const Self* curve; 404 typedef size_type argument_type; 405 typedef Status_line_1 result_type; 406 result_type operator() (argument_type index) const { 407 return curve->status_line_at_event(index); 408 } 409 }; 410 411 struct Intermediate_functor { 412 Intermediate_functor(const Self* curve) : curve(curve) {} 413 const Self* curve; 414 typedef size_type argument_type; 415 typedef Status_line_1 result_type; 416 result_type operator() (argument_type index) const { 417 return curve->status_line_of_interval(index); 418 } 419 }; 420 421 struct Stha_functor { 422 Stha_functor(const Self* curve) : curve(curve) {} 423 const Self* curve; 424 typedef size_type argument_type; 425 typedef Polynomial_1 result_type; 426 result_type operator() (argument_type index) const { 427 return curve->principal_sturm_habicht_of_primitive(index); 428 } 429 }; 430 431 //! @} 432 433 public: 434 435 //! \name Iterators 436 //! @{ 437 438 //! Iterator type for status lines at events 439 typedef boost::transform_iterator<Event_functor, 440 boost::counting_iterator<size_type> > 441 Event_line_iterator; 442 443 //! Iterator type for status lines of intervals 444 typedef boost::transform_iterator<Intermediate_functor, 445 boost::counting_iterator<size_type> > 446 Intermediate_line_iterator; 447 448 //! Iterator type for the principal sturm habicht coefficients of the curve 449 typedef boost::transform_iterator<Stha_functor, 450 boost::counting_iterator<size_type> > 451 Principal_sturm_habicht_iterator; 452 453 //! @} 454 455 public: 456 457 //!\name Constructors 458 //!@{ 459 460 //! Default constructor, constructs an empty and invalid curve analysis 461 Curve_analysis_2() :Base(Rep()) { 462 } 463 464 /*! 465 * \brief Constructs the curve analysis for the given polynomial 466 * 467 * Analyses the curve that is defined by the vanishing set of the 468 * polynomial \c f. 469 * \pre \c f is square free. 470 * \param strategy The default strategy 471 * (\c SHEAR_ONLY_AT_IRRATIONAL_STRATEGY) 472 * is to \c shear the curve 473 * if a degenerate situation is detected during the analysis, 474 * except at rational x-coordinates where the curve can be analysed 475 * more directly. The analysis 476 * is then performed in the sheared system, and finally translated back 477 * into the original system. 478 * Using \c SHEAR_STRATEGY, a shear is triggered also for degeneracies 479 * at rational x-coordinate. With both strategies, it is guaranteed that 480 * the analysis works successfully for any square free input curve. 481 * On the other hand, the EXCEPTION_STRATEGY throws an exception of type 482 * \c internal::Zero_resultant_exception<Polynomial_2>, 483 * instead of performing a shear. 484 * 485 * \Todo Currently the defualt strategy has been changed to SHEAR_STRATEGY 486 * because there exist a problem if vertical asymtotes are present at 487 * the rational x-coordinate. 488 */ 489 explicit Curve_analysis_2(Algebraic_kernel_with_analysis_2 *kernel, 490 const Polynomial_2& f, 491 CGAL::Degeneracy_strategy strategy 492 = CGAL_ACK_DEFAULT_DEGENERACY_STRATEGY) 493 : Base(Rep(kernel,f,strategy)) 494 { 495 496 } 497 498 //! \brief Copy constructor 499 #ifdef DOXYGEN_RUNNING 500 Curve_analysis_2(const Self& alg_curve) 501 : Base(static_cast<const Base&>(alg_curve)) 502 { 503 } 504 #endif 505 506 //!@} 507 508 509 //! \name Members 510 //!@{ 511 512 private: 513 514 /* 515 * \brief sets all status lines at events and of intervals 516 * 517 * Writes the status lines of events and interval into the object. 518 * The value type of both \c InputIterator1 and \c InputIterator2 519 * is \c Status_line_1. 520 */ 521 template<typename InputIterator1,typename InputIterator2> 522 void set_event_lines(InputIterator1 event_begin, 523 InputIterator1 event_end, 524 InputIterator2 intermediate_begin, 525 InputIterator2 CGAL_precondition_code(intermediate_end)) const { 526 527 if(! this->ptr()->event_coordinates) { 528 529 std::vector<Event_coordinate_1> event_coordinate_vector; 530 531 for(InputIterator1 it = event_begin; it != event_end; it++) { 532 Event_coordinate_1 curr_event; 533 curr_event.val = it->x(); 534 event_coordinate_vector.push_back(curr_event); 535 } 536 this->ptr()->event_coordinates = event_coordinate_vector; 537 538 } 539 540 InputIterator1 it1 = event_begin; 541 for(size_type i = 0; i < number_of_status_lines_with_event() ; i++ ) { 542 this->ptr()->vert_line_map[event_coordinates()[i].val] = *it1; 543 event_coordinates()[i].stack = *it1; 544 545 it1++; 546 } 547 CGAL_assertion(it1 == event_end); 548 549 if(! this->ptr()->intermediate_values) { 550 this->ptr()->intermediate_values 551 = std::vector<boost::optional<Bound> > 552 (number_of_status_lines_with_event()+1); 553 } 554 555 InputIterator2 it2 = intermediate_begin; 556 for(size_type i = 0; 557 i < static_cast<int>(intermediate_values().size()); 558 i++,it2++) { 559 560 CGAL_assertion(it2->x().is_rational()); 561 Bound q = it2->x().rational(); 562 563 intermediate_values()[i] = q; 564 this->ptr()->vert_line_map[it2->x()] = *it2; 565 this->ptr()->vert_line_at_rational_map[q] = *it2; 566 567 } 568 CGAL_assertion(it2 == intermediate_end); 569 570 } 571 572 public: 573 574 /*! \brief returns whether the curve has a valid defining polynomial 575 */ 576 bool has_defining_polynomial() const { 577 return bool(this->ptr()->f); 578 } 579 580 public: 581 582 /*! \brief sets the defining polynomial. 583 * 584 * \pre The object has no defining polynomial yet. 585 */ 586 void set_f(Polynomial_2 f) { 587 CGAL_precondition(! has_defining_polynomial()); 588 if((! this->ptr()->f) || f!=this->ptr()->f.get()) { 589 this->copy_on_write(); 590 this->ptr()->f=f; 591 } 592 } 593 594 595 public: 596 597 /*! 598 * \brief returns whether the curve is y-regular 599 * 600 * A curve is called y-regular if the leading coefficient of its defining 601 * polynomial wrt y is a constant, i.e., contains no x 602 */ 603 bool is_y_regular() const { 604 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 605 if(CGAL::degree(polynomial_2(),1)==2) { 606 return this->conic_is_y_regular(); 607 } 608 #endif 609 return CGAL::degree(CGAL::leading_coefficient(polynomial_2())) == 0; 610 } 611 612 public: 613 614 /*! 615 * \brief returns whether the curve contains a vertical line as a component 616 * 617 * In algebraic terms, this methods computes whether the content 618 * of its defining polynomial has a real root. 619 */ 620 bool has_vertical_component() const { 621 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 622 if(CGAL::degree(polynomial_2(),1)==2) { 623 return this->conic_has_vertical_components(); 624 } 625 #endif 626 if(is_y_regular()) { 627 this->ptr()->has_vertical_component = false; 628 } 629 if(! this->ptr()->has_vertical_component) { 630 // This is computed as side effect 631 // when the event coordinates are computed 632 event_coordinates(); 633 CGAL_assertion(this->ptr()->has_vertical_component); 634 } 635 return this->ptr()->has_vertical_component.get(); 636 } 637 638 public: 639 640 //! Returns the defining polynomial 641 Polynomial_2 polynomial_2() const { 642 CGAL_precondition(bool(this->ptr()->f)); 643 return this->ptr()->f.get(); 644 } 645 646 public: 647 648 /*! 649 * \brief returns the number of event lines of the curve 650 * 651 * Algebraically, the number of real roots of the discriminant of 652 * the curve's defining equation is returned. 653 */ 654 size_type number_of_status_lines_with_event() const { 655 CGAL_precondition(bool(this->ptr()->f)); 656 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 657 if(CGAL::degree(polynomial_2(),1)==2) { 658 return this->conic_number_of_status_lines_with_event(); 659 } 660 #endif 661 return static_cast<size_type>(event_coordinates().size()); 662 } 663 664 public: 665 666 /*! 667 * \brief returns whether the given x-coordinate is critical for the curve 668 * and which event or interval index the x-coordinate belongs to. 669 * 670 * \param is_event is set to \c true if the curve has an event 671 * at this x-coordinate, or in other words, if the discriminant of its 672 * defining polynomial vanishes at \c x 673 * \param i is set to the index of the event if \c x is an event. Otherwise 674 * \c i is set to the index of the interval \c x is contained in. 675 */ 676 void x_to_index(Algebraic_real_1 x,size_type& i,bool& is_event) const { 677 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 678 if(CGAL::degree(polynomial_2(),1)==2) { 679 return this->conic_x_to_index(x,i,is_event); 680 } 681 #endif 682 CGAL_precondition(has_defining_polynomial()); 683 typename Rep::Val_functor xval; 684 i = static_cast<size_type>(std::lower_bound( 685 ::boost::make_transform_iterator(event_coordinates().begin(), 686 xval), 687 ::boost::make_transform_iterator(event_coordinates().end(), 688 xval), 689 x 690 ) - ::boost::make_transform_iterator(event_coordinates().begin(), 691 xval)); 692 is_event = (i < static_cast<size_type>(event_coordinates().size()) && 693 (event_coordinates()[i].val == x) ); 694 } 695 696 public: 697 698 //! Returns the status line at the <tt>i</tt>-th event of the curve. 699 Status_line_1& status_line_at_event(size_type i) const { 700 701 CGAL_precondition(has_defining_polynomial()); 702 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 703 if(CGAL::degree(polynomial_2(),1)==2) { 704 return this->conic_status_line_at_event(i); 705 } 706 #endif 707 CGAL_precondition_code( 708 size_type n = 709 static_cast<size_type>(event_coordinates().size()); 710 ); 711 CGAL_precondition(i>=0 && i<n); 712 if(! event_coordinates()[i].stack) { 713 Status_line_1 event_line = create_status_line_at_event(i); 714 this->ptr()->vert_line_map[event_coordinates()[i].val] 715 = event_line; 716 event_coordinates()[i].stack = event_line; 717 } 718 CGAL_postcondition(event_coordinates()[i].stack.get().is_event()); 719 return event_coordinates()[i].stack.get(); 720 } 721 722 public: 723 724 //! Returns a status line at the rational <tt>x</tt>-coordinate \c b 725 Status_line_1& status_line_at_exact_x(Bound b) const { 726 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 727 if(CGAL::degree(polynomial_2(),1)==2) { 728 return this->conic_status_line_at_exact_x(b); 729 } 730 #endif 731 return status_line_at_exact_x(Algebraic_real_1(b)); 732 } 733 734 private: 735 736 /* 737 * \brief returns a status line for an exact value \c alpha that 738 * is not an event of the curve 739 * 740 * This function controls the internal cache that stores already created 741 * status line at non-events. 742 */ 743 Status_line_1& status_line_at_exact_non_event_x(Algebraic_real_1 alpha) 744 const { 745 746 if(alpha.is_rational()) { 747 748 typename Rep::Vert_line_at_rational_map::iterator it = 749 this->ptr()->vert_line_at_rational_map.find 750 (alpha.rational()); 751 752 if (it != this->ptr()->vert_line_at_rational_map.end()) { 753 CGAL_assertion(!it->second.is_event()); 754 return it->second; 755 } 756 } 757 758 typename Rep::Vert_line_map::iterator it = 759 this->ptr()->vert_line_map.find(alpha); 760 761 if (it != this->ptr()->vert_line_map.end()) { 762 CGAL_assertion(!it->second.is_event()); 763 return it->second; 764 } 765 766 767 // Not stored yet, so create it and store it 768 Status_line_1 cvl 769 = create_status_line_at_non_event(alpha); 770 CGAL_assertion(!cvl.is_event()); 771 this->ptr()->vert_line_map[alpha] = cvl; 772 773 if(alpha.is_rational()) { 774 this->ptr()->vert_line_at_rational_map[alpha.rational()] = cvl; 775 } 776 return this->ptr()->vert_line_map[alpha]; 777 } 778 779 public: 780 781 //! Returns a vert line for the <tt>x</tt>-coordinate alpha 782 Status_line_1& status_line_at_exact_x(Algebraic_real_1 alpha) const { 783 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 784 if(CGAL::degree(polynomial_2(),1)==2) { 785 return this->conic_status_line_at_exact_x(alpha); 786 } 787 #endif 788 bool is_event_value; 789 size_type index; 790 this->x_to_index(alpha,index,is_event_value); 791 if(is_event_value) { 792 return status_line_at_event(index); 793 } 794 else { 795 return status_line_at_exact_non_event_x(alpha); 796 } 797 } 798 799 private: 800 801 // Creates a status line for the curve's <tt>index</tt>th critical point 802 Status_line_1 create_status_line_at_event(size_type index) const 803 { 804 805 Event_coordinate_1& event = event_coordinates()[index]; 806 807 Algebraic_real_1 x = event.val; 808 809 try { 810 811 Event_coordinate_1& event = event_coordinates()[index]; 812 813 Algebraic_real_1 x = event.val; 814 815 #if CGAL_ACK_SHEAR_ALL_NOT_Y_REGULAR_CURVES 816 if(event.mult_of_prim_lcoeff_root > 0) { 817 throw CGAL::internal::Non_generic_position_exception(); 818 } 819 #else 820 if(event.mult_of_prim_lcoeff_root > 0) { 821 if(event.mult_of_prim_lcoeff_root > 1 || 822 event.mult_of_prim_res_root > 1) { 823 throw CGAL::internal::Non_generic_position_exception(); 824 } 825 } 826 827 #endif 828 829 #if CGAL_ACK_DEBUG_FLAG 830 double ev_approx = CGAL::to_double(x); 831 CGAL_ACK_DEBUG_PRINT << (index+1) << "th line: " 832 << std::setw(6) << std::setprecision(3) 833 << ev_approx 834 << ".." 835 << std::flush; 836 #endif 837 size_type left_arcs 838 = status_line_for_x(x,CGAL::NEGATIVE).number_of_events(); 839 size_type right_arcs 840 = status_line_for_x(x,CGAL::POSITIVE).number_of_events(); 841 842 bool root_of_resultant=(event.mult_of_prim_res_root>0); 843 bool root_of_content=(event.mult_of_content_root>0); 844 845 size_type mult_of_resultant = event.mult_of_prim_res_root; 846 847 /* 848 #if CGAL_ACK_DEBUG_FLAG 849 CGAL_ACK_DEBUG_PRINT << "Event line for " << index << " " 850 << root_of_resultant << " " 851 << root_of_content << " " 852 << mult_of_resultant << " " 853 << left_arcs << " " << right_arcs 854 << std::endl; 855 #endif 856 */ 857 858 Status_line_1 ev_line 859 = event_line_builder().create_event_line(index, 860 x, 861 left_arcs, 862 right_arcs, 863 root_of_resultant, 864 root_of_content, 865 mult_of_resultant); 866 867 event.stack = ev_line; 868 869 #if CGAL_ACK_DEBUG_FLAG 870 CGAL_ACK_DEBUG_PRINT << "done" << std::endl; 871 #endif 872 873 return ev_line; 874 } catch(CGAL::internal::Non_generic_position_exception /* exc */) { 875 switch(this->ptr()->degeneracy_strategy) { 876 case(CGAL::EXCEPTION_STRATEGY): { 877 throw CGAL::internal::Non_generic_position_exception(); 878 break; 879 } 880 // Feature does not working atm 881 case(CGAL::SHEAR_ONLY_AT_IRRATIONAL_STRATEGY): { 882 CGAL_error_msg("Currently not supported"); 883 /* 884 if(x.is_rational()) { 885 return create_non_generic_event_at_rational(x,index); 886 } 887 // FALL INTO NEXT CASE 888 */ 889 } 890 case(CGAL::SHEAR_STRATEGY): { 891 return create_non_generic_event_with_shear(index); 892 break; 893 } 894 default:{ 895 CGAL_assertion(false); // !!! Never reached 896 } 897 } 898 } 899 // !!! Never reached 900 return Status_line_1(); 901 } 902 903 private: 904 905 /* 906 * \brief Method to create a status line using shear and backshear 907 * 908 * Note that this methods creates <b>all</b> event lines of the object 909 * at once, and stores them in the object. 910 */ 911 Status_line_1 create_non_generic_event_with_shear(size_type index) const { 912 913 #if CGAL_ACK_DEBUG_FLAG 914 CGAL_ACK_DEBUG_PRINT << "Use sheared technique..." << std::endl; 915 #endif 916 internal::Shear_controller<Integer> shear_controller; 917 Integer s(0); 918 while(true) { 919 try { 920 s = shear_controller.get_shear_factor(); 921 #if CGAL_ACK_DEBUG_FLAG 922 CGAL_ACK_DEBUG_PRINT << "Trying shear factor " 923 << s << std::endl; 924 #endif 925 // TODO: Move shear somewhere else 926 Self D(kernel(), 927 CGAL::internal::shear 928 (primitive_polynomial_2(),Coefficient(s)), 929 CGAL::EXCEPTION_STRATEGY); 930 Shear_transformation< Algebraic_kernel_with_analysis_2 > 931 shear_transformation(kernel()); 932 shear_transformation.report_sheared_disc_roots 933 (boost::make_transform_iterator( 934 event_coordinates().begin(), 935 typename Rep::Val_functor()), 936 boost::make_transform_iterator( 937 event_coordinates().end(), 938 typename Rep::Val_functor()) 939 ); 940 941 // Store the sheared curve for later use 942 this->ptr()->sheared_curves.insert(std::make_pair(s,D)); 943 shear_transformation(D,-s,(Self&)*this,false); 944 set_vertical_line_components(); 945 946 break; 947 } 948 catch(CGAL::internal::Non_generic_position_exception /* err */) { 949 950 shear_controller.report_failure(s); 951 #if CGAL_ACK_DEBUG_FLAG 952 CGAL_ACK_DEBUG_PRINT << "Bad shear factor, retrying..." 953 << std::endl; 954 #endif 955 } 956 } 957 958 return status_line_at_event(index); 959 } 960 961 private: 962 963 /* 964 * \brief creates a status line for a rational event x-coordinate 965 * 966 * If an event coordinate is rational, a shear can be prevented 967 * by plugging in the x-coordinate for x and explicitly computing 968 * the square free part of the defining polynomial at this position. 969 * 970 * COMMENTED OUT 971 972 Status_line_1 create_non_generic_event_at_rational(Algebraic_real_1 x, 973 size_type index) const { 974 975 976 #if CGAL_ACK_DEBUG_FLAG 977 CGAL_ACK_DEBUG_PRINT << "Non-generic, rational position x = " 978 << CGAL::to_double(x) 979 << std::flush; 980 #endif 981 982 CGAL_precondition(x.is_rational()); 983 Bound r = x.rational(); 984 985 Polynomial_1 f_at_x = kernel()->evaluate_utcf_2_object() 986 (typename Polynomial_traits_2::Swap() 987 (primitive_polynomial_2(),0, 1), 988 r); 989 990 f_at_x_sq_free 991 = typename CGAL::Polynomial_traits_d<typename FT::Numerator_type> 992 ::Make_square_free() (f_at_x); 993 994 Bitstream_coefficient_kernel coeff_kernel(kernel(),x); 995 Bitstream_traits traits(coeff_kernel); 996 997 // We need to make an artificial bivariate polynomial 998 typedef typename 999 CGAL::Polynomial_traits_d<typename FT::Numerator_type> 1000 ::template Rebind<typename FT::Numerator_type,1>::Other::Type 1001 Poly_coer_num_2; 1002 1003 std::vector<typename FT::Numerator_type> coeffs; 1004 for(int i = 0; i <= CGAL::degree(f_at_x_sq_free); i++) { 1005 coeffs.push_back(typename FT::Numerator_type(f_at_x_sq_free[i])); 1006 } 1007 Poly_coer_num_2 f_at_x_ext(coeffs.begin(), coeffs.end()); 1008 1009 Bitstream_descartes isolator(CGAL::internal::Square_free_descartes_tag(), 1010 f_at_x_ext, 1011 traits); 1012 1013 // Now adjacencies 1014 std::vector<Bound> bucket_borders; 1015 1016 int n = isolator.number_of_real_roots(); 1017 1018 if(n==0) { 1019 bucket_borders.push_back(0); 1020 } else { 1021 bucket_borders.push_back( 1022 CGAL::internal::bound_left_of 1023 (kernel(),Algebraic_real_1(isolator.left_bound(0)))); 1024 for(int i = 1; i < n; i++) { 1025 while(Algebraic_real_1(isolator.right_bound(i-1))== 1026 Algebraic_real_1(isolator.left_bound(i))) { 1027 isolator.refine_interval(i-1); 1028 isolator.refine_interval(i); 1029 } 1030 bucket_borders.push_back( 1031 kernel()->bound_between_1_object() 1032 (Algebraic_real_1(isolator.right_bound(i-1)), 1033 Algebraic_real_1(isolator.left_bound(i))) 1034 ); 1035 } 1036 1037 bucket_borders.push_back( 1038 CGAL::internal::bound_right_of 1039 (kernel(), 1040 Algebraic_real_1(isolator.right_bound(n-1)))); 1041 } 1042 1043 Bound left = bound_value_in_interval(index); 1044 Bound right = bound_value_in_interval(index+1); 1045 1046 typedef boost::numeric::interval<Coercion_type> Coercion_interval; 1047 1048 typename Coercion::Cast cast; 1049 1050 for(int i = 0; i < static_cast<int>(bucket_borders.size()); i++) { 1051 1052 Poly_coer_1 curr_pol 1053 = primitive_polynomial_2().evaluate(bucket_borders[i]); 1054 1055 CGAL::internal::Interval_evaluate_1 1056 <Poly_coer_1,Bound> 1057 interval_evaluate_1; 1058 1059 while(true) { 1060 std::pair<Bound,Bound> curr_interval_pair 1061 = interval_evaluate_1(curr_pol,std::make_pair(left,right)); 1062 Coercion_interval curr_interval(curr_interval_pair.first, 1063 curr_interval_pair.second); 1064 1065 if(boost::numeric::in_zero(curr_interval)) { 1066 // "refine" 1067 Bound middle = (left+right)/2; 1068 if(middle==r) { 1069 left=(left+middle)/2; 1070 right = (right+middle)/2; 1071 } else if(middle>r) { 1072 right=middle; 1073 } else { 1074 left=middle; 1075 } 1076 } else { 1077 break; 1078 } 1079 } 1080 } 1081 1082 Status_line_1 left_line 1083 = status_line_at_exact_non_event_x(Algebraic_real_1(left)), 1084 right_line 1085 = status_line_at_exact_non_event_x(Algebraic_real_1(right)); 1086 1087 int n_left = left_line.number_of_events(); 1088 int n_right = right_line.number_of_events(); 1089 1090 std::vector<int> left_arcs(bucket_borders.size()+1), 1091 right_arcs(bucket_borders.size()+1); 1092 1093 for(unsigned int i=0;i<left_arcs.size();i++) { 1094 left_arcs[i]=0; 1095 } 1096 for(unsigned int i=0;i<right_arcs.size();i++) { 1097 right_arcs[i]=0; 1098 } 1099 1100 int curr_index=0; 1101 for(int i=0; i < n_left; i++) { 1102 1103 while(true) { 1104 if(curr_index==static_cast<int>(bucket_borders.size())) { 1105 left_arcs[curr_index]++; 1106 break; 1107 } else if(left_line.lower_bound(i)> 1108 bucket_borders[curr_index]) { 1109 curr_index++; 1110 } else if(left_line.upper_bound(i)< 1111 bucket_borders[curr_index]) { 1112 left_arcs[curr_index]++; 1113 break; 1114 } else { 1115 left_line.refine(i); 1116 } 1117 } 1118 } 1119 curr_index=0; 1120 for(int i=0; i < n_right; i++) { 1121 1122 while(true) { 1123 if(curr_index==static_cast<int>(bucket_borders.size())) { 1124 right_arcs[curr_index]++; 1125 break; 1126 } else if(right_line.lower_bound(i)> 1127 bucket_borders[curr_index]) { 1128 curr_index++; 1129 } else if(right_line.upper_bound(i)< 1130 bucket_borders[curr_index]) { 1131 right_arcs[curr_index]++; 1132 break; 1133 } else { 1134 right_line.refine(i); 1135 } 1136 } 1137 1138 } 1139 1140 typename Status_line_1::Arc_container arc_container; 1141 1142 for(int i = 0; i < n; i++) { 1143 arc_container.push_back(std::make_pair(left_arcs[i+1], 1144 right_arcs[i+1])); 1145 } 1146 1147 Status_line_1 status_line(x,index,*this,n_left,n_right,arc_container); 1148 1149 status_line._set_number_of_branches_approaching_infinity 1150 (std::make_pair(left_arcs[0],right_arcs[0]), 1151 std::make_pair(left_arcs[n+1],right_arcs[n+1])); 1152 1153 status_line.set_isolator(isolator); 1154 1155 if(event_coordinates()[index].mult_of_content_root > 0) { 1156 status_line._set_v_line(); 1157 } 1158 1159 #if CGAL_ACK_DEBUG_FLAG 1160 CGAL_ACK_DEBUG_PRINT << "done" << std::endl; 1161 #endif 1162 1163 return status_line; 1164 } 1165 */ 1166 1167 public: 1168 1169 /*! 1170 * \brief returns the status line for the interval 1171 * preceeding the <tt>i</tt>th event 1172 * 1173 * Returns a status line for a reference x-coordinate of the <tt>i</tt>th 1174 * interval of the curve. If called multiple times for the same <tt>i</tt>, 1175 * the same status line is returned. 1176 */ 1177 Status_line_1 status_line_of_interval(size_type i) const 1178 { 1179 CGAL_precondition(i >= 0 && i <= number_of_status_lines_with_event()); 1180 1181 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 1182 if(CGAL::degree(polynomial_2(),1)==2) { 1183 return this->conic_status_line_of_interval(i); 1184 } 1185 #endif 1186 1187 Bound b = bound_value_in_interval(i); 1188 1189 Status_line_1 intermediate_line 1190 = status_line_at_exact_non_event_x(Algebraic_real_1(b)); 1191 1192 CGAL_postcondition(! intermediate_line.is_event()); 1193 1194 return intermediate_line; 1195 } 1196 1197 1198 public: 1199 1200 /*! 1201 * \brief returns a status line at position \c x 1202 * 1203 * If \c x is not an event of the curve, and lies in the <tt>i</tt>th 1204 * interval, the result is equal to <tt>status_line_of_interval(i)</tt>. 1205 * Different from <tt>status_line_at_exact_x(x)</tt> 1206 * the status line \c s returned does not satisft <tt>s.x()==x</tt>. 1207 * If \c x is an event, and \c perturb is set to \c CGAL::ZERO, 1208 * the status line for the event is returned. Otherwise, the status line 1209 * for the left or right neighboring interval is returned, depending 1210 * on whether \c perturb is set to \c CGAL::NEGATIVE or \c CGAL::POSITIVE. 1211 * If \c x is not an event, \c perturb has no effect. 1212 */ 1213 Status_line_1 status_line_for_x(Algebraic_real_1 x, 1214 CGAL::Sign perturb = CGAL::ZERO) const 1215 { 1216 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 1217 if(CGAL::degree(polynomial_2(),1)==2) { 1218 return this->conic_status_line_for_x(x,perturb); 1219 } 1220 #endif 1221 1222 size_type i; 1223 bool is_evt; 1224 x_to_index(x, i, is_evt); 1225 if(is_evt) { 1226 if(perturb == CGAL::ZERO) 1227 return status_line_at_event(i); 1228 if(perturb == CGAL::POSITIVE) 1229 i++; 1230 } 1231 return status_line_of_interval(i); 1232 } 1233 1234 1235 private: 1236 1237 /* 1238 * \brief creates an intermediate line at position \c ar. 1239 * 1240 * It is required that none of the following situations occurs at position 1241 * <tt>ar</tt>: singularity, vertical tangent line, vertical asymptote.\n 1242 * Otherwise, the method might run into an infinite loop. 1243 * 1244 * \param index if set to -1, the interval containing \c ar is computed 1245 * within the method, and the index of the status line is set accordingly. 1246 */ 1247 Status_line_1 1248 create_status_line_at_non_event(Algebraic_real_1 ar, int index = -1) 1249 const { 1250 1251 if(index==-1) { 1252 bool event; 1253 x_to_index(ar,index,event); 1254 CGAL_assertion(!event); 1255 } 1256 CGAL_assertion(index>=0); 1257 1258 // TODO .. delay creation of refinement object 1259 // especially when ar is rational 1260 1261 Bitstream_coefficient_kernel coeff_kernel(kernel(),ar); 1262 Bitstream_traits traits(coeff_kernel); 1263 1264 Bitstream_descartes 1265 bitstream_descartes(CGAL::internal::Square_free_descartes_tag(), 1266 primitive_polynomial_2(), 1267 traits); 1268 1269 size_type root_number=bitstream_descartes.number_of_real_roots(); 1270 1271 Status_line_1 status_line(ar, index, *this, root_number); 1272 status_line.set_isolator(bitstream_descartes); 1273 1274 CGAL_assertion(! status_line.is_event()); 1275 1276 return status_line; 1277 } 1278 1279 private: 1280 1281 /* 1282 * \brief returns an Event_line_builder instance 1283 * 1284 * Note: So far, a new instance is created each time the function is called 1285 */ 1286 Event_line_builder event_line_builder() const { 1287 1288 return Event_line_builder(kernel(), *this, primitive_polynomial_2()); 1289 } 1290 1291 public: 1292 1293 /*! 1294 * \brief Number of arcs over the given interval 1295 * 1296 * Shortcut for <tt>status_line_of_interval(i).number_of_events()</tt> 1297 */ 1298 size_type arcs_over_interval(size_type i) const { 1299 CGAL_precondition(has_defining_polynomial()); 1300 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 1301 if(CGAL::degree(polynomial_2(),1)==2) { 1302 return this->conic_arcs_over_interval(i); 1303 } 1304 #endif 1305 CGAL_assertion_code( 1306 size_type n 1307 = static_cast<size_type>(intermediate_values().size()); 1308 ); 1309 CGAL_precondition(i>=0 && i<=n); 1310 return status_line_of_interval(i).number_of_events(); 1311 } 1312 1313 public: 1314 1315 /*! 1316 * \brief Rational number in the <tt>i</tt>th interval between events 1317 * 1318 * The result of this method is taken as the reference x-coordinate 1319 * for the status lines of intervals. 1320 */ 1321 Bound bound_value_in_interval(size_type i) const { 1322 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 1323 if(CGAL::degree(polynomial_2(),1)==2) { 1324 return this->conic_bound_value_in_interval(i); 1325 } 1326 #endif 1327 CGAL_assertion(i>=0 && 1328 i < static_cast<size_type> 1329 (intermediate_values().size())); 1330 if(! intermediate_values()[i]) { 1331 // Create it 1332 if(event_coordinates().size()==0) { 1333 CGAL_assertion(i==0); 1334 intermediate_values()[0]=Bound(0); 1335 } else { 1336 if(i==0) { 1337 intermediate_values()[i] 1338 = bound_left_of(kernel(),event_coordinates()[i].val); 1339 } else if(i == static_cast<size_type> 1340 (event_coordinates().size())) { 1341 intermediate_values()[i] 1342 = bound_right_of 1343 (kernel(),event_coordinates()[i-1].val); 1344 1345 } else { 1346 intermediate_values()[i] 1347 = kernel()->bound_between_1_object() 1348 (event_coordinates()[i-1].val, 1349 event_coordinates()[i].val); 1350 } 1351 } 1352 } 1353 return intermediate_values()[i].get(); 1354 } 1355 1356 1357 public: 1358 1359 /*! 1360 * Returns the content of the defining polynomial 1361 * 1362 * The content is the gcd of its coefficients (the polynomial is considered 1363 * as polynomial in \c y) 1364 */ 1365 Polynomial_1 content() const { 1366 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 1367 if(CGAL::degree(polynomial_2(),1)==2) { 1368 return this->conic_content(); 1369 } 1370 #endif 1371 if(! this->ptr()->content) { 1372 compute_content_and_primitive_part(); 1373 } 1374 return this->ptr()->content.get(); 1375 } 1376 1377 public: 1378 1379 /*! 1380 * Returns the primitive part of the defining polynomial 1381 * 1382 * The primitive part of \c f is the \c f divided by its content. 1383 */ 1384 Polynomial_2 primitive_polynomial_2() const { 1385 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 1386 if(CGAL::degree(polynomial_2(),1)==2) { 1387 return this->conic_primitive_polynomial_2(); 1388 } 1389 #endif 1390 if(! this->ptr()->f_primitive) { 1391 compute_content_and_primitive_part(); 1392 } 1393 return this->ptr()->f_primitive.get(); 1394 } 1395 1396 Algebraic_kernel_with_analysis_2* kernel() const { 1397 return this->ptr()->_m_kernel; 1398 } 1399 1400 private: 1401 1402 1403 // computes and sets the content and the primitive part for the curve 1404 void compute_content_and_primitive_part() const { 1405 1406 CGAL_assertion(has_defining_polynomial()); 1407 #if CGAL_ACK_DEBUG_FLAG 1408 CGAL_ACK_DEBUG_PRINT << "Computing the content..." << std::flush; 1409 #endif 1410 this->ptr()->content 1411 = typename CGAL::Polynomial_traits_d< Polynomial_2 >:: 1412 Univariate_content_up_to_constant_factor()( polynomial_2() ); 1413 if(CGAL::degree(content())==0) { 1414 #if CGAL_ACK_DEBUG_FLAG 1415 CGAL_ACK_DEBUG_PRINT << "no vertical lines as components" 1416 << std::endl; 1417 #endif 1418 this->ptr()->f_primitive=polynomial_2(); 1419 } 1420 else { 1421 #if CGAL_ACK_DEBUG_FLAG 1422 CGAL_ACK_DEBUG_PRINT << "non-trivial content found" << std::endl; 1423 #endif 1424 // Content must be square free, because the curve is square free 1425 CGAL_assertion( typename CGAL::Polynomial_traits_d< Polynomial_1 > 1426 ::Is_square_free()(content())); 1427 this->ptr()->f_primitive=polynomial_2() / content(); 1428 1429 } 1430 1431 } 1432 1433 private: 1434 1435 //! Returns the Sturm-Habicht sequence of the primitive part of f 1436 std::vector<Polynomial_2>& sturm_habicht_of_primitive() const 1437 { 1438 if(! this->ptr()->sturm_habicht_of_primitive) { 1439 compute_sturm_habicht_of_primitive(); 1440 } 1441 return this->ptr()->sturm_habicht_of_primitive.get(); 1442 } 1443 1444 public: 1445 1446 /*! 1447 * \brief returns the <tt>i</tt>th Sturm-Habicht polynomial 1448 * of the primitive part of the defining polynomial 1449 */ 1450 Polynomial_2 sturm_habicht_of_primitive(size_type i) const 1451 { 1452 CGAL_assertion(i>=0 && 1453 i < static_cast<size_type> 1454 (sturm_habicht_of_primitive().size())); 1455 return sturm_habicht_of_primitive()[i]; 1456 } 1457 1458 public: 1459 1460 /*! 1461 * \brief returns the <tt>i</tt>th principal Sturm-Habicht coefficient 1462 * of the primitive part of the defining polynomial 1463 */ 1464 Polynomial_1 principal_sturm_habicht_of_primitive(size_type i) const 1465 { 1466 CGAL_assertion(i>=0 && 1467 i < static_cast<size_type> 1468 (sturm_habicht_of_primitive().size())); 1469 1470 CGAL_assertion(CGAL::degree(sturm_habicht_of_primitive()[i])<=i); 1471 if(CGAL::degree(sturm_habicht_of_primitive()[i]) < i) { 1472 return Polynomial_1(0); 1473 } // else: 1474 return sturm_habicht_of_primitive()[i][i]; 1475 } 1476 1477 public: 1478 1479 /*! 1480 * \brief returns the <tt>i</tt>th coprincipal Sturm-Habicht coefficient 1481 * of the primitive part of the defining polynomial 1482 * 1483 * The coprincipal Sturm-Habicht coefficient is the coefficient 1484 * of <tt>y^{i-1}</tt> of the <tt>i</tt>th Sturm-Habicht polynomial 1485 */ 1486 Polynomial_1 coprincipal_sturm_habicht_of_primitive(size_type i) const 1487 { 1488 CGAL_assertion(i>=1 && 1489 i < static_cast<size_type> 1490 (sturm_habicht_of_primitive().size())); 1491 CGAL_assertion(CGAL::degree(sturm_habicht_of_primitive()[i])<=i); 1492 if(CGAL::degree(sturm_habicht_of_primitive()[i]) < i-1) { 1493 return Polynomial_1(0); 1494 } // else: 1495 return sturm_habicht_of_primitive()[i][i-1]; 1496 } 1497 1498 public: 1499 1500 /*! 1501 * \brief returns an iterator to the principal Sturm-Habicht coefficients, 1502 * starting with the <tt>0</tt>th one (the resultant) 1503 */ 1504 Principal_sturm_habicht_iterator principal_sturm_habicht_begin() const { 1505 return boost::make_transform_iterator 1506 (boost::counting_iterator<size_type>(0), 1507 Stha_functor(this)); 1508 } 1509 1510 //! Returns an iterator to the end of principal Sturm-Habicht coefficients 1511 Principal_sturm_habicht_iterator principal_sturm_habicht_end() const { 1512 return boost::make_transform_iterator 1513 (boost::counting_iterator<size_type> 1514 (static_cast<int>(sturm_habicht_of_primitive().size())), 1515 Stha_functor(this)); 1516 } 1517 1518 private: 1519 1520 // Internal method to compute the Sturm-Habicht sequence 1521 void compute_sturm_habicht_of_primitive() const 1522 { 1523 1524 #if CGAL_ACK_DEBUG_FLAG 1525 CGAL_ACK_DEBUG_PRINT << "Compute Sturm-Habicht.." << std::flush; 1526 #endif 1527 std::vector<Polynomial_2> stha; 1528 1529 // Fix a problem for constant primitive part. 1530 // In this case, the St.-Ha. sequence is never needed 1531 if(CGAL::degree(primitive_polynomial_2()) == 0) { 1532 // Set the resultant 1533 stha.push_back(primitive_polynomial_2()); 1534 } else { 1535 1536 #if CGAL_ACK_USE_BEZOUT_MATRIX_FOR_SUBRESULTANTS 1537 #warning USES BEZOUT MATRIX FOR SUBRESULTANTS 1538 CGAL::internal::bezout_polynomial_subresultants<Polynomial_traits_2> 1539 (primitive_polynomial_2(), 1540 CGAL::differentiate(primitive_polynomial_2()), 1541 std::back_inserter(stha)); 1542 stha.push_back(primitive_polynomial_2()); 1543 size_type p = CGAL::degree(primitive_polynomial_2()); 1544 CGAL_assertion(static_cast<size_type>(stha.size()) == p+1); 1545 for(size_type i=0;i<p; i++) { 1546 if((p-i)%4==0 || (p-i)%4==1) { 1547 stha[i] = stha[i]; 1548 } else { 1549 stha[i] = -stha[i]; 1550 } 1551 } 1552 1553 #else 1554 typename Polynomial_traits_2::Sturm_habicht_sequence() 1555 (primitive_polynomial_2(),std::back_inserter(stha)); 1556 #endif 1557 } 1558 // Also set the resultant, if not yet set 1559 if(! this->ptr()->resultant_of_primitive_and_derivative_y) { 1560 this->ptr()->resultant_of_primitive_and_derivative_y = stha[0][0]; 1561 if(this->ptr()->resultant_of_primitive_and_derivative_y. 1562 get().is_zero()) { 1563 throw internal::Zero_resultant_exception<Polynomial_2> 1564 (polynomial_2()); 1565 } 1566 } 1567 1568 this->ptr()->sturm_habicht_of_primitive = stha; 1569 CGAL_assertion(CGAL::canonicalize 1570 (resultant_of_primitive_and_derivative_y()) == 1571 CGAL::canonicalize 1572 (principal_sturm_habicht_of_primitive(0))); 1573 #if CGAL_ACK_DEBUG_FLAG 1574 CGAL_ACK_DEBUG_PRINT << "done" << std::endl; 1575 #endif 1576 } 1577 1578 private: 1579 1580 //! Returns the resultant of the primitive part of f and its y-derivative 1581 Polynomial_1 resultant_of_primitive_and_derivative_y() const 1582 { 1583 if(! this->ptr()->resultant_of_primitive_and_derivative_y) { 1584 compute_resultant_of_primitive_and_derivative_y(); 1585 } 1586 return this->ptr()->resultant_of_primitive_and_derivative_y.get(); 1587 } 1588 1589 private: 1590 1591 //! Returns the resultant of the primitive part of f with its x-derivative 1592 Polynomial_1 resultant_of_primitive_and_derivative_x() const 1593 { 1594 if(! this->ptr()->resultant_of_primitive_and_derivative_x) { 1595 compute_resultant_of_primitive_and_derivative_x(); 1596 } 1597 return this->ptr()->resultant_of_primitive_and_derivative_x.get(); 1598 } 1599 1600 private: 1601 // Computes <tt>res_y(f,f_y)</tt>, where \c f is the defining polynomial 1602 void compute_resultant_of_primitive_and_derivative_y() const 1603 { 1604 1605 #if CGAL_ACK_DEBUG_FLAG 1606 CGAL_ACK_DEBUG_PRINT << "Compute resultant.." << std::flush; 1607 #endif 1608 1609 CGAL_assertion(has_defining_polynomial()); 1610 1611 #if CGAL_ACK_RESULTANT_FIRST_STRATEGY 1612 #ifndef CGAL_ACK_RESULTANT_FIRST_STRATEGY_DEGREE_THRESHOLD 1613 bool speed_up = true; 1614 #else 1615 bool speed_up=CGAL::degree(polynomial_2()) >= 1616 CGAL_ACK_RESULTANT_FIRST_STRATEGY_DEGREE_THRESHOLD; 1617 #endif 1618 #else 1619 bool speed_up=false; 1620 #endif 1621 1622 if(CGAL::degree(polynomial_2()) == 0) { 1623 this->ptr()->resultant_of_primitive_and_derivative_y 1624 = Polynomial_1(1); 1625 } else { 1626 1627 if(! speed_up) { 1628 1629 // Compute resultant using the Sturm-Habicht sequence 1630 this->ptr()->resultant_of_primitive_and_derivative_y 1631 = principal_sturm_habicht_of_primitive(0); 1632 1633 } else { 1634 typename Polynomial_traits_2::Differentiate diff; 1635 this->ptr()->resultant_of_primitive_and_derivative_y 1636 = CGAL::resultant 1637 (primitive_polynomial_2(), 1638 diff(primitive_polynomial_2(),1)); 1639 } 1640 1641 } 1642 1643 #if CGAL_ACK_DEBUG_FLAG 1644 CGAL_ACK_DEBUG_PRINT << "done" << std::endl; 1645 #endif 1646 1647 if(resultant_of_primitive_and_derivative_y().is_zero()) { 1648 throw internal::Zero_resultant_exception<Polynomial_2> 1649 (polynomial_2()); 1650 } 1651 } 1652 1653 // Computes <tt>res_y(f,f_x)</tt>, where \c f is the defining polynomial 1654 void compute_resultant_of_primitive_and_derivative_x() const 1655 { 1656 1657 #if CGAL_ACK_DEBUG_FLAG 1658 CGAL_ACK_DEBUG_PRINT << "Compute x-resultant.." << std::flush; 1659 #endif 1660 1661 CGAL_assertion(has_defining_polynomial()); 1662 1663 // Transpose the polynomial 1664 Polynomial_2 f_yx = typename Polynomial_traits_2::Swap() 1665 (polynomial_2(),0,1); 1666 1667 if( CGAL::degree(f_yx) == 0 ) { 1668 // Polynomial only consists of horizontal lines 1669 // primitive resultant is set to 1 1670 this->ptr()->resultant_of_primitive_and_derivative_x 1671 = Polynomial_1(1); 1672 } else { 1673 1674 Polynomial_2 f_yx_primitive; 1675 1676 Polynomial_1 content_yx 1677 = typename CGAL::Polynomial_traits_d< Polynomial_2 >:: 1678 Univariate_content_up_to_constant_factor()( f_yx ); 1679 if(CGAL::degree(content_yx)==0) { 1680 f_yx_primitive=f_yx; 1681 } 1682 else { 1683 CGAL_assertion 1684 (typename CGAL::Polynomial_traits_d< Polynomial_1 >:: 1685 Is_square_free()(content_yx)); 1686 f_yx_primitive=f_yx / content_yx; 1687 1688 } 1689 1690 this->ptr()->resultant_of_primitive_and_derivative_x 1691 = CGAL::resultant 1692 (typename Polynomial_traits_2::Swap() (f_yx_primitive,0,1), 1693 typename Polynomial_traits_2::Swap() 1694 (CGAL::differentiate(f_yx_primitive),0,1) ); 1695 } 1696 1697 #if CGAL_ACK_DEBUG_FLAG 1698 CGAL_ACK_DEBUG_PRINT << "done" << std::endl; 1699 #endif 1700 1701 if(resultant_of_primitive_and_derivative_x().is_zero()) { 1702 throw internal::Zero_resultant_exception<Polynomial_2> 1703 (polynomial_2()); 1704 } 1705 } 1706 1707 1708 1709 1710 private: 1711 1712 // Returns the critical event coordinates 1713 std::vector<Event_coordinate_1>& event_coordinates() const 1714 { 1715 if(! this->ptr()->event_coordinates) { 1716 compute_event_coordinates(); 1717 } 1718 return this->ptr()->event_coordinates.get(); 1719 } 1720 1721 private: 1722 1723 // Returns the intermediate values for intervals between events 1724 std::vector<boost::optional<Bound> >& intermediate_values() const 1725 { 1726 if(! this->ptr()->intermediate_values) { 1727 // This is created during event_coordiantes() 1728 event_coordinates(); 1729 CGAL_assertion(bool(this->ptr()->intermediate_values)); 1730 } 1731 return this->ptr()->intermediate_values.get(); 1732 } 1733 1734 1735 private: 1736 1737 /* 1738 * \brief Computes the event coordinates of the curve. 1739 * 1740 * This function computes the content of the defining polynomial, 1741 * and the roots of its discriminant. These two sets form the critical 1742 * x-coordinates of the curve. 1743 */ 1744 void compute_event_coordinates() const 1745 { 1746 1747 #if CGAL_ACK_DEBUG_FLAG 1748 CGAL_ACK_DEBUG_PRINT << "compute events..." << std::flush; 1749 #endif 1750 1751 Solve_1 solve_1; 1752 1753 std::vector<std::pair<Algebraic_real_1,size_type> > content_pairs; 1754 std::vector<Algebraic_real_1> content_roots; 1755 std::vector<size_type> content_mults; 1756 solve_1(content(), 1757 std::back_inserter(content_pairs)); 1758 1759 for(int i=0; i < static_cast<int>(content_pairs.size()); i++ ) { 1760 content_roots.push_back(content_pairs[i].first); 1761 content_mults.push_back(content_pairs[i].second); 1762 } 1763 1764 // Set the vertical_line_components flag as side effect 1765 this->ptr()->has_vertical_component = (content_roots.size() > 0); 1766 1767 std::vector<std::pair<Algebraic_real_1,size_type> > res_pairs; 1768 std::vector<Algebraic_real_1> res_roots; 1769 std::vector<size_type> res_mults; 1770 Polynomial_1 R = resultant_of_primitive_and_derivative_y(); 1771 solve_1(R,std::back_inserter(res_pairs)); 1772 1773 for(int i=0; i < static_cast<int>(res_pairs.size()); i++ ) { 1774 res_roots.push_back(res_pairs[i].first); 1775 res_mults.push_back(res_pairs[i].second); 1776 } 1777 1778 std::vector<std::pair<Algebraic_real_1,size_type> > lcoeff_pairs; 1779 std::vector<Algebraic_real_1> lcoeff_roots; 1780 std::vector<size_type> lcoeff_mults; 1781 solve_1(CGAL::leading_coefficient(primitive_polynomial_2()), 1782 std::back_inserter(lcoeff_pairs)); 1783 1784 for(int i=0; i < static_cast<int>(lcoeff_pairs.size()); i++ ) { 1785 lcoeff_roots.push_back(lcoeff_pairs[i].first); 1786 lcoeff_mults.push_back(lcoeff_pairs[i].second); 1787 } 1788 1789 1790 //Now, merge the vertical line positions with the resultant roots 1791 typename 1792 CGAL::Real_embeddable_traits<Algebraic_real_1>::Compare compare; 1793 1794 std::vector<Algebraic_real_1> event_values; 1795 std::vector<CGAL::internal::Three_valued> event_values_info; 1796 1797 CGAL::internal::set_union_with_source 1798 (res_roots.begin(), 1799 res_roots.end(), 1800 content_roots.begin(), 1801 content_roots.end(), 1802 std::back_inserter(event_values), 1803 std::back_inserter(event_values_info), 1804 compare); 1805 1806 // Now, build the Event_coordinate_1 entries 1807 // for each element of event_values 1808 size_type curr_res_index = 0, curr_content_index = 0, 1809 curr_lcoeff_index = 0; 1810 std::vector<Event_coordinate_1> event_coordinate_vector; 1811 1812 for(size_type i = 0; 1813 i < static_cast<size_type>(event_values.size()); 1814 i++ ) { 1815 1816 Event_coordinate_1 curr_event; 1817 curr_event.val = event_values[i]; 1818 switch(event_values_info[i]) { 1819 1820 case(CGAL::internal::ROOT_OF_FIRST_SET): { 1821 curr_event.index_of_prim_res_root = curr_res_index; 1822 CGAL_expensive_assertion(res_roots[curr_res_index] == 1823 event_values[i]); 1824 curr_event.mult_of_prim_res_root 1825 = res_mults[curr_res_index]; 1826 curr_res_index++; 1827 if(curr_lcoeff_index < 1828 static_cast<size_type>(lcoeff_roots.size()) && 1829 event_values[i]==lcoeff_roots[curr_lcoeff_index]) { 1830 // We have a root of the leading coefficient 1831 // of the primitve polynomial 1832 curr_event.index_of_prim_lcoeff_root = curr_lcoeff_index; 1833 curr_event.mult_of_prim_lcoeff_root 1834 = lcoeff_mults[curr_lcoeff_index]; 1835 curr_lcoeff_index++; 1836 } else { 1837 curr_event.index_of_prim_lcoeff_root = -1; 1838 curr_event.mult_of_prim_lcoeff_root = 0; 1839 } 1840 1841 curr_event.index_of_content_root = -1; 1842 curr_event.mult_of_content_root = 0; 1843 break; 1844 } 1845 case(CGAL::internal::ROOT_OF_SECOND_SET): { 1846 curr_event.index_of_content_root = curr_content_index; 1847 CGAL_expensive_assertion(content_roots[curr_content_index] == 1848 event_values[i]); 1849 curr_event.mult_of_content_root 1850 = content_mults[curr_content_index]; 1851 curr_content_index++; 1852 curr_event.index_of_prim_res_root = -1; 1853 curr_event.mult_of_prim_res_root = 0; 1854 CGAL_expensive_assertion(event_values[i]!= 1855 lcoeff_roots[curr_lcoeff_index]); 1856 curr_event.index_of_prim_lcoeff_root = -1; 1857 curr_event.mult_of_prim_lcoeff_root = 0; 1858 break; 1859 } 1860 case(CGAL::internal::ROOT_OF_BOTH_SETS): { 1861 curr_event.index_of_prim_res_root = curr_res_index; 1862 CGAL_expensive_assertion(res_roots[curr_res_index] == 1863 event_values[i]); 1864 curr_event.mult_of_prim_res_root 1865 = res_mults[curr_res_index]; 1866 curr_res_index++; 1867 if(curr_lcoeff_index < 1868 static_cast<size_type>(lcoeff_roots.size()) && 1869 event_values[i]==lcoeff_roots[curr_lcoeff_index]) { 1870 // We have a root of the leading coefficient 1871 // of the primitve polynomial 1872 curr_event.index_of_prim_lcoeff_root = curr_lcoeff_index; 1873 curr_event.mult_of_prim_lcoeff_root 1874 = lcoeff_mults[curr_lcoeff_index]; 1875 curr_lcoeff_index++; 1876 } else { 1877 curr_event.index_of_prim_lcoeff_root = -1; 1878 curr_event.mult_of_prim_lcoeff_root = 0; 1879 } 1880 curr_event.index_of_content_root = curr_content_index; 1881 CGAL_expensive_assertion(content_roots[curr_content_index] == 1882 event_values[i]); 1883 curr_event.mult_of_content_root 1884 = content_mults[curr_content_index]; 1885 curr_content_index++; 1886 break; 1887 } 1888 } // of switch 1889 /* 1890 #if CGAL_ACK_DEBUG_FLAG 1891 CGAL_ACK_DEBUG_PRINT << "Constructed event_coordinate: " 1892 << CGAL::to_double(curr_event.val) << " " 1893 << "\nmult_of_prim_res_root : " 1894 << curr_event.mult_of_prim_res_root 1895 << "\nindex_of_prim_res_root : " 1896 << curr_event.index_of_prim_res_root 1897 << "\nmult_of_content_root : " 1898 << curr_event.mult_of_content_root 1899 << "\nindex_of_content_root : " 1900 << curr_event.index_of_content_root 1901 << "\nmult_of_lcoeff_root : " 1902 << curr_event.mult_of_prim_lcoeff_root 1903 << "\nindex_of_lcoeff_root : " 1904 << curr_event.index_of_prim_lcoeff_root 1905 << std::endl; 1906 #endif 1907 */ 1908 event_coordinate_vector.push_back(curr_event); 1909 } 1910 1911 1912 CGAL_assertion(curr_lcoeff_index == 1913 static_cast<size_type>(lcoeff_roots.size())); 1914 CGAL_assertion(curr_res_index == 1915 static_cast<size_type>(res_roots.size())); 1916 CGAL_assertion(curr_content_index == 1917 static_cast<size_type>(content_roots.size())); 1918 1919 this->ptr()->intermediate_values 1920 = std::vector<boost::optional<Bound> > 1921 (event_coordinate_vector.size()+1); 1922 this->ptr()->event_coordinates = event_coordinate_vector; 1923 1924 #if CGAL_ACK_DEBUG_FLAG 1925 CGAL_ACK_DEBUG_PRINT << "done" << std::endl; 1926 #endif 1927 1928 } 1929 1930 public: 1931 1932 /*! 1933 * \brief returns a \c Curve_analysis_2 object for a sheared curve. 1934 * 1935 * The shear factor is given by the integer \c s. 1936 * This functions only shears the primitive part of the defining equation. 1937 * Internal caching is used to avoid repeated shears. 1938 * 1939 * \todo The sheared curves are not inserted into the curve_cache 1940 * of the Algebraic_curve_kernel_2 yet. 1941 */ 1942 Self& shear_primitive_part(Integer s) const 1943 { 1944 CGAL_assertion(s!=0); 1945 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 1946 if(CGAL::degree(polynomial_2(),1)==2) { 1947 return this->conic_shear_primitive_part(); 1948 } 1949 #endif 1950 if(this->ptr()->bad_shears.find(s) != 1951 this->ptr()->bad_shears.end()) { 1952 throw CGAL::internal::Non_generic_position_exception(); 1953 } 1954 typedef typename std::map<Integer,Self>::iterator 1955 Map_iterator; 1956 Map_iterator it = this->ptr()->sheared_curves.find(s); 1957 if(it != this->ptr()->sheared_curves.end()) { 1958 return it->second; 1959 } 1960 try { 1961 Shear_transformation<Algebraic_kernel_with_analysis_2> 1962 shear_transformation(kernel()); 1963 Self D=shear_transformation((Self&)*this, s); 1964 std::pair<Map_iterator,bool> insertion = 1965 this->ptr()->sheared_curves.insert(std::make_pair(s,D)); 1966 CGAL_assertion(insertion.second); 1967 return insertion.first->second; 1968 } 1969 catch(CGAL::internal::Non_generic_position_exception /* err */) { 1970 this->ptr()->bad_shears.insert(s); 1971 throw CGAL::internal::Non_generic_position_exception(); 1972 } 1973 } 1974 1975 public: 1976 1977 //! Iterator for sheared curves 1978 typename std::map<Coefficient,Self>::const_iterator shear_begin() { 1979 return this->ptr()->sheared_curves.begin(); 1980 } 1981 1982 //! Iterator for sheared curves 1983 typename std::map<Coefficient,Self>::const_iterator shear_end() { 1984 return this->ptr()->sheared_curves.end(); 1985 } 1986 1987 private: 1988 1989 // Sets the flag for vertical lines in all status lines that need it 1990 void set_vertical_line_components() const { 1991 for(size_type i = 0; 1992 i < static_cast<size_type>(event_coordinates().size()); 1993 i++ ) { 1994 1995 if(event_coordinates()[i].mult_of_content_root > 0) { 1996 status_line_at_event(i)._set_v_line(); 1997 } 1998 } 1999 2000 } 2001 2002 2003 public: 2004 2005 /*! 2006 * \brief Increases the precision of all status lines 2007 * 2008 * For each status line at an event and each status line that represents 2009 * an interval, all y-coordinates are approximated such that their 2010 * isolating interval has absolute size smaller than \c precision. 2011 */ 2012 void refine_all(Bound precision) { 2013 2014 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 2015 if(CGAL::degree(polynomial_2(),1)==2) { 2016 return this->conic_refine_all(precision); 2017 } 2018 #endif 2019 2020 for(size_type i=0; 2021 i<static_cast<size_type>(event_coordinates().size()); 2022 i++) { 2023 /* 2024 #if CGAL_ACK_DEBUG_FLAG 2025 CGAL_ACK_DEBUG_PRINT << i << ": " << std::flush; 2026 #endif 2027 */ 2028 Status_line_1& el = status_line_at_event(i); 2029 2030 for(size_type j=0;j<el.number_of_events();j++) { 2031 /* 2032 #if CGAL_ACK_DEBUG_FLAG 2033 CGAL_ACK_DEBUG_PRINT << j << " " << std::flush; 2034 #endif 2035 */ 2036 el.refine_to(j,precision); 2037 } 2038 } 2039 for(size_type i=0; 2040 i<static_cast<size_type>(intermediate_values().size()); 2041 i++) { 2042 Status_line_1 il = status_line_of_interval(i); 2043 for(size_type j=0;j<il.number_of_events();j++) { 2044 il.refine_to(j,precision); 2045 } 2046 } 2047 } 2048 2049 public: 2050 2051 //! \brief Iterator for the status lines at events 2052 Event_line_iterator event_begin() const { 2053 return boost::make_transform_iterator 2054 (boost::counting_iterator<size_type>(0), 2055 Event_functor(this)); 2056 } 2057 2058 //! \brief Iterator for the status lines at events 2059 Event_line_iterator event_end() const { 2060 return boost::make_transform_iterator 2061 (boost::counting_iterator<size_type> 2062 (number_of_status_lines_with_event()), 2063 Event_functor(this)); 2064 } 2065 2066 public: 2067 2068 //! \brief Iterator for the status lines for intervals 2069 Intermediate_line_iterator intermediate_begin() const { 2070 return boost::make_transform_iterator 2071 (boost::counting_iterator<size_type>(0), 2072 Intermediate_functor(this)); 2073 } 2074 2075 //! \brief Iterator for the status lines for intervals 2076 Intermediate_line_iterator intermediate_end() const { 2077 return boost::make_transform_iterator 2078 (boost::counting_iterator<size_type>(intermediate_values().size()), 2079 Intermediate_functor(this)); 2080 } 2081 2082 public: 2083 2084 /*! 2085 * \brief returns the limit an infinite arc converges to 2086 * 2087 * \pre <tt>loc==CGAL::LEFT_BOUNDARY || 2088 * loc==CGAL::RIGHT_BOUNDARY</tt> 2089 * 2090 * This method returns for the <tt>arcno</tt>th arc that goes to -infinity 2091 * or +infinity (depending on \c loc) the y-coordinate it converges to. 2092 * Possible values are either a \c Algebraic_real_1 object, or one of the 2093 * values \c CGAL::TOP_BOUNDARY, \c CGAL::BOTTOM_BOUNDARY 2094 * that denote that the arc is unbounded in y-direction. 2095 * The result is wrapped into a \c CGAL::Object object. 2096 */ 2097 Asymptote_y asymptotic_value_of_arc(CGAL::Box_parameter_space_2 loc, 2098 size_type arcno) const { 2099 2100 CGAL_precondition(loc == CGAL::LEFT_BOUNDARY || 2101 loc == CGAL::RIGHT_BOUNDARY); 2102 2103 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 2104 if(CGAL::degree(polynomial_2(),1)==2) { 2105 return this->conic_asymptotic_value_of_arc(loc,arcno); 2106 } 2107 #endif 2108 2109 if(loc == CGAL::LEFT_BOUNDARY) { 2110 2111 if(! this->ptr()->horizontal_asymptotes_left) { 2112 compute_horizontal_asymptotes(); 2113 } 2114 std::vector<Asymptote_y>& asym_info 2115 = this->ptr()->horizontal_asymptotes_left.get(); 2116 CGAL_precondition(arcno>=0 && 2117 arcno<static_cast<size_type>(asym_info.size())); 2118 return asym_info[arcno]; 2119 } // else loc == CGAL::RIGHT_BOUNDARY 2120 2121 if(! this->ptr()->horizontal_asymptotes_right) { 2122 compute_horizontal_asymptotes(); 2123 } 2124 std::vector<Asymptote_y>& asym_info 2125 = this->ptr()->horizontal_asymptotes_right.get(); 2126 CGAL_precondition(arcno>=0 && 2127 arcno<static_cast<size_type>(asym_info.size())); 2128 return asym_info[arcno]; 2129 2130 } 2131 2132 2133 private: 2134 2135 // Internal method to compute horizontal asymptotes 2136 void compute_horizontal_asymptotes() const { 2137 2138 // TODO: Filter out curves with no arc to +/- infty 2139 2140 Solve_1 solve_1 = kernel()->solve_1_object(); 2141 2142 Polynomial_1 leading_coefficient_in_x 2143 = CGAL::leading_coefficient(typename Polynomial_traits_2::Swap() 2144 (polynomial_2(),0,1)); 2145 std::vector<Algebraic_real_1> roots_of_lcoeff; 2146 2147 solve_1(leading_coefficient_in_x, 2148 std::back_inserter(roots_of_lcoeff), 2149 false); 2150 2151 2152 std::vector<Bound> stripe_bounds; 2153 find_intermediate_values(kernel(), 2154 roots_of_lcoeff.begin(), 2155 roots_of_lcoeff.end(), 2156 std::back_inserter(stripe_bounds)); 2157 Bound leftmost_bound = bound_value_in_interval(0), 2158 rightmost_bound = bound_value_in_interval 2159 (this->number_of_status_lines_with_event()); 2160 for(size_type i=0;i<static_cast<size_type>(stripe_bounds.size());i++) { 2161 Bound& beta = stripe_bounds[i]; 2162 Polynomial_1 poly_at_beta 2163 = kernel()->evaluate_utcf_2_object()(this->polynomial_2(),beta); 2164 std::vector<Algebraic_real_1> x_coordinates_at_beta; 2165 solve_1(poly_at_beta,std::back_inserter(x_coordinates_at_beta), 2166 false); 2167 size_type number_of_roots 2168 = static_cast<size_type>(x_coordinates_at_beta.size()); 2169 if(number_of_roots>0) { 2170 if(leftmost_bound > x_coordinates_at_beta[0].low()) { 2171 leftmost_bound = x_coordinates_at_beta[0].low(); 2172 } 2173 if(rightmost_bound 2174 < x_coordinates_at_beta[number_of_roots-1].high()) { 2175 rightmost_bound 2176 = x_coordinates_at_beta[number_of_roots-1].high(); 2177 } 2178 } 2179 } 2180 2181 // Just to be sure... 2182 leftmost_bound = leftmost_bound - 1; 2183 rightmost_bound = rightmost_bound + 1; 2184 2185 Polynomial_1 curve_at_left_end 2186 = kernel()->evaluate_utcf_2_object() 2187 (typename Polynomial_traits_2::Swap() (this->polynomial_2(),0,1), 2188 leftmost_bound); 2189 std::vector<Algebraic_real_1> roots_at_left_end; 2190 solve_1(curve_at_left_end,std::back_inserter(roots_at_left_end),false); 2191 size_type number_of_roots_at_left_end 2192 = static_cast<size_type>(roots_at_left_end.size()); 2193 std::vector<Asymptote_y> asym_left_info; 2194 size_type current_stripe=0,i=0; 2195 while(i<number_of_roots_at_left_end) { 2196 if(current_stripe==static_cast<size_type>(stripe_bounds.size())) { 2197 asym_left_info.push_back( CGAL::make_object 2198 (CGAL::TOP_BOUNDARY) ); 2199 i++; 2200 continue; 2201 } 2202 if(roots_at_left_end[i].low() > stripe_bounds[current_stripe]) { 2203 current_stripe++; 2204 continue; 2205 } 2206 if(roots_at_left_end[i].high() < stripe_bounds[current_stripe]) { 2207 if(current_stripe==0) { 2208 asym_left_info.push_back(CGAL::make_object 2209 (CGAL::BOTTOM_BOUNDARY)); 2210 i++; 2211 continue; 2212 } else { 2213 asym_left_info.push_back(CGAL::make_object 2214 (roots_of_lcoeff[current_stripe-1])); 2215 i++; 2216 continue; 2217 } 2218 } 2219 roots_at_left_end[i].refine(); 2220 } 2221 this->ptr()->horizontal_asymptotes_left = asym_left_info; 2222 2223 Polynomial_1 curve_at_right_end 2224 = kernel()->evaluate_utcf_2_object() 2225 (typename Polynomial_traits_2::Swap() (this->polynomial_2(),0,1), 2226 rightmost_bound); 2227 std::vector<Algebraic_real_1> roots_at_right_end; 2228 solve_1(curve_at_right_end,std::back_inserter(roots_at_right_end),false); 2229 size_type number_of_roots_at_right_end 2230 = static_cast<size_type>(roots_at_right_end.size()); 2231 std::vector<Asymptote_y> asym_right_info; 2232 current_stripe=0; 2233 i=0; 2234 while(i<number_of_roots_at_right_end) { 2235 if(current_stripe==static_cast<size_type>(stripe_bounds.size())) { 2236 asym_right_info.push_back(CGAL::make_object 2237 (CGAL::TOP_BOUNDARY) ); 2238 i++; 2239 continue; 2240 } 2241 if(roots_at_right_end[i].low() > stripe_bounds[current_stripe]) { 2242 current_stripe++; 2243 continue; 2244 } 2245 if(roots_at_right_end[i].high() < stripe_bounds[current_stripe]) { 2246 if(current_stripe==0) { 2247 asym_right_info.push_back(CGAL::make_object 2248 (CGAL::BOTTOM_BOUNDARY)); 2249 i++; 2250 continue; 2251 } else { 2252 asym_right_info.push_back 2253 (CGAL::make_object(roots_of_lcoeff[current_stripe-1])); 2254 i++; 2255 continue; 2256 } 2257 } 2258 roots_at_right_end[i].refine(); 2259 } 2260 this->ptr()->horizontal_asymptotes_right = asym_right_info; 2261 2262 } 2263 2264 //! @} 2265 2266 public: 2267 2268 template<typename OutputIterator> void get_roots_at_rational 2269 (Bound r, OutputIterator it) const { 2270 2271 typename Rep::Intermediate_cache::Find_result find_result 2272 = this->ptr()->intermediate_cache.find(r); 2273 2274 std::vector<Algebraic_real_1> p_roots; 2275 2276 if(find_result.second) { 2277 p_roots = find_result.first->second; 2278 } else { 2279 Polynomial_2 swapped = typename Polynomial_traits_2::Swap() 2280 (this->polynomial_2(), 0, 1); 2281 Polynomial_1 p = kernel()->evaluate_utcf_2_object()(swapped,r); 2282 kernel()->solve_1_object()(p,std::back_inserter(p_roots),false); 2283 2284 this->ptr()->intermediate_cache.insert(std::make_pair(r,p_roots)); 2285 2286 } 2287 std::copy(p_roots.begin(),p_roots.end(),it); 2288 } 2289 2290 2291 // \name Internal functions for Conic optimization 2292 //! @{ 2293 2294 #if CGAL_ACK_USE_SPECIAL_TREATMENT_FOR_CONIX 2295 2296 private: 2297 2298 bool conic_is_y_regular() const { 2299 CGAL_error_msg("Implement me"); 2300 return false; 2301 } 2302 2303 bool conic_has_vertical_component() const { 2304 CGAL_error_msg("Implement me"); 2305 return false; 2306 } 2307 2308 size_type conic_number_of_status_lines_with_event() const { 2309 CGAL_error_msg("Implement me"); 2310 return 0; 2311 } 2312 2313 void conic_x_to_index(Algebraic_real_1 x,size_type& i,bool& is_event) const 2314 { 2315 CGAL_error_msg("Implement me"); 2316 } 2317 2318 Status_line_1& conic_status_line_at_event(size_type i) const { 2319 CGAL_error_msg("Implement me"); 2320 // Just a random status line to make compiler happy 2321 return this->ptr()->vert_line_at_rational_map[Bound(0)]; 2322 } 2323 2324 Status_line_1& conic_status_line_at_exact_x(Bound b) const { 2325 CGAL_error_msg("Implement me"); 2326 return this->ptr()->vert_line_at_rational_map[Bound(0)]; 2327 } 2328 2329 Status_line_1& conic_status_line_at_exact_x(Algebraic_real_1 alpha) const { 2330 CGAL_error_msg("Implement me"); 2331 return this->ptr()->vert_line_at_rational_map[Bound(0)]; 2332 } 2333 2334 Status_line_1 conic_status_line_of_interval(size_type i) const { 2335 CGAL_error_msg("Implement me"); 2336 return this->ptr()->vert_line_at_rational_map[Bound(0)]; 2337 } 2338 2339 Status_line_1 conic_status_line_for_x 2340 (Algebraic_real_1 x, 2341 CGAL::Sign perturb = CGAL::ZERO) const { 2342 CGAL_error_msg("Implement me"); 2343 return this->ptr()->vert_line_at_rational_map[Bound(0)]; 2344 } 2345 2346 size_type conic_arcs_over_interval(size_type i) const { 2347 CGAL_error_msg("Implement me"); 2348 return -1; 2349 } 2350 2351 Bound conic_bound_value_in_interval(size_type i) const { 2352 CGAL_error_msg("Implement me"); 2353 return Bound(0); 2354 } 2355 2356 Polynomial_1 conic_content() const { 2357 CGAL_error_msg("Implement me"); 2358 return Polynomial_1(); 2359 } 2360 2361 Polynomial_2 conic_primitive_polynomial_2() const { 2362 CGAL_error_msg("Implement me"); 2363 return Polynomial_2(); 2364 } 2365 2366 Self& conic_shear_primitive_part(Integer s) const { 2367 CGAL_error_msg("Implement me"); 2368 return Self(); 2369 } 2370 2371 void conic_refine_all(Bound precision) { 2372 CGAL_error_msg("Implement me"); 2373 } 2374 2375 Asymptote_y conic_asymptotic_value_of_arc(CGAL::Box_parameter_space_2 loc, 2376 size_type arcno) const { 2377 CGAL_error_msg("Implement me"); 2378 return Asymptote_y(); 2379 } 2380 2381 #endif 2382 2383 2384 //! @} 2385 2386 //! \name friends 2387 //! @{ 2388 2389 // friend function for id-based hashing 2390 friend std::size_t hash_value(const Self& x) { 2391 return static_cast<std::size_t>(x.id()); 2392 } 2393 2394 // another friend 2395 friend class Shear_transformation<Algebraic_kernel_with_analysis_2>; 2396 2397 //! @} 2398 2399 }; // class Algebraic_curve_2_2 2400 2401 2402 //! \brief prints the objects. 2403 template<typename AlgebraicKernelWithAnalysis_2, 2404 typename Rep_> 2405 std::ostream& operator<< ( 2406 std::ostream& out, 2407 const Curve_analysis_2< AlgebraicKernelWithAnalysis_2, 2408 Rep_ >& curve) { 2409 2410 typedef AlgebraicKernelWithAnalysis_2 Algebraic_kernel_with_analysis_2; 2411 2412 typedef Rep_ Rep; 2413 2414 typedef Curve_analysis_2< Algebraic_kernel_with_analysis_2, Rep > Curve; 2415 2416 typedef typename Curve::size_type size_type; 2417 typedef typename Curve::Asymptote_y Asymptote_y; 2418 2419 2420 switch (::CGAL::IO::get_mode(out)) { 2421 case ::CGAL::IO::PRETTY: { 2422 2423 out << "--------------- Analysis results ---------------" << std::endl; 2424 out << "Number of constructed event lines: " 2425 << curve.number_of_status_lines_with_event() 2426 << std::endl; 2427 out << "(Horizontal) asymptotes at -infty: " << std::flush; 2428 for (size_type i = 0; i < curve.arcs_over_interval(0); i++) { 2429 2430 const Asymptote_y& curr_asym_info_obj 2431 = curve.asymptotic_value_of_arc(CGAL::LEFT_BOUNDARY,i); 2432 typename Curve::Algebraic_real_1 curr_asym_info; 2433 bool is_finite = CGAL::assign(curr_asym_info,curr_asym_info_obj); 2434 if (!is_finite) { 2435 // Assignment to prevent compiler warning 2436 CGAL::Box_parameter_space_2 loc = CGAL::LEFT_BOUNDARY; 2437 CGAL_assertion_code(bool is_valid = ) 2438 CGAL::assign(loc, curr_asym_info_obj); 2439 CGAL_assertion(is_valid); 2440 if (loc == CGAL::TOP_BOUNDARY) { 2441 out << "+infty " << std::flush; 2442 } else { 2443 CGAL_assertion(loc == CGAL::BOTTOM_BOUNDARY); 2444 out << "-infty " << std::flush; 2445 } 2446 } else { // is_finite 2447 out << CGAL::to_double(curr_asym_info) 2448 << " " << std::flush; 2449 } 2450 } 2451 2452 out << std::endl; 2453 2454 out << "Intermediate line at " 2455 << CGAL::to_double(curve.bound_value_in_interval(0)) 2456 << ": " << curve.arcs_over_interval(0) << " passing arcs" 2457 << std::endl 2458 << std::endl; 2459 for (size_type i = 0; i < curve.number_of_status_lines_with_event(); 2460 i++) { 2461 out << curve.status_line_at_event(i) << std::endl; 2462 out << "Intermediate line at " 2463 << CGAL::to_double(curve.bound_value_in_interval(i+1)) 2464 << ": " << curve.arcs_over_interval(i+1) 2465 << " passing arcs" << std::endl 2466 << std::endl; 2467 } 2468 out << "(Horizontal) asymptotes at +infty: " << std::flush; 2469 size_type no_events = curve.number_of_status_lines_with_event(); 2470 for (size_type i = 0; i < curve.arcs_over_interval(no_events); i++) { 2471 2472 const Asymptote_y& curr_asym_info_obj 2473 = curve.asymptotic_value_of_arc(CGAL::RIGHT_BOUNDARY,i); 2474 typename Curve::Algebraic_real_1 curr_asym_info; 2475 bool is_finite = CGAL::assign(curr_asym_info,curr_asym_info_obj); 2476 if(! is_finite) { 2477 // Assignment to prevent compiler warning 2478 CGAL::Box_parameter_space_2 loc = CGAL::LEFT_BOUNDARY; 2479 CGAL_assertion_code(bool is_valid = ) 2480 CGAL::assign(loc, curr_asym_info_obj); 2481 CGAL_assertion(is_valid); 2482 if(loc == CGAL::TOP_BOUNDARY) { 2483 out << "+infty " << std::flush; 2484 } else { 2485 CGAL_assertion(loc == CGAL::BOTTOM_BOUNDARY); 2486 out << "-infty " << std::flush; 2487 } 2488 } else { // is_finite 2489 out << CGAL::to_double(curr_asym_info) 2490 << " " << std::flush; 2491 } 2492 } 2493 2494 out << std::endl; 2495 2496 out << "------------------------------------------------" << std::endl; 2497 break; 2498 } 2499 case ::CGAL::IO::BINARY: 2500 std::cerr << "BINARY format not yet implemented" << std::endl; 2501 break; 2502 default: 2503 // ASCII 2504 out << curve.polynomial_2(); 2505 } 2506 2507 return out; 2508 } 2509 2510 //! \brief reads the objects from stream 2511 template<typename AlgebraicKernelWithAnalysis_2, 2512 typename Rep_> 2513 std::istream& operator>> ( 2514 std::istream& is, 2515 Curve_analysis_2< AlgebraicKernelWithAnalysis_2, Rep_ >& curve) { 2516 2517 CGAL_precondition(CGAL::IO::is_ascii(is)); 2518 2519 typedef AlgebraicKernelWithAnalysis_2 Algebraic_kernel_with_analysis_2; 2520 2521 typedef Rep_ Rep; 2522 2523 typename Curve_analysis_2< Algebraic_kernel_with_analysis_2, Rep >:: 2524 Polynomial_2 f; 2525 2526 is >> f; 2527 2528 // TODO is get_static_instance the right way? 2529 curve = Algebraic_kernel_with_analysis_2::get_static_instance(). 2530 construct_curve_2_object()(f); 2531 2532 return is; 2533 } 2534 2535 2536 } //namespace CGAL 2537 2538 2539 #include <CGAL/enable_warnings.h> 2540 2541 #endif // ALGEBRAIC_CURVE_2_H 2542