1 // Copyright (c) 2006-2009 Max-Planck-Institute Saarbruecken (Germany).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org)
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/bound_between_1.h $
7 // $Id: bound_between_1.h 26355e2 2020-06-25T12:31:21+02:00 Mael Rouxel-Labbé
8 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 //
12 // Author(s) : Michael Kerber <mkerber@mpi-inf.mpg.de>
13 //
14 // ============================================================================
15
16
17 #ifndef CGAL_BOUND_BETWEEN_1_H
18 #define CGAL_BOUND_BETWEEN_1_H 1
19
20 #include <iterator>
21
22 #include <CGAL/basic.h>
23 #include <CGAL/Algebraic_kernel_d/enums.h>
24 #include <CGAL/Arithmetic_kernel.h>
25 #include <CGAL/Algebraic_kernel_d/Float_traits.h>
26 #include <CGAL/convert_to_bfi.h>
27 #include <CGAL/Algebraic_kernel_d/Real_embeddable_extension.h>
28 #include <CGAL/Polynomial_traits_d.h>
29 #include <CGAL/Algebraic_kernel_d/Bitstream_coefficient_kernel.h>
30 #include <CGAL/Bigfloat_interval_traits.h>
31 #include <boost/numeric/interval.hpp>
32
33 #include <CGAL/Coercion_traits.h>
34
35 namespace CGAL {
36
37 namespace internal {
38
39 /*! \brief tries to find a SIMPLE rational q with a<q<b.
40 *
41 * In this context, simple means that the denominator of <tt>q</tt>
42 * is a power of two, and is not too big. There is no guarantee to find
43 * the rational value between <tt>a</tt> and <tt>b</tt> of minimal
44 * bit size.
45 */
46 template<typename Algebraic_real>
47 typename Algebraic_real::Rational
simple_bound_between(const Algebraic_real & a,const Algebraic_real & b)48 simple_bound_between(const Algebraic_real& a,
49 const Algebraic_real&b) {
50
51 //srb.start();
52 typedef typename Algebraic_real::Rational Rational;
53 typename CGAL::Fraction_traits<Rational>::Compose compose;
54 typedef typename
55 CGAL::Get_arithmetic_kernel<Rational>::Arithmetic_kernel AK;
56 typedef typename AK::Bigfloat_interval Bigfloat_interval;
57 typedef typename CGAL::Bigfloat_interval_traits<Bigfloat_interval>
58 ::Bound Bigfloat;
59 typedef typename AK::Integer Integer;
60
61 long old_prec = CGAL::get_precision(Bigfloat_interval());
62
63 CGAL_assertion(a!=b);
64 if(a>b) {
65 return simple_bound_between(b,a);
66 }
67
68 //std::cout << "Intermediate1: " << CGAL::to_double(a) << " " << CGAL::to_double(b) << std::endl;
69 /*
70 * First, refine a and b until their isolating intervals are disjoint
71 * Therefore, the bigger interval is refined in each substep
72 */
73 //srb_a.start();
74 if(a.high() >= b.low()) {
75 Rational size_a=a.high()-a.low(),
76 size_b=b.high() - b.low();
77 while(a.high() >= b.low()) {
78 if(size_a < size_b) {
79 b.refine();
80 size_b=b.high() - b.low();
81 } else {
82 a.refine();
83 size_a=a.high()-a.low();
84 }
85 }
86 }
87 //srb_a.stop();
88
89 //srb_b.start();
90 Bigfloat x=CGAL::upper(CGAL::convert_to_bfi(a.high()));
91 Bigfloat y=CGAL::lower(CGAL::convert_to_bfi(b.low()));
92
93 if(x>=y) {
94 Rational size_a=a.high() - a.low(),
95 size_b=b.high() - b.low(),
96 size_max = size_a>size_b ? size_a : size_b,
97 size_int = b.low()-a.high();
98 while(x>=y) {
99 //std::cout << "x and y: " << x << " and " << y << std::endl;
100 //std::cout << "sizes: " << CGAL::to_double(size_int) << " " << CGAL::to_double(size_max) << std::endl;
101 if(size_int>size_max) {
102 CGAL::set_precision(Bigfloat_interval(),
103 2*CGAL::get_precision(Bigfloat_interval()));
104 x=CGAL::upper(CGAL::convert_to_bfi(a.high()));
105 y=CGAL::lower(CGAL::convert_to_bfi(b.low()));
106 } else {
107 if(size_a < size_b) {
108 b.refine();
109 size_b=b.high() - b.low();
110 y=CGAL::lower(CGAL::convert_to_bfi(b.low()));
111 } else {
112 a.refine();
113 size_a=a.high()-a.low();
114 x=CGAL::upper(CGAL::convert_to_bfi(a.high()));
115 }
116 size_max = size_a>size_b ? size_a : size_b;
117 size_int = b.low()-a.high();
118 }
119 }
120 }
121 CGAL_assertion(x<y);
122
123 //srb_b.stop();
124 //std::cout << "Intermediate2: " << x << " " << y << std::endl;
125 typename CGAL::internal::Float_traits<Bigfloat>::Get_mantissa mantissa;
126 typename CGAL::internal::Float_traits<Bigfloat>::Get_exponent exponent;
127
128 // std::cout << CGAL::to_double(x) << " < " << CGAL::to_double(y) << std::endl;
129
130 Integer x_m = mantissa(x), y_m=mantissa(y);
131 long x_e = exponent(x), y_e = exponent(y);
132 //std::cout << "Floats1: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
133
134
135 if (((x_m > 0) && (y_m < 0)) || ((x_m < 0) && (y_m > 0))) {
136 //srb.stop();
137 return Rational(0);
138 }
139 bool negative=false;
140 if(x_m<=0 && y_m <=0) {
141 x_m=-x_m;
142 y_m=-y_m;
143 std::swap(x_m,y_m);
144 std::swap(x_e,y_e);
145 negative=true;
146 }
147 // Now, we have that (x_m,x_e) represents a number smaller than (y_m,y_e)
148 //srb_c.start();
149 //std::cout << "Floats2: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
150
151 // As long as the mantissa is even, simplify
152 while(x_m != 0 && (x_m & 1)==0 ) {
153 x_m=x_m >> 1;
154 x_e++;
155 }
156 while(y_m != 0 && (y_m & 1)==0 ) {
157 y_m=y_m >> 1;
158 y_e++;
159 }
160 //srb_c.stop();
161 //std::cout << "Floats3: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
162
163 // Bring both numbers to a common exponent
164 //srb_d.start();
165 long min_e = x_e < y_e ? x_e : y_e;
166 while(x_e > min_e) {
167 x_m=x_m << 1;
168 x_e--;
169 }
170 while(y_e > min_e) {
171 y_m=y_m << 1;
172 y_e--;
173 }
174 //srb_d.stop();
175 CGAL_assertion(y_e==x_e && x_e==min_e);
176 CGAL_assertion(x_m < y_m);
177 //std::cout << "Floats4: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
178
179 // Avoid mantissas to have difference one
180 if(y_m-x_m==Integer(1)) {
181 x_m=x_m << 1;
182 y_m=y_m << 1;
183 x_e--;
184 y_e--;
185 min_e--;
186 }
187 //std::cout << "Floats5: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
188 Integer final_mantissa(0);
189 //srb_e.start();
190 long x_log = x_m==Integer(0) ? -1 : CGAL::internal::floor_log2_abs(x_m),
191 y_log = y_m==Integer(0) ? -1 : CGAL::internal::floor_log2_abs(y_m),
192 old_log = y_log;
193 //std::cout << x_log << " < " << y_log << std::endl;
194 while(x_log==y_log) {
195 //std::cout << "here" << std::endl;
196 while(old_log > y_log) {
197 final_mantissa = final_mantissa << 1;
198 old_log--;
199 }
200 CGAL_assertion((x_m & ((Integer(1) << x_log) - 1)) == x_m - CGAL::ipower(Integer(2),x_log));
201 x_m = x_m & ((Integer(1) << x_log) - 1); // x_m - CGAL::ipower(Integer(2),x_log);
202 y_m = y_m & ((Integer(1) << y_log) - 1); // y_m - CGAL::ipower(Integer(2),y_log);
203
204 final_mantissa++;
205 old_log=y_log;
206 x_log = x_m==0 ? -1 : CGAL::internal::floor_log2_abs(x_m);
207 y_log = y_m==0 ? -1 : CGAL::internal::floor_log2_abs(y_m);
208 }
209 //srb_e.stop();
210 // Now, x_log != y_log, in fact, y_log is greater
211 CGAL_assertion(x_log<y_log);
212 //srb_f.start();
213 while(old_log > y_log) {
214 final_mantissa = final_mantissa << 1;
215 old_log--;
216 }
217 if((y_m & ((Integer(1) << y_log) - 1 ))==0) { // y_m - CGAL::ipower(Integer(2),y_log)==0) {
218 // Now, the constructed value would be equal to
219 while(y_log!=0 && x_log==y_log-1) {
220 final_mantissa = final_mantissa << 1;
221 final_mantissa++;
222 y_log--;
223 x_m = x_m==0 ? 0 : x_m & ((Integer(1) << x_log) - 1); //x_m - CGAL::ipower(Integer(2),x_log);
224 x_log = x_m==0 ? -1 : CGAL::internal::floor_log2_abs(x_m);
225 }
226 final_mantissa = final_mantissa << 1;
227 final_mantissa++;
228 y_log--;
229 } else {
230 final_mantissa++;
231 }
232 //srb_f.stop();
233 min_e += y_log;
234 Rational rat_between;
235 //std::cout << "Min_e: " << min_e << std::endl;
236 if(min_e > 0) {
237 rat_between = compose(final_mantissa << min_e,
238 Integer(1));
239 } else {
240 rat_between = compose(final_mantissa, Integer(1) << -min_e);
241 }
242 if(negative) {
243 rat_between = -rat_between;
244 }
245 //std::cout << "Result: " << a.high() << " " << rat_between << " " << b.low() << std::endl;
246 CGAL_assertion(a.high() < rat_between);
247 CGAL_assertion(b.low() > rat_between);
248 CGAL::set_precision(Bigfloat_interval(),old_prec);
249 //srb.stop();
250 return rat_between;
251 }
252
253
254 } // namespace internal
255
256
257 } //namespace CGAL
258
259 #endif // CGAL_BOUND_BETWEEN_1_H
260