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_pair_analysis_2.h $
7 // $Id: Curve_pair_analysis_2.h 26355e2 2020-06-25T12:31:21+02:00 Mael Rouxel-Labbé
8 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Eric Berberich <eric@mpi-inf.mpg.de>
12 // Michael Kerber <mkerber@mpi-inf.mpg.de>
13 //
14 // ============================================================================
15
16
17 #ifndef CGAL_ACK_CURVE_PAIR_ANALYSIS_H
18 #define CGAL_ACK_CURVE_PAIR_ANALYSIS_H 1
19
20 #include <CGAL/disable_warnings.h>
21
22 #include <vector>
23 #include <algorithm>
24
25 #include <boost/optional.hpp>
26
27 #include <CGAL/Handle_with_policy.h>
28 #include <CGAL/boost/iterator/transform_iterator.hpp>
29
30 #include <CGAL/Polynomial.h>
31 #include <CGAL/Polynomial_traits_d.h>
32
33 #include <CGAL/Algebraic_kernel_d/Shear_controller.h>
34 #include <CGAL/Algebraic_kernel_d/Shear_transformation.h>
35 #include <CGAL/Algebraic_kernel_d/enums.h>
36 #include <CGAL/Algebraic_kernel_d/exceptions.h>
37 #include <CGAL/Algebraic_kernel_d/Status_line_CPA_1.h>
38
39
40
41 namespace CGAL {
42
43 namespace internal {
44
45 template<class AlgebraicReal_1>
46 class Distinct_compare {
47
48 public:
49
50 typedef AlgebraicReal_1 Algebraic_real_1;
51
52 typedef ::CGAL::Comparison_result result_type;
53 typedef Algebraic_real_1 first_argument_type;
54 typedef Algebraic_real_1 second_argument_type;
55
operator()56 ::CGAL::Comparison_result operator()
57 (Algebraic_real_1 a,Algebraic_real_1 b) {
58 return a.compare_distinct(b);
59 }
60
61 };
62
63 }// namespace internal
64
65 //////////////////////////////////////////////////////////////////////////////
66 // Curve_pair_2
67
68 // Forwards
69 template < typename AlgebraicKernelWithAnalysis_2 >
70 class Curve_pair_analysis_2;
71
72 template<typename AlgebraicKernelWithAnalysis_2>
73 std::ostream& operator<<
74 (std::ostream&,const Curve_pair_analysis_2
75 <AlgebraicKernelWithAnalysis_2>&);
76
77 namespace internal {
78
79 // Internally used enums and structs
80
81 enum Slice_type {
82 FIRST_CURVE = 0,
83 SECOND_CURVE = 1,
84 INTERSECTION = 2,
85 CANDIDATE = 3
86 };
87
88
89 /*!
90 * An x-event of the curve pair is either a root of a dicriminant of a single
91 * curve, or a root of the resultant of both curves, or both.
92 * The \c Event_indices vector stores a triple <tt>(fg,ffy,ggy)</tt> denoting
93 * that some event is the <tt>fg</tt> root of <tt>res(f,g,y)</tt>,
94 * the <tt>ffy</tt>th root of <tt>disc(f,y)</tt> and
95 * the <tt>ggy</tt>th root of <tt>disc(g,y)</tt>.
96 */
97 template<typename size_type>
98 struct Event_indices {
99
100 size_type fg;
101 size_type ffy;
102 size_type ggy;
Event_indicesEvent_indices103 Event_indices(size_type fg,size_type ffy, size_type ggy)
104 : fg(fg), ffy(ffy), ggy(ggy) {}
105 };
106
107 // Representation class for curve pairs
108 template < class AlgebraicKernelWithAnalysis_2 >
109 class Curve_pair_analysis_2_rep {
110 public:
111
112 //! \name public typedefs
113 //! @{
114 typedef AlgebraicKernelWithAnalysis_2 Algebraic_kernel_with_analysis_2;
115
116 typedef Curve_pair_analysis_2_rep<Algebraic_kernel_with_analysis_2> Self;
117
118 typedef Curve_pair_analysis_2<Algebraic_kernel_with_analysis_2> Handle;
119
120 typedef typename Algebraic_kernel_with_analysis_2::Curve_analysis_2
121 Curve_analysis_2;
122
123 typedef typename Curve_analysis_2::size_type size_type;
124
125 typedef typename Curve_analysis_2::Polynomial_2 Polynomial_2;
126
127 typedef typename Curve_analysis_2::Algebraic_real_1 Algebraic_real_1;
128
129 typedef typename Polynomial_2::NT Polynomial_1;
130
131 typedef typename Curve_analysis_2::Bound Bound;
132
133 typedef CGAL::internal::Status_line_CPA_1<Handle> Status_line_CPA_1;
134
135 typedef std::pair<Slice_type,size_type> Slice_element;
136
137 typedef std::vector<Slice_element> Slice_info;
138
139 typedef boost::optional<Slice_info> Lazy_slice_info;
140
141 typedef boost::optional<Bound> Lazy_bound;
142
143 typedef CGAL::internal::Event_indices<size_type> Event_indices;
144
145 struct Intersection_info {
146 typename Curve_analysis_2::Status_line_1 ev;
147 size_type index;
148 size_type mult;
149 };
150
151 typedef std::vector<std::vector<Intersection_info> >
152 Intersection_info_container;
153
154 typedef boost::optional<Intersection_info_container>
155 Lazy_intersection_info_container;
156
157 // For lazy evaluation of Status_line_CPA_1s.
158 typedef boost::optional<Status_line_CPA_1> Lazy_status_line_CPA_1;
159
160 //! @}
161
162 //! \name Constructors
163 //! @{
164
165 // DefaultConstructible
Curve_pair_analysis_2_rep()166 Curve_pair_analysis_2_rep() :
167 c1_(), c2_() {
168 }
169
Curve_pair_analysis_2_rep(Algebraic_kernel_with_analysis_2 * kernel,Curve_analysis_2 c1,Curve_analysis_2 c2,CGAL::Degeneracy_strategy strategy)170 Curve_pair_analysis_2_rep(Algebraic_kernel_with_analysis_2 *kernel,
171 Curve_analysis_2 c1, Curve_analysis_2 c2,
172 CGAL::Degeneracy_strategy strategy) :
173 _m_kernel(kernel),
174 c1_(c1), c2_(c2), f(c1.polynomial_2()), g(c2.polynomial_2()),
175 degeneracy_strategy(strategy) {
176 }
177
178 //! @}
179
180 private:
181
182 //! \name members
183 //! @{
184
185 Algebraic_kernel_with_analysis_2* _m_kernel;
186
187 Curve_analysis_2 c1_;
188 Curve_analysis_2 c2_;
189
190 Polynomial_2 f;
191 Polynomial_2 g;
192
193
194 mutable boost::optional<std::vector<Polynomial_2> > subresultants;
195
196 mutable boost::optional<std::vector<Polynomial_1> >
197 principal_subresultants;
198 mutable boost::optional<std::vector<Polynomial_1> >
199 coprincipal_subresultants;
200
201 mutable boost::optional<Polynomial_1> resultant;
202
203 mutable boost::optional<std::vector<Algebraic_real_1> > resultant_roots;
204 mutable boost::optional<std::vector<Algebraic_real_1> >
205 event_x_coordinates;
206 mutable boost::optional<std::vector<size_type> >
207 multiplicities_of_resultant_roots;
208
209 mutable boost::optional<std::vector<Bound> > stripe_values;
210
211 mutable std::vector< Lazy_status_line_CPA_1 > event_slices;
212
213 mutable boost::optional<std::vector< Lazy_bound > > intermediate_values;
214
215 mutable boost::optional< std::vector< Lazy_status_line_CPA_1 > >
216 intermediate_slices;
217
218 mutable boost::optional<std::vector<Event_indices> > event_indices;
219
220 mutable Lazy_intersection_info_container intersection_info_container;
221
222 typedef typename Curve_analysis_2::Integer Integer;
223
224 CGAL::Degeneracy_strategy degeneracy_strategy;
225
226 mutable CGAL::internal::Shear_controller<Integer> shear_controller;
227
228 //! @}
229
230 //! \name friends
231 //! @{
232
233 friend class Curve_pair_analysis_2<Algebraic_kernel_with_analysis_2>;
234
235 //!@}
236
237 };
238
239 } // namespace internal
240
241 /*!
242 * A model for <tt>AlgebraicKernelWithAnalysis_2::CurvePairAnalysis_2</tt>
243 * It provides topological-geometric information about the intersection
244 * points, and the vertical order of arcs of two algebraic plane curves.
245 *
246 * The curve pair is passed by two \c Curve_analysis_2 instances.
247 * It is required that they do not share a component, i.e., the number
248 * of common points must be finite. Note that overlapping curves are handled
249 * by \c Algebraic_curve_kernel_2::Construct_curve_pair_2.
250 * Also for caching reasons, it is recommended to construct curve pairs
251 * always with this method.
252 *
253 * As for the single-curve analysis, the curve pair analysis is implemented
254 * in a "lazy" fashion. That means, any computation is triggered when
255 * the result is actually queried by the user. This prevents
256 * expensive symbolic computations in some cases.
257 *
258 * For all algorithmic details of the curve pair analysis, we refer to
259 * Arno Eigenwillig, Michael Kerber: Exact and Efficient 2D-Arrangements
260 * of Arbitrary Algebraic Curves. Proceedings of the Nineteenth Annual
261 * ACM-SIAM Symposium on Discrete Algorithms (SODA 2008), pp. 122-131
262 */
263 template < typename AlgebraicKernelWithAnalysis_2 >
264 class Curve_pair_analysis_2 :
265 public ::CGAL::Handle_with_policy
266 < CGAL::internal::Curve_pair_analysis_2_rep
267 < AlgebraicKernelWithAnalysis_2 > > {
268
269
270 public:
271
272 //! \name typedefs
273 //! @{
274
275 //! The algebraic kernel that uses the curve pair analysis
276 typedef AlgebraicKernelWithAnalysis_2 Algebraic_kernel_with_analysis_2;
277
278 private:
279
280 //! Representation class
281 typedef CGAL::internal::Curve_pair_analysis_2_rep
282 < Algebraic_kernel_with_analysis_2 > Rep;
283
284 //! Base class
285 typedef ::CGAL::Handle_with_policy< Rep > Base;
286
287 public:
288 //! The Curve_pair_analysis_2 itself
289 typedef Curve_pair_analysis_2<Algebraic_kernel_with_analysis_2> Self;
290
291 //! The corresponding Curve_analysis_2 class
292 typedef typename Rep::Curve_analysis_2 Curve_analysis_2;
293
294 //! Index type
295 typedef typename Rep::size_type size_type;
296
297 //! Univariate polynomials
298 typedef typename Rep::Polynomial_1 Polynomial_1;
299
300 //! Bivariate polynomials
301 typedef typename Rep::Polynomial_2 Polynomial_2;
302
303 //! Type for algebraic numbers (one dimension)
304 typedef typename Rep::Algebraic_real_1 Algebraic_real_1;
305
306 //! Type for points with algebraic coordinates
307 typedef typename Algebraic_kernel_with_analysis_2::Algebraic_real_2
308 Algebraic_real_2;
309
310 //! Bound type (for rational numbers)
311 typedef typename Rep::Bound Bound;
312
313 private:
314 // Optional for boundaries
315 typedef typename Rep::Lazy_bound Lazy_bound;
316
317 // Object to store information about intersection points
318 typedef typename Rep::Intersection_info_container
319 Intersection_info_container;
320
321 // Its lazy version
322 typedef typename Rep::Lazy_intersection_info_container
323 Lazy_intersection_info_container;
324
325 // Type for indices of events.
326 typedef typename Rep::Event_indices Event_indices;
327
328 // Integer type
329 typedef typename Curve_analysis_2::Integer Integer;
330
331 // Status line of single curve analysis
332 typedef typename Curve_analysis_2::Status_line_1 Status_line_CA_1;
333
334 // Coefficient type
335 typedef typename Curve_analysis_2::Coefficient Coefficient;
336
337 // Polynomial traits class
338 typedef CGAL::Polynomial_traits_d<Polynomial_2> Polynomial_traits_2;
339
340 // Polynomial traits class
341 typedef CGAL::Polynomial_traits_d<Polynomial_1> Polynomial_traits_1;
342
343 public:
344
345 //! The event slice object type
346 typedef typename Rep::Status_line_CPA_1 Status_line_CPA_1;
347
348 /*!
349 * Required by the concept. The name is not used internally
350 * to distinguish from one curve status_lines syntactically
351 */
352 typedef Status_line_CPA_1 Status_line_1;
353
354 private:
355
356 // Lazy version of status lines
357 typedef typename Rep::Lazy_status_line_CPA_1 Lazy_status_line_CPA_1;
358
359 // Coercion between Bound and Coefficient type
360 typedef CGAL::Coercion_traits<Bound, Coefficient> Coercion;
361
362 // The common supertype
363 typedef typename Coercion::Type Coercion_type;
364
365 // Polynomials over that supertype
366 typedef typename CGAL::Polynomial_traits_d<Polynomial_2>
367 ::template Rebind<Coercion_type,1>::Other::Type Poly_coer_1;
368
369 // Functor to isolate real roots of univariate polynomials
370 typedef typename Algebraic_kernel_with_analysis_2::Solve_1 Solve_1;
371
372 // Slice info objects
373 typedef typename Rep::Slice_info Slice_info;
374
375 // Lazy version
376 typedef typename Rep::Lazy_slice_info Lazy_slice_info;
377
378 //! @}
379
380 private:
381
382 //! \name Internal structs
383 //! @{
384
385 struct Curves_at_event_functor {
386
387 typedef size_type argument_type;
388 typedef CGAL::internal::Slice_type result_type;
389
Curves_at_event_functorCurves_at_event_functor390 Curves_at_event_functor(const Status_line_CPA_1& status_line)
391 : status_line(status_line)
392 {}
393
operatorCurves_at_event_functor394 CGAL::internal::Slice_type operator() (size_type i) const {
395 typedef typename Status_line_CPA_1::size_type
396 Status_line_size_type;
397 std::pair<Status_line_size_type,Status_line_size_type> pair =
398 status_line.curves_at_event(i);
399 CGAL_assertion(pair.first>=0 || pair.second >=0);
400 if(pair.first==-1) {
401 return CGAL::internal::SECOND_CURVE;
402 }
403 if(pair.second==-1) {
404 return CGAL::internal::FIRST_CURVE;
405 }
406 return CGAL::internal::INTERSECTION;
407 }
408 private:
409
410 const Status_line_CPA_1& status_line;
411
412 };
413
414 typedef boost::transform_iterator<Curves_at_event_functor,
415 boost::counting_iterator<size_type> >
416 Status_line_CPA_iterator;
417
418 struct Xval_of_status_line_CA_1 {
419 typedef Status_line_CA_1 argument_type;
420 typedef Algebraic_real_1 result_type;
operatorXval_of_status_line_CA_1421 Algebraic_real_1 operator() (const Status_line_CA_1& status_line)
422 const {
423 return status_line.x();
424 }
425 };
426
427 // @}
428
429 //! \name Constructors
430 //! @{
431
432 public:
433
434 //! DefaultConstructible
Curve_pair_analysis_2()435 Curve_pair_analysis_2() :
436 Base(Rep()) {
437 };
438
439 //! \brief Copy constructor
440 #ifdef DOXYGEN_RUNNING
Curve_pair_analysis_2(const Self & alg_curve_pair)441 Curve_pair_analysis_2(const Self& alg_curve_pair)
442 : Base(static_cast<const Base&>(alg_curve_pair))
443 {
444 }
445 #endif
446 // Assignable
447
448 /*!
449 * \brief Constructable from two curves
450 *
451 * Create a curve pair object for the two curves \c c1 and \c c2,
452 * given by their curve analysis object. The two curves are checked
453 * to have no common vertical line component (if they have, an
454 * exception of type \c CGAL::internal::Non_generic_position_exception
455 * is thrown), no further computation is performed.
456 *
457 * \param strategy If a degenerate situation (e.g., two covertical
458 * intersection at the same x-coordinate) occurs during the analysis,
459 * this value controls the strategy to handle it. If set to
460 * CGAL::EXCEPTION_STRATEGY, an exception of type
461 * \c CGAL::internal::Non_generic_position_exception is thrown whenever
462 * such a degeneracy occurs. If set to \c CGAL::SHEAR_STRATEGY, a shear
463 * transformation is performed, and the sheared curve pair is used
464 * to handle degenerate situations. Finally, if set to
465 * CGAL::SHEAR_ONLY_AT_IRRATIONAL_STRATEGY, degeneracies at rational
466 * x-ccordinates are handled directly, and a shear is only applied
467 * in other situations. The default argument for \c strategy is
468 * \c CGAL::SHEAR_ONLY_AT_IRRATIONAL_STRATEGY.
469 */
470 Curve_pair_analysis_2(Algebraic_kernel_with_analysis_2* kernel,
471 Curve_analysis_2 c1,
472 Curve_analysis_2 c2,
473 CGAL::Degeneracy_strategy strategy
474 = CGAL_ACK_DEFAULT_DEGENERACY_STRATEGY)
Base(Rep (kernel,c1,c2,strategy))475 : Base(Rep(kernel,c1, c2, strategy))
476 {
477
478 #if CGAL_ACK_DEBUG_FLAG
479 CGAL::set_pretty_mode(CGAL_ACK_DEBUG_PRINT);
480 #endif
481
482 #if CGAL_ACK_DEBUG_FLAG
483 CGAL_ACK_DEBUG_PRINT << "Check content for squarefreeness.."
484 << std::flush;
485 #endif
486 if(CGAL::degree(this->ptr()->c1_.content())>0 &&
487 CGAL::degree(this->ptr()->c2_.content())>0) {
488 typename Polynomial_traits_1::Gcd_up_to_constant_factor gcd_utcf;
489 if(CGAL::degree(gcd_utcf
490 (this->ptr()->c1_.content(),
491 this->ptr()->c2_.content())) >= 1) {
492
493 #if CGAL_ACK_DEBUG_FLAG
494 CGAL_ACK_DEBUG_PRINT << "Common vertical line discovered"
495 << std::endl;
496 #endif
497 throw CGAL::internal::Non_generic_position_exception();
498 } else {
499 #if CGAL_ACK_DEBUG_FLAG
500 CGAL_ACK_DEBUG_PRINT << "done" << std::endl;
501 #endif
502 }
503 }
504
505 }
506
507 //! @}
508
509
510 private:
511
512 // Computes the resultant of the defining polynomials wrt \c y
513 void compute_resultant() const;
514
515 // Computes the subresultant coefficients of the defining polynomials
516 void compute_subresultants() const;
517
518 /*
519 * Computes the roots of the resultants (via isolation) and their
520 * multiplicities
521 */
522 void compute_resultant_roots_with_multiplicities() const;
523
524 /*
525 * Computes all x-events of the curve pair,
526 * together with their event indices
527 */
528 void compute_event_x_coordinates_with_event_indices() const;
529
530 /*
531 * \brief Computes the intermediate x-coordinates and their status lines
532 *
533 * In fact, it only fills the data fields with boost::none instances,
534 * according to the lazy philosophy of the whole class.
535 */
536 void compute_intermediate_values_and_slices() const;
537
538 public:
539
kernel()540 Algebraic_kernel_with_analysis_2* kernel() const {
541 return this->ptr()->_m_kernel;
542 }
543
544 //! Returns the resultant of the defing polynomials wrt \c y
resultant()545 Polynomial_1 resultant() const {
546 if(! this->ptr()->resultant) {
547 compute_resultant();
548 }
549 CGAL_assertion(bool(this->ptr()->resultant));
550 return this->ptr()->resultant.get();
551 }
552
resultant_roots()553 std::vector<Algebraic_real_1>& resultant_roots() const {
554 if(! this->ptr()->resultant_roots) {
555 compute_resultant_roots_with_multiplicities();
556 }
557 CGAL_assertion(bool(this->ptr()->resultant_roots));
558 return this->ptr()->resultant_roots.get();
559 }
560
resultant_roots(size_type i)561 Algebraic_real_1& resultant_roots(size_type i) const {
562 CGAL_assertion(i>=0 &&
563 i < static_cast<size_type>(resultant_roots().size()));
564 return resultant_roots()[i];
565 }
566
multiplicities_of_resultant_roots()567 std::vector<size_type>& multiplicities_of_resultant_roots() const {
568 if(! this->ptr()->multiplicities_of_resultant_roots) {
569 compute_resultant_roots_with_multiplicities();
570 }
571 CGAL_assertion(bool(this->ptr()->multiplicities_of_resultant_roots));
572 return this->ptr()->multiplicities_of_resultant_roots.get();
573 }
574
multiplicities_of_resultant_roots(size_type i)575 size_type multiplicities_of_resultant_roots(size_type i) const {
576 CGAL_assertion(i>=0 &&
577 i < static_cast<size_type>
578 (multiplicities_of_resultant_roots().size()));
579 return multiplicities_of_resultant_roots()[i];
580 }
581
stripe_values()582 std::vector<Bound>& stripe_values() const {
583 if(! this->ptr()->stripe_values) {
584 this->ptr()->stripe_values = std::vector<Bound>();
585 find_intermediate_values
586 (kernel(),
587 resultant_roots().begin(),
588 resultant_roots().end(),
589 std::back_inserter(this->ptr()->stripe_values.get()));
590 }
591 CGAL_assertion(bool(this->ptr()->stripe_values));
592 return this->ptr()->stripe_values.get();
593 }
594
event_x_coordinates()595 std::vector<Algebraic_real_1>& event_x_coordinates() const {
596 if(! this->ptr()->event_x_coordinates) {
597 compute_event_x_coordinates_with_event_indices();
598 }
599 CGAL_assertion(bool(this->ptr()->event_x_coordinates));
600 return this->ptr()->event_x_coordinates.get();
601 }
602
event_indices()603 std::vector<Event_indices>& event_indices() const {
604 if(! this->ptr()->event_indices) {
605 compute_event_x_coordinates_with_event_indices();
606 }
607 CGAL_assertion(bool(this->ptr()->event_indices));
608 return this->ptr()->event_indices.get();
609 }
610
611 public:
612
613 /*
614 * \brief returns the indices of the <tt>i</tt>th event value
615 *
616 * Returns a Event_indices <tt>(fg,ffy,ggy)</tt> such that
617 * the <tt>i</tt>th event root is the <tt>fg</tt>th root of the
618 * resultant of \c f and \c g, the <tt>ffy</tt>th root of the
619 * discriminant of \c f, and the <tt>ggy</tt>th root of the
620 * discriminant of \c g.
621 */
event_indices(size_type i)622 Event_indices event_indices(size_type i) const {
623 CGAL_assertion(i>=0 &&
624 i < static_cast<size_type>
625 (event_indices().size()));
626 return event_indices()[i];
627 }
628
629 private:
630
intermediate_values()631 std::vector<Lazy_bound>& intermediate_values() const {
632 if(! this->ptr()->intermediate_values) {
633 compute_intermediate_values_and_slices();
634 }
635 CGAL_assertion(bool(this->ptr()->intermediate_values));
636 return this->ptr()->intermediate_values.get();
637 }
638
intermediate_slices()639 std::vector<Lazy_status_line_CPA_1>& intermediate_slices() const {
640 if(! this->ptr()->intermediate_slices) {
641 compute_intermediate_values_and_slices();
642 }
643 CGAL_assertion(bool(this->ptr()->intermediate_slices));
644 return this->ptr()->intermediate_slices.get();
645 }
646
647
648 private:
649
subresultants()650 std::vector<Polynomial_2>& subresultants() const {
651 if(! this->ptr()->subresultants) {
652 compute_subresultants();
653 }
654 CGAL_assertion(bool(this->ptr()->subresultants));
655 return this->ptr()->subresultants.get();
656 }
657
subresultants(size_type i)658 Polynomial_2& subresultants(size_type i) const {
659 CGAL_assertion(i>=0 &&
660 i < static_cast<size_type>(subresultants().size()));
661 return subresultants()[i];
662 }
663
principal_subresultants()664 std::vector<Polynomial_1>& principal_subresultants() const {
665 if(! this->ptr()->principal_subresultants) {
666 compute_subresultants();
667 }
668 CGAL_assertion(bool(this->ptr()->principal_subresultants));
669 return this->ptr()->principal_subresultants.get();
670 }
671
principal_subresultants(size_type i)672 Polynomial_1& principal_subresultants(size_type i) const {
673 CGAL_assertion(i>=0 &&
674 i < static_cast<size_type>
675 (principal_subresultants().size()));
676 return principal_subresultants()[i];
677 }
678
coprincipal_subresultants()679 std::vector<Polynomial_1>& coprincipal_subresultants() const {
680 if(! this->ptr()->coprincipal_subresultants) {
681 compute_subresultants();
682 }
683 CGAL_assertion(bool(this->ptr()->coprincipal_subresultants));
684 return this->ptr()->coprincipal_subresultants.get();
685 }
686
coprincipal_subresultants(size_type i)687 Polynomial_1& coprincipal_subresultants(size_type i) const {
688 CGAL_assertion(i>=0 &&
689 i < static_cast<size_type>
690 (coprincipal_subresultants().size()));
691 return coprincipal_subresultants()[i];
692 }
693
694
695
696 private:
697
698 /*
699 * Refines the isolating intervals until they are disjoint
700 * Returns CGAL::SMALLER, if the y-coordinate defined by <tt>(e1,i1)</tt>
701 * is smaller than the y-coordinate <tt>(e2,i2)</tt>,
702 * and CGAL::GREATER otherwise
703 *
704 * If both y-coordinates are equal, this method does not terminate
705 */
split_compare(Status_line_CA_1 & e1,size_type i1,Status_line_CA_1 & e2,size_type i2)706 CGAL::Sign split_compare(Status_line_CA_1& e1, size_type i1,
707 Status_line_CA_1& e2, size_type i2) const {
708 while(overlap(e1,i1,e2,i2)) {
709 if(e1.interval_length(i1)<e2.interval_length(i2)) {
710 e2.refine(i2);
711 }
712 else {
713 e1.refine(i1);
714 }
715 }
716 return (e1.lower_bound(i1) < e2.lower_bound(i2))
717 ? CGAL::SMALLER
718 : CGAL::LARGER;
719 }
720
721 private:
722
723 /*!
724 * TODO doc
725 */
create_event_slice(size_type i)726 Status_line_CPA_1 create_event_slice(size_type i)
727 const {
728 #if !CGAL_ACK_NO_ARC_FLIP
729 size_type index_in_fg = event_indices(i).fg;
730 if(index_in_fg == -1 ) {
731 return create_slice_with_multiplicity_zero_or_one(i);
732 } else {
733 size_type mult_of_alpha
734 = multiplicities_of_resultant_roots(index_in_fg);
735 if(mult_of_alpha == 1) {
736 return create_slice_with_multiplicity_zero_or_one(i);
737 } else {
738 #endif
739 return create_slice_of_higher_multiplicity(i);
740 #if !CGAL_ACK_NO_ARC_FLIP
741 }
742 }
743 #endif
744 }
745
create_slice_of_higher_multiplicity(size_type i)746 Status_line_CPA_1 create_slice_of_higher_multiplicity(size_type i)
747 const {
748 bool is_resultant_root = event_indices(i).fg >=0;
749 if(is_resultant_root &&
750 this->ptr()->intersection_info_container) {
751 return create_event_slice_with_shear(i);
752 }
753 try {
754 Status_line_CPA_1 slice = construct_generic_case(i);
755
756 return slice;
757 } catch(CGAL::internal::Non_generic_position_exception ex) {
758 switch(this->ptr()->degeneracy_strategy) {
759 case(CGAL::EXCEPTION_STRATEGY): {
760 throw ex;
761 break;
762 }
763 case(CGAL::SHEAR_ONLY_AT_IRRATIONAL_STRATEGY): {
764 if(event_x(i).is_rational()) {
765 return create_event_slice_at_rational(i);
766 }
767
768 CGAL_FALLTHROUGH;
769 }
770 case(CGAL::SHEAR_STRATEGY): {
771 return create_event_slice_with_shear(i);
772 }
773 }
774
775
776 // NEVER HAPPENS
777 return Status_line_CPA_1();
778
779 }
780 }
781
782 private:
create_event_slice_at_rational(size_type i)783 Status_line_CPA_1 create_event_slice_at_rational(size_type i) const {
784
785 Algebraic_real_1& x = event_x(i);
786
787 CGAL_precondition(x.is_rational());
788 Bound r = x.rational();
789
790 int k = degree_of_local_gcd(event_indices(i).fg,x);
791 Polynomial_2 sres = subresultants(k);
792
793 Polynomial_1 gcd = kernel()->evaluate_utcf_2_object()
794 (typename Polynomial_traits_2::Swap() (sres,0,1),r);
795 std::vector<Algebraic_real_1> gcd_roots;
796 kernel()->solve_1_object()(gcd,std::back_inserter(gcd_roots),false);
797 size_type m = static_cast<size_type>(gcd_roots.size());
798
799 Slice_info slice_info = construct_slice_info(x);
800 reduce_number_of_candidates_and_intersections_to
801 (m,
802 this->ptr()->c1_.status_line_at_exact_x(x),
803 this->ptr()->c2_.status_line_at_exact_x(x),
804 slice_info);
805 for(typename Slice_info::iterator it=slice_info.begin();
806 it!=slice_info.end();
807 it++) {
808
809 if(it->first==CGAL::internal::CANDIDATE) {
810 it->first=CGAL::internal::INTERSECTION;
811 }
812 }
813
814 return create_slice_from_slice_info(i,slice_info,true);
815 }
816
817 private:
818
819 /*!
820 * TODO doc
821 */
822 Status_line_CPA_1 create_slice_with_multiplicity_zero_or_one(size_type i)
823 const;
824
825 private:
826
827 // Creates an intermediate slice at a rational value
828 Status_line_CPA_1 create_intermediate_slice_at(int i) const;
829
830 private:
831
832 // Create a slice with id \c id from the Slice_info object
833 Status_line_CPA_1 create_slice_from_slice_info(size_type id,
834 const Slice_info& slice,
835 bool event_flag) const;
836
837 private:
838
839 // Computes a slice_info object at Algebraic_real_1 \c alpha
840 Slice_info construct_slice_info(Algebraic_real_1 alpha) const;
841
842 private:
843
844 Status_line_CPA_1 construct_generic_case(size_type i) const;
845
846 private:
847
848 bool check_candidate_by_arc_pattern(size_type index,
849 Status_line_CA_1& e1,
850 size_type i1,
851 Status_line_CA_1& e2,
852 size_type i2) const;
853 private:
854
855 /*
856 * TODO update doc
857 * Checks the point on e1 with index i1, and
858 * the point on e2 with index i2 really intersect. The \c slice_info
859 * is updated accordingly: If not intersecting, the corresponding
860 * points are refined until they can be arranged in the correct order.
861 * If intersecting, the corresponding Slice_info element is set to
862 * INTERSECTION.
863 */
864 template<typename InputIterator>
865 void check_candidate(Status_line_CA_1& e1,size_type i1,
866 Status_line_CA_1& e2,size_type i2,
867 size_type k,
868 Slice_info& slice_info,
869 InputIterator slice_it,
870 size_type root_index) const;
871
872 private:
873
874 /*
875 * Checks intersection with symbolic methods
876 */
check_candidate_symbolically(Status_line_CA_1 & e1,size_type,Status_line_CA_1 & CGAL_precondition_code (e2),size_type,size_type k)877 bool check_candidate_symbolically(Status_line_CA_1& e1,size_type ,
878 Status_line_CA_1& CGAL_precondition_code(e2),size_type ,
879 size_type k) const {
880 Polynomial_1 p = -coprincipal_subresultants(k-1);
881 Polynomial_1 q = principal_subresultants(k)*Coefficient(k);
882 Algebraic_real_1 alpha = e1.x();
883 CGAL_assertion(alpha==e2.x());
884 if(CGAL::internal::zero_test_bivariate
885 <Algebraic_kernel_with_analysis_2>
886 (kernel(),alpha,this->ptr()->f,p,q) &&
887 CGAL::internal::zero_test_bivariate
888 <Algebraic_kernel_with_analysis_2>
889 (kernel(),alpha,this->ptr()->g,p,q)) {
890 return true;
891 }
892 else {
893 throw CGAL::internal::Non_generic_position_exception();
894 }
895 return false; // never happens
896 }
897
898 private:
899
900 /*
901 * Checks whether the isolting intervals for the point on \c e1 with
902 * index \c index1, and for the point on \c e2 with index \c index2
903 * overlap
904 */
overlap(Status_line_CA_1 & e1,size_type index1,Status_line_CA_1 & e2,size_type index2)905 bool overlap(Status_line_CA_1& e1,
906 size_type index1,
907 Status_line_CA_1& e2,
908 size_type index2) const {
909 if(e1.lower_bound(index1) > e2.upper_bound(index2)) {
910 return false;
911 }
912 else if(e1.upper_bound(index1) < e2.lower_bound(index2)) {
913 return false;
914 }
915 else {
916 return true;
917 }
918 }
919
920 /*
921 * For the point \c p on \c e1 with index \c index1, find the
922 * unique point on \c e2 which might be equal to \c p. If no point
923 * can be equal, -1 is returned.
924 */
925 size_type find_possible_matching(Status_line_CA_1& e1,
926 size_type index1,
927 Status_line_CA_1& e2) const;
928
929
degree_of_local_gcd(size_type index_of_fg,Algebraic_real_1 alpha)930 size_type degree_of_local_gcd(size_type index_of_fg,
931 Algebraic_real_1 alpha) const {
932
933 if(multiplicities_of_resultant_roots(index_of_fg) == 1) {
934 return 1;
935 } else {
936 size_type k=1;
937 while(kernel()->is_zero_at_1_object()
938 (principal_subresultants(k),alpha)) {
939 k++;
940 }
941 return k;
942 }
943 }
944
945 public:
946
947 //! Returns curve analysis for the cth curve
curve_analysis(bool c)948 Curve_analysis_2 curve_analysis(bool c) const {
949 return c ? this->ptr()->c2_ : this->ptr()->c1_;
950 }
951
event_of_curve_analysis(size_type i,bool c)952 size_type event_of_curve_analysis(size_type i, bool c) const {
953 Event_indices& ev_ind = event_indices(i);
954 return c ? ev_ind.ggy : ev_ind.ffy;
955 }
956
event_of_curve_analysis(size_type i,const Curve_analysis_2 & c)957 size_type event_of_curve_analysis(size_type i,
958 const Curve_analysis_2& c) const {
959 CGAL_assertion(c.id()==curve_analysis(false).id() ||
960 c.id()==curve_analysis(true).id());
961 Event_indices& ev_ind = event_indices(i);
962 return (c.id()==curve_analysis(false).id()) ? ev_ind.ffy : ev_ind.ggy;
963 }
964
965 /*!
966 * \brief returns the number of event slices
967 *
968 * Precisely, this is the number of points which are either root of
969 * the resultant of the two curves, or root of discriminant of one
970 * of the curves
971 */
number_of_status_lines_with_event()972 size_type number_of_status_lines_with_event() const {
973 return static_cast<size_type>(event_x_coordinates().size());
974 }
975
976 //! Returns the x-coordinate of the <tt>i</tt>th event
event_x(size_type i)977 Algebraic_real_1& event_x(size_type i) const {
978 CGAL_assertion(i>=0 &&
979 i<static_cast<size_type>(event_x_coordinates().size()));
980 return event_x_coordinates()[i];
981 }
982
983 /*!
984 * \brief The index of the x-coordinate
985 *
986 * For x-value \c x, the index of the suitable slice is computed. For
987 * event value, the \c event flag is set to true, otherwise to false
988 * and the slice of the interval to which \c x belongs is returned
989 */
x_to_index(Algebraic_real_1 x,size_type & idx,bool & event)990 void x_to_index(Algebraic_real_1 x,
991 size_type& idx, bool& event) const {
992 const std::vector<Algebraic_real_1>& sl = event_x_coordinates();
993 idx = static_cast<size_type>(std::lower_bound(sl.begin(),
994 sl.end(),
995 x) - sl.begin());
996 event = (idx < static_cast<size_type>(sl.size()) && (sl[idx] == x));
997
998 }
999
1000
1001 Status_line_CPA_1 status_line_for_x(Algebraic_real_1 x,
1002 CGAL::Sign perturb = CGAL::ZERO)
1003 const {
1004 size_type index;
1005 bool evt;
1006 x_to_index(x,index,evt);
1007 if(evt) {
1008 switch(perturb) {
1009 case(CGAL::ZERO): return status_line_at_event(index);
1010 case(CGAL::NEGATIVE): return status_line_of_interval(index);
1011 case(CGAL::POSITIVE): return status_line_of_interval(index+1);
1012 }
1013 } // else:
1014 return status_line_of_interval(index);
1015
1016
1017 }
1018
1019
status_line_at_exact_x(Algebraic_real_1 x)1020 Status_line_CPA_1 status_line_at_exact_x(Algebraic_real_1 x) {
1021 return status_line_for_x(x);
1022 }
1023
1024
1025 public:
1026
1027 //! Returns the Status_line_CPA_1 at the <tt>i</tt>th event
status_line_at_event(size_type i)1028 const Status_line_CPA_1& status_line_at_event(size_type i) const {
1029 if(! this->ptr()->event_slices[i]) {
1030 this->ptr()->event_slices[i] = create_event_slice(i);
1031 }
1032 CGAL_assertion(bool(this->ptr()->event_slices[i]));
1033 return this->ptr()->event_slices[i].get();
1034 }
1035
1036
1037
1038 //! Returns the Status_line_CPA_1 at the <tt>i</tt>th interval
status_line_of_interval(size_type i)1039 const Status_line_CPA_1& status_line_of_interval(size_type i) const {
1040
1041 if(! intermediate_slices()[i]) {
1042
1043 intermediate_slices()[i]
1044 = create_intermediate_slice_at(i);
1045
1046 }
1047
1048 return intermediate_slices()[i].get();
1049 }
1050
1051 //! Returns bound representative value at the <tt>i</tt>th interval
bound_value_in_interval(size_type i)1052 const Bound bound_value_in_interval(size_type i) const {
1053
1054 const std::vector<Algebraic_real_1>& events = event_x_coordinates();
1055
1056 if(! intermediate_values()[i]) {
1057 // Create the intermediate x-coordinate first
1058 if(events.size()==0) {
1059 CGAL_assertion(i==0);
1060 intermediate_values()[0]=Bound(0);
1061 } else {
1062 if(i==0) {
1063 intermediate_values()[i]
1064 = bound_left_of(kernel(),events[i]);
1065 } else if(i == static_cast<size_type>(events.size())) {
1066 intermediate_values()[i]
1067 = bound_right_of(kernel(),events[i-1]);
1068
1069 } else {
1070 intermediate_values()[i]
1071 = kernel()->bound_between_1_object()
1072 (events[i-1],events[i]);
1073 }
1074 }
1075 }
1076 CGAL_assertion(bool(intermediate_values()[i]));
1077 return intermediate_values()[i].get();
1078
1079 }
1080
1081 private:
1082
1083 struct Bound_to_coercion_functor {
1084
1085 typedef Bound argument_type;
1086 typedef Coercion_type result_type;
1087
operatorBound_to_coercion_functor1088 result_type operator() (argument_type x) const {
1089 typename CGAL::Coercion_traits<Bound,Coefficient>::Cast cast;
1090 return cast(x);
1091 }
1092 };
1093
1094 struct Coefficient_to_coercion_functor {
1095
1096 typedef Coefficient argument_type;
1097 typedef Coercion_type result_type;
1098
operatorCoefficient_to_coercion_functor1099 result_type operator() (argument_type x) const {
1100 typename CGAL::Coercion_traits<Bound,Coefficient>::Cast cast;
1101 return cast(x);
1102 }
1103 };
1104
1105
1106 // If a new shear was used, update intersection multiplicities
merge_new_intersection_info(const Intersection_info_container & new_info_container)1107 void merge_new_intersection_info
1108 (const Intersection_info_container& new_info_container) const {
1109 if(! this->ptr()->intersection_info_container) {
1110 // ok, nothing existed, so take the new intersection info
1111 this->ptr()->intersection_info_container
1112 = new_info_container;
1113 return;
1114 }
1115 Intersection_info_container& old_info_container
1116 = *(this->ptr()->intersection_info_container);
1117 std::size_t n = old_info_container.size();
1118 CGAL_assertion(n == new_info_container.size());
1119 //iterate through the vector and update
1120 // (-1 stands for "multiplicity unknown")
1121 for(std::size_t i=0;i<n;i++) {
1122 std::size_t m = old_info_container[i].size();
1123 CGAL_assertion(m == new_info_container[i].size());
1124 for(std::size_t j=0;j<m;j++) {
1125 old_info_container[i][j].mult
1126 = (std::max)(new_info_container[i][j].mult,
1127 old_info_container[i][j].mult);
1128 }
1129
1130 }
1131 }
1132
1133 void new_shear_for_intersection_info
1134 (Intersection_info_container& info_container) const;
1135
create_event_slice_with_shear(size_type i)1136 Status_line_CPA_1 create_event_slice_with_shear(size_type i) const {
1137 while(true) { // we know that it works at some point
1138 try {
1139 if(! this->ptr()->intersection_info_container) {
1140 Intersection_info_container info_container;
1141 new_shear_for_intersection_info(info_container);
1142 merge_new_intersection_info(info_container);
1143 }
1144 Status_line_CPA_1 slice
1145 = create_event_slice_from_current_intersection_info(i);
1146
1147 return slice;
1148 } catch(CGAL::internal::Non_generic_position_exception /* ex */) {
1149 // just try the next one
1150 Intersection_info_container info_container;
1151 new_shear_for_intersection_info(info_container);
1152 merge_new_intersection_info(info_container);
1153 }
1154 }
1155 }
1156
1157
1158 Status_line_CPA_1
1159 create_event_slice_from_current_intersection_info (size_type i) const;
1160
x_sheared(Bound x,Bound y,Integer sh)1161 Bound x_sheared(Bound x, Bound y,Integer sh) const {
1162 return x-sh*y;
1163 }
1164
1165 void update_intersection_info(Intersection_info_container&
1166 info_container,
1167 Self& sh_pair,
1168 Status_line_CPA_1 slice,
1169 size_type i,
1170 size_type j,
1171 Integer s) const;
1172
1173 /*
1174 * \brief reduces the number of possible intersections
1175 *
1176 * At the position given by the event lins \c e1 and \c e2 and the slice
1177 * info object \c slice, the points on the event lines are further refined
1178 * until there are only \c n possible intersection points. The method can
1179 * be interrupted if all possible intersection points are known to have
1180 * a maximal intersection mulipicity smaller \c k, and a
1181 * Non_generic_position_exception is thrown then.
1182 */
1183 size_type reduce_number_of_candidates_and_intersections_to
1184 (size_type n,
1185 Status_line_CA_1& e1,
1186 Status_line_CA_1& e2,
1187 Slice_info& slice,
1188 size_type k=-1) const;
1189
1190 // Handle provides
1191 // .id()
1192 // .is_identical
1193
1194 friend std::ostream& operator<< <>
1195 (std::ostream& out,
1196 const Self& curve_pair);
1197
1198 }; // end of Curve_pair_analysis_2
1199
1200 //! \brief prints the objects.
1201 template<typename AlgebraicKernelWithAnalysis_2>
1202 std::ostream& operator<<
1203 (std::ostream& out,
1204 const Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>& curve_pair) {
1205 typedef Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>
1206 Curve_pair_analysis_2;
1207 typedef typename Curve_pair_analysis_2::size_type size_type;
1208 typedef typename Curve_pair_analysis_2::Event_indices Event_indices;
1209 typedef typename Curve_pair_analysis_2::Status_line_CPA_1 Slice;
1210 out << "--------------- Analysis results ---------------" << std::endl;
1211 out << "Number of constructed event lines: "
1212 << curve_pair.number_of_status_lines_with_event()
1213 << std::endl;
1214
1215 out << "Intermediate line: " << std::flush;
1216 Slice slice=curve_pair.status_line_of_interval(0);
1217
1218 out << slice.number_of_events() << " passing arcs" << std::endl ;
1219 out << "in order: " << std::flush;
1220 for(size_type i=0;i<slice.number_of_events();i++) {
1221 CGAL_assertion(slice.curves_at_event(i).first==-1 ||
1222 slice.curves_at_event(i).second==-1 );
1223 if(slice.curves_at_event(i).second==-1) {
1224 out << "First " <<std::flush;
1225 } else {
1226 out << "Second " <<std::flush;
1227 }
1228 }
1229 out << std::endl << std::endl;
1230 for(size_type j = 0;
1231 j < curve_pair.number_of_status_lines_with_event();
1232 j++) {
1233
1234 out << "Event line at " << CGAL::to_double(curve_pair.event_x(j))
1235 << ": " << std::endl;
1236 out << "Indices: ";
1237 Event_indices ev_ind = curve_pair.event_indices(j);
1238 out << "fg: " << ev_ind.fg << ", ffy: "
1239 << ev_ind.ffy <<", ggy: " << ev_ind.ggy
1240 << std::endl;
1241 slice = curve_pair.status_line_at_event(j);
1242 out << slice.number_of_events() << " passing arcs" << std::endl ;
1243 out << "in order: " << std::flush;
1244 for(size_type i=0;i<slice.number_of_events();i++) {
1245 if(slice.curves_at_event(i).second==-1) {
1246 out << "First " <<std::flush;
1247 } else if(slice.curves_at_event(i).first==-1) {
1248 out << "Second " <<std::flush;
1249 } else {
1250 out << "Inter," << slice.multiplicity_of_intersection(i)
1251 << " " << std::flush;
1252 }
1253
1254 }
1255 out << std::endl << std::endl;
1256
1257 out << "Intermediate line:" << std::flush;
1258 Slice slice=curve_pair.status_line_of_interval(j+1);
1259 out << slice.number_of_events() << " passing arcs" << std::endl ;
1260 out << "in order: " << std::flush;
1261 for(size_type i=0;i<slice.number_of_events();i++) {
1262 CGAL_assertion(slice.curves_at_event(i).first==-1 ||
1263 slice.curves_at_event(i).second==-1 );
1264 if(slice.curves_at_event(i).second==-1) {
1265 out << "First " <<std::flush;
1266 } else {
1267 out << "Second " <<std::flush;
1268 }
1269
1270 }
1271 out << std::endl << std::endl;
1272 }
1273 out << "------------------------------------------------" << std::endl;
1274
1275 return out;
1276 }
1277
1278
1279
1280 // Implementation of functions from Curve_pair_analysis class
1281
1282 //////////////////// compute_resultant()
1283
1284 template <typename AlgebraicKernelWithAnalysis_2>
compute_resultant()1285 void Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::compute_resultant()
1286 const {
1287
1288 #if CGAL_ACK_RESULTANT_FIRST_STRATEGY
1289 #ifndef CGAL_ACK_RESULTANT_FIRST_STRATEGY_DEGREE_THRESHOLD
1290 bool speed_up = true;
1291 #else
1292 bool speed_up = (std::min)
1293 (CGAL::degree(curve_analysis(false).polynomial_2(),1),
1294 CGAL::degree(curve_analysis(true).polynomial_2(),1)) >=
1295 CGAL_ACK_RESULTANT_FIRST_STRATEGY_DEGREE_THRESHOLD;
1296 #endif
1297 #else
1298 bool speed_up=false;
1299 #endif
1300 if(speed_up) {
1301 #if CGAL_ACK_DEBUG_FLAG
1302 CGAL_ACK_DEBUG_PRINT << "Compute the resultant of f and g..."
1303 << std::flush;
1304 #endif
1305 this->ptr()->resultant
1306 = CGAL::resultant(this->ptr()->f,this->ptr()->g);
1307 } else {
1308 #if CGAL_ACK_DEBUG_FLAG
1309 CGAL_ACK_DEBUG_PRINT << "Compute the subres-seq of f and g..."
1310 << std::flush;
1311 #endif
1312 compute_subresultants();
1313
1314 this->ptr()->resultant
1315 = this->ptr()->principal_subresultants.get()[0];
1316 }
1317
1318
1319 if(this->ptr()->resultant.get().is_zero()) {
1320 throw CGAL::internal::Zero_resultant_exception<Polynomial_2>
1321 (this->ptr()->f,
1322 this->ptr()->g);
1323 }
1324 #if CGAL_ACK_DEBUG_FLAG
1325 CGAL_ACK_DEBUG_PRINT << "done" << std::endl;
1326 #endif
1327
1328 }
1329
1330 //////////////////// compute_resultant_roots_with_multiplicities()
1331
1332 template<typename AlgebraicKernelWithAnalysis_2>
1333 void Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
compute_resultant_roots_with_multiplicities()1334 compute_resultant_roots_with_multiplicities() const {
1335
1336 #if CGAL_ACK_DEBUG_FLAG
1337 CGAL_ACK_DEBUG_PRINT << "Isolate the real roots of resultant..."
1338 << std::flush;
1339 #endif
1340 Solve_1 solve_1;
1341 this->ptr()->resultant_roots = std::vector<Algebraic_real_1>();
1342 this->ptr()->multiplicities_of_resultant_roots
1343 = std::vector<size_type>();
1344 std::vector<std::pair<Algebraic_real_1, size_type> > res_pairs;
1345 solve_1(resultant(), std::back_inserter(res_pairs));
1346
1347 for(int i=0; i < static_cast<int>(res_pairs.size()); i++ ) {
1348 this->ptr()->resultant_roots.get().push_back(res_pairs[i].first);
1349 this->ptr()->multiplicities_of_resultant_roots.get()
1350 .push_back(res_pairs[i].second);
1351 }
1352
1353 #if CGAL_ACK_DEBUG_FLAG
1354 CGAL_ACK_DEBUG_PRINT << "done" << std::endl;
1355 #endif
1356
1357 #if CGAL_ACK_DEBUG_FLAG
1358 for(size_type i = 0;
1359 i<static_cast<size_type>
1360 (this->ptr()->resultant_roots.get().size());
1361 i++) {
1362 CGAL_ACK_DEBUG_PRINT
1363 << "Root at "
1364 << CGAL::to_double(this->ptr()->resultant_roots.get()[i])
1365 << " with multiplicity "
1366 << this->ptr()->multiplicities_of_resultant_roots.get()[i]
1367 << std::endl;
1368 }
1369 #endif
1370
1371 }
1372
1373 //////////////////// compute_event_x_coordinates_with_event_indices
1374
1375 template<typename AlgebraicKernelWithAnalysis_2>
1376 void Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
compute_event_x_coordinates_with_event_indices()1377 compute_event_x_coordinates_with_event_indices() const {
1378
1379 Xval_of_status_line_CA_1 xval;
1380 const Curve_analysis_2& c1=this->ptr()->c1_, c2=this->ptr()->c2_;
1381
1382 std::vector<Algebraic_real_1> one_curve_events;
1383
1384 std::vector<CGAL::internal::Three_valued> one_curve_events_type;
1385
1386 typename CGAL::Real_embeddable_traits<Algebraic_real_1>::Compare compare;
1387
1388 CGAL::internal::set_union_with_source
1389 (::boost::make_transform_iterator(c1.event_begin(),xval),
1390 ::boost::make_transform_iterator(c1.event_end(),xval),
1391 ::boost::make_transform_iterator(c2.event_begin(),xval),
1392 ::boost::make_transform_iterator(c2.event_end(),xval),
1393 std::back_inserter(one_curve_events),
1394 std::back_inserter(one_curve_events_type),
1395 compare);
1396
1397 this->ptr()->event_x_coordinates = std::vector<Algebraic_real_1>();
1398 std::vector<CGAL::internal::Three_valued> events_type;
1399 CGAL::internal::set_union_with_source
1400 (one_curve_events.begin(),
1401 one_curve_events.end(),
1402 resultant_roots().begin(),
1403 resultant_roots().end(),
1404 std::back_inserter(this->ptr()->event_x_coordinates.get()),
1405 std::back_inserter(events_type),
1406 compare);
1407 std::vector<Algebraic_real_1>& events
1408 = this->ptr()->event_x_coordinates.get();
1409
1410 typename std::vector<CGAL::internal::Three_valued>::iterator one_curve_it
1411 =one_curve_events_type.begin();
1412 size_type inter_count=0, f_count=0,g_count=0;
1413 this->ptr()->event_indices = std::vector<Event_indices>();
1414 std::vector<Event_indices>& event_indices
1415 = this->ptr()->event_indices.get();
1416 for(size_type i=0;i<static_cast<size_type>(events.size());i++) {
1417 /*
1418 #if CGAL_ACK_DEBUG_FLAG
1419 CGAL_ACK_DEBUG_PRINT << CGAL::to_double(events[i]) << std::flush;
1420 #endif
1421 */
1422 switch(events_type[i]) {
1423 case(CGAL::internal::ROOT_OF_FIRST_SET): {
1424 /*
1425 #if CGAL_ACK_DEBUG_FLAG
1426 CGAL_ACK_DEBUG_PRINT << " one curve event" << std::endl;
1427 #endif
1428 */
1429 // Fix a warning by using `emplace_back()` instead of
1430 // copying a non-initialized `optional
1431 this->ptr()->event_slices.emplace_back();
1432 switch(*(one_curve_it++)) {
1433 case(CGAL::internal::ROOT_OF_FIRST_SET): {
1434 event_indices.push_back(Event_indices(-1,f_count,-1));
1435 f_count++;
1436 break;
1437 }
1438 case(CGAL::internal::ROOT_OF_SECOND_SET): {
1439 event_indices.push_back(Event_indices(-1,-1,g_count));
1440 g_count++;
1441 break;
1442 }
1443 case(CGAL::internal::ROOT_OF_BOTH_SETS): {
1444 event_indices.push_back(Event_indices(-1,f_count,g_count));
1445 f_count++;
1446 g_count++;
1447 break;
1448 }
1449 }
1450 break;
1451 }
1452 case(CGAL::internal::ROOT_OF_SECOND_SET): {
1453 /*
1454 #if CGAL_ACK_DEBUG_FLAG
1455 CGAL_ACK_DEBUG_PRINT << " two curve event" << std::endl;
1456 #endif
1457 */
1458 this->ptr()->event_slices.emplace_back();
1459
1460 event_indices.push_back
1461 (Event_indices(inter_count,-1,-1));
1462 inter_count++;
1463 break;
1464 }
1465 case(CGAL::internal::ROOT_OF_BOTH_SETS): {
1466 /*
1467 #if CGAL_ACK_DEBUG_FLAG
1468 CGAL_ACK_DEBUG_PRINT << " one and two curve event"
1469 << std::endl;
1470 #endif
1471 */
1472 this->ptr()->event_slices.emplace_back();
1473
1474 switch(*(one_curve_it++)) {
1475 case(CGAL::internal::ROOT_OF_FIRST_SET): {
1476 event_indices.push_back
1477 (Event_indices(inter_count,f_count,-1));
1478 f_count++;
1479 break;
1480 }
1481 case(CGAL::internal::ROOT_OF_SECOND_SET): {
1482 event_indices.push_back
1483 (Event_indices(inter_count,-1,g_count));
1484 g_count++;
1485 break;
1486 }
1487 case(CGAL::internal::ROOT_OF_BOTH_SETS): {
1488 event_indices.push_back
1489 (Event_indices(inter_count,f_count,g_count));
1490 f_count++;
1491 g_count++;
1492 break;
1493 }
1494 }
1495 inter_count++;
1496 break;
1497 }
1498 }
1499 }
1500 CGAL_assertion(inter_count
1501 == static_cast<size_type>
1502 (resultant_roots().size()));
1503 CGAL_assertion(one_curve_it==one_curve_events_type.end());
1504 #if CGAL_ACK_DEBUG_FLAG
1505 CGAL_ACK_DEBUG_PRINT << "done" << std::endl;
1506 #endif
1507
1508 }
1509
1510 //////////////////// compute_intermediate_values_and_slices()
1511
1512 template<typename AlgebraicKernelWithAnalysis_2>
1513 void Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
compute_intermediate_values_and_slices()1514 compute_intermediate_values_and_slices() const {
1515
1516 #if CGAL_ACK_DEBUG_FLAG
1517 CGAL_ACK_DEBUG_PRINT << "Prepare intermediate slices.." << std::flush;
1518 #endif
1519 std::size_t size = event_x_coordinates().size()+1;
1520 this->ptr()->intermediate_values=std::vector<Lazy_bound>();
1521 this->ptr()->intermediate_slices=std::vector<Lazy_status_line_CPA_1>();
1522 this->ptr()->intermediate_values.get().resize(size);
1523 this->ptr()->intermediate_slices.get().resize(size);
1524 #if CGAL_ACK_DEBUG_FLAG
1525 CGAL_ACK_DEBUG_PRINT << "done" << std::endl;
1526 #endif
1527 }
1528
1529 //////////////////// compute_subresultants
1530
1531 template<typename AlgebraicKernelWithAnalysis_2>
1532 void Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
compute_subresultants()1533 compute_subresultants() const {
1534 typedef std::vector<Polynomial_1> Polynomial_container;
1535 this->ptr()->principal_subresultants = Polynomial_container();
1536 this->ptr()->coprincipal_subresultants = Polynomial_container();
1537 const Polynomial_2& f = this->ptr()->f, g=this->ptr()->g;
1538 this->ptr()->subresultants = std::vector<Polynomial_2>();
1539 if(CGAL::degree(f,1)<CGAL::degree(g,1)) {
1540 #if CGAL_ACK_USE_BEZOUT_MATRIX_FOR_SUBRESULTANTS
1541 CGAL::internal::bezout_polynomial_subresultants
1542 (g,f,std::back_inserter(this->ptr()->subresultants.get()));
1543 #else
1544 typename CGAL::Polynomial_traits_d<Polynomial_2>
1545 ::Polynomial_subresultants()
1546 (g,f,std::back_inserter(this->ptr()->subresultants.get()));
1547 #endif
1548 } else {
1549 #if CGAL_ACK_USE_BEZOUT_MATRIX_FOR_SUBRESULTANTS
1550 CGAL::internal::bezout_polynomial_subresultants
1551 (f,g,std::back_inserter(this->ptr()->subresultants.get()));
1552 #else
1553 typename CGAL::Polynomial_traits_d<Polynomial_2>
1554 ::Polynomial_subresultants()
1555 (f,g,std::back_inserter(this->ptr()->subresultants.get()));
1556 #endif
1557 }
1558
1559 std::vector<Polynomial_2>& subresultants
1560 = this->ptr()->subresultants.get();
1561
1562 size_type n = static_cast<size_type>(subresultants.size());
1563
1564 for(size_type i=0;i<n;i++) {
1565 if(CGAL::degree(subresultants[i]) < i) {
1566 this->ptr()->principal_subresultants->
1567 push_back(Polynomial_1(0));
1568 }
1569 else {
1570 this->ptr()->principal_subresultants->
1571 push_back(subresultants[i][i]);
1572 }
1573 }
1574 for(size_type i=1;i<n;i++) {
1575 if(CGAL::degree(subresultants[i]) < i-1) {
1576 this->ptr()->coprincipal_subresultants->
1577 push_back(Polynomial_1(0));
1578 }
1579 else {
1580 this->ptr()->coprincipal_subresultants->
1581 push_back(subresultants[i][i-1]);
1582 }
1583 }
1584 // This must be corrected, if f and g have same degree:
1585 if(CGAL::degree(f,1) == CGAL::degree(g,1)) {
1586 if(n>=1) {
1587 this->ptr()->principal_subresultants.get()[n-1]
1588 = Polynomial_1(CGAL::leading_coefficient(g));
1589 }
1590 if(n>=2) {
1591 this->ptr()->coprincipal_subresultants.get()[n-2]
1592 = Polynomial_1(g[CGAL::degree(g,1)-1]);
1593 }
1594 }
1595
1596 }
1597
1598 //////////////////// create_slice_with_multiplicity_zero_or_one
1599
1600 template<typename AlgebraicKernelWithAnalysis_2>
1601 typename Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>
1602 ::Status_line_CPA_1
1603 Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
create_slice_with_multiplicity_zero_or_one(size_type i)1604 create_slice_with_multiplicity_zero_or_one(size_type i) const {
1605
1606 const std::vector<Algebraic_real_1>& events
1607 = event_x_coordinates();
1608 Algebraic_real_1 alpha = events[i];
1609 const Curve_analysis_2& c1=curve_analysis(false), c2=curve_analysis(true);
1610 size_type i1,i2;
1611 bool flag1,flag2;
1612 c1.x_to_index(alpha,i1,flag1);
1613 c2.x_to_index(alpha,i2,flag2);
1614
1615 bool exactly_at_alpha_1 = flag1, exactly_at_alpha_2 = flag2;
1616 Status_line_CA_1 e1=flag1 ? c1.status_line_at_event(i1)
1617 : c1.status_line_of_interval(i1);
1618 Status_line_CA_1 e2=flag2 ? c2.status_line_at_event(i2)
1619 : c2.status_line_of_interval(i2);
1620
1621 Status_line_CPA_1 left_slice = this->status_line_of_interval(i),
1622 right_slice = this->status_line_of_interval(i+1);
1623
1624 Status_line_CPA_iterator left_it
1625 = ::boost::make_transform_iterator
1626 (::boost::counting_iterator<size_type>(0),
1627 Curves_at_event_functor(left_slice));
1628 Status_line_CPA_iterator right_it
1629 = ::boost::make_transform_iterator
1630 (::boost::counting_iterator<size_type>(0),
1631 Curves_at_event_functor(right_slice));
1632
1633 Status_line_CPA_iterator left_end
1634 = ::boost::make_transform_iterator
1635 (::boost::counting_iterator<size_type>
1636 (left_slice.number_of_events()),
1637 Curves_at_event_functor(left_slice));
1638 Status_line_CPA_iterator right_end
1639 = ::boost::make_transform_iterator
1640 (::boost::counting_iterator<size_type>
1641 (right_slice.number_of_events()),
1642 Curves_at_event_functor(right_slice));
1643
1644 // Take out asymptotes
1645 size_type asym_lm_1
1646 = e1.number_of_branches_approaching_minus_infinity().first;
1647 size_type asym_rm_1
1648 = e1.number_of_branches_approaching_minus_infinity().second;
1649 size_type asym_lp_1
1650 = e1.number_of_branches_approaching_plus_infinity().first;
1651 size_type asym_rp_1
1652 = e1.number_of_branches_approaching_plus_infinity().second;
1653 size_type asym_lm_2
1654 = e2.number_of_branches_approaching_minus_infinity().first;
1655 size_type asym_rm_2
1656 = e2.number_of_branches_approaching_minus_infinity().second;
1657 size_type asym_lp_2
1658 = e2.number_of_branches_approaching_plus_infinity().first;
1659 size_type asym_rp_2
1660 = e2.number_of_branches_approaching_plus_infinity().second;
1661
1662 while(asym_lm_1 != 0 || asym_lm_2 != 0) {
1663 CGAL_assertion(*left_it != CGAL::internal::INTERSECTION);
1664 if(*left_it == CGAL::internal::FIRST_CURVE) {
1665 CGAL_assertion(asym_lm_1!=0);
1666 asym_lm_1--;
1667 }
1668 if(*left_it == CGAL::internal::SECOND_CURVE) {
1669 CGAL_assertion(asym_lm_2!=0);
1670 asym_lm_2--;
1671 }
1672 left_it++;
1673 }
1674 while(asym_rm_1 != 0 || asym_rm_2 != 0) {
1675 if(*right_it == CGAL::internal::FIRST_CURVE) {
1676 CGAL_assertion(asym_rm_1!=0);
1677 asym_rm_1--;
1678 }
1679 if(*right_it == CGAL::internal::SECOND_CURVE) {
1680 CGAL_assertion(asym_rm_2!=0);
1681 asym_rm_2--;
1682 }
1683 right_it++;
1684 }
1685 while(asym_lp_1 != 0 || asym_lp_2 != 0) {
1686 left_end--;
1687 if(*left_end == CGAL::internal::FIRST_CURVE) {
1688 CGAL_assertion(asym_lp_1!=0);
1689 asym_lp_1--;
1690 }
1691 if(*left_end == CGAL::internal::SECOND_CURVE) {
1692 CGAL_assertion(asym_lp_2!=0);
1693 asym_lp_2--;
1694 }
1695 }
1696 while(asym_rp_1 != 0 || asym_rp_2 != 0) {
1697 right_end--;
1698 if(*right_end == CGAL::internal::FIRST_CURVE) {
1699 CGAL_assertion(asym_rp_1!=0);
1700 asym_rp_1--;
1701 }
1702 if(*right_end == CGAL::internal::SECOND_CURVE) {
1703 CGAL_assertion(asym_rp_2!=0);
1704 asym_rp_2--;
1705 }
1706 }
1707 // Now, the iterator ranges [left_it,left_end)
1708 // and [right_it,right_end) give the arcs really
1709 // going into the event line
1710
1711 Slice_info slice_info;
1712 CGAL::internal::Slice_type curr_lowest_arc;
1713 size_type curr_multiplicity;
1714
1715 size_type event_index_1=0, event_index_2=0;
1716
1717 while(event_index_1 != e1.number_of_events() ||
1718 event_index_2 != e2.number_of_events()) {
1719
1720 CGAL_assertion(event_index_1 != e1.number_of_events() ||
1721 event_index_2 != e2.number_of_events());
1722 if(event_index_1==e1.number_of_events()) {
1723 curr_lowest_arc=CGAL::internal::SECOND_CURVE;
1724 }
1725 else if(event_index_2==e2.number_of_events()) {
1726 curr_lowest_arc=CGAL::internal::FIRST_CURVE;
1727 }
1728 else if((e1.number_of_incident_branches(event_index_1).first>0 &&
1729 e2.number_of_incident_branches(event_index_2).first>0)) {
1730 // The next arc on the left must come as next:
1731 curr_lowest_arc=*left_it;
1732 }
1733 else if((e1.number_of_incident_branches(event_index_1).second>0 &&
1734 e2.number_of_incident_branches(event_index_2).second>0)) {
1735 // The next arc on the right must come as next:
1736 curr_lowest_arc=*right_it;
1737 }
1738 else {
1739 // We cannot decide it from the arcs, so we have to compare
1740 // isolating intervals
1741 if(! exactly_at_alpha_1) {
1742 e1 = c1.status_line_at_exact_x(alpha);
1743 CGAL_assertion(e1.number_of_events()>event_index_1);
1744 }
1745 if(! exactly_at_alpha_2) {
1746 e2 = c2.status_line_at_exact_x(alpha);
1747 CGAL_assertion(e2.number_of_events()>event_index_2);
1748 }
1749 CGAL::Sign e1_smaller
1750 = split_compare(e1,event_index_1,e2,event_index_2);
1751 curr_lowest_arc
1752 = (e1_smaller==CGAL::SMALLER)
1753 ? CGAL::internal::FIRST_CURVE : CGAL::internal::SECOND_CURVE;
1754 }
1755
1756 curr_multiplicity = -1;
1757
1758 // Move the iterators
1759 size_type arcs_of_other_curve_left=0, arcs_of_other_curve_right=0;
1760 if(curr_lowest_arc==CGAL::internal::FIRST_CURVE) {
1761 size_type j=0;
1762 while(j<e1.number_of_incident_branches(event_index_1).first) {
1763 if(*left_it==CGAL::internal::FIRST_CURVE) {
1764 j++;
1765 } else {
1766 CGAL_assertion(event_index_2 < e2.number_of_events());
1767 arcs_of_other_curve_left++;
1768 }
1769 left_it++;
1770 }
1771
1772 j=0;
1773 while(j<e1.number_of_incident_branches(event_index_1).second) {
1774 if(*right_it==CGAL::internal::FIRST_CURVE) {
1775 j++;
1776 } else {
1777 CGAL_assertion(event_index_2 < e2.number_of_events());
1778 arcs_of_other_curve_right++;
1779 }
1780 right_it++;
1781 }
1782 event_index_1++;
1783 if(arcs_of_other_curve_left+arcs_of_other_curve_right>0) {
1784 // Intersection! Iterate over the remaining arcs
1785 // on both sides belonging to this intersection
1786 for(size_type j=arcs_of_other_curve_left;
1787 j<e2.number_of_incident_branches(event_index_2).first;
1788 j++) {
1789 CGAL_assertion(*left_it==CGAL::internal::SECOND_CURVE);
1790 left_it++;
1791 }
1792 for(size_type j=arcs_of_other_curve_right;
1793 j<e2.number_of_incident_branches(event_index_2).second;
1794 j++) {
1795 CGAL_assertion(*right_it==CGAL::internal::SECOND_CURVE);
1796 right_it++;
1797 }
1798 event_index_2++;
1799 curr_lowest_arc=CGAL::internal::INTERSECTION;
1800 curr_multiplicity=1;
1801 }
1802 } else { // curr_lowest_arc=CGAL::internal::SECOND_CURVE
1803 size_type j=0;
1804 while(j<e2.number_of_incident_branches(event_index_2).first) {
1805 if(*left_it==CGAL::internal::SECOND_CURVE) {
1806 j++;
1807 } else {
1808 CGAL_assertion(event_index_1 < e1.number_of_events());
1809 arcs_of_other_curve_left++;
1810 }
1811 left_it++;
1812 }
1813
1814 j=0;
1815 while(j<e2.number_of_incident_branches(event_index_2).second) {
1816 if(*right_it==CGAL::internal::SECOND_CURVE) {
1817 j++;
1818 } else {
1819 CGAL_assertion(event_index_1 < e1.number_of_events());
1820 arcs_of_other_curve_right++;
1821 }
1822 right_it++;
1823 }
1824 event_index_2++;
1825 if(arcs_of_other_curve_left+arcs_of_other_curve_right>0) {
1826 // Intersection! Iterate over the remaining arcs
1827 // on both sides belonging to this intersection
1828 for(size_type j=arcs_of_other_curve_left;
1829 j<e1.number_of_incident_branches(event_index_1).first;
1830 j++) {
1831 CGAL_assertion(*left_it==CGAL::internal::FIRST_CURVE);
1832 left_it++;
1833 }
1834 for(size_type j=arcs_of_other_curve_right;
1835 j<e1.number_of_incident_branches(event_index_1).second;
1836 j++) {
1837 CGAL_assertion(*right_it==CGAL::internal::FIRST_CURVE);
1838 right_it++;
1839 }
1840 event_index_1++;
1841 curr_lowest_arc=CGAL::internal::INTERSECTION;
1842 curr_multiplicity=1;
1843 }
1844
1845 }
1846 slice_info.push_back(std::make_pair(curr_lowest_arc,
1847 curr_multiplicity));
1848 }
1849 CGAL_assertion(left_it == left_end &&
1850 right_it == right_end);
1851
1852 return create_slice_from_slice_info(i,slice_info,true);
1853 }
1854
1855 //////////////////// create_intermediate_slice_at
1856
1857 template<typename AlgebraicKernelWithAnalysis_2>
1858 typename Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>
1859 ::Status_line_CPA_1
1860 Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
create_intermediate_slice_at(int i)1861 create_intermediate_slice_at(int i) const {
1862
1863 Bound r = bound_value_in_interval(i);
1864
1865 std::vector<Algebraic_real_1> p1_roots,p2_roots;
1866
1867 this->ptr()->c1_.get_roots_at_rational(r,std::back_inserter(p1_roots));
1868 this->ptr()->c2_.get_roots_at_rational(r,std::back_inserter(p2_roots));
1869
1870 size_type number_of_roots
1871 = static_cast<size_type>(p1_roots.size() + p2_roots.size());
1872 std::vector<Algebraic_real_1> p12_roots;
1873 p12_roots.reserve(number_of_roots);
1874 std::vector<CGAL::internal::Three_valued> p12_order;
1875 p12_order.reserve(number_of_roots);
1876
1877 CGAL::internal::Distinct_compare<Algebraic_real_1> distinct_compare;
1878 set_union_with_source(p1_roots.begin(),
1879 p1_roots.end(),
1880 p2_roots.begin(),
1881 p2_roots.end(),
1882 std::back_inserter(p12_roots),
1883 std::back_inserter(p12_order),
1884 distinct_compare);
1885
1886 Slice_info slice_info;
1887
1888 for(typename std::vector<CGAL::internal::Three_valued>::const_iterator
1889 it = p12_order.begin();
1890 it!=p12_order.end();
1891 it++) {
1892 switch(*it){
1893 case(CGAL::internal::ROOT_OF_FIRST_SET): {
1894 slice_info.push_back
1895 (std::make_pair(CGAL::internal::FIRST_CURVE,-1));
1896 break;
1897 }
1898 case(CGAL::internal::ROOT_OF_SECOND_SET): {
1899 slice_info.push_back
1900 (std::make_pair(CGAL::internal::SECOND_CURVE,-1));
1901 break;
1902 }
1903 case(CGAL::internal::ROOT_OF_BOTH_SETS): {
1904 CGAL_assertion(false);
1905 break;
1906 }
1907 }
1908 }
1909
1910 Status_line_CPA_1 new_slice
1911 = create_slice_from_slice_info(i,slice_info,false);
1912
1913 return new_slice;
1914 }
1915
1916 //////////////////// create_slice_from_slice_info
1917
1918 template<typename AlgebraicKernelWithAnalysis_2>
1919 typename Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>
1920 ::Status_line_CPA_1
1921 Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
create_slice_from_slice_info(size_type id,const Slice_info & slice,bool event_flag)1922 create_slice_from_slice_info(size_type id,
1923 const Slice_info& slice,
1924 bool event_flag) const {
1925 typedef typename Status_line_CPA_1::Arc_pair Arc_pair;
1926 typedef typename Status_line_CPA_1::Arc_container Arc_container;
1927 typedef typename Status_line_CPA_1::Int_container Int_container;
1928 Arc_container arc_container;
1929 Int_container int_container;
1930
1931 for(typename Slice_info::const_iterator it = slice.begin();
1932 it!=slice.end();
1933 it++) {
1934 CGAL_assertion(it->first != CGAL::internal::CANDIDATE);
1935 switch(it->first) {
1936 case(CGAL::internal::FIRST_CURVE): {
1937 if(event_flag) {
1938 arc_container.push_back(Arc_pair(0,it->second));
1939 } else {
1940 int_container.push_back(0);
1941 }
1942 break;
1943 }
1944 case(CGAL::internal::SECOND_CURVE): {
1945 if(event_flag) {
1946 arc_container.push_back(Arc_pair(1,it->second));
1947 } else {
1948 int_container.push_back(1);
1949 }
1950 break;
1951 }
1952 case(CGAL::internal::INTERSECTION): {
1953 CGAL_assertion(event_flag);
1954 arc_container.push_back(Arc_pair(2,it->second));
1955 break;
1956 }
1957 case(CGAL::internal::CANDIDATE): {
1958 CGAL_assertion(false);
1959 break;
1960 }
1961 }
1962 }
1963
1964 return event_flag
1965 ? Status_line_CPA_1(id,arc_container,*this)
1966 : Status_line_CPA_1(id,int_container,*this);
1967 }
1968
1969 //////////////////// construct_slice_info
1970
1971 template<typename AlgebraicKernelWithAnalysis_2>
1972 typename Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::Slice_info
1973 Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
construct_slice_info(Algebraic_real_1 alpha)1974 construct_slice_info(Algebraic_real_1 alpha) const {
1975
1976 /*
1977 #if CGAL_ACK_DEBUG_FLAG
1978 CGAL_ACK_DEBUG_PRINT << "Consider alpha=" << CGAL::to_double(alpha)
1979 << std::endl;
1980 #endif
1981 */
1982
1983 Status_line_CA_1 e1 = this->ptr()->c1_.status_line_at_exact_x(alpha);
1984
1985 Status_line_CA_1 e2 = this->ptr()->c2_.status_line_at_exact_x(alpha);
1986
1987 std::vector<std::pair<size_type,size_type> > matchings;
1988 for(size_type i=0;i<e1.number_of_events();i++) {
1989 size_type match=find_possible_matching(e1,i,e2);
1990 if(match==-1) {
1991 continue;
1992 }
1993 if(find_possible_matching(e2,match,e1) != i) {
1994 continue;
1995 }
1996 /*
1997 #if CGAL_ACK_DEBUG_FLAG
1998 CGAL_ACK_DEBUG_PRINT << "New matching: (" << i
1999 << "," << match << ")" << std::endl;
2000 #endif
2001 */
2002 matchings.push_back(std::make_pair(i,match));
2003 }
2004 size_type i1=0, i2=0,
2005 n1=e1.number_of_events(), n2=e2.number_of_events();
2006
2007 typename std::vector<std::pair<size_type,size_type> >::const_iterator
2008 match = matchings.begin();
2009 Slice_info slice_info;
2010 while(i1<n1 || i2<n2) {
2011 if(i1==n1) {
2012 slice_info.push_back
2013 (std::make_pair(CGAL::internal::SECOND_CURVE,-1));
2014 i2++;
2015 continue;
2016 }
2017 if(i2==n2) {
2018 slice_info.push_back
2019 (std::make_pair(CGAL::internal::FIRST_CURVE,-1));
2020 i1++;
2021 continue;
2022 }
2023 if(match!=matchings.end() &&
2024 i1==match->first &&
2025 i2==match->second) {
2026 slice_info.push_back(std::make_pair(CGAL::internal::CANDIDATE,1));
2027 i1++;
2028 i2++;
2029 match++;
2030 continue;
2031 }
2032 CGAL_assertion(!overlap(e1,i1,e2,i2));
2033 if(e1.lower_bound(i1) < e2.lower_bound(i2)) {
2034 slice_info.push_back
2035 (std::make_pair(CGAL::internal::FIRST_CURVE,-1));
2036 i1++;
2037 continue;
2038 } else {
2039 slice_info.push_back
2040 (std::make_pair(CGAL::internal::SECOND_CURVE,-1));
2041 i2++;
2042 continue;
2043 }
2044 }
2045 CGAL_assertion(match==matchings.end());
2046 return slice_info;
2047 }
2048
2049 //////////////////// construct_generic_case
2050
2051 template<typename AlgebraicKernelWithAnalysis_2>
2052 typename Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>
2053 ::Status_line_CPA_1
2054 Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
construct_generic_case(size_type i)2055 construct_generic_case(size_type i) const {
2056
2057 Algebraic_real_1 alpha = event_x(i);
2058
2059 Slice_info slice_info;
2060
2061 size_type index_of_fg = event_indices(i).fg;
2062 size_type index_of_ffy =event_indices(i).ffy;
2063 size_type index_of_ggy =event_indices(i).ggy;
2064 if(index_of_fg>=0) {
2065 if(kernel()->is_zero_at_1_object()
2066 (CGAL::leading_coefficient
2067 (this->ptr()->c1_.polynomial_2()),alpha)
2068 ||
2069 kernel()->is_zero_at_1_object()
2070 (CGAL::leading_coefficient
2071 (this->ptr()->c2_.polynomial_2()),alpha)) {
2072 throw CGAL::internal::Non_generic_position_exception();
2073 }
2074 size_type k = -1; // not yet computed
2075 if(index_of_ffy==-1 && index_of_ggy==-1) {
2076 // this means, we need the multiplicity of the intersections
2077 if(kernel()->is_zero_at_1_object()
2078 (principal_subresultants(1),alpha)) {
2079 // multiplicity cannot be determined, throw exception
2080 throw CGAL::internal::Non_generic_position_exception();
2081 } else {
2082 k=1;
2083 }
2084 } else {
2085 k = degree_of_local_gcd(index_of_fg,alpha);
2086 }
2087 Status_line_CA_1 e1
2088 = this->ptr()->c1_.status_line_at_exact_x(alpha);
2089 Status_line_CA_1 e2
2090 = this->ptr()->c2_.status_line_at_exact_x(alpha);
2091 slice_info = construct_slice_info(alpha);
2092 size_type no_candidates=
2093 reduce_number_of_candidates_and_intersections_to
2094 (1,e1,e2,slice_info,k);
2095 CGAL_assertion(no_candidates==0 || no_candidates==1);
2096 if(no_candidates==1) {
2097 typename Slice_info::iterator slice_it
2098 = slice_info.begin();
2099 size_type i1=0,i2=0;
2100 while(slice_it->first!=CGAL::internal::CANDIDATE) {
2101 if(slice_it->first==CGAL::internal::FIRST_CURVE) {
2102 i1++;
2103 }
2104 if(slice_it->first==CGAL::internal::SECOND_CURVE) {
2105 i2++;
2106 }
2107 if(slice_it->first==CGAL::internal::INTERSECTION) {
2108 i1++;
2109 i2++;
2110 }
2111 slice_it++;
2112 }
2113 check_candidate(e1,i1,e2,i2,k,
2114 slice_info,
2115 slice_it,i);
2116 }
2117 } else {
2118 Status_line_CA_1 e1
2119 = this->ptr()->c1_.status_line_at_exact_x(alpha);
2120 Status_line_CA_1 e2
2121 = this->ptr()->c2_.status_line_at_exact_x(alpha);
2122 slice_info = construct_slice_info(alpha);
2123 reduce_number_of_candidates_and_intersections_to
2124 (0,e1,e2,slice_info,0);
2125 }
2126 return create_slice_from_slice_info(i,slice_info,true);
2127 }
2128
2129 //////////////////// check_candidate_by_arc_pattern
2130
2131 template<typename AlgebraicKernelWithAnalysis_2>
2132
2133 bool Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
check_candidate_by_arc_pattern(size_type index,Status_line_CA_1 & e1,size_type i1,Status_line_CA_1 & e2,size_type i2)2134 check_candidate_by_arc_pattern(size_type index,
2135 Status_line_CA_1& e1,
2136 size_type i1,
2137 Status_line_CA_1& e2,
2138 size_type i2) const {
2139
2140 Status_line_CPA_1 left_slice = status_line_of_interval(index),
2141 right_slice = status_line_of_interval(index+1);
2142 size_type left_index=0,right_index=0;
2143 for(size_type i=0;i<i1;i++) {
2144 left_index+=e1.number_of_incident_branches(i).first;
2145 right_index+=e1.number_of_incident_branches(i).second;
2146 }
2147 for(size_type i=0;i<i2;i++) {
2148 left_index+=e2.number_of_incident_branches(i).first;
2149 right_index+=e2.number_of_incident_branches(i).second;
2150 }
2151 // left_index and right_index now point at the position
2152 // of the first arc into the candidate
2153 size_type num_of_arcs_to_candidate_left
2154 = e1.number_of_incident_branches(i1).first
2155 + e2.number_of_incident_branches(i2).first,
2156 num_of_arcs_to_candidate_right
2157 = e1.number_of_incident_branches(i1).second
2158 + e2.number_of_incident_branches(i2).second;
2159 CGAL_assertion(left_index + num_of_arcs_to_candidate_left <=
2160 left_slice.number_of_events());
2161 CGAL_assertion(right_index + num_of_arcs_to_candidate_right <=
2162 right_slice.number_of_events());
2163
2164 CGAL::internal::Slice_type curr;
2165 Curves_at_event_functor left_functor(left_slice);
2166 size_type number_of_changes;
2167 if(left_index < left_slice.number_of_events()) {
2168 curr = left_functor(left_index);
2169 number_of_changes=0;
2170 /*
2171 #if CGAL_ACK_DEBUG_FLAG
2172 CGAL_ACK_DEBUG_PRINT << num_of_arcs_to_candidate_left
2173 << num_of_arcs_to_candidate_right
2174 << left_index << right_index << std::endl;
2175 #endif
2176 */
2177 for(size_type i=1;i<num_of_arcs_to_candidate_left;i++) {
2178 if(curr != left_functor(left_index+i)) {
2179 curr = left_functor(left_index+i);
2180 number_of_changes++;
2181 if(number_of_changes>=2) {
2182 return true;
2183 }
2184 }
2185 }
2186 }
2187
2188 Curves_at_event_functor right_functor(right_slice);
2189
2190 if(right_index < right_slice.number_of_events()) {
2191 curr = right_functor(right_index);
2192 number_of_changes=0;
2193 for(size_type i=1;i<num_of_arcs_to_candidate_right;i++) {
2194 if(curr != right_functor(right_index+i)) {
2195 curr = right_functor(right_index+i);
2196 number_of_changes++;
2197 if(number_of_changes>=2) {
2198 return true;
2199 }
2200 }
2201 }
2202 }
2203
2204 return false;
2205 }
2206
2207 //////////////////// check_candidate
2208
2209 template<typename AlgebraicKernelWithAnalysis_2>
2210 template<typename InputIterator>
2211 void Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
check_candidate(Status_line_CA_1 & e1,size_type i1,Status_line_CA_1 & e2,size_type i2,size_type k,Slice_info & slice_info,InputIterator slice_it,size_type root_index)2212 check_candidate(Status_line_CA_1& e1,size_type i1,
2213 Status_line_CA_1& e2,size_type i2,
2214 size_type k,
2215 Slice_info& slice_info,
2216 InputIterator slice_it,
2217 size_type root_index) const {
2218
2219 Algebraic_real_1 xval = e1.x();
2220 CGAL_assertion(xval==e2.x());
2221 bool is_intersection=false;
2222 size_type index_in_fg = event_indices(root_index).fg;
2223 size_type mult_of_resultant
2224 = multiplicities_of_resultant_roots(index_in_fg);
2225
2226 if(k%2==1 || mult_of_resultant%2==1) {
2227 is_intersection=true;
2228 } else {
2229 if(check_candidate_by_arc_pattern(root_index,e1,i1,e2,i2)) {
2230 is_intersection=true;
2231 } else {
2232 is_intersection=check_candidate_symbolically(e1,i1,e2,i2,k);
2233 }
2234 }
2235 if(is_intersection) {
2236 slice_it=slice_info.erase(slice_it);
2237 size_type mult_of_intersection;
2238 if(k==1) {
2239 mult_of_intersection = mult_of_resultant;
2240 } else {
2241 mult_of_intersection = -1;
2242 }
2243 slice_info.insert(slice_it,
2244 std::make_pair(CGAL::internal::INTERSECTION,
2245 mult_of_intersection));
2246 }
2247 else {
2248 CGAL::Sign e1_smaller=split_compare(e1,i1,e2,i2);
2249
2250 slice_it=slice_info.erase(slice_it);
2251 if(e1_smaller==CGAL::SMALLER) {
2252 slice_it = slice_info.insert
2253 (slice_it,std::make_pair(CGAL::internal::FIRST_CURVE,-1));
2254 slice_it++;
2255 slice_it = slice_info.insert
2256 (slice_it,std::make_pair(CGAL::internal::SECOND_CURVE,-1));
2257 } else {
2258 slice_it = slice_info.insert
2259 (slice_it,std::make_pair(CGAL::internal::SECOND_CURVE,-1));
2260 slice_it++;
2261 slice_it = slice_info.insert
2262 (slice_it,std::make_pair(CGAL::internal::FIRST_CURVE,-1));
2263 }
2264 }
2265 }
2266
2267 //////////////////// find_possible_matching
2268
2269 template<typename AlgebraicKernelWithAnalysis_2>
2270 typename Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::size_type
2271 Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
find_possible_matching(Status_line_CA_1 & e1,size_type index1,Status_line_CA_1 & e2)2272 find_possible_matching(Status_line_CA_1& e1,
2273 size_type index1,
2274 Status_line_CA_1& e2) const {
2275
2276 std::vector<size_type> possible_overlaps;
2277 for(size_type i=0;i<e2.number_of_events();i++) {
2278 if(overlap(e1,index1,e2,i)) {
2279 possible_overlaps.push_back(i);
2280 }
2281 }
2282 while(possible_overlaps.size()>1) {
2283 if(possible_overlaps.size()==2) {
2284 // Prevent that both intervals touch in a bound
2285 while(overlap(e2,possible_overlaps[0],
2286 e2,possible_overlaps[1])) {
2287 e2.refine(possible_overlaps[0]);
2288 e2.refine(possible_overlaps[1]);
2289 }
2290 }
2291 e1.refine(index1);
2292
2293 typename std::vector<size_type>::iterator it
2294 = possible_overlaps.begin();
2295 while(it!=possible_overlaps.end()) {
2296 if(!overlap(e1,index1,e2,*it)) {
2297 it=possible_overlaps.erase(it);
2298 }
2299 else {
2300 it++;
2301 }
2302 }
2303 }
2304 if(possible_overlaps.size()==0) {
2305 return -1;
2306 }
2307 else {
2308 return possible_overlaps[0];
2309 }
2310 }
2311
2312
2313 //////////////////// new_shear_for_intersection_info
2314
2315 template<typename AlgebraicKernelWithAnalysis_2>
2316 void Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
new_shear_for_intersection_info(Intersection_info_container & info_container)2317 new_shear_for_intersection_info(Intersection_info_container& info_container)
2318 const {
2319 #if CGAL_ACK_DEBUG_FLAG
2320 CGAL_ACK_DEBUG_PRINT << "Use shear for intersections.." << std::endl;
2321 #endif
2322 bool good_direction_found=false;
2323 Integer s;
2324
2325 while(! good_direction_found) {
2326 try {
2327 info_container.clear();
2328 info_container.resize(resultant_roots().size());
2329 s = this->ptr()->shear_controller.get_shear_factor();
2330 #if CGAL_ACK_DEBUG_FLAG
2331 CGAL_ACK_DEBUG_PRINT << "Try shear factor " << s << std::endl;
2332 CGAL_ACK_DEBUG_PRINT
2333 << ">>>>>>>>>>> Transform first curve" << std::endl;
2334 #endif
2335
2336 Curve_analysis_2 sh1
2337 = this->ptr()->c1_.shear_primitive_part(s);
2338 #if CGAL_ACK_DEBUG_FLAG
2339 CGAL_ACK_DEBUG_PRINT
2340 << "<<<<<<<<<<< End of transform first curve" << std::endl;
2341 CGAL_ACK_DEBUG_PRINT << ">>>>>>>>>>> Transform second curve"
2342 << std::endl;
2343 #endif
2344 Curve_analysis_2 sh2 = this->ptr()->c2_.shear_primitive_part(s);
2345 #if CGAL_ACK_DEBUG_FLAG
2346 CGAL_ACK_DEBUG_PRINT
2347 << "<<<<<<<<<<< End of transform second curve"
2348 << std::endl;
2349 #endif
2350 Self sh_pair(kernel(),sh1,sh2,CGAL::EXCEPTION_STRATEGY);
2351
2352 #if CGAL_ACK_DEBUG_FLAG
2353 CGAL_ACK_DEBUG_PRINT << "Shear back intersection points..."
2354 << std::flush;
2355 #endif
2356 for(size_type i=0;
2357 i<static_cast<size_type>
2358 (sh_pair.event_x_coordinates().size());
2359 i++) {
2360 if(sh_pair.event_indices(i).fg==-1) {
2361 continue;
2362 }
2363 Status_line_CPA_1 slice
2364 = sh_pair.status_line_at_event(i);
2365 Curves_at_event_functor functor(slice);
2366 for(size_type j=0;j<slice.number_of_events();j++) {
2367 if(functor(j) == CGAL::internal::INTERSECTION) {
2368 this->update_intersection_info(info_container,
2369 sh_pair,
2370 slice,
2371 i,j,s);
2372 }
2373 }
2374 }
2375 good_direction_found=true;
2376 }
2377 catch(CGAL::internal::Non_generic_position_exception /* ex */) {
2378 this->ptr()->shear_controller.report_failure(s);
2379 }
2380 }
2381
2382 #if CGAL_ACK_DEBUG_FLAG
2383 CGAL_ACK_DEBUG_PRINT << "done" << std::endl;
2384 #endif
2385 return;
2386 }
2387
2388 //////////////////// create_event_slice_from_current_intersection_info
2389
2390 template<typename AlgebraicKernelWithAnalysis_2>
2391 typename Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>
2392 ::Status_line_CPA_1
2393 Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
create_event_slice_from_current_intersection_info(size_type i)2394 create_event_slice_from_current_intersection_info (size_type i) const{
2395 #if CGAL_ACK_DEBUG_FLAG
2396 CGAL_ACK_DEBUG_PRINT << "Reduce the candidates.." << std::flush;
2397 #endif
2398 Event_indices ev_ind = event_indices(i);
2399 size_type index_of_fg = ev_ind.fg;
2400 Intersection_info_container& intersection_info_container
2401 = *(this->ptr()->intersection_info_container);
2402 Algebraic_real_1 alpha = event_x(i);
2403 CGAL_assertion(index_of_fg <
2404 static_cast<size_type>
2405 (intersection_info_container.size()));
2406 #if CGAL_ACK_DEBUG_FLAG
2407 CGAL_ACK_DEBUG_PRINT << i << "th slice has "
2408 << intersection_info_container[index_of_fg].size()
2409 << " intersections" << std::endl;
2410 #endif
2411 Status_line_CA_1 e1=this->ptr()->c1_.
2412 status_line_at_exact_x(resultant_roots(index_of_fg)),
2413 e2=this->ptr()->c2_.
2414 status_line_at_exact_x(resultant_roots(index_of_fg));
2415 Slice_info slice=construct_slice_info(alpha);
2416 CGAL_assertion_code(size_type no_intersections=)
2417 reduce_number_of_candidates_and_intersections_to
2418 (static_cast<size_type>
2419 (intersection_info_container[index_of_fg].size()),
2420 e1,
2421 e2,
2422 slice,
2423 -1);
2424 CGAL_assertion(no_intersections==static_cast<size_type>
2425 (intersection_info_container[index_of_fg].size()));
2426 typename std::vector<typename Rep::Intersection_info>::iterator
2427 inter_info_it
2428 = intersection_info_container[index_of_fg].begin();
2429 for(size_type j=0;j<static_cast<size_type>(slice.size());j++) {
2430 if(slice[j].first==CGAL::internal::INTERSECTION) {
2431 inter_info_it++;
2432 }
2433 if(slice[j].first==CGAL::internal::CANDIDATE) {
2434 slice[j].first=CGAL::internal::INTERSECTION;
2435 if(ev_ind.ffy==-1 && ev_ind.ggy==-1 && inter_info_it->mult==-1) {
2436 // Multiplicity unknown for case where we need it
2437 throw CGAL::internal::Non_generic_position_exception();
2438 }
2439 slice[j].second=inter_info_it->mult;
2440 inter_info_it++;
2441 }
2442 }
2443
2444 #if CGAL_ACK_DEBUG_FLAG
2445 CGAL_ACK_DEBUG_PRINT << "done" << std::endl;
2446 #endif
2447 return create_slice_from_slice_info(i,slice,true);
2448
2449 }
2450
2451 //////////////////// update_intersection_info
2452
2453 template<typename AlgebraicKernelWithAnalysis_2>
2454 void Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
update_intersection_info(Intersection_info_container & info_container,Self & sh_pair,Status_line_CPA_1 slice,size_type i,size_type j,Integer s)2455 update_intersection_info(Intersection_info_container&
2456 info_container,
2457 Self& sh_pair,
2458 Status_line_CPA_1 slice,
2459 size_type i,
2460 size_type j,
2461 Integer s) const {
2462 typedef typename Rep::Intersection_info Intersection_info;
2463 const Algebraic_real_1& xval = sh_pair.event_x(i);
2464 CGAL_assertion(Curves_at_event_functor(slice)(j)
2465 ==CGAL::internal::INTERSECTION);
2466 Status_line_CA_1 ev = sh_pair.ptr()->c1_.status_line_at_exact_x(xval);
2467 // x_coordinate is given by xval
2468 // y_coordinate by ev[index]
2469 Intersection_info intersection_info;
2470 intersection_info.ev=ev;
2471 int index = slice.curves_at_event(j).first;
2472 intersection_info.index=index;
2473 intersection_info.mult=slice.multiplicity_of_intersection(j);
2474 // Find the right position to insert the object
2475 // first the "x-coordiante"
2476 size_type left_index = -1,
2477 right_index = static_cast<size_type>(stripe_values().size()-1);
2478 Algebraic_real_1 xv = ev.x();
2479 Bound lx = xv.low(), rx=xv.high(),
2480 x_iv_size = rx-lx;
2481 Bound ly = ev.lower_bound(index),
2482 ry = ev.upper_bound(index);;
2483 while(left_index < right_index) {
2484 if(x_iv_size > ry-ly) {
2485 xv.refine();
2486 lx = xv.low();
2487 rx=xv.high();
2488 x_iv_size=rx-lx;
2489 continue;
2490 }
2491 ev.refine(index);
2492 ly = ev.lower_bound(index);
2493 ry = ev.upper_bound(index);
2494 Bound left=(s<0) ? x_sheared(lx,ry,-s): x_sheared(lx,ly,-s);
2495 Bound right = (s<0) ? x_sheared(rx,ly,-s) : x_sheared(rx,ry,-s);
2496 CGAL_assertion(left<right);
2497 while(left_index<right_index &&
2498 stripe_values()[left_index+1]<left) {
2499 ++left_index;
2500 }
2501 while(left_index<right_index &&
2502 right<stripe_values()[right_index]) {
2503 --right_index;
2504 }
2505 }
2506 CGAL_assertion(left_index==right_index);
2507
2508 // Now, the "y-coordinate"
2509 typedef std::vector<Intersection_info> Intersection_info_vector;
2510 Intersection_info_vector& info_vec
2511 = info_container[left_index];
2512 typename Intersection_info_vector::iterator info_it=info_vec.begin();
2513 while(info_it!=info_vec.end()) {
2514 Status_line_CA_1& comp_ev=info_it->ev;
2515 size_type comp_index = info_it->index;
2516 CGAL::Sign index_smaller
2517 = split_compare(ev,index,comp_ev,comp_index);
2518 if(index_smaller==CGAL::LARGER) {
2519 info_it++;
2520 } else {
2521 break;
2522 }
2523 }
2524 info_vec.insert(info_it,intersection_info);
2525 }
2526
2527 //////////////////// reduce_number_of_candidates_and_intersections_to
2528
2529 template<typename AlgebraicKernelWithAnalysis_2>
2530 typename Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::size_type
2531 Curve_pair_analysis_2<AlgebraicKernelWithAnalysis_2>::
reduce_number_of_candidates_and_intersections_to(size_type n,Status_line_CA_1 & e1,Status_line_CA_1 & e2,Slice_info & slice,size_type k)2532 reduce_number_of_candidates_and_intersections_to(size_type n,
2533 Status_line_CA_1& e1,
2534 Status_line_CA_1& e2,
2535 Slice_info& slice,
2536 size_type k) const {
2537 /*
2538 #if CGAL_ACK_DEBUG_FLAG
2539 CGAL_ACK_DEBUG_PRINT << "Reduce: " << n << " "
2540 << CGAL::to_double(e1.x()) << " " << k
2541 << std::endl;
2542 #endif
2543 */
2544 size_type number_of_intersections=0;
2545 size_type number_of_candidates=0;
2546 for(size_type i=0;i<static_cast<size_type>(slice.size());i++) {
2547 if(slice[i].first==CGAL::internal::CANDIDATE) {
2548 number_of_candidates++;
2549 }
2550 if(slice[i].first==CGAL::internal::INTERSECTION) {
2551 number_of_intersections++;
2552 }
2553 }
2554 CGAL_assertion(number_of_intersections<=n);
2555
2556 typename Slice_info::iterator slice_it=slice.begin();
2557 size_type i1=0,i2=0;
2558 size_type max_candidate_mult=0;
2559 while(n<number_of_candidates+number_of_intersections) {
2560 if(slice_it==slice.end()) {
2561 CGAL_assertion(e1.number_of_events()==i1 &&
2562 e2.number_of_events()==i2);
2563 if(max_candidate_mult<k) {
2564 throw CGAL::internal::Non_generic_position_exception();
2565 } else {
2566 slice_it=slice.begin();
2567 max_candidate_mult=0;
2568 i1=i2=0;
2569 }
2570 }
2571 switch(slice_it->first) {
2572 case(CGAL::internal::FIRST_CURVE): {
2573 i1++;
2574 break;
2575 }
2576 case(CGAL::internal::SECOND_CURVE): {
2577 i2++;
2578 break;
2579 }
2580 case(CGAL::internal::CANDIDATE): {
2581 if(e1.interval_length(i1)<e2.interval_length(i2)) {
2582 e2.refine(i2);
2583 }
2584 else {
2585 e1.refine(i1);
2586 }
2587 if(! overlap(e1,i1,e2,i2)) {
2588 number_of_candidates--;
2589 slice_it=slice.erase(slice_it);
2590 if(e1.lower_bound(i1)<e2.lower_bound(i2)) {
2591 slice_it=slice.insert
2592 (slice_it,std::make_pair(CGAL::internal::FIRST_CURVE,-1));
2593 slice_it++;
2594 slice_it=slice.insert
2595 (slice_it,std::make_pair
2596 (CGAL::internal::SECOND_CURVE,-1));
2597 } else {
2598 slice_it=slice.
2599 insert(slice_it,std::make_pair
2600 (CGAL::internal::SECOND_CURVE,-1));
2601 slice_it++;
2602 slice_it=slice.
2603 insert(slice_it,std::make_pair
2604 (CGAL::internal::FIRST_CURVE,-1));
2605 }
2606 } else {
2607 size_type m1 = e1.get_upper_bound_for_multiplicity(i1),
2608 m2 = e2.get_upper_bound_for_multiplicity(i2),
2609 min_m = m1<m2 ? m1 : m2;
2610 max_candidate_mult = max_candidate_mult>min_m
2611 ? max_candidate_mult : min_m;
2612 }
2613 i1++;
2614 i2++;
2615 break;
2616 }
2617 case(CGAL::internal::INTERSECTION): {
2618 i1++;
2619 i2++;
2620 break;
2621 }
2622 }
2623 slice_it++;
2624 }
2625 return number_of_intersections+number_of_candidates;
2626 }
2627
2628 } //namespace CGAL
2629
2630 #include <CGAL/enable_warnings.h>
2631
2632 #endif
2633