← Back to branch summary

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

« back to all changes in this revision

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