← Back to branch summary

~ubuntu-branches/debian/stretch/cgal/stretch

« back to all changes in this revision

Viewing changes to include/CGAL/RS/functors_1.h

  • Committer: Package Import Robot
  • Author(s): Joachim Reichel
  • Date: 2014-04-05 10:56:43 UTC
  • mfrom: (1.2.4)
  • Revision ID: package-import@ubuntu.com-20140405105643-jgnrpu2thtx23zfs
Tags: 4.4-1
* New upstream release.
* Remove patches do-not-link-example-with-qt4-support-library.patch and
  fix_jet_fitting_3.patch (applied upstream).

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/CGAL/RS/functors_1.h
1
 
// Copyright (c) 2006-2010 Inria Lorraine (France). All rights reserved.
 
1
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). All rights reserved.
2
2
//
3
3
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
4
4
// modify it under the terms of the GNU Lesser General Public License as
5
5
// published by the Free Software Foundation; either version 3 of the License,
6
6
// or (at your option) any later version.
 
7
 
 
8
// See the file LICENSE.LGPL distributed with CGAL.
7
9
//
8
10
// Licensees holding a valid commercial license may use this file in
9
11
// accordance with the commercial license agreement provided with the software.
16
18
//
17
19
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
18
20
 
19
 
#ifndef CGAL_RS_FUNCTORS_H
20
 
#define CGAL_RS_FUNCTORS_H
21
 
 
22
 
#include <CGAL/basic.h>
23
 
#include <mpfi.h>
24
 
#include <CGAL/RS/polynomial_1.h>
25
 
#include <CGAL/RS/algebraic_1.h>
26
 
#include <CGAL/RS/polynomial_1_utils.h>
27
 
#include <CGAL/RS/solve_1.h>
28
 
#include <CGAL/RS/sign_1.h>
29
 
#include <CGAL/RS/sign_1_rs.h>
30
 
#include <CGAL/RS/refine_1_rs.h>
31
 
#include <CGAL/RS/compare_1.h>
32
 
#include <CGAL/RS/polynomial_converter.h>
33
 
#include <CGAL/Gmpfr.h>
34
 
 
35
 
#ifdef IEEE_DBL_MANT_DIG
36
 
#  define CGAL_RS_FUNCTORS_DBL_PREC IEEE_DBL_MANT_DIG
37
 
#else
38
 
#  define CGAL_RS_FUNCTORS_DBL_PREC 53
39
 
#endif
40
 
 
41
 
namespace RSFunctors{
42
 
 
43
 
typedef CGAL::RS_polynomial_1           Polynomial;
44
 
typedef CGAL::Algebraic_1               Algebraic;
45
 
typedef CGAL::Gmpfr                     Bound;
46
 
typedef int                             Multiplicity;
47
 
 
48
 
template <class PolynomialType>
 
21
#ifndef CGAL_RS_FUNCTORS_1_H
 
22
#define CGAL_RS_FUNCTORS_1_H
 
23
 
 
24
#include <vector>
 
25
#include <CGAL/Gmpfi.h>
 
26
 
 
27
namespace CGAL{
 
28
namespace RS_AK1{
 
29
 
 
30
template <class Polynomial_,
 
31
          class Algebraic_,
 
32
          class Bound_,
 
33
          class Coefficient_,
 
34
          class Isolator_>
 
35
struct Construct_algebraic_real_1{
 
36
        typedef Polynomial_                                     Polynomial;
 
37
        typedef Algebraic_                                      Algebraic;
 
38
        typedef Bound_                                          Bound;
 
39
        typedef Coefficient_                                    Coefficient;
 
40
        typedef Isolator_                                       Isolator;
 
41
        typedef Algebraic                                       result_type;
 
42
 
 
43
        template <class T>
 
44
        Algebraic operator()(const T &a)const{
 
45
                return Algebraic(a);
 
46
        }
 
47
 
 
48
        Algebraic operator()(const Polynomial &p,size_t i)const{
 
49
                Isolator isol(p);
 
50
                return Algebraic(p,isol.left_bound(i),isol.right_bound(i));
 
51
        }
 
52
 
 
53
        Algebraic operator()(const Polynomial &p,
 
54
                             const Bound &l,
 
55
                             const Bound &r)const{
 
56
                return Algebraic(p,l,r);
 
57
        }
 
58
 
 
59
}; // struct Construct_algebraic_1
 
60
 
 
61
template <class Polynomial_,class Algebraic_>
49
62
struct Compute_polynomial_1:
50
 
public std::unary_function<Algebraic,PolynomialType>{
51
 
        typedef PolynomialType          P;
52
 
        typedef CGAL::from_rs_poly<P>   back;
53
 
        P operator()(const Algebraic &a)const{
54
 
                return back()(a.pol());
55
 
        }
56
 
};      // Compute_polynomial_1
57
 
 
58
 
template <class PolynomialType,class GcdPolicy>
59
 
struct Is_square_free_1:
60
 
public std::unary_function<PolynomialType,bool>{
61
 
        typedef PolynomialType          P;
62
 
        typedef CGAL::to_rs_poly<P>     convert;
63
 
        typedef GcdPolicy               Gcd;
64
 
        bool operator()(const P &p)const{
65
 
                Polynomial rsp=convert()(p);
66
 
                return(!(Gcd()(rsp,rsp.derive()).get_degree_static()));
67
 
        }
68
 
};      // Is_square_free_1
69
 
 
70
 
template <class PolynomialType,class GcdPolicy>
71
 
struct Make_square_free_1:
72
 
public std::unary_function<PolynomialType,PolynomialType>{
73
 
        typedef PolynomialType          P;
74
 
        typedef CGAL::to_rs_poly<P>     convert;
75
 
        typedef CGAL::from_rs_poly<P>   back;
76
 
        typedef GcdPolicy               Gcd;
77
 
        P operator()(const P &p)const{
78
 
                return back()(CGAL::sfpart_1<Gcd>()(convert()(p)));
79
 
        }
80
 
};      // Make_square_free_1
81
 
 
82
 
template <class PolynomialType,class GcdPolicy>
83
 
struct Square_free_factorize_1{
84
 
        typedef PolynomialType          P;
85
 
        typedef CGAL::to_rs_poly<P>     convert;
86
 
        typedef CGAL::from_rs_poly<P>   back;
87
 
        typedef GcdPolicy               Gcd;
88
 
        template <class OutputIterator>
89
 
