1 /*++
2 Copyright (c) 2012 Microsoft Corporation
3
4 Module Name:
5
6 interval.h
7
8 Abstract:
9
10 Goodies/Templates for interval arithmetic
11
12 Author:
13
14 Leonardo de Moura (leonardo) 2012-07-19.
15
16 Revision History:
17
18 --*/
19 #pragma once
20
21 #include "util/mpq.h"
22 #include "util/ext_numeral.h"
23 #include "util/rlimit.h"
24
25 /**
26 \brief Default configuration for interval manager.
27 It is used for documenting the required interface.
28 */
29 class im_default_config {
30 unsynch_mpq_manager & m_manager;
31 public:
32 typedef unsynch_mpq_manager numeral_manager;
33 typedef mpq numeral;
34
35 // Every configuration object must provide an interval type.
36 // The actual fields are irrelevant, the interval manager
37 // accesses interval data using the following API.
38 struct interval {
39 numeral m_lower;
40 numeral m_upper;
41 unsigned m_lower_open:1;
42 unsigned m_upper_open:1;
43 unsigned m_lower_inf:1;
44 unsigned m_upper_inf:1;
intervalinterval45 interval():
46 m_lower_open(false),
47 m_upper_open(false),
48 m_lower_inf(true),
49 m_upper_inf(true) {}
50 };
51
52 // Should be NOOPs for precise numeral types.
53 // For imprecise types (e.g., floats) it should set the rounding mode.
round_to_minus_inf()54 void round_to_minus_inf() {}
round_to_plus_inf()55 void round_to_plus_inf() {}
set_rounding(bool to_plus_inf)56 void set_rounding(bool to_plus_inf) {}
57
58 // Getters
lower(interval const & a)59 numeral const & lower(interval const & a) const { return a.m_lower; }
upper(interval const & a)60 numeral const & upper(interval const & a) const { return a.m_upper; }
lower(interval & a)61 numeral & lower(interval & a) { return a.m_lower; }
upper(interval & a)62 numeral & upper(interval & a) { return a.m_upper; }
lower_is_open(interval const & a)63 bool lower_is_open(interval const & a) const { return a.m_lower_open; }
upper_is_open(interval const & a)64 bool upper_is_open(interval const & a) const { return a.m_upper_open; }
lower_is_inf(interval const & a)65 bool lower_is_inf(interval const & a) const { return a.m_lower_inf; }
upper_is_inf(interval const & a)66 bool upper_is_inf(interval const & a) const { return a.m_upper_inf; }
67
68 // Setters
set_lower(interval & a,numeral const & n)69 void set_lower(interval & a, numeral const & n) { m_manager.set(a.m_lower, n); }
set_upper(interval & a,numeral const & n)70 void set_upper(interval & a, numeral const & n) { m_manager.set(a.m_upper, n); }
set_lower_is_open(interval & a,bool v)71 void set_lower_is_open(interval & a, bool v) { a.m_lower_open = v; }
set_upper_is_open(interval & a,bool v)72 void set_upper_is_open(interval & a, bool v) { a.m_upper_open = v; }
set_lower_is_inf(interval & a,bool v)73 void set_lower_is_inf(interval & a, bool v) { a.m_lower_inf = v; }
set_upper_is_inf(interval & a,bool v)74 void set_upper_is_inf(interval & a, bool v) { a.m_upper_inf = v; }
75
76 // Reference to numeral manager
m()77 numeral_manager & m() const { return m_manager; }
78
im_default_config(numeral_manager & m)79 im_default_config(numeral_manager & m):m_manager(m) {}
80 };
81
82 #define DEP_IN_LOWER1 1
83 #define DEP_IN_UPPER1 2
84 #define DEP_IN_LOWER2 4
85 #define DEP_IN_UPPER2 8
86
87 typedef short deps_combine_rule;
dep_in_lower1(deps_combine_rule d)88 inline bool dep_in_lower1(deps_combine_rule d) { return (d & DEP_IN_LOWER1) != 0; }
dep_in_lower2(deps_combine_rule d)89 inline bool dep_in_lower2(deps_combine_rule d) { return (d & DEP_IN_LOWER2) != 0; }
dep_in_upper1(deps_combine_rule d)90 inline bool dep_in_upper1(deps_combine_rule d) { return (d & DEP_IN_UPPER1) != 0; }
dep_in_upper2(deps_combine_rule d)91 inline bool dep_in_upper2(deps_combine_rule d) { return (d & DEP_IN_UPPER2) != 0; }
dep1_to_dep2(deps_combine_rule d)92 inline deps_combine_rule dep1_to_dep2(deps_combine_rule d) {
93 SASSERT(!dep_in_lower2(d) && !dep_in_upper2(d));
94 deps_combine_rule r = d << 2;
95 SASSERT(dep_in_lower1(d) == dep_in_lower2(r));
96 SASSERT(dep_in_upper1(d) == dep_in_upper2(r));
97 SASSERT(!dep_in_lower1(r) && !dep_in_upper1(r));
98 return r;
99 }
100
101 /**
102 \brief Interval dependencies for unary and binary operations on intervals.
103 It contains the dependencies for the output lower and upper bounds
104 for the resultant interval.
105 */
106 struct interval_deps_combine_rule {
107 deps_combine_rule m_lower_combine;
108 deps_combine_rule m_upper_combine;
resetinterval_deps_combine_rule109 void reset() {
110 m_lower_combine = m_upper_combine = 0;
111 }
112 };
113
114 template<typename C>
115 class interval_manager {
116 public:
117 typedef typename C::numeral_manager numeral_manager;
118 typedef typename numeral_manager::numeral numeral;
119 typedef typename C::interval interval;
120 private:
121 reslimit& m_limit;
122 C m_c;
123 numeral m_result_lower;
124 numeral m_result_upper;
125 numeral m_mul_ad;
126 numeral m_mul_bc;
127 numeral m_mul_ac;
128 numeral m_mul_bd;
129 numeral m_one;
130 numeral m_minus_one;
131 numeral m_inv_k;
132
133 unsigned m_pi_n;
134 interval m_pi_div_2;
135 interval m_pi;
136 interval m_3_pi_div_2;
137 interval m_2_pi;
138
139
round_to_minus_inf()140 void round_to_minus_inf() { m_c.round_to_minus_inf(); }
round_to_plus_inf()141 void round_to_plus_inf() { m_c.round_to_plus_inf(); }
set_rounding(bool to_plus_inf)142 void set_rounding(bool to_plus_inf) { m_c.set_rounding(to_plus_inf); }
143
lower_kind(interval const & a)144 ext_numeral_kind lower_kind(interval const & a) const { return m_c.lower_is_inf(a) ? EN_MINUS_INFINITY : EN_NUMERAL; }
upper_kind(interval const & a)145 ext_numeral_kind upper_kind(interval const & a) const { return m_c.upper_is_inf(a) ? EN_PLUS_INFINITY : EN_NUMERAL; }
146
set_lower(interval & a,numeral const & n)147 void set_lower(interval & a, numeral const & n) { m_c.set_lower(a, n); }
set_upper(interval & a,numeral const & n)148 void set_upper(interval & a, numeral const & n) { m_c.set_upper(a, n); }
set_lower_is_open(interval & a,bool v)149 void set_lower_is_open(interval & a, bool v) { m_c.set_lower_is_open(a, v); }
set_upper_is_open(interval & a,bool v)150 void set_upper_is_open(interval & a, bool v) { m_c.set_upper_is_open(a, v); }
set_lower_is_inf(interval & a,bool v)151 void set_lower_is_inf(interval & a, bool v) { m_c.set_lower_is_inf(a, v); }
set_upper_is_inf(interval & a,bool v)152 void set_upper_is_inf(interval & a, bool v) { m_c.set_upper_is_inf(a, v); }
153
154 void nth_root_slow(numeral const & a, unsigned n, numeral const & p, numeral & lo, numeral & hi);
155 void A_div_x_n(numeral const & A, numeral const & x, unsigned n, bool to_plus_inf, numeral & r);
156 void rough_approx_nth_root(numeral const & a, unsigned n, numeral & o);
157 void approx_nth_root(numeral const & a, unsigned n, numeral const & p, numeral & o);
158 void nth_root_pos(numeral const & A, unsigned n, numeral const & p, numeral & lo, numeral & hi);
159 void nth_root(numeral const & a, unsigned n, numeral const & p, numeral & lo, numeral & hi);
160
161 void pi_series(int x, numeral & r, bool to_plus_inf);
162 void fact(unsigned n, numeral & o);
163 void sine_series(numeral const & a, unsigned k, bool upper, numeral & o);
164 void cosine_series(numeral const & a, unsigned k, bool upper, numeral & o);
165 void e_series(unsigned k, bool upper, numeral & o);
166
167 void div_mul(numeral const & k, interval const & a, interval & b, bool inv_k);
168
169 void checkpoint();
170
171 public:
172 interval_manager(reslimit& lim, C && c);
173 ~interval_manager();
174
m()175 numeral_manager & m() const { return m_c.m(); }
176
177 void del(interval & a);
178
lower(interval const & a)179 numeral const & lower(interval const & a) const { return m_c.lower(a); }
upper(interval const & a)180 numeral const & upper(interval const & a) const { return m_c.upper(a); }
lower(interval & a)181 numeral & lower(interval & a) { return m_c.lower(a); }
upper(interval & a)182 numeral & upper(interval & a) { return m_c.upper(a); }
lower_is_open(interval const & a)183 bool lower_is_open(interval const & a) const { return m_c.lower_is_open(a); }
upper_is_open(interval const & a)184 bool upper_is_open(interval const & a) const { return m_c.upper_is_open(a); }
lower_is_inf(interval const & a)185 bool lower_is_inf(interval const & a) const { return m_c.lower_is_inf(a); }
upper_is_inf(interval const & a)186 bool upper_is_inf(interval const & a) const { return m_c.upper_is_inf(a); }
187
lower_is_neg(interval const & a)188 bool lower_is_neg(interval const & a) const { return ::is_neg(m(), lower(a), lower_kind(a)); }
lower_is_pos(interval const & a)189 bool lower_is_pos(interval const & a) const { return ::is_pos(m(), lower(a), lower_kind(a)); }
lower_is_zero(interval const & a)190 bool lower_is_zero(interval const & a) const { return ::is_zero(m(), lower(a), lower_kind(a)); }
191
upper_is_neg(interval const & a)192 bool upper_is_neg(interval const & a) const { return ::is_neg(m(), upper(a), upper_kind(a)); }
upper_is_pos(interval const & a)193 bool upper_is_pos(interval const & a) const { return ::is_pos(m(), upper(a), upper_kind(a)); }
upper_is_zero(interval const & a)194 bool upper_is_zero(interval const & a) const { return ::is_zero(m(), upper(a), upper_kind(a)); }
195
is_P(interval const & n)196 bool is_P(interval const & n) const { return lower_is_pos(n) || lower_is_zero(n); }
is_P0(interval const & n)197 bool is_P0(interval const & n) const { return lower_is_zero(n) && !lower_is_open(n); }
is_P1(interval const & n)198 bool is_P1(interval const & n) const { return lower_is_pos(n) || (lower_is_zero(n) && lower_is_open(n)); }
is_N(interval const & n)199 bool is_N(interval const & n) const { return upper_is_neg(n) || upper_is_zero(n); }
is_N0(interval const & n)200 bool is_N0(interval const & n) const { return upper_is_zero(n) && !upper_is_open(n); }
is_N1(interval const & n)201 bool is_N1(interval const & n) const { return upper_is_neg(n) || (upper_is_zero(n) && upper_is_open(n)); }
is_M(interval const & n)202 bool is_M(interval const & n) const { return lower_is_neg(n) && upper_is_pos(n); }
is_zero(interval const & n)203 bool is_zero(interval const & n) const { return lower_is_zero(n) && upper_is_zero(n); }
204
205 void set(interval & t, interval const & s);
206
207 void set(interval & t, numeral const& n);
208
209 bool eq(interval const & a, interval const & b) const;
210
211 /**
212 \brief Return true if all values in 'a' are less than all values in 'b'.
213 */
214 bool before(interval const & a, interval const & b) const;
215
216 /**
217 \brief Set lower bound to -oo.
218 */
219 void reset_lower(interval & a);
220
221 /**
222 \brief Set upper bound to +oo.
223 */
224 void reset_upper(interval & a);
225
226 /**
227 \brief Set interval to (-oo, oo)
228 */
229 void reset(interval & a);
230
231 /**
232 \brief Return true if the given interval contains 0.
233 */
234 bool contains_zero(interval const & n) const;
235
236 /**
237 \brief Return true if n contains v.
238 */
239 bool contains(interval const & n, numeral const & v) const;
240
241 void display(std::ostream & out, interval const & n) const;
242 void display_pp(std::ostream & out, interval const & n) const;
243
244 bool check_invariant(interval const & n) const;
245
246 /**
247 \brief b <- -a
248 */
249 void neg(interval const & a, interval & b, interval_deps_combine_rule & b_deps);
250 void neg(interval const & a, interval & b);
251 void neg_jst(interval const & a, interval_deps_combine_rule & b_deps);
252
253 /**
254 \brief c <- a + b
255 */
256 void add(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps);
257 void add(interval const & a, interval const & b, interval & c);
258 void add_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps);
259
260 /**
261 \brief c <- a - b
262 */
263 void sub(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps);
264 void sub(interval const & a, interval const & b, interval & c);
265 void sub_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps);
266
267 /**
268 \brief b <- k * a
269 */
270 void mul(numeral const & k, interval const & a, interval & b, interval_deps_combine_rule & b_deps);
mul(numeral const & k,interval const & a,interval & b)271 void mul(numeral const & k, interval const & a, interval & b) { div_mul(k, a, b, false); }
272 void mul_jst(numeral const & k, interval const & a, interval_deps_combine_rule & b_deps);
273 /**
274 \brief b <- (n/d) * a
275 */
276 void mul(int n, int d, interval const & a, interval & b);
277
278 /**
279 \brief b <- a/k
280
281 \remark For imprecise numerals, this is not equivalent to
282 m().inv(k)
283 mul(k, a, b)
284
285 That is, we must invert k rounding towards +oo or -oo depending whether we
286 are computing a lower or upper bound.
287 */
288 void div(interval const & a, numeral const & k, interval & b, interval_deps_combine_rule & b_deps);
div(interval const & a,numeral const & k,interval & b)289 void div(interval const & a, numeral const & k, interval & b) { div_mul(k, a, b, true); }
div_jst(interval const & a,numeral const & k,interval_deps_combine_rule & b_deps)290 void div_jst(interval const & a, numeral const & k, interval_deps_combine_rule & b_deps) { mul_jst(k, a, b_deps); }
291
292 /**
293 \brief c <- a * b
294 */
295 void mul(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps);
296 void mul(interval const & a, interval const & b, interval & c);
297 void mul_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps);
298
299 /**
300 \brief b <- a^n
301 */
302 void power(interval const & a, unsigned n, interval & b, interval_deps_combine_rule & b_deps);
303 void power(interval const & a, unsigned n, interval & b);
304 void power_jst(interval const & a, unsigned n, interval_deps_combine_rule & b_deps);
305
306 /**
307 \brief b <- a^(1/n) with precision p.
308
309 \pre if n is even, then a must not contain negative numbers.
310 */
311 void nth_root(interval const & a, unsigned n, numeral const & p, interval & b, interval_deps_combine_rule & b_deps);
312 void nth_root(interval const & a, unsigned n, numeral const & p, interval & b);
313 void nth_root_jst(interval const & a, unsigned n, numeral const & p, interval_deps_combine_rule & b_deps);
314
315 /**
316 \brief Given an equation x^n = y and an interval for y, compute the solution set for x with precision p.
317
318 \pre if n is even, then !lower_is_neg(y)
319 */
320 void xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x, interval_deps_combine_rule & x_deps);
321 void xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x);
322 void xn_eq_y_jst(interval const & y, unsigned n, numeral const & p, interval_deps_combine_rule & x_deps);
323
324 /**
325 \brief b <- 1/a
326 \pre !contains_zero(a)
327 */
328 void inv(interval const & a, interval & b, interval_deps_combine_rule & b_deps);
329 void inv(interval const & a, interval & b);
330 void inv_jst(interval const & a, interval_deps_combine_rule & b_deps);
331
332 /**
333 \brief c <- a/b
334 \pre !contains_zero(b)
335 \pre &a == &c (that is, c should not be an alias for a)
336 */
337 void div(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps);
338 void div(interval const & a, interval const & b, interval & c);
339 void div_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps);
340
341 /**
342 \brief Store in r an interval that contains the number pi.
343 The size of the interval is (1/15)*(1/16^n)
344 */
345 void pi(unsigned n, interval & r);
346
347 /**
348 \brief Set the precision of the internal interval representing pi.
349 */
350 void set_pi_prec(unsigned n);
351
352 /**
353 \brief Set the precision of the internal interval representing pi to a precision of at least n.
354 */
355 void set_pi_at_least_prec(unsigned n);
356
357 void sine(numeral const & a, unsigned k, numeral & lo, numeral & hi);
358
359 void cosine(numeral const & a, unsigned k, numeral & lo, numeral & hi);
360
361 /**
362 \brief Store in r the Euler's constant e.
363 The size of the interval is 4/(k+1)!
364 */
365 void e(unsigned k, interval & r);
366 };
367
368 template<typename Manager>
369 class _scoped_interval {
370 public:
371 typedef typename Manager::interval interval;
372 private:
373 Manager & m_manager;
374 interval m_interval;
375 public:
_scoped_interval(Manager & m)376 _scoped_interval(Manager & m):m_manager(m) {}
~_scoped_interval()377 ~_scoped_interval() { m_manager.del(m_interval); }
378
m()379 Manager & m() const { return m_manager; }
380
381 operator interval const &() const { return m_interval; }
382 operator interval&() { return m_interval; }
get()383 interval const & get() const { return m_interval; }
get()384 interval & get() { return m_interval; }
385 interval * operator->() {
386 return &m_interval;
387 }
388 interval const * operator->() const {
389 return &m_interval;
390 }
391 };
392
393