1 // Copyright (c) 2011 Tel-Aviv University (Israel), INRIA Sophia-Antipolis (France).
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/Arrangement_on_surface_2/include/CGAL/Arr_rat_arc/Rational_function.h $
7 // $Id: Rational_function.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 // Author(s)     : Oren Salzman <orenzalz@post.tau.ac.il >
11 //                 Michael Hemmer <Michael.Hemmer@sophia.inria.fr>
12 
13 #ifndef CGAL_RATIONAL_FUNCTION_H
14 #define CGAL_RATIONAL_FUNCTION_H
15 
16 #include <CGAL/license/Arrangement_on_surface_2.h>
17 
18 
19 #include <CGAL/tss.h>
20 #include <CGAL/Arr_rat_arc/Base_rational_arc_ds_1.h>
21 #include <CGAL/Handle_with_policy.h>
22 namespace CGAL {
23 namespace Arr_rational_arc {
24 
25 template <class AlgebraicKernel_d_1 >
26 class Rational_function_rep : public Base_rational_arc_ds_1<AlgebraicKernel_d_1>
27 {
28 public:
29   typedef AlgebraicKernel_d_1                          Algebraic_kernel_d_1;
30   typedef Base_rational_arc_ds_1<Algebraic_kernel_d_1> Base;
31 
32   typedef typename Base::Polynomial_1                  Polynomial_1;
33   typedef typename Base::Algebraic_real_1              Algebraic_real_1;
34   typedef typename Base::Algebraic_vector              Algebraic_vector;
35   typedef typename Base::Multiplicity                  Multiplicity;
36   typedef typename Base::Multiplicity_vector           Multiplicity_vector;
37   typedef typename Base::Root_multiplicity_vector      Root_multiplicity_vector;
38   typedef typename Base::Solve_1                       Solve_1;
39 
40 public:
Rational_function_rep()41   Rational_function_rep() : _ak_ptr(nullptr){}
Rational_function_rep(const Polynomial_1 & numer,const Polynomial_1 & denom,Algebraic_kernel_d_1 * ak_ptr)42   Rational_function_rep(const Polynomial_1& numer,
43                         const Polynomial_1& denom,
44                         Algebraic_kernel_d_1* ak_ptr):
45     _numer(numer), _denom(denom),_ak_ptr(ak_ptr)
46   {
47     initialize();
48   }
49 
50   CGAL::Sign sign_at (const Algebraic_real_1& x,
51                       CGAL::Sign epsilon = CGAL::ZERO) const
52   {
53     //find interval
54     typename Algebraic_vector::const_iterator iter =
55       std::lower_bound(_event_roots.begin(), _event_roots.end(),x);
56 
57     //case of a value larger than largest root
58     if (iter == _event_roots.end())
59       return (_sign.back());
60 
61     typename Algebraic_vector::iterator::difference_type dist =
62       iter - _event_roots.begin();
63 
64     //if x is not a root, ignore epsilons
65     if (*iter != x)
66       return (_sign[dist]);
67 
68     //x is a root
69     if (epsilon == CGAL::ZERO)
70       return (CGAL::EQUAL);
71     else if (epsilon == CGAL::NEGATIVE)
72       return (_sign[dist] );
73     else // CGAL::POSITIVE
74       return (_sign[dist+1]);
75   }
76 
sign_near_minus_infinity()77   CGAL::Sign sign_near_minus_infinity() const
78   {
79     return _sign.front();
80   }
81 
82   bool operator==(const Rational_function_rep& other) const
83   {
84     return ((this->_numer == other.numer()) &&
85             (this->_denom == other.denom()) );
86   }
87 
numer()88   const Polynomial_1& numer() const
89   {
90     return _numer;
91   }
92 
denom()93   const Polynomial_1& denom() const
94   {
95     return _denom;
96   }
97 
poles()98   const Algebraic_vector& poles() const
99   {
100     return _poles;
101   }
102 
pole_multiplicities()103   const Multiplicity_vector& pole_multiplicities() const
104   {
105     return _pole_multiplicities;
106   }
107 
108 private:
initialize()109   void initialize()
110   {
111     CGAL_precondition(_ak_ptr != nullptr);
112     CGAL_precondition(CGAL::is_zero(_denom) == false);
113     if (CGAL::is_zero(_numer))
114     {
115       //function does not change sign
116       _sign.push_back(CGAL::ZERO);
117       return;
118     }
119 
120     Solve_1 solve_1 (_ak_ptr->solve_1_object());
121     Root_multiplicity_vector rm_poles_vec,rm_intersctions_vec;
122     solve_1(_denom, std::back_inserter(rm_poles_vec));  //poles
123     solve_1(_numer, std::back_inserter(rm_intersctions_vec)); //intersections with zero
124 
125     //reserve memory
126     typename Root_multiplicity_vector::size_type num_of_poles =
127       rm_poles_vec.size();
128     typename Root_multiplicity_vector::size_type num_of_intersections =
129       rm_intersctions_vec.size();
130 
131     _poles.reserve(num_of_poles);
132     _pole_multiplicities.reserve(num_of_poles);
133 
134     _event_roots.reserve(num_of_poles + num_of_intersections);
135     _pole_multiplicities.reserve(num_of_poles + num_of_intersections);
136 
137     //initialize poles
138     for ( typename Root_multiplicity_vector::iterator it = rm_poles_vec.begin();
139           it != rm_poles_vec.end() ;
140           ++it)
141     {
142       _poles.push_back(it->first);
143       _pole_multiplicities.push_back(it->second);
144     }
145 
146     //initialize events
147     Root_multiplicity_vector events_mult_vec;
148     std::merge( rm_poles_vec.begin()  , rm_poles_vec.end()  ,
149                 rm_intersctions_vec.begin() , rm_intersctions_vec.end() ,
150                 std::back_inserter(events_mult_vec));
151 
152     typename Root_multiplicity_vector::iterator it;
153     for (it = events_mult_vec.begin(); it != events_mult_vec.end(); ++it)
154     {
155       _event_roots.push_back(it->first);
156       _event_multiplicities.push_back(it->second);
157     }
158 
159     //initialize left most interval (at -oo)
160     //for (ax^n+.../x^m+...) the sign is:  sign(a) * (-1)^(n +m)
161     CGAL::Sign curr_sign = CGAL::sign(CGAL::leading_coefficient(_numer));
162     if ((CGAL::degree(_numer) + CGAL::degree(_denom)) % 2 == 1)
163       curr_sign = curr_sign * CGAL::NEGATIVE;
164     _sign.push_back(curr_sign);
165 
166     typename Multiplicity_vector::iterator it2;
167     for (it2 = _event_multiplicities.begin();
168          it2 != _event_multiplicities.end();
169          ++it2)
170     {
171       if (*it2 % 2 == 1)
172         curr_sign = curr_sign * CGAL::NEGATIVE;
173       _sign.push_back(curr_sign);
174     }
175   }
176 
177 private:
178   Polynomial_1        _numer;
179   Polynomial_1        _denom;
180   Algebraic_vector    _poles;     //roots of the denominator
181   Multiplicity_vector _pole_multiplicities; //multiplicities of the poles
182   Algebraic_vector    _event_roots;   //poles and intersection points with y=0, function can change signs in these events
183   Multiplicity_vector _event_multiplicities; //multiplicities of the events
184   std::vector<CGAL::Sign> _sign;    //function's sign in the corresponding interval induced by _event_roots (if no roots then only one value)
185   mutable Algebraic_kernel_d_1*   _ak_ptr;
186 
187 };//Rational_function_rep
188 
189 template < class Algebraic_kernel_ >
190 class Rational_function:
191     public Handle_with_policy<Rational_function_rep<Algebraic_kernel_> >
192 {
193 public:
194   typedef Algebraic_kernel_                             Algebraic_kernel_d_1;
195   typedef Handle_with_policy<Rational_function_rep<Algebraic_kernel_> >
196                                                         Base;
197   typedef Rational_function<Algebraic_kernel_d_1>       Self;
198   typedef Rational_function_rep<Algebraic_kernel_>      Rep;
199   typedef typename Rep::Algebraic_real_1                Algebraic_real_1;
200   typedef typename Rep::Polynomial_1                    Polynomial_1;
201   typedef typename Rep::Algebraic_vector                Algebraic_vector;
202   typedef typename Rep::Multiplicity_vector             Multiplicity_vector;
203 
204   typedef typename Base::Id_type                        Id_type;
205 private:
get_default_instance()206   static Self& get_default_instance()
207   {
208     CGAL_STATIC_THREAD_LOCAL_VARIABLE_0(Algebraic_kernel_d_1, kernel);
209     CGAL_STATIC_THREAD_LOCAL_VARIABLE_3(Self, x, Polynomial_1(0), Polynomial_1(1), &kernel);
210     return x;
211   }
212 public:
Rational_function(const Polynomial_1 & numer,const Polynomial_1 & denom,Algebraic_kernel_d_1 * ak_ptr)213   Rational_function(const Polynomial_1& numer,
214                     const Polynomial_1& denom,
215                     Algebraic_kernel_d_1* ak_ptr) :
216     Base(numer,denom,ak_ptr) {}
217 
218   //used to solve VS bug...
Rational_function()219   Rational_function () :
220     Base(static_cast<const Base &> (get_default_instance())) {}
221 
222   // explicit copy-constructor, required by VC9
Rational_function(const Self & r)223   Rational_function (const Self & r)
224     : Base(static_cast<const Base &> (r)) {}
225 
226   Self& operator=(const Self&)=default;
227 
228   CGAL::Sign sign_at(const Algebraic_real_1& x,
229                      CGAL::Sign epsilon = CGAL::ZERO) const
230   {
231     return this->ptr()->sign_at(x,epsilon);
232   }
233 
sign_near_minus_infinity()234   CGAL::Sign sign_near_minus_infinity() const
235   {
236     return this->ptr()->sign_near_minus_infinity ();
237   }
238   bool operator== (const Self & other) const
239   {
240     if (this->is_identical (other))
241       return true;
242     return (*(this->ptr()) == *(other.ptr()));
243   }
244 
numer()245   const Polynomial_1& numer () const
246   {
247     return this->ptr()->numer();
248   }
249 
denom()250   const Polynomial_1& denom () const
251   {
252     return this->ptr()->denom();
253   }
254 
poles()255   const Algebraic_vector& poles () const
256   {
257     return this->ptr()->poles();
258   }
259 
pole_multiplicities()260   const Multiplicity_vector& pole_multiplicities () const
261   {
262     return this->ptr()->pole_multiplicities();
263   }
264 };  //Rational_function
265 
266 }   //namespace Arr_rational_arc
267 }   //namespace CGAL {
268 
269 #endif //CGAL_RATIONAL_FUNCTION_H
270