        OutputIterator operator()(const P &p,OutputIterator oi)const{
90
 
                Polynomial rsp=convert()(p);
91
 
                CGAL::sqfrvec factorization(CGAL::sqfr_1<Gcd>()(rsp));
92
 
                for(CGAL::sqfrvec::iterator i=factorization.begin();
93
 
                    i!=factorization.end();
94
 
                    ++i){
95
 
                        *oi++=std::make_pair(back()((*i).first),(*i).second);
96
 
                    }
97
 
                return oi;
98
 
        }
99
 
};      // Square_free_factorize_1
100
 
 
101
 
template <class PolynomialType,class GcdPolicy>
 
63
public std::unary_function<Algebraic_,Polynomial_>{
 
64
        typedef Polynomial_                                     Polynomial;
 
65
        typedef Algebraic_                                      Algebraic;
 
66
        Polynomial operator()(const Algebraic &x)const{
 
67
                return x.get_pol();
 
68
        }
 
69
}; // struct Compute_polynomial_1
 
70
 
 
71
template <class Polynomial_,class Ptraits_>
102
72
struct Is_coprime_1:
103
 
public std::binary_function<PolynomialType,PolynomialType,bool>{
104
 
        typedef PolynomialType          P;
105
 
        typedef CGAL::to_rs_poly<P>     convert;
106
 
        typedef GcdPolicy               Gcd;
107
 
        bool operator()(const P &p1,const P &p2)const{
108
 
                CGAL::RS_polynomial_1 rsp1=convert()(p1);
109
 
                CGAL::RS_polynomial_1 rsp2=convert()(p2);
110
 
                return(!Gcd()(rsp1,rsp2).get_degree_static());
 
73
public std::binary_function<Polynomial_,Polynomial_,bool>{
 
74
        typedef Polynomial_                                     Polynomial;
 
75
        typedef Ptraits_                                        Ptraits;
 
76
        typedef typename Ptraits::Gcd_up_to_constant_factor     Gcd;
 
77
        typedef typename Ptraits::Degree                        Degree;
 
78
        inline bool operator()(const Polynomial &p1,const Polynomial &p2)const{
 
79
                return Degree()(Gcd()(p1,p2))==0;
111
80
        }
112
 
};      // Is_coprime_1
 
81
}; // struct Is_coprime_1
113
82
 
114
 
template <class PolynomialType,class GcdPolicy>
 
83
template <class Polynomial_,class Ptraits_>
115
84
struct Make_coprime_1{
116
 
        typedef PolynomialType          P;
117
 
        typedef CGAL::to_rs_poly<P>     convert;
118
 
        typedef CGAL::from_rs_poly<P>   back;
119
 
        typedef GcdPolicy               Gcd;
120
 
        bool operator()(const P &p1,const P &p2,P &g,P &q1,P &q2)const{
121
 
                typedef GcdPolicy       Gcd;
122
 
                CGAL::RS_polynomial_1 rsp1=convert()(p1);
123
 
                CGAL::RS_polynomial_1 rsp2=convert()(p2);
124
 
                CGAL::RS_polynomial_1 rsg=convert()(g);
125
 
                rsg=Gcd()(rsp1,rsp2);
126
 
                g=back()(rsg);
127
 
                // even when g==1, we calculate q1 and q2
128
 
                q1=back()(*CGAL::Ediv_1()(rsp1,rsg));
129
 
                q2=back()(*CGAL::Ediv_1()(rsp2,rsg));
130
 
                return rsg.get_degree_static()?false:true;
131
 
        }
132
 
};      // Make_coprime_1
133
 
 
134
 
template <class GcdPolicy>
135
 
struct Solve_RS_1{
136
 
    typedef GcdPolicy           Gcd;
137
 
    typedef CGAL::sfpart_1<Gcd> Sfpart;
138
 
    typedef CGAL::sqfr_1<Gcd>   Sqfr;
139
 
 
140
 
    template <class OutputIterator>
141
 
    OutputIterator operator()(const Polynomial &p,OutputIterator res)const{
142
 
        int nr,*m;
143
 
        mpfi_ptr *x;
144
 
        CGAL::sqfrvec sfv=Sqfr()(p);
145
 
        x=(mpfi_ptr*)malloc(Sfpart()(p).get_degree()*sizeof(mpfi_ptr));
146
 
        m=(int*)calloc(Sfpart()(p).get_degree(),sizeof(int));
147
 
        nr=solve_1(x,Sfpart()(p));
148
 
        CGAL_assertion_msg(nr>=0,"error in resolution");
149
 
        for(CGAL::sqfrvec::size_type i=0;i<sfv.size();++i){
150
 
            int k=sfv[i].first.get_degree();
151
 
            for(int j=0;k&&j<nr;++j){
152
 
                if(!m[j]){
153
 
                    CGAL::Sign sg_l=
154
 
                        CGAL::RSSign::signat(sfv[i].first,&(x[j]->left));
155
 
                    CGAL::Sign sg_r=
156
 
                        CGAL::RSSign::signat(sfv[i].first,&(x[j]->right));
157
 
                    if((sg_l!=sg_r)||((sg_l==CGAL::ZERO)&&(sg_r==CGAL::ZERO))){
158
 
                        m[j]=sfv[i].second;
159
 
                        --k;
160
 
                    }
161
 
                }
162
 
            }
163
 
        }
164
 
        for(int i=0;i<nr;++i){
165
 
            *res++=std::make_pair(*new Algebraic(x[i],p,i,m[i]
166
 
                                                 //,i?x[i-1]:NULL,
167
 
                                                 //i==nr-1?NULL:x[i+1]
168
 
                                                 ),
169
 
                                  m[i]);
170
 
        }
171
 
        free(m);
172
 
        free(x);
173
 
        return res;
174
 
    }
175
 
 
176
 
    template <class OutputIterator>
177
 
    OutputIterator operator()(const Polynomial &p,
178
 
                              bool known_to_be_square_free,
179
 
                              OutputIterator res)const{
180
 
        int nr,m;
181
 
        mpfi_ptr *x;
182
 
        if(known_to_be_square_free){
183
 
            p.set_sf();
184
 
            x=(mpfi_ptr*)malloc(p.get_degree()*sizeof(mpfi_ptr));
185
 
            nr=solve_1(x,p);
186
 
            CGAL_assertion_msg(nr>=0,"error in resolution");
187
 
            m=1;    // we know in this case that multiplicity is 1
188
 
        }else{
189
 
            x=(mpfi_ptr*)malloc(Sfpart()(p).get_degree()*sizeof(mpfi_ptr));
190
 
            nr=solve_1(x,Sfpart()(p));
191
 
            CGAL_assertion_msg(nr>=0,"error in resolution");
192
 
            m=0;    // we won't calculate multiplicities
193
 
        }
194
 
        for(int i=0;i<nr;++i)
195
 
            *res++=*new Algebraic(x[i],p,i,m
196
 
                            //,i?x[i-1]:NULL,
197
 
                            //i==nr-1?NULL:x[i+1]
198
 
                            );
199
 
        free(x);
200
 
        return res;
201
 
    }
202
 
};  // Solve_RS_1
203
 
 
204
 
template <class PolynomialType,class GcdPolicy>
 
85
        typedef Polynomial_                                     Polynomial;
 
86
        typedef Ptraits_                                        Ptraits;
 
87
        typedef typename Ptraits::Gcd_up_to_constant_factor     Gcd;
 
88
        typedef typename Ptraits::Degree                        Degree;
 
89
        typedef typename Ptraits::Integral_division_up_to_constant_factor
 
90
                                                                IDiv;
 
91
        bool operator()(const Polynomial &p1,
 
92
                        const Polynomial &p2,
 
93
                        Polynomial &g,
 
94
                        Polynomial &q1,
 
95
                        Polynomial &q2)const{
 
96
                g=Gcd()(p1,p2);
 
97
                q1=IDiv()(p1,g);
 
98
                q2=IDiv()(p2,g);
 
99
                return Degree()(Gcd()(p1,p2))==0;
 
100
        }
 
101
}; // struct Make_coprime_1
 
102
 
 
103
template <class Polynomial_,
 
104
          class Bound_,
 
105
          class Algebraic_,
 
106
          class Isolator_,
 
107
          class Signat_,
 
108
          class Ptraits_>
205
109
struct Solve_1{
206
 
