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