← Back to branch summary

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

« back to all changes in this revision

Viewing changes to include/CGAL/RS/bisection_refiner_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/bisection_refiner_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
// This file contains the simplest refiner, that bisects the interval a given
 
22
// number of times.
 
23
 
 
24
#ifndef CGAL_RS_BISECTION_REFINER_1_H
 
25
#define CGAL_RS_BISECTION_REFINER_1_H
 
26
 
 
27
#include <CGAL/Polynomial_traits_d.h>
 
28
#include "signat_1.h"
 
29
#include "Gmpfr_make_unique.h"
 
30
 
 
31
namespace CGAL{
 
32
 
 
33
template <class Polynomial_,class Bound_>
 
34
struct Bisection_refiner_1{
 
35
        typedef CGAL::RS_AK1::Signat_1<Polynomial_,Bound_>           Signat;
 
36
        void operator()(const Polynomial_&,Bound_&,Bound_&,int);
 
37
}; // class Bisection_refiner_1
 
38
 
 
39
// TODO: Write in a generic way, if possible (see next function).
 
40
template <class Polynomial_,class Bound_>
 
41
void
 
42
Bisection_refiner_1<Polynomial_,Bound_>::
 
43
operator()(const Polynomial_&,Bound_&,Bound_&,int){
 
44
        CGAL_error_msg("bisection refiner not implemented for these types");
 
45
        return;
 
46
}
 
47
 
 
48
// This works with any type of polynomial, but only for Gmpfr bounds.
 
49
// TODO: Beyond writing generically, optimize this function. This would
 
50
// remove the need for the next function, which essentially the same.
 
51
template<>
 
52
void
 
53
Bisection_refiner_1<Polynomial<Gmpz>,Gmpfr>::
 
54
operator()(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec){
 
55
        typedef Polynomial<Gmpz>                        Polynomial;
 
56
        typedef Polynomial_traits_d<Polynomial>         Ptraits;
 
57
        typedef Ptraits::Make_square_free               Sfpart;
 
58
        typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
 
59
                                                        Signat;
 
60
        CGAL_precondition(left<=right);
 
61
        // TODO: add precondition to check whether the interval is a point
 
62
        // or the evaluations on its endpoints have different signs
 
63
        //std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
 
64
 
 
65
        CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
 
66
        CGAL_RS_GMPFR_MAKE_UNIQUE(right,temp_right);
 
67
 
 
68
        Polynomial sfpp=Sfpart()(pol);
 
69
        Signat signof(sfpp);
 
70
        CGAL::Sign sl,sc;
 
71
        mp_prec_t pl,pc;
 
72
        mpfr_t center;
 
73
 
 
74
        sl=signof(left);
 
75
        if(sl==ZERO)
 
76
                return;
 
77
        pl=left.get_precision();
 
78
        pc=right.get_precision();
 
79
        pc=(pl>pc?pl:pc)+(mp_prec_t)prec;
 
80
        mpfr_init2(center,pc);
 
81
        CGAL_assertion_code(int round=)
 
82
        mpfr_prec_round(left.fr(),pc,GMP_RNDN);
 
83
        CGAL_assertion(!round);
 
84
        CGAL_assertion_code(round=)
 
85
        mpfr_prec_round(right.fr(),pc,GMP_RNDN);
 
86
        CGAL_assertion(!round);
 
87
        for(int i=0;i<prec;++i){
 
88
                CGAL_assertion_code(round=)
 
89
                mpfr_add(center,left.fr(),right.fr(),GMP_RNDN);
 
90
                CGAL_assertion(!round);
 
91
                CGAL_assertion_code(round=)
 
92
                mpfr_div_2ui(center,center,1,GMP_RNDN);
 
93
                CGAL_assertion(!round);
 
94
                sc=signof(Gmpfr(center));
 
95
                if(sc==ZERO){   // we have a root
 
96
                        CGAL_assertion_code(round=)
 
97
                        mpfr_set(left.fr(),center,GMP_RNDN);
 
98
                        CGAL_assertion(!round);
 
99
                        mpfr_swap(right.fr(),center);
 
100
                        break;
 
101
                }
 
102
                if(sc==sl)
 
103
                        mpfr_swap(left.fr(),center);
 
104
                else
 
105
                        mpfr_swap(right.fr(),center);
 
106
        }
 
107
        mpfr_clear(center);
 
108
        CGAL_postcondition(left<=right);
 
109
        //std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
 
110
        return;
 
111
}
 
112
 
 
113
template<>
 
114
void
 
115
Bisection_refiner_1<Polynomial<Gmpq>,Gmpfr>::
 
116
operator()(const Polynomial<Gmpq> &pol,Gmpfr &left,Gmpfr &right,int prec){
 
117
        typedef Polynomial<Gmpq>                        Polynomial;
 
118
        typedef Polynomial_traits_d<Polynomial>         Ptraits;
 
119
        typedef Ptraits::Make_square_free               Sfpart;
 
120
        typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
 
121
                                                        Signat;
 
122
        CGAL_precondition(left<=right);
 
123
        // TODO: add precondition to check whether the interval is a point
 
124
        // or the evaluations on its endpoints have different signs
 
125
        //std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
 
126
 
 
127
        CGAL_RS_GMPFR_MAKE_UNIQUE(left,temp_left);
 
128
        CGAL_RS_GMPFR_MAKE_UNIQUE(right,temp_right);
 
129
 
 
130
        Polynomial sfpp=Sfpart()(pol);
 
131
        Signat signof(sfpp);
 
132
        CGAL::Sign sl,sc;
 
133
        mp_prec_t pl,pc;
 
134
        mpfr_t center;
 
135
 
 
136
        sl=signof(left);
 
137
        if(sl==ZERO)
 
138
                return;
 
139
        pl=left.get_precision();
 
140
        pc=right.get_precision();
 
141
        pc=(pl>pc?pl:pc)+(mp_prec_t)prec;
 
142
        mpfr_init2(center,pc);
 
143
        CGAL_assertion_code(int round=)
 
144
        mpfr_prec_round(left.fr(),pc,GMP_RNDN);
 
145
        CGAL_assertion(!round);
 
146
        CGAL_assertion_code(round=)
 
147
        mpfr_prec_round(right.fr(),pc,GMP_RNDN);
 
148
        CGAL_assertion(!round);
 
149
        for(int i=0;i<prec;++i){
 
150
                CGAL_assertion_code(round=)
 
151
                mpfr_add(center,left.fr(),right.fr(),GMP_RNDN);
 
152
                CGAL_assertion(!round);
 
153
                CGAL_assertion_code(round=)
 
154
                mpfr_div_2ui(center,center,1,GMP_RNDN);
 
155
                CGAL_assertion(!round);
 
156
                sc=signof(Gmpfr(center));
 
157
                if(sc==ZERO){   // we have a root
 
158
                        CGAL_assertion_code(round=)
 
159
                        mpfr_set(left.fr(),center,GMP_RNDN);
 
160
                        CGAL_assertion(!round);
 
161
                        mpfr_swap(right.fr(),center);
 
162
                        break;
 
163
                }
 
164
                if(sc==sl)
 
165
                        mpfr_swap(left.fr(),center);
 
166
                else
 
167
                        mpfr_swap(right.fr(),center);
 
168
        }
 
169
        mpfr_clear(center);
 
170
        CGAL_postcondition(left<=right);
 
171
        //std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
 
172
        return;
 
173
}
 
174
 
 
175
} // namespace CGAL
 
176
 
 
177
#endif // CGAL_RS_BISECTION_REFINER_1_H