    typedef PolynomialType      P;
207
 
    typedef CGAL::to_rs_poly<P> convert;
208
 
    typedef GcdPolicy           Gcd;
209
 
    typedef Solve_RS_1<Gcd>     Solve_RS;
210
 
 
211
 
    template <class OutputIterator>
212
 
    OutputIterator operator()(const P &p,OutputIterator res)const{
213
 
        return Solve_RS()(convert()(p),res);
214
 
    }
215
 
 
216
 
    template <class OutputIterator>
217
 
    OutputIterator operator()(const P &p,
218
 
                              bool known_to_be_square_free,
219
 
                              OutputIterator res)const{
220
 
        return Solve_RS()(convert()(p),known_to_be_square_free,res);
221
 
    }
222
 
 
223
 
  template <class OutputIterator>
224
 
  OutputIterator operator()(const P &p,
225
 
                            const Bound& lower,
226
 
                            const Bound& upper,
227
 
                            OutputIterator res)const{
228
 
    typedef std::vector<std::pair<Algebraic,Multiplicity> > RMV;
229
 
    RMV roots;
230
 
    this->operator()(p,std::back_inserter(roots));
231
 
    for(typename RMV::iterator it = roots.begin(); it != roots.end();it++){
232
 
      if(lower <= it->first && it->first <= upper)
233
 
        *res++=*it;
234
 
    }
235
 
    return res;
236
 
  }
237
 
 
238
 
  template <class OutputIterator>
239
 
  OutputIterator operator()(const P &p,
240
 
                            bool known_to_be_square_free,
241
 
                            const Bound& lower,
242
 
                            const Bound& upper,
243
 
                            OutputIterator res)const{
244
 
    typedef std::vector< Algebraic > RV;
245
 
    RV roots;
246
 
    this->operator()(p,known_to_be_square_free,std::back_inserter(roots));
247
 
    for(typename RV::iterator it = roots.begin(); it != roots.end();it++){
248
 
      if(lower <= *it && *it <= upper)
249
 
        *res++=*it;
250
 
    }
251
 
    return res;
252
 
  }
253
 
};  // Solve_1
254
 
 
255
 
template <class PolynomialType,class CoeffType,class GcdFunction>
256
 
struct Construct_alg_1{
257
 
    typedef PolynomialType              Poly;
258
 
    typedef CoeffType                   Coeff;
259
 
    typedef GcdFunction                 Gcd;
260
 
    typedef CGAL::to_rs_poly<Poly>      convert;
261
 
    typedef Solve_1<Poly,Gcd>           Solve;
262
 
 
263
 
    Algebraic operator()(int a) const {
264
 
        return Algebraic(a);
265
 
    }
266
 
 
267
 
    Algebraic operator()(const Bound a) const {
268
 
        return Algebraic(a);
269
 
    }
270
 
 
271
 
    Algebraic operator()(const Coeff a) const {
272
 
        return Algebraic(a);
273
 
    }
274
 
 
275
 
    Algebraic operator()(const Poly &p,int i) const {
276
 
      CGAL_precondition(CGAL::is_square_free(p));
277
 
      std::vector<Algebraic> roots;
278
 
      std::back_insert_iterator<std::vector<Algebraic> > rootsit(roots);
279
 
      Solve()(p,true,rootsit);
280
 
      return roots[i];
281
 
    }
282
 
 
283
 
    Algebraic operator()(const Poly &p,Bound l,Bound u) const {
284
 
        mpfi_t i;
285
 
        mpfi_init(i);
286
 
        mpfr_set(&i->left,l.fr(),GMP_RNDD);
287
 
        mpfr_set(&i->right,u.fr(),GMP_RNDU);
288
 
        return Algebraic(i,convert()(p),0,0
289
 
                         //,NULL,NULL
290
 
                        );
291
 
    }
292
 
};  // Construct_alg_1
293
 
 
294
 
template <class PolynomialType>
 
110
        typedef Polynomial_                                     Polynomial_1;
 
111
        typedef Bound_                                          Bound;
 
112
        typedef Algebraic_                                      Algebraic;
 
113
        typedef Isolator_                                       Isolator;
 
114
        typedef Signat_                                         Signat;
 
115
        typedef Ptraits_                                        Ptraits;
 
116
        typedef typename Ptraits::Gcd_up_to_constant_factor     Gcd;
 
117
        typedef typename Ptraits::Square_free_factorize_up_to_constant_factor
 
118
                                                                Sqfr;
 
119
        typedef typename Ptraits::Degree                        Degree;
 
120
        typedef typename Ptraits::Make_square_free              Sfpart;
 
121
 
 
122
        template <class OutputIterator>
 
123
        OutputIterator operator()(const Polynomial_1 &p,
 
124
                                  OutputIterator res)const{
 
125
                typedef std::pair<Polynomial_1,int>     polmult;
 
126
                typedef std::vector<polmult>            sqvec;
 
127
 
 
128
                Polynomial_1 sfp=Sfpart()(p);
 
129
                sqvec sfv;
 
130
                Sqfr()(p,std::back_inserter(sfv));
 
131
                Isolator isol(sfp);
 
132
                int *m=(int*)calloc(isol.number_of_real_roots(),sizeof(int));
 
133
                for(typename sqvec::iterator i=sfv.begin();i!=sfv.end();++i){
 
134
                        int k=Degree()(i->first);
 
135
                        Signat signof(i->first);
 
136
                        for(int j=0;k&&j<isol.number_of_real_roots();++j){
 
137
                                if(!m[j]){
 
138
                                        CGAL::Sign sg_l=
 
139
                                                signof(isol.left_bound(j));
 
140
                                        CGAL::Sign sg_r=
 
141
                                                signof(isol.right_bound(j));
 
142
                                        if((sg_l!=sg_r)||
 
143
                                           ((sg_l==CGAL::ZERO)&&
 
144
                                            (sg_r==CGAL::ZERO))){
 
145
                                                m[j]=i->second;
 
146
                                                --k;
 
147
                                        }
 
148
                                }
 
149
                        }
 
150
                }
 
151
                for(int l=0;l<isol.number_of_real_roots();++l)
 
152
                        *res++=std::make_pair(Algebraic(p,
 
153
                                                        isol.left_bound(l),
 
154
                                                        isol.right_bound(l)),
 
155
                                              m[l]);
 
156
                free(m);
 
157
                return res;
 
158
        }
 
159
 
 
160
        template <class OutputIterator>
 
161
        OutputIterator operator()(const Polynomial_1 &p,
 
162
                                  bool known_to_be_square_free,
 
163
                                  OutputIterator res)const{
 
164
                Isolator isol(p);
 
165
                for(int l=0;l<isol.number_of_real_roots();++l)
 
166
                        *res++=Algebraic(p,
 
167
                                         isol.left_bound(l),
 
168
                                         isol.right_bound(l));
 
169
                return res;
 
170
        }
 
171
 
 
172
        template <class OutputIterator>
 
173
        OutputIterator operator()(const Polynomial_1 &p,
 
174
                                  const Bound &l,
 
175
                                  const Bound &u,
 
176
                                  OutputIterator res)const{
 
177
                typedef std::vector<Algebraic>                  RV;
 
178
                typedef std::pair<Polynomial_1,int>             PM;
 
179
                typedef std::vector<PM>                         PMV;
 
180
                typedef typename PMV::iterator                  PMVI;
 
181
                CGAL_precondition_msg(l<=u,
 
182
                                      "left bound must be <= right bound");
 
183
                RV roots; // all roots of the polynomial
 
184
                this->operator()(p,false,std::back_inserter(roots));
 
185
                size_t nb_roots=roots.size();
 
186
                // indices of the first and last roots to be reported:
 
187
                size_t index_l=0,index_u;
 
188
                while(index_l<nb_roots&&roots[index_l]<l)
 
189
                        ++index_l;
 
190
                CGAL_assertion(index_l<=nb_roots);
 
191
                if(index_l==nb_roots)
 
192
                        return res;
 
193
                index_u=index_l;
 
194
                while(index_u<nb_roots&&roots[index_u]<u)
 
195
                        ++index_u;
 
196
                CGAL_assertion(index_u<=nb_roots);
 
197
                if(index_u==index_l)
 
198
                        return res;
 
199
                // now, we have to return roots in [index_l,index_u)
 
200
                PMV sfv;
 
201
                Sqfr()(p,std::back_inserter(sfv)); // square-free fact. of p
 
202
                // array to store the multiplicities
 
203
                int *m=(int*)calloc(nb_roots,sizeof(int));
 
204
                // we iterate over all the pairs <root,mult> and match the
 
205
                // roots in the interval [index_l,index_u)
 
206
                for(PMVI i=sfv.begin();i!=sfv.end();++i){
 
207
                        int k=Degree()(i->first);
 
208
                        Signat signof(i->first);
 
209
                        for(size_t j=index_l;k&&j<index_u;++j){
 
210
                                if(!m[j]){
 
211
                                        CGAL::Sign sg_l=
 
212
                                                signof(roots[j].get_left());
 
213
                                        CGAL::Sign sg_r=
 
214
                                                signof(roots[j].get_right());
 
215
                                        if((sg_l!=sg_r)||
 
216
                                           ((sg_l==CGAL::ZERO)&&
 
217
                                            (sg_r==CGAL::ZERO))){
 
218
                                                m[j]=i->second;
 
219
                                                --k;
 
220
                                        }
 
221
                                }
 
222
                        }
 
223
                }
 
224
                for(size_t l=index_l;l<index_u;++l)
 
225
                        *res++=std::make_pair(roots[l],m[l]);
 
226
                free(m);
 
227
                return res;
 
228
        }
 
229
 
 
230
        template <class OutputIterator>
 
231
        OutputIterator operator()(const Polynomial_1 &p,
 
232
                                  bool known_to_be_square_free,
 
233
                                  const Bound &l,
 
234
                                  const Bound &u,
 
235
                                  OutputIterator res)const{
 
236
                typedef std::vector<Algebraic>                  RV;
 
237
                typedef typename RV::iterator                   RVI;
 
238
                CGAL_precondition_msg(l<=u,
 
239
                                      "left bound must be <= right bound");
 
240
                RV roots;
 
241
                this->operator()(p,
 
242
                                 known_to_be_square_free,
 
243
                                 std::back_inserter(roots));
 
244
                for(RVI it=roots.begin();it!=roots.end();it++)
 
245
                        if(*it>=l&&*it<=u)
 
246
                                *res++=*it;
 
247
                return res;
 
248
        }
 
249
 
 
250
}; // struct Solve_1
 
251
 
 
252
template <class Polynomial_,
 
253
          class Bound_,
 
254
          class Algebraic_,
 
255
          class Refiner_,
 
256
          class Signat_,
 
257
          class Ptraits_>
 
258
class Sign_at_1:
 
259
public std::binary_function<Polynomial_,Algebraic_,CGAL::Sign>{
 
260
        // This implementation will work with any polynomial type whose
 
261
        // coefficient type is explicit interoperable with Gmpfi.
 
262
        // TODO: Make this function generic.
 
263
        public:
 
264
        typedef Polynomial_                                     Polynomial_1;
 
265
        typedef Bound_                                          Bound;
 
266
        typedef Algebraic_                                      Algebraic;
 
267
        typedef Refiner_                                        Refiner;
 
268
        typedef Signat_                                         Signat;
 
269
        typedef Ptraits_                                        Ptraits;
 
270
 
 
271
        private:
 
272
        CGAL::Uncertain<CGAL::Sign> eval_interv(const Polynomial_1 &p,
 
273
                                                const Bound &l,
 
274
                                                const Bound &r)const{
 
275
                typedef typename Ptraits::Substitute                    Subst;
 
276
                std::vector<CGAL::Gmpfi> substitutions;
 
277
                substitutions.push_back(CGAL::Gmpfi(l,r));
 
278
                CGAL::Gmpfi eval=Subst()(p,
 
279
                                         substitutions.begin(),
 
280
                                         substitutions.end());
 
281
                return eval.sign();
 
282
        }
 
283
 
 
284
        // This function assumes that the sign of the evaluation is not zero,
 
285
        // it just refines x until having the correct sign.
 
286
        CGAL::Sign refine_and_return(const Polynomial_1 &p,Algebraic x)const{
 
287
                CGAL::Gmpfr xl(x.get_left());
 
288
                CGAL::Gmpfr xr(x.get_right());
 
289
                CGAL::Uncertain<CGAL::Sign> s;
 
290
                for(;;){
 
291
                        Refiner()(x.get_pol(),
 
292
                                  xl,
 
293
                                  xr,
 
294
                                  2*CGAL::max(xl.get_precision(),
 
295
                                              xr.get_precision()));
 
296
                        s=eval_interv(p,xl,xr);
 
297
                        if(!s.is_same(Uncertain<CGAL::Sign>::indeterminate())){
 
298
                                x.set_left(xl);
 
299
                                x.set_right(xr);
 
300
                                return s;
 
301
                        }
 
302
                }
 
303
        }
 
304
 
 
305
        public:
 
306
        CGAL::Sign operator()(const Polynomial_1 &p,Algebraic x)const{
 
307
                typedef typename Ptraits::Gcd_up_to_constant_factor     Gcd;
 
308
                typedef typename Ptraits::Make_square_free              Sfpart;
 
309
                typedef typename Ptraits::Degree                        Degree;
 
310
                typedef typename Ptraits::Differentiate                 Deriv;
 
311
                CGAL::Uncertain<CGAL::Sign> unknown=
 
312
                                        Uncertain<CGAL::Sign>::indeterminate();
 
313
                CGAL::Uncertain<CGAL::Sign> s=eval_interv(p,
 
314
                                                          x.get_left(),
 
315
                                                          x.get_right());
 
316
                if(!s.is_same(unknown))
 
317
                        return s;
 
318
                // We are not sure about the sign. We calculate the gcd in
 
319
                // order to know if both polynomials have common roots.
 
320
                Polynomial_1 sfpp=Sfpart()(p);
 
321
                Polynomial_1 gcd=Gcd()(sfpp,Sfpart()(x.get_pol()));
 
322
                if(Degree()(gcd)==0)
 
323
                        return refine_and_return(p,x);
 
324
 
 
325
                // At this point, gcd is not 1; we proceed as follows:
 
326
                // -we refine x until having p monotonic in x's interval (to be
 
327
                // sure that p has at most one root on that interval),
 
328
                // -if the gcd has a root on this interval, both roots are
 
329
                // equal (we return 0), otherwise, we refine until having a
 
330
                // result.
 
331
 
 
332
                // How to assure that p is monotonic in an interval: when its
 
333
                // derivative is never zero in that interval.
 
334
                Polynomial_1 dsfpp=Deriv()(sfpp);
 
335
                CGAL::Gmpfr xl(x.get_left());
 
336
                CGAL::Gmpfr xr(x.get_right());
 
337
                while(eval_interv(dsfpp,xl,xr).is_same(unknown)){
 
338
                        Refiner()(x.get_pol(),
 
339
                                  xl,
 
340
                                  xr,
 
341
                                  2*CGAL::max(xl.get_precision(),
 
342
                                              xr.get_precision()));
 
343
                }
 
344
                x.set_left(xl);
 
345
                x.set_right(xr);
 
346
 
 
347
                // How to know that the gcd has a root: evaluate endpoints.
 
348
                CGAL::Sign sleft,sright;
 
349
                Signat sign_at_gcd(gcd);
 
350
                if((sleft=sign_at_gcd(x.get_left()))==CGAL::ZERO||
 
351
                   (sright=sign_at_gcd(x.get_right()))==CGAL::ZERO||
 
352
                   (sleft!=sright))
 
353
                        return CGAL::ZERO;
 
354
                return refine_and_return(p,x);
 
355
        }
 
356
}; // struct Sign_at_1
 
357
 
 
358
template <class Polynomial_,
 
359
          class Bound_,
 
360
          class Algebraic_,
 
361
          class Refiner_,
 
362
          class Signat_,
 
363
          class Ptraits_>
 
364
class Is_zero_at_1:
 
365
public std::binary_function<Polynomial_,Algebraic_,bool>{
 
366
        // This implementation will work with any polynomial type whose
 
367
        // coefficient type is explicit interoperable with Gmpfi.
 
368
        // TODO: Make this function generic.
 
369
        public:
 
370
        typedef Polynomial_                                     Polynomial_1;
 
371
        typedef Bound_                                          Bound;
 
372
        typedef Algebraic_                                      Algebraic;
 
373
        typedef Refiner_                                        Refiner;
 
374
        typedef Signat_                                         Signat;
 
375
        typedef Ptraits_                                        Ptraits;
 
376
 
 
377
        private:
 
378
        CGAL::Uncertain<CGAL::Sign> eval_interv(const Polynomial_1 &p,
 
379
                                                const Bound &l,
 
380
                                                const Bound &r)const{
 
381
                typedef typename Ptraits::Substitute                    Subst;
 
382
                std::vector<CGAL::Gmpfi> substitutions;
 
383
                substitutions.push_back(CGAL::Gmpfi(l,r));
 
384
                CGAL::Gmpfi eval=Subst()(p,     
 
385
                                         substitutions.begin(),
 
386
                                         substitutions.end());
 
387
                return eval.sign();
 
388
        }
 
389
 
 
390
        public:
 
391
        bool operator()(const Polynomial_1 &p,Algebraic x)const{
 
392
                typedef typename Ptraits::Gcd_up_to_constant_factor     Gcd;
 
393
                typedef typename Ptraits::Make_square_free              Sfpart;
 
394
                typedef typename Ptraits::Degree                        Degree;
 
395
                typedef typename Ptraits::Differentiate                 Deriv;
 
396
                CGAL::Uncertain<CGAL::Sign> unknown=
 
397
                                        Uncertain<CGAL::Sign>::indeterminate();
 
398
                CGAL::Uncertain<CGAL::Sign> s=eval_interv(p,
 
399
                                                          x.get_left(),
 
400
                                                          x.get_right());
 
401
                if(!s.is_same(unknown))
 
402
                        return (s==CGAL::ZERO);
 
403
                // We are not sure about the sign. We calculate the gcd in
 
404
                // order to know if both polynomials have common roots.
 
405
                Polynomial_1 sfpp=Sfpart()(p);
 
406
                Polynomial_1 gcd=Gcd()(sfpp,Sfpart()(x.get_pol()));
 
407
                if(Degree()(gcd)==0)
 
408
                        return false;
 
409
 
 
410
                // At this point, gcd is not 1; we proceed as follows:
 
411
                // -we refine x until having p monotonic in x's interval (to be
 
412
                // sure that p has at most one root on that interval),
 
413
                // -if the gcd has a root on this interval, both roots are
 
414
                // equal (we return 0), otherwise, we refine until having a
 
415
                // result.
 
416
 
 
417
                // How to assure that p is monotonic in an interval: when its
 
418
                // derivative is never zero in that interval.
 
419
                Polynomial_1 dsfpp=Deriv()(sfpp);
 
420
                CGAL::Gmpfr xl(x.get_left());
 
421
                CGAL::Gmpfr xr(x.get_right());
 
422
                while(eval_interv(dsfpp,xl,xr).is_same(unknown)){
 
423
                        Refiner()(x.get_pol(),
 
424
                                  xl,
 
425
                                  xr,
 
426
                                  2*CGAL::max(xl.get_precision(),
 
427
                                              xr.get_precision()));
 
428
                }
 
429
                x.set_left(xl);
 
430
                x.set_right(xr);
 
431
 
 
432
                // How to know that the gcd has a root: evaluate endpoints.
 
433
                CGAL::Sign sleft,sright;
 
434
                Signat sign_at_gcd(gcd);
 
435
                return((sleft=sign_at_gcd(x.get_left()))==CGAL::ZERO||
 
436
                       (sright=sign_at_gcd(x.get_right()))==CGAL::ZERO||
 
437
                       (sleft!=sright));
 
438
        }
 
439
}; // class Is_zero_at_1
 
440
 
 
441
// TODO: it says in the manual that this should return a size_type, but test
 
442
// programs assume that this is equal to int
 
443
template <class Polynomial_,class Isolator_>
295
444
struct Number_of_solutions_1:
296
 
public std::unary_function<PolynomialType,int>{
297
 
    typedef PolynomialType      P;
298
 
    typedef CGAL::to_rs_poly<P> convert;
299
 
 
300
 
    int operator()(const P &p)const{
301
 
        int nr;
302
 
        mpfi_ptr *x;
303
 
        CGAL::RS_polynomial_1 rspoly=convert()(p);
304
 
        x=(mpfi_ptr*)malloc(rspoly.get_degree()*sizeof(mpfi_ptr));
305
 
        nr=solve_1(x,rspoly);
306
 
        CGAL_assertion_msg(nr>=0,"error in resolution");
307
 
        free(x);
308
 
        return nr;
309
 
    }
310
 
};  // Number_of_solutions_1
311
 
 
312
 
template <class PolynomialType,class GcdPolicy>
313
 
struct Sign_at_1:
314
 
public std::binary_function<PolynomialType,Algebraic,CGAL::Sign>{
315
 
    typedef PolynomialType      P;
316
 
    typedef GcdPolicy           Gcd;
317
 
    typedef CGAL::to_rs_poly<P> convert;
318
 
 
319
 
    CGAL::Sign operator()(const P &p,const Algebraic &a)const{
320
 
        return RS3::sign_1(convert()(p),a);
321
 
    }
322
 
};  // Sign_at_1
323
 
 
324
 
template <class PolynomialType,class GcdPolicy>
325
 
struct Is_zero_at_1:
326
 
public std::binary_function<PolynomialType,Algebraic,bool>{
327
 
    typedef PolynomialType      P;
328
 
    typedef GcdPolicy           Gcd;
329
 
    typedef Sign_at_1<P,Gcd>    Sign_at;
330
 
 
331
 
    bool operator()(const P &p,const Algebraic &a)const{
332
 
        return (Sign_at()(p,a)==CGAL::ZERO);
333
 
    }
334
 
};  // Is_zero_at_1
335
 
 
336
 
template <class GcdPolicy>
 
445
public std::unary_function<Polynomial_,int>{
 
446
        typedef Polynomial_                                     Polynomial_1;
 
447
        typedef Isolator_                                       Isolator;
 
448
        size_t operator()(const Polynomial_1 &p)const{
 
449
                // TODO: make sure that p is square free (precondition)
 
450
                Isolator isol(p);
 
451
                return isol.number_of_real_roots();
 
452
        }
 
453
}; // struct Number_of_solutions_1
 
454
 
 
455
// This functor not only compares two algebraic numbers. In case they are
 
456
// different, it refines them until they do not overlap.
 
457
template <class Algebraic_,
 
458
          class Bound_,
 
459
          class Comparator_>
337
460
struct Compare_1:
338
 
    public std::binary_function<Algebraic,Algebraic,CGAL::Comparison_result>{
339
 
  typedef GcdPolicy                     Gcd;
340
 
  typedef CGAL::Comparison_result       Comparison_result;
341
 
  typedef CGAL::Gmpz                    Gmpz;
342
 
  typedef CGAL::Gmpq                    Gmpq;
343
 
 
344
 
  Comparison_result operator()(const Algebraic &r1,const Algebraic &r2)const{
345
 
    return CGAL::RS_COMPARE::compare_1<Gcd>(r1,r2);
346
 
  }
347
 
 
348
 
  Comparison_result operator()(const int &r1,  const Algebraic &r2)const{
349
 
    return this->operator()(Algebraic(r1),r2);}
350
 
  Comparison_result operator()(const Bound &r1,const Algebraic &r2)const{
351
 
    return this->operator()(Algebraic(r1),r2);}
352
 
  Comparison_result operator()(const Gmpz &r1, const Algebraic &r2)const{
353
 
    return this->operator()(Algebraic(r1),r2);}
354
 
  Comparison_result operator()(const Gmpq &r1, const Algebraic &r2)const{
355
 
    return this->operator()(Algebraic(r1),r2);}
356
 
  Comparison_result operator()(const Algebraic &r1,const int   &r2)const{
357
 
    return this->operator()(r1,Algebraic(r2));}
358
 
  Comparison_result operator()(const Algebraic &r1,const Bound &r2)const{
359
 
    return this->operator()(r1,Algebraic(r2));}
360
 
  Comparison_result operator()(const Algebraic &r1,const Gmpz  &r2)const{
361
 
    return this->operator()(r1,Algebraic(r2));}
362
 
  Comparison_result operator()(const Algebraic &r1,const Gmpq  &r2)const{
363
 
    return this->operator()(r1,Algebraic(r2));}
364
 
};  // Compare_1
365
 
 
366
 
template <class PolynomialType,class GcdFunction>
 
461
public std::binary_function<Algebraic_,Algebraic_,CGAL::Comparison_result>{
 
462
        typedef Algebraic_                                      Algebraic;
 
463
        typedef Bound_                                          Bound;
 
464
        typedef Comparator_                                     Comparator;
 
465
 
 
466
        CGAL::Comparison_result operator()(Algebraic a,Algebraic b)const{
 
467
                Bound al=a.get_left();
 
468
                Bound ar=a.get_right();
 
469
                Bound bl=b.get_left();
 
470
                Bound br=b.get_right();
 
471
                CGAL::Comparison_result c=Comparator()(a.get_pol(),al,ar,
 
472
                                                       b.get_pol(),bl,br);
 
473
                a.set_left(al);
 
474
                a.set_right(ar);
 
475
                b.set_left(bl);
 
476
                b.set_right(br);
 
477
                return c;
 
478
        }
 
479
 
 
480
        CGAL::Comparison_result operator()(Algebraic a,const Bound &b)const{
 
481
                Bound al=a.get_left();
 
482
                Bound ar=a.get_right();
 
483
                Algebraic balg(b);
 
484
                CGAL::Comparison_result c=Comparator()(a.get_pol(),al,ar,
 
485
                                                       balg.get_pol(),b,b);
 
486
                a.set_left(al);
 
487
                a.set_right(ar);
 
488
                return c;
 
489
        }
 
490
 
 
491
        template <class T>
 
492
        CGAL::Comparison_result operator()(Algebraic a,const T &b)const{
 
493
                Bound al=a.get_left();
 
494
                Bound ar=a.get_right();
 
495
                Algebraic balg(b);
 
496
                CGAL::Comparison_result c=Comparator()(a.get_pol(),
 
497
                                                       al,
 
498
                                                       ar,
 
499
                                                       balg.get_pol(),
 
500
                                                       balg.get_left(),
 
501
                                                       balg.get_right());
 
502
                a.set_left(al);
 
503
                a.set_right(ar);
 
504
                return c;
 
505
        }
 
506
 
 
507
}; // class Compare_1
 
508
 
 
509
template <class Algebraic_,
 
510
          class Bound_,
 
511
          class Comparator_>
 
512
struct Bound_between_1:
 
513
public std::binary_function<Algebraic_,Algebraic_,Bound_>{
 
514
        typedef Algebraic_                                      Algebraic;
 
515
        typedef Bound_                                          Bound;
 
516
        typedef Comparator_                                     Comparator;
 
517
 
 
518
        Bound operator()(Algebraic a,Algebraic b)const{
 
519
                typedef Compare_1<Algebraic,Bound,Comparator>   Compare;
 
520
                typename Bound::Precision_type prec;
 
521
                switch(Compare()(a,b)){
 
522
                        case CGAL::LARGER:
 
523
                                CGAL_assertion(b.get_right()<a.get_left());
 
524
                                prec=CGAL::max(b.get_right().get_precision(),
 
525
                                               a.get_left().get_precision());
 
526
                                        return Bound::add(b.get_right(),
 
527
                                                          a.get_left(),
 
528
                                                          1+prec)/2;
 
529
                                break;
 
530
                        case CGAL::SMALLER:
 
531
                                CGAL_assertion(a.get_right()<b.get_left());
 
532
                                prec=CGAL::max(a.get_right().get_precision(),
 
533
                                               b.get_left().get_precision());
 
534
                                        return Bound::add(a.get_right(),
 
535
                                                          b.get_left(),
 
536
                                                          1+prec)/2;
 
537
                                break;
 
538
                        default:
 
539
                                CGAL_error_msg(
 
540
                                        "bound between two equal numbers");
 
541
                }
 
542
        }
 
543
}; // struct Bound_between_1
 
544
 
 
545
template <class Polynomial_,
 
546
          class Bound_,
 
547
          class Algebraic_,
 
548
          class Isolator_,
 
549
          class Comparator_,
 
550
          class Signat_,
 
551
          class Ptraits_>
367
552
struct Isolate_1:
368
 
public std::binary_function<Algebraic,PolynomialType,std::pair<Bound,Bound> >{
369
 
    typedef PolynomialType              Poly;
370
 
    typedef GcdFunction                 Gcd;
371
 
    typedef CGAL::to_rs_poly<Poly>      convert;
372
 
    typedef Solve_1<Poly,Gcd>           Solve;
373
 
    typedef Compare_1<Gcd>              Compare;
374
 
 
375
 
    std::pair<Bound,Bound> operator()(const Algebraic &a,const Poly &p)const{
376
 
        std::vector<Algebraic> roots;
377
 
        std::back_insert_iterator<std::vector<Algebraic> > rootsit(roots);
378
 
        Solve()(p,true,rootsit);
379
 
        for(std::vector<Algebraic>::size_type i=0;i<roots.size();++i)
380
 
                Compare()(a,roots[i]);
381
 
        return std::make_pair(Bound(a.left()),Bound(a.right()));
382
 
    }
383
 
};  // Isolate_1
384
 
 
385
 
template <class GcdPolicy>
386
 
struct Bound_between_1:
387
 
    public std::binary_function<Algebraic,Algebraic,Bound>{
388
 
        typedef GcdPolicy       Gcd;
389
 
        Bound operator()(const Algebraic &x1,const Algebraic &x2)const{
390
 
            double l,r,m;
391
 
            switch(CGAL::RS_COMPARE::compare_1<Gcd>(x1,x2)){
392
 
                case CGAL::LARGER:
393
 
                    CGAL_assertion(x2.sup()<x1.inf());
394
 
                    l=x2.sup().to_double(std::round_toward_infinity);
395
 
                    r=x1.inf().to_double(std::round_toward_neg_infinity);
396
 
                    m=(l+r)/2;
397
 
                    if(l<m&&m<r){
398
 
                        return Bound(m,CGAL_RS_FUNCTORS_DBL_PREC);
399
 
                    }
400
 
                    return Bound::add(x2.sup(),
401
 
                                      x1.inf(),
402
 
                                      (x2.sup().get_precision()>
403
 
                                                x1.inf().get_precision()?
404
 
                                       1+x2.sup().get_precision():
405
 
                                       1+x1.inf().get_precision()))/2;
406
 
                    break;
407
 
                case CGAL::SMALLER:
408
 
                    CGAL_assertion(x1.sup()<x2.inf());
409
 
                    l=x1.sup().to_double(std::round_toward_infinity);
410
 
                    r=x2.inf().to_double(std::round_toward_neg_infinity);
411
 
                    m=(l+r)/2;
412
 
                    if(l<m&&m<r){
413
 
                        return Bound(m,CGAL_RS_FUNCTORS_DBL_PREC);
414
 
                    }
415
 
                    return Bound::add(x1.sup(),
416
 
                                      x2.inf(),
417
 
                                      (x1.sup().get_precision()>
418
 
                                                x2.inf().get_precision()?
419
 
                                       1+x1.sup().get_precision():
420
 
                                       1+x2.inf().get_precision()))/2;
421
 
                    break;
422
 
                default:
423
 
                    CGAL_error_msg("bound between two equal numbers");
424
 
            }
 
553
public std::binary_function<Algebraic_,Polynomial_,std::pair<Bound_,Bound_> >{
 
554
        typedef Polynomial_                                     Polynomial_1;
 
555
        typedef Bound_                                          Bound;
 
556
        typedef Algebraic_                                      Algebraic;
 
557
        typedef Isolator_                                       Isolator;
 
558
        typedef Comparator_                                     Comparator;
 
559
        typedef Signat_                                         Signat;
 
560
        typedef Ptraits_                                        Ptraits;
 
561
 
 
562
        std::pair<Bound,Bound>
 
563
        operator()(const Algebraic &a,const Polynomial_1 &p)const{
 
564
                std::vector<Algebraic> roots;
 
565
                std::back_insert_iterator<std::vector<Algebraic> > rit(roots);
 
566
                typedef Solve_1<Polynomial_1,
 
567
                                Bound,
 
568
                                Algebraic,
 
569
                                Isolator,
 
570
                                Signat,
 
571
                                Ptraits>                        Solve;
 
572
                typedef Compare_1<Algebraic,Bound,Comparator>   Compare;
 
573
                Solve()(p,false,rit);
 
574
                for(typename std::vector<Algebraic>::size_type i=0;
 
575
                    i<roots.size();
 
576
                    ++i){
 
577
                        // we use the comparison functor, that makes both
 
578
                        // intervals disjoint iff the algebraic numbers they
 
579
                        // represent are not equal
 
580
                        Compare()(a,roots[i]);
 
581
                }
 
582
                return std::make_pair(a.left(),a.right());
425
583
        }
426
 
    };  // Bound_between_1
 
584
}; // Isolate_1
427
585
 
 
586
template <class Polynomial_,
 
587
          class Bound_,
 
588
          class Algebraic_,
 
589
          class Refiner_>
428
590
struct Approximate_absolute_1:
429
 
    public std::binary_function<Algebraic,int,std::pair<Bound,Bound> >{
430
 
 
431
 
  std::pair<Bound,Bound>
432
 
  operator()(const Algebraic& x, int prec) const {
433
 
    RS3::refine_1(x,CGAL::abs(prec));
434
 
 
435
 
    CGAL_postcondition_code(
436
 
     CGAL::Gmpfr::Precision_type
437
 
        subprec=1+
438
 
            std::max<CGAL::Gmpfr::Precision_type>(x.inf().get_precision(),
439
 
                                                  x.sup().get_precision());
440
 
        CGAL::Gmpfr width=CGAL::Gmpfr::sub(x.sup(),x.inf(),subprec);
441
 
    )
442
 
    CGAL_postcondition_code(if(prec>0))
443
 
    CGAL_postcondition(width*CGAL::ipower(Bound(2),prec)<=Bound(1));
444
 
    CGAL_postcondition_code(else)
445
 
    CGAL_postcondition(width<=CGAL::ipower(Bound(2),-prec));
446
 
 
447
 
    return std::make_pair(x.inf(),x.sup());
448
 
  }
449
 
};
450
 
 
451
 
struct Approximate_relative_1
452
 
  :public std::binary_function<Algebraic,int,std::pair<Bound,Bound> >{
453
 
 
454
 
  std::pair<Bound,Bound> operator()(const Algebraic &x, int prec) const {
455
 
    if(CGAL::is_zero(x))
456
 
      return std::make_pair(Bound(0),Bound(0));
457
 
 
458
 
    Bound error=CGAL::ipower(Bound(2),CGAL::abs(prec));
459
 
    Bound max_b=(CGAL::max)(CGAL::abs(x.sup()),CGAL::abs(x.inf()));
460
 
    while(prec>0?
461
 
          (x.sup()-x.inf())*error>max_b:
462
 
          (x.sup()-x.inf())>error*max_b){
463
 
      RS3::refine_1(x,mpfi_get_prec(x.mpfi())+CGAL::abs(prec));
464
 
      max_b=(CGAL::max)(CGAL::abs(x.sup()),CGAL::abs(x.inf()));
465
 
    }
466
 
 
467
 
    CGAL_postcondition_code(if(prec>0))
468
 
    CGAL_postcondition((x.sup()-x.inf())*CGAL::ipower(Bound(2),prec)<=max_b);
469
 
    CGAL_postcondition_code(else)
470
 
    CGAL_postcondition((x.sup()-x.inf())<=CGAL::ipower(Bound(2),-prec)*max_b);
471
 
 
472
 
    return std::make_pair(x.inf(),x.sup());
473
 
  }
474
 
};
475
 
 
476
 
} // namespace RSFunctors
477
 
 
478
 
#undef CGAL_RS_FUNCTORS_DBL_PREC
479
 
 
480
 
#endif  // CGAL_RS_FUNCTORS_H
 
591
public std::binary_function<Algebraic_,int,std::pair<Bound_,Bound_> >{
 
592
        typedef Polynomial_                                     Polynomial_1;
 
593
        typedef Bound_                                          Bound;
 
594
        typedef Algebraic_                                      Algebraic;
 
595
        typedef Refiner_                                        Refiner;
 
596
 
 
597
        // This implementation assumes that Bound is Gmpfr.
 
598
        // TODO: make generic.
 
599
        std::pair<Bound,Bound> operator()(const Algebraic &x,const int &a)const{
 
600
                Bound xl(x.get_left()),xr(x.get_right());
 
601
                // refsteps=log2(xl-xr)
 
602
                mpfr_t temp;
 
603
                mpfr_init(temp);
 
604
                mpfr_sub(temp,xr.fr(),xl.fr(),GMP_RNDU);
 
605
                mpfr_log2(temp,temp,GMP_RNDU);
 
606
                long refsteps=mpfr_get_si(temp,GMP_RNDU);
 
607
                mpfr_clear(temp);
 
608
                Refiner()(x.get_pol(),xl,xr,CGAL::abs(refsteps+a));
 
609
                x.set_left(xl);
 
610
                x.set_right(xr);
 
611
                CGAL_assertion(a>0?
 
612
                               (xr-xl)*CGAL::ipower(Bound(2),a)<=Bound(1):
 
613
                               (xr-xl)<=CGAL::ipower(Bound(2),-a));
 
614
                return std::make_pair(xl,xr);
 
615
        }
 
616
}; // struct Approximate_absolute_1
 
617
 
 
618
template <class Polynomial_,
 
619
          class Bound_,
 
620
          class Algebraic_,
 
621
          class Refiner_>
 
622
struct Approximate_relative_1:
 
623
public std::binary_function<Algebraic_,int,std::pair<Bound_,Bound_> >{
 
624
        typedef Polynomial_                                     Polynomial_1;
 
625
        typedef Bound_                                          Bound;
 
626
        typedef Algebraic_                                      Algebraic;
 
627
        typedef Refiner_                                        Refiner;
 
628
 
 
629
        std::pair<Bound,Bound> operator()(const Algebraic &x,const int &a)const{
 
630
                if(CGAL::is_zero(x))
 
631
                        return std::make_pair(Bound(0),Bound(0));
 
632
                Bound error=CGAL::ipower(Bound(2),CGAL::abs(a));
 
633
                Bound xl(x.get_left()),xr(x.get_right());
 
634
                Bound max_b=(CGAL::max)(CGAL::abs(xr),CGAL::abs(xl));
 
635
                while(a>0?(xr-xl)*error>max_b:(xr-xl)>error*max_b){
 
636
                        Refiner()(x.get_pol(),
 
637
                                  xl,
 
638
                                  xr,
 
639
                                  std::max<unsigned>(
 
640
                                        CGAL::abs(a),
 
641
                                        CGAL::max(xl.get_precision(),
 
642
                                                  xr.get_precision())));
 
643
                        max_b=(CGAL::max)(CGAL::abs(xr),CGAL::abs(xl));
 
644
                }
 
645
                x.set_left(xl);
 
646
                x.set_right(xr);
 
647
                CGAL_assertion(
 
648
                                a>0?
 
649
                                (xr-xl)*CGAL::ipower(Bound(2),a)<=max_b:
 
650
                                (xr-xl)<=CGAL::ipower(Bound(2),-a)*max_b);
 
651
                return std::make_pair(xl,xr);
 
652
        }
 
653
}; // struct Approximate_relative_1
 
654
 
 
655
} // namespace RS_AK1
 
656
} // namespace CGAL
 
657
 
 
658
#endif // CGAL_RS_FUNCTORS_1_H