1 /*++
2 Copyright (c) 2012 Microsoft Corporation
3
4 Module Name:
5
6 interval_def.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 "math/interval/interval.h"
22 #include "util/debug.h"
23 #include "util/trace.h"
24 #include "util/scoped_numeral.h"
25 #include "util/common_msgs.h"
26
27 #define DEFAULT_PI_PRECISION 2
28
29 // #define TRACE_NTH_ROOT
30
31 template<typename C>
interval_manager(reslimit & lim,C && c)32 interval_manager<C>::interval_manager(reslimit& lim, C && c): m_limit(lim), m_c(std::move(c)) {
33 m().set(m_minus_one, -1);
34 m().set(m_one, 1);
35 m_pi_n = 0;
36 }
37
38 template<typename C>
~interval_manager()39 interval_manager<C>::~interval_manager() {
40 del(m_pi_div_2);
41 del(m_pi);
42 del(m_3_pi_div_2);
43 del(m_2_pi);
44 m().del(m_result_lower);
45 m().del(m_result_upper);
46 m().del(m_mul_ad);
47 m().del(m_mul_bc);
48 m().del(m_mul_ac);
49 m().del(m_mul_bd);
50 m().del(m_minus_one);
51 m().del(m_one);
52 m().del(m_inv_k);
53 }
54
55 template<typename C>
del(interval & a)56 void interval_manager<C>::del(interval & a) {
57 m().del(lower(a));
58 m().del(upper(a));
59 }
60
61
62 template<typename C>
checkpoint()63 void interval_manager<C>::checkpoint() {
64 if (!m_limit.inc())
65 throw default_exception(Z3_CANCELED_MSG);
66 }
67
68 /*
69 Compute the n-th root of a with precision p. The result hi - lo <= p
70 lo and hi are lower/upper bounds for the value of the n-th root of a.
71 That is, the n-th root is in the interval [lo, hi]
72
73 If n is even, then a is assumed to be nonnegative.
74
75 If numeral_manager is not precise, the procedure does not guarantee the precision p.
76 */
77 template<typename C>
nth_root_slow(numeral const & a,unsigned n,numeral const & p,numeral & lo,numeral & hi)78 void interval_manager<C>::nth_root_slow(numeral const & a, unsigned n, numeral const & p, numeral & lo, numeral & hi) {
79 #ifdef TRACE_NTH_ROOT
80 static unsigned counter = 0;
81 static unsigned loop_counter = 0;
82 counter++;
83 if (counter % 1000 == 0)
84 std::cerr << "[nth-root] " << counter << " " << loop_counter << " " << ((double)loop_counter)/((double)counter) << std::endl;
85 #endif
86
87 bool n_is_even = (n % 2 == 0);
88 SASSERT(!n_is_even || m().is_nonneg(a));
89 if (m().is_zero(a) || m().is_one(a) || (!n_is_even && m().eq(a, m_minus_one))) {
90 m().set(lo, a);
91 m().set(hi, a);
92 return;
93 }
94 if (m().lt(a, m_minus_one)) {
95 m().set(lo, a);
96 m().set(hi, -1);
97 }
98 else if (m().is_neg(a)) {
99 m().set(lo, -1);
100 m().set(hi, 0);
101 }
102 else if (m().lt(a, m_one)) {
103 m().set(lo, 0);
104 m().set(hi, 1);
105 }
106 else {
107 m().set(lo, 1);
108 m().set(hi, a);
109 }
110 SASSERT(m().le(lo, hi));
111 _scoped_numeral<numeral_manager> c(m()), cn(m());
112 _scoped_numeral<numeral_manager> two(m());
113 m().set(two, 2);
114 while (true) {
115 checkpoint();
116 #ifdef TRACE_NTH_ROOT
117 loop_counter++;
118 #endif
119 m().add(hi, lo, c);
120 m().div(c, two, c);
121 if (m().precise()) {
122 m().power(c, n, cn);
123 if (m().gt(cn, a)) {
124 m().set(hi, c);
125 }
126 else if (m().eq(cn, a)) {
127 // root is precise
128 m().set(lo, c);
129 m().set(hi, c);
130 return;
131 }
132 else {
133 m().set(lo, c);
134 }
135 }
136 else {
137 round_to_minus_inf();
138 m().power(c, n, cn);
139 if (m().gt(cn, a)) {
140 m().set(hi, c);
141 }
142 else {
143 round_to_plus_inf();
144 m().power(c, n, cn);
145 if (m().lt(cn, a)) {
146 m().set(lo, c);
147 }
148 else {
149 // can't improve, numeral_manager is not precise enough,
150 // a is between round-to-minus-inf(c^n) and round-to-plus-inf(c^n)
151 return;
152 }
153 }
154 }
155 round_to_plus_inf();
156 m().sub(hi, lo, c);
157 if (m().le(c, p))
158 return; // result is precise enough
159 }
160 }
161
162 /**
163 \brief Store in o a rough approximation of a^1/n.
164
165 It uses 2^Floor[Floor(Log2(a))/n]
166
167 \pre is_pos(a)
168 */
169 template<typename C>
rough_approx_nth_root(numeral const & a,unsigned n,numeral & o)170 void interval_manager<C>::rough_approx_nth_root(numeral const & a, unsigned n, numeral & o) {
171 SASSERT(m().is_pos(a));
172 SASSERT(n > 0);
173 round_to_minus_inf();
174 unsigned k = m().prev_power_of_two(a);
175 m().set(o, 2);
176 m().power(o, k/n, o);
177 }
178
179 /*
180 Compute the n-th root of \c a with (suggested) precision p.
181 The only guarantee provided by this method is that a^(1/n) is in [lo, hi].
182
183 If n is even, then a is assumed to be nonnegative.
184 */
185 template<typename C>
nth_root(numeral const & a,unsigned n,numeral const & p,numeral & lo,numeral & hi)186 void interval_manager<C>::nth_root(numeral const & a, unsigned n, numeral const & p, numeral & lo, numeral & hi) {
187 // nth_root_slow(a, n, p, lo, hi);
188 // return;
189 SASSERT(n > 0);
190 SASSERT(n % 2 != 0 || m().is_nonneg(a));
191 if (n == 1 || m().is_zero(a) || m().is_one(a) || m().is_minus_one(a)) {
192 // easy cases: 1, -1, 0
193 m().set(lo, a);
194 m().set(hi, a);
195 return;
196 }
197 bool is_neg = m().is_neg(a);
198 _scoped_numeral<numeral_manager> A(m());
199 m().set(A, a);
200 m().abs(A);
201
202 nth_root_pos(A, n, p, lo, hi);
203 STRACE("nth_root_trace",
204 tout << "[nth-root] ("; m().display(tout, A); tout << ")^(1/" << n << ") >= "; m().display(tout, lo); tout << "\n";
205 tout << "[nth-root] ("; m().display(tout, A); tout << ")^(1/" << n << ") <= "; m().display(tout, hi); tout << "\n";);
206 if (is_neg) {
207 m().swap(lo, hi);
208 m().neg(lo);
209 m().neg(hi);
210 }
211 }
212
213 /**
214 r <- A/(x^n)
215
216 If to_plus_inf, then r >= A/(x^n)
217 If not to_plus_inf, then r <= A/(x^n)
218
219 */
220 template<typename C>
A_div_x_n(numeral const & A,numeral const & x,unsigned n,bool to_plus_inf,numeral & r)221 void interval_manager<C>::A_div_x_n(numeral const & A, numeral const & x, unsigned n, bool to_plus_inf, numeral & r) {
222 if (n == 1) {
223 if (m().precise()) {
224 m().div(A, x, r);
225 }
226 else {
227 set_rounding(to_plus_inf);
228 m().div(A, x, r);
229 }
230 }
231 else {
232 if (m().precise()) {
233 m().power(x, n, r);
234 m().div(A, r, r);
235 }
236 else {
237 set_rounding(!to_plus_inf);
238 m().power(x, n, r);
239 set_rounding(to_plus_inf);
240 m().div(A, r, r);
241 }
242 }
243 }
244
245 /**
246 \brief Compute an approximation of A^(1/n) using the sequence
247
248 x' = 1/n((n-1)*x + A/(x^(n-1)))
249
250 The computation stops when the difference between current and new x is less than p.
251 The procedure may not terminate if m() is not precise and p is very small.
252
253 */
254 template<typename C>
approx_nth_root(numeral const & A,unsigned n,numeral const & p,numeral & x)255 void interval_manager<C>::approx_nth_root(numeral const & A, unsigned n, numeral const & p, numeral & x) {
256 SASSERT(m().is_pos(A));
257 SASSERT(n > 1);
258 #ifdef TRACE_NTH_ROOT
259 static unsigned counter = 0;
260 static unsigned loop_counter = 0;
261 counter++;
262 if (counter % 1000 == 0)
263 std::cerr << "[nth-root] " << counter << " " << loop_counter << " " << ((double)loop_counter)/((double)counter) << std::endl;
264 #endif
265
266 _scoped_numeral<numeral_manager> x_prime(m()), d(m());
267
268 m().set(d, 1);
269 if (m().lt(A, d))
270 m().set(x, A);
271 else
272 rough_approx_nth_root(A, n, x);
273
274 round_to_minus_inf();
275
276 if (n == 2) {
277 _scoped_numeral<numeral_manager> two(m());
278 m().set(two, 2);
279 while (true) {
280 checkpoint();
281 #ifdef TRACE_NTH_ROOT
282 loop_counter++;
283 #endif
284 m().div(A, x, x_prime);
285 m().add(x, x_prime, x_prime);
286 m().div(x_prime, two, x_prime);
287 m().sub(x_prime, x, d);
288 m().abs(d);
289 m().swap(x, x_prime);
290 if (m().lt(d, p))
291 return;
292 }
293 }
294 else {
295 _scoped_numeral<numeral_manager> _n(m()), _n_1(m());
296 m().set(_n, n); // _n contains n
297 m().set(_n_1, n);
298 m().dec(_n_1); // _n_1 contains n-1
299
300 while (true) {
301 checkpoint();
302 #ifdef TRACE_NTH_ROOT
303 loop_counter++;
304 #endif
305 m().power(x, n-1, x_prime);
306 m().div(A, x_prime, x_prime);
307 m().mul(_n_1, x, d);
308 m().add(d, x_prime, x_prime);
309 m().div(x_prime, _n, x_prime);
310 m().sub(x_prime, x, d);
311 m().abs(d);
312 TRACE("nth_root",
313 tout << "A: "; m().display(tout, A); tout << "\n";
314 tout << "x: "; m().display(tout, x); tout << "\n";
315 tout << "x_prime: "; m().display(tout, x_prime); tout << "\n";
316 tout << "d: "; m().display(tout, d); tout << "\n";
317 );
318 m().swap(x, x_prime);
319 if (m().lt(d, p))
320 return;
321 }
322 }
323 }
324
325 template<typename C>
nth_root_pos(numeral const & A,unsigned n,numeral const & p,numeral & lo,numeral & hi)326 void interval_manager<C>::nth_root_pos(numeral const & A, unsigned n, numeral const & p, numeral & lo, numeral & hi) {
327 approx_nth_root(A, n, p, hi);
328 if (m().precise()) {
329 // Assuming hi has a upper bound for A^(n-1)
330 // Then, A/(x^(n-1)) must be lower bound
331 A_div_x_n(A, hi, n-1, false, lo);
332 // Check if we were wrong
333 if (m().lt(hi, lo)) {
334 // swap if wrong
335 m().swap(lo, hi);
336 }
337 }
338 else {
339 // Check if hi is really a upper bound for A^(n-1)
340 A_div_x_n(A, hi, n-1, true /* lo will be greater than the actual lower bound */, lo);
341 TRACE("nth_root_bug",
342 tout << "Assuming upper\n";
343 tout << "A: "; m().display(tout, A); tout << "\n";
344 tout << "hi: "; m().display(tout, hi); tout << "\n";
345 tout << "lo: "; m().display(tout, hi); tout << "\n";);
346 if (m().le(lo, hi)) {
347 // hi is really the upper bound
348 // Must compute lo again but approximating to -oo
349 A_div_x_n(A, hi, n-1, false, lo);
350 }
351 else {
352 // hi should be lower bound
353 m().swap(lo, hi);
354 // check if lo is lower bound
355 A_div_x_n(A, lo, n-1, false /* hi will less than the actual upper bound */, hi);
356 if (m().le(lo, hi)) {
357 // lo is really the lower bound
358 // Must compute hi again but approximating to +oo
359 A_div_x_n(A, lo, n-1, true, hi);
360 }
361 else {
362 // we don't have anything due to rounding errors
363 // Be supper conservative
364 // This should not really happen very often.
365 _scoped_numeral<numeral_manager> one(m());
366 if (m().lt(A, one)) {
367 m().set(lo, 0);
368 m().set(hi, 1);
369 }
370 else {
371 m().set(lo, 1);
372 m().set(hi, A);
373 }
374 }
375 }
376 }
377 }
378
379
380 /**
381 \brief o <- n!
382 */
383 template<typename C>
fact(unsigned n,numeral & o)384 void interval_manager<C>::fact(unsigned n, numeral & o) {
385 _scoped_numeral<numeral_manager> aux(m());
386 m().set(o, 1);
387 for (unsigned i = 2; i <= n; i++) {
388 m().set(aux, static_cast<int>(i));
389 m().mul(aux, o, o);
390 TRACE("fact_bug", tout << "i: " << i << ", o: " << m().to_rational_string(o) << "\n";);
391 }
392 }
393
394 template<typename C>
sine_series(numeral const & a,unsigned k,bool upper,numeral & o)395 void interval_manager<C>::sine_series(numeral const & a, unsigned k, bool upper, numeral & o) {
396 SASSERT(k % 2 == 1);
397 // Compute sine using taylor series up to k
398 // x - x^3/3! + x^5/5! - x^7/7! + ...
399 // The result should be greater than or equal to the actual value if upper == true
400 // Otherwise it must be less than or equal to the actual value.
401 // The argument upper only matter if the numeral_manager is not precise.
402
403 // Taylor series up to k with rounding to
404 _scoped_numeral<numeral_manager> f(m());
405 _scoped_numeral<numeral_manager> aux(m());
406 m().set(o, a);
407 bool sign = true;
408 bool upper_factor = !upper; // since the first sign is negative, we must minimize factor to maximize result
409 for (unsigned i = 3; i <= k; i+=2) {
410 TRACE("sine_bug", tout << "[begin-loop] o: " << m().to_rational_string(o) << "\ni: " << i << "\n";
411 tout << "upper: " << upper << ", upper_factor: " << upper_factor << "\n";
412 tout << "o (default): " << m().to_string(o) << "\n";);
413 set_rounding(upper_factor);
414 m().power(a, i, f);
415 TRACE("sine_bug", tout << "a^i " << m().to_rational_string(f) << "\n";);
416 set_rounding(!upper_factor);
417 fact(i, aux);
418 TRACE("sine_bug", tout << "i! " << m().to_rational_string(aux) << "\n";);
419 set_rounding(upper_factor);
420 m().div(f, aux, f);
421 TRACE("sine_bug", tout << "a^i/i! " << m().to_rational_string(f) << "\n";);
422 set_rounding(upper);
423 if (sign)
424 m().sub(o, f, o);
425 else
426 m().add(o, f, o);
427 TRACE("sine_bug", tout << "o: " << m().to_rational_string(o) << "\n";);
428 sign = !sign;
429 upper_factor = !upper_factor;
430 }
431 }
432
433 template<typename C>
sine(numeral const & a,unsigned k,numeral & lo,numeral & hi)434 void interval_manager<C>::sine(numeral const & a, unsigned k, numeral & lo, numeral & hi) {
435 TRACE("sine", tout << "sine(a), a: " << m().to_rational_string(a) << "\na: " << m().to_string(a) << "\n";);
436 SASSERT(&lo != &hi);
437 if (m().is_zero(a)) {
438 m().reset(lo);
439 m().reset(hi);
440 return;
441 }
442
443 // Compute sine using taylor series
444 // x - x^3/3! + x^5/5! - x^7/7! + ...
445 //
446 // Note that, the coefficient of even terms is 0.
447 // So, we force k to be odd to make sure the error is minimized.
448 if (k % 2 == 0)
449 k++;
450
451 // Taylor series error = |x|^(k+1)/(k+1)!
452 _scoped_numeral<numeral_manager> error(m());
453 _scoped_numeral<numeral_manager> aux(m());
454 round_to_plus_inf();
455 m().set(error, a);
456 if (m().is_neg(error))
457 m().neg(error);
458 m().power(error, k+1, error);
459 TRACE("sine", tout << "a^(k+1): " << m().to_rational_string(error) << "\nk : " << k << "\n";);
460 round_to_minus_inf();
461 fact(k+1, aux);
462 TRACE("sine", tout << "(k+1)!: " << m().to_rational_string(aux) << "\n";);
463 round_to_plus_inf();
464 m().div(error, aux, error);
465 TRACE("sine", tout << "error: " << m().to_rational_string(error) << "\n";);
466
467 // Taylor series up to k with rounding to -oo
468 sine_series(a, k, false, lo);
469
470 if (m().precise()) {
471 m().set(hi, lo);
472 m().sub(lo, error, lo);
473 if (m().lt(lo, m_minus_one)) {
474 m().set(lo, -1);
475 m().set(hi, 1);
476 }
477 else {
478 m().add(hi, error, hi);
479 }
480 }
481 else {
482 // We must recompute the series with rounding to +oo
483 TRACE("sine", tout << "lo before -error: " << m().to_rational_string(lo) << "\n";);
484 round_to_minus_inf();
485 m().sub(lo, error, lo);
486 TRACE("sine", tout << "lo: " << m().to_rational_string(lo) << "\n";);
487 if (m().lt(lo, m_minus_one)) {
488 m().set(lo, -1);
489 m().set(hi, 1);
490 return;
491 }
492 sine_series(a, k, true, hi);
493 round_to_plus_inf();
494 m().add(hi, error, hi);
495 TRACE("sine", tout << "hi: " << m().to_rational_string(hi) << "\n";);
496 }
497 }
498
499 template<typename C>
cosine_series(numeral const & a,unsigned k,bool upper,numeral & o)500 void interval_manager<C>::cosine_series(numeral const & a, unsigned k, bool upper, numeral & o) {
501 SASSERT(k % 2 == 0);
502 // Compute cosine using taylor series up to k
503 // 1 - x^2/2! + x^4/4! - x^6/6! + ...
504 // The result should be greater than or equal to the actual value if upper == true
505 // Otherwise it must be less than or equal to the actual value.
506 // The argument upper only matter if the numeral_manager is not precise.
507
508
509 // Taylor series up to k with rounding to -oo
510 _scoped_numeral<numeral_manager> f(m());
511 _scoped_numeral<numeral_manager> aux(m());
512 m().set(o, 1);
513 bool sign = true;
514 bool upper_factor = !upper; // since the first sign is negative, we must minimize factor to maximize result
515 for (unsigned i = 2; i <= k; i+=2) {
516 set_rounding(upper_factor);
517 m().power(a, i, f);
518 set_rounding(!upper_factor);
519 fact(i, aux);
520 set_rounding(upper_factor);
521 m().div(f, aux, f);
522 set_rounding(upper);
523 if (sign)
524 m().sub(o, f, o);
525 else
526 m().add(o, f, o);
527 sign = !sign;
528 upper_factor = !upper_factor;
529 }
530 }
531
532 template<typename C>
cosine(numeral const & a,unsigned k,numeral & lo,numeral & hi)533 void interval_manager<C>::cosine(numeral const & a, unsigned k, numeral & lo, numeral & hi) {
534 TRACE("cosine", tout << "cosine(a): "; m().display_decimal(tout, a, 32); tout << "\n";);
535 SASSERT(&lo != &hi);
536 if (m().is_zero(a)) {
537 m().set(lo, 1);
538 m().set(hi, 1);
539 return;
540 }
541
542 // Compute cosine using taylor series
543 // 1 - x^2/2! + x^4/4! - x^6/6! + ...
544 //
545 // Note that, the coefficient of odd terms is 0.
546 // So, we force k to be even to make sure the error is minimized.
547 if (k % 2 == 1)
548 k++;
549
550 // Taylor series error = |x|^(k+1)/(k+1)!
551 _scoped_numeral<numeral_manager> error(m());
552 _scoped_numeral<numeral_manager> aux(m());
553 round_to_plus_inf();
554 m().set(error, a);
555 if (m().is_neg(error))
556 m().neg(error);
557 m().power(error, k+1, error);
558 round_to_minus_inf();
559 fact(k+1, aux);
560 round_to_plus_inf();
561 m().div(error, aux, error);
562 TRACE("sine", tout << "error: "; m().display_decimal(tout, error, 32); tout << "\n";);
563
564 // Taylor series up to k with rounding to -oo
565 cosine_series(a, k, false, lo);
566
567 if (m().precise()) {
568 m().set(hi, lo);
569 m().sub(lo, error, lo);
570 if (m().lt(lo, m_minus_one)) {
571 m().set(lo, -1);
572 m().set(hi, 1);
573 }
574 else {
575 m().add(hi, error, hi);
576 }
577 }
578 else {
579 // We must recompute the series with rounding to +oo
580 round_to_minus_inf();
581 m().sub(lo, error, lo);
582 if (m().lt(lo, m_minus_one)) {
583 m().set(lo, -1);
584 m().set(hi, 1);
585 return;
586 }
587 cosine_series(a, k, true, hi);
588 round_to_plus_inf();
589 m().add(hi, error, hi);
590 }
591 }
592
593 template<typename C>
reset_lower(interval & a)594 void interval_manager<C>::reset_lower(interval & a) {
595 m().reset(lower(a));
596 set_lower_is_open(a, true);
597 set_lower_is_inf(a, true);
598 }
599
600 template<typename C>
reset_upper(interval & a)601 void interval_manager<C>::reset_upper(interval & a) {
602 m().reset(upper(a));
603 set_upper_is_open(a, true);
604 set_upper_is_inf(a, true);
605 }
606
607 template<typename C>
reset(interval & a)608 void interval_manager<C>::reset(interval & a) {
609 reset_lower(a);
610 reset_upper(a);
611 }
612
613 template<typename C>
contains_zero(interval const & n)614 bool interval_manager<C>::contains_zero(interval const & n) const {
615 return
616 (lower_is_neg(n) || (lower_is_zero(n) && !lower_is_open(n))) &&
617 (upper_is_pos(n) || (upper_is_zero(n) && !upper_is_open(n)));
618 }
619
620
621 template<typename C>
contains(interval const & n,numeral const & v)622 bool interval_manager<C>::contains(interval const & n, numeral const & v) const {
623 if (!lower_is_inf(n)) {
624 if (m().lt(v, lower(n))) return false;
625 if (m().eq(v, lower(n)) && lower_is_open(n)) return false;
626 }
627 if (!upper_is_inf(n)) {
628 if (m().gt(v, upper(n))) return false;
629 if (m().eq(v, upper(n)) && upper_is_open(n)) return false;
630 }
631 return true;
632 }
633
634 template<typename C>
display(std::ostream & out,interval const & n)635 void interval_manager<C>::display(std::ostream & out, interval const & n) const {
636 out << (lower_is_open(n) ? "(" : "[");
637 ::display(out, m(), lower(n), lower_kind(n));
638 out << ", ";
639 ::display(out, m(), upper(n), upper_kind(n));
640 out << (upper_is_open(n) ? ")" : "]");
641 }
642
643 template<typename C>
display_pp(std::ostream & out,interval const & n)644 void interval_manager<C>::display_pp(std::ostream & out, interval const & n) const {
645 out << (lower_is_open(n) ? "(" : "[");
646 ::display_pp(out, m(), lower(n), lower_kind(n));
647 out << ", ";
648 ::display_pp(out, m(), upper(n), upper_kind(n));
649 out << (upper_is_open(n) ? ")" : "]");
650 }
651
652 template<typename C>
check_invariant(interval const & n)653 bool interval_manager<C>::check_invariant(interval const & n) const {
654 if (::eq(m(), lower(n), lower_kind(n), upper(n), upper_kind(n))) {
655 SASSERT(!lower_is_open(n));
656 SASSERT(!upper_is_open(n));
657 }
658 else {
659 SASSERT(lt(m(), lower(n), lower_kind(n), upper(n), upper_kind(n)));
660 }
661 return true;
662 }
663
664 template<typename C>
set(interval & t,interval const & s)665 void interval_manager<C>::set(interval & t, interval const & s) {
666 if (&t == &const_cast<interval&>(s))
667 return;
668 if (lower_is_inf(s)) {
669 set_lower_is_inf(t, true);
670 }
671 else {
672 m().set(lower(t), lower(s));
673 set_lower_is_inf(t, false);
674 }
675 if (upper_is_inf(s)) {
676 set_upper_is_inf(t, true);
677 }
678 else {
679 m().set(upper(t), upper(s));
680 set_upper_is_inf(t, false);
681 }
682 set_lower_is_open(t, lower_is_open(s));
683 set_upper_is_open(t, upper_is_open(s));
684 SASSERT(check_invariant(t));
685 }
686
687 template<typename C>
set(interval & t,numeral const & n)688 void interval_manager<C>::set(interval & t, numeral const& n) {
689 m().set(lower(t), n);
690 set_lower_is_inf(t, false);
691 m().set(upper(t), n);
692 set_upper_is_inf(t, false);
693 set_lower_is_open(t, false);
694 set_upper_is_open(t, false);
695 SASSERT(check_invariant(t));
696 }
697
698 template<typename C>
eq(interval const & a,interval const & b)699 bool interval_manager<C>::eq(interval const & a, interval const & b) const {
700 return
701 ::eq(m(), lower(a), lower_kind(a), lower(b), lower_kind(b)) &&
702 ::eq(m(), upper(a), upper_kind(a), upper(b), upper_kind(b)) &&
703 lower_is_open(a) == lower_is_open(b) &&
704 upper_is_open(a) == upper_is_open(b);
705 }
706
707 template<typename C>
before(interval const & a,interval const & b)708 bool interval_manager<C>::before(interval const & a, interval const & b) const {
709 if (upper_is_inf(a) || lower_is_inf(b))
710 return false;
711 return m().lt(upper(a), lower(b)) || (upper_is_open(a) && m().eq(upper(a), lower(b)));
712 }
713
714 template<typename C>
neg_jst(interval const & a,interval_deps_combine_rule & b_deps)715 void interval_manager<C>::neg_jst(interval const & a, interval_deps_combine_rule & b_deps) {
716 if (lower_is_inf(a)) {
717 if (upper_is_inf(a)) {
718 b_deps.m_lower_combine = 0;
719 b_deps.m_upper_combine = 0;
720 }
721 else {
722 b_deps.m_lower_combine = DEP_IN_UPPER1;
723 b_deps.m_upper_combine = 0;
724 }
725 }
726 else {
727 if (upper_is_inf(a)) {
728 b_deps.m_lower_combine = 0;
729 b_deps.m_upper_combine = DEP_IN_LOWER1;
730 }
731 else {
732 b_deps.m_lower_combine = DEP_IN_UPPER1;
733 b_deps.m_upper_combine = DEP_IN_LOWER1;
734 }
735 }
736 }
737
738 template<typename C>
neg(interval const & a,interval & b,interval_deps_combine_rule & b_deps)739 void interval_manager<C>::neg(interval const & a, interval & b, interval_deps_combine_rule & b_deps) {
740 neg_jst(a, b_deps);
741 neg(a, b);
742 }
743
744 template<typename C>
neg(interval const & a,interval & b)745 void interval_manager<C>::neg(interval const & a, interval & b) {
746 if (lower_is_inf(a)) {
747 if (upper_is_inf(a)) {
748 reset(b);
749 }
750 else {
751 m().set(lower(b), upper(a));
752 m().neg(lower(b));
753 set_lower_is_inf(b, false);
754 set_lower_is_open(b, upper_is_open(a));
755
756 m().reset(upper(b));
757 set_upper_is_inf(b, true);
758 set_upper_is_open(b, true);
759 }
760 }
761 else {
762 if (upper_is_inf(a)) {
763 m().set(upper(b), lower(a));
764 m().neg(upper(b));
765 set_upper_is_inf(b, false);
766 set_upper_is_open(b, lower_is_open(a));
767
768 m().reset(lower(b));
769 set_lower_is_inf(b, true);
770 set_lower_is_open(b, true);
771 }
772 else {
773 if (&a == &b) {
774 m().swap(lower(b), upper(b));
775 }
776 else {
777 m().set(lower(b), upper(a));
778 m().set(upper(b), lower(a));
779 }
780 m().neg(lower(b));
781 m().neg(upper(b));
782 set_lower_is_inf(b, false);
783 set_upper_is_inf(b, false);
784 bool l_o = lower_is_open(a);
785 bool u_o = upper_is_open(a);
786 set_lower_is_open(b, u_o);
787 set_upper_is_open(b, l_o);
788 }
789 }
790 SASSERT(check_invariant(b));
791 }
792
793 template<typename C>
add_jst(interval const & a,interval const & b,interval_deps_combine_rule & c_deps)794 void interval_manager<C>::add_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps) {
795 c_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2;
796 c_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2;
797 }
798
799 template<typename C>
add(interval const & a,interval const & b,interval & c,interval_deps_combine_rule & c_deps)800 void interval_manager<C>::add(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps) {
801 add_jst(a, b, c_deps);
802 add(a, b, c);
803 }
804
805 template<typename C>
add(interval const & a,interval const & b,interval & c)806 void interval_manager<C>::add(interval const & a, interval const & b, interval & c) {
807 ext_numeral_kind new_l_kind, new_u_kind;
808 round_to_minus_inf();
809 ::add(m(), lower(a), lower_kind(a), lower(b), lower_kind(b), lower(c), new_l_kind);
810 round_to_plus_inf();
811 ::add(m(), upper(a), upper_kind(a), upper(b), upper_kind(b), upper(c), new_u_kind);
812 set_lower_is_inf(c, new_l_kind == EN_MINUS_INFINITY);
813 set_upper_is_inf(c, new_u_kind == EN_PLUS_INFINITY);
814 set_lower_is_open(c, lower_is_open(a) || lower_is_open(b));
815 set_upper_is_open(c, upper_is_open(a) || upper_is_open(b));
816 SASSERT(check_invariant(c));
817 }
818
819 template<typename C>
sub_jst(interval const & a,interval const & b,interval_deps_combine_rule & c_deps)820 void interval_manager<C>::sub_jst(interval const & a, interval const & b, interval_deps_combine_rule & c_deps) {
821 c_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2;
822 c_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2;
823 }
824
825 template<typename C>
sub(interval const & a,interval const & b,interval & c,interval_deps_combine_rule & c_deps)826 void interval_manager<C>::sub(interval const & a, interval const & b, interval & c, interval_deps_combine_rule & c_deps) {
827 sub_jst(a, b, c_deps);
828 sub(a, b, c);
829 }
830
831 template<typename C>
sub(interval const & a,interval const & b,interval & c)832 void interval_manager<C>::sub(interval const & a, interval const & b, interval & c) {
833 ext_numeral_kind new_l_kind, new_u_kind;
834 round_to_minus_inf();
835 ::sub(m(), lower(a), lower_kind(a), upper(b), upper_kind(b), lower(c), new_l_kind);
836 round_to_plus_inf();
837 ::sub(m(), upper(a), upper_kind(a), lower(b), lower_kind(b), upper(c), new_u_kind);
838 set_lower_is_inf(c, new_l_kind == EN_MINUS_INFINITY);
839 set_upper_is_inf(c, new_u_kind == EN_PLUS_INFINITY);
840 set_lower_is_open(c, lower_is_open(a) || upper_is_open(b));
841 set_upper_is_open(c, upper_is_open(a) || lower_is_open(b));
842 SASSERT(check_invariant(c));
843 }
844
845 template<typename C>
mul_jst(numeral const & k,interval const & a,interval_deps_combine_rule & b_deps)846 void interval_manager<C>::mul_jst(numeral const & k, interval const & a, interval_deps_combine_rule & b_deps) {
847 if (m().is_zero(k)) {
848 b_deps.m_lower_combine = 0;
849 b_deps.m_upper_combine = 0;
850 }
851 else if (m().is_neg(k)) {
852 b_deps.m_lower_combine = DEP_IN_UPPER1;
853 b_deps.m_upper_combine = DEP_IN_LOWER1;
854 }
855 else {
856 b_deps.m_lower_combine = DEP_IN_LOWER1;
857 b_deps.m_upper_combine = DEP_IN_UPPER1;
858 }
859 }
860
861 template<typename C>
div_mul(numeral const & k,interval const & a,interval & b,bool inv_k)862 void interval_manager<C>::div_mul(numeral const & k, interval const & a, interval & b, bool inv_k) {
863 if (m().is_zero(k)) {
864 reset(b);
865 }
866 else {
867 numeral const & l = lower(a); ext_numeral_kind l_k = lower_kind(a);
868 numeral const & u = upper(a); ext_numeral_kind u_k = upper_kind(a);
869 numeral & new_l_val = m_result_lower;
870 numeral & new_u_val = m_result_upper;
871 ext_numeral_kind new_l_kind, new_u_kind;
872 bool l_o = lower_is_open(a);
873 bool u_o = upper_is_open(a);
874 if (m().is_pos(k)) {
875 set_lower_is_open(b, l_o);
876 set_upper_is_open(b, u_o);
877 if (inv_k) {
878 round_to_minus_inf();
879 m().inv(k, m_inv_k);
880 ::mul(m(), l, l_k, m_inv_k, EN_NUMERAL, new_l_val, new_l_kind);
881
882 round_to_plus_inf();
883 m().inv(k, m_inv_k);
884 ::mul(m(), u, u_k, m_inv_k, EN_NUMERAL, new_u_val, new_u_kind);
885 }
886 else {
887 round_to_minus_inf();
888 ::mul(m(), l, l_k, k, EN_NUMERAL, new_l_val, new_l_kind);
889 round_to_plus_inf();
890 ::mul(m(), u, u_k, k, EN_NUMERAL, new_u_val, new_u_kind);
891 }
892 }
893 else {
894 set_lower_is_open(b, u_o);
895 set_upper_is_open(b, l_o);
896 if (inv_k) {
897 round_to_minus_inf();
898 m().inv(k, m_inv_k);
899 ::mul(m(), u, u_k, m_inv_k, EN_NUMERAL, new_l_val, new_l_kind);
900
901 round_to_plus_inf();
902 m().inv(k, m_inv_k);
903 ::mul(m(), l, l_k, m_inv_k, EN_NUMERAL, new_u_val, new_u_kind);
904 }
905 else {
906 round_to_minus_inf();
907 ::mul(m(), u, u_k, k, EN_NUMERAL, new_l_val, new_l_kind);
908 round_to_plus_inf();
909 ::mul(m(), l, l_k, k, EN_NUMERAL, new_u_val, new_u_kind);
910 }
911 }
912 m().swap(lower(b), new_l_val);
913 m().swap(upper(b), new_u_val);
914 set_lower_is_inf(b, new_l_kind == EN_MINUS_INFINITY);
915 set_upper_is_inf(b, new_u_kind == EN_PLUS_INFINITY);
916 }
917 }
918
919 template<typename C>
mul(numeral const & k,interval const & a,interval & b,interval_deps_combine_rule & b_deps)920 void interval_manager<C>::mul(numeral const & k, interval const & a, interval & b, interval_deps_combine_rule & b_deps) {
921 mul_jst(k, a, b_deps);
922 mul(k, a, b);
923 }
924
925 template<typename C>
mul(int n,int d,interval const & a,interval & b)926 void interval_manager<C>::mul(int n, int d, interval const & a, interval & b) {
927 _scoped_numeral<numeral_manager> aux(m());
928 m().set(aux, n, d);
929 mul(aux, a, b);
930 }
931
932 template<typename C>
div(interval const & a,numeral const & k,interval & b,interval_deps_combine_rule & b_deps)933 void interval_manager<C>::div(interval const & a, numeral const & k, interval & b, interval_deps_combine_rule & b_deps) {
934 div_jst(a, k, b_deps);
935 div(a, k, b);
936 }
937
938 template<typename C>
mul_jst(interval const & i1,interval const & i2,interval_deps_combine_rule & r_deps)939 void interval_manager<C>::mul_jst(interval const & i1, interval const & i2, interval_deps_combine_rule & r_deps) {
940 if (is_zero(i1)) {
941 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
942 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
943 }
944 else if (is_zero(i2)) {
945 r_deps.m_lower_combine = DEP_IN_LOWER2 | DEP_IN_UPPER2;
946 r_deps.m_upper_combine = DEP_IN_LOWER2 | DEP_IN_UPPER2;
947 }
948 else if (is_N(i1)) {
949 if (is_N(i2)) {
950 // x <= b <= 0, y <= d <= 0 --> b*d <= x*y
951 // a <= x <= b <= 0, c <= y <= d <= 0 --> x*y <= a*c (we can use the fact that x or y is always negative (i.e., b is neg or d is neg))
952 r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2;
953 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER1; // we can replace DEP_IN_UPPER1 with DEP_IN_UPPER2
954 }
955 else if (is_M(i2)) {
956 // a <= x <= b <= 0, y <= d, d > 0 --> a*d <= x*y (uses the fact that b is not positive)
957 // a <= x <= b <= 0, c <= y, c < 0 --> x*y <= a*c (uses the fact that b is not positive)
958 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_UPPER1;
959 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER1;
960 }
961 else {
962 // a <= x <= b <= 0, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that x is neg (b is not positive) or y is pos (c is not negative))
963 // x <= b <= 0, 0 <= c <= y --> x*y <= b*c
964 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_UPPER1; // we can replace DEP_IN_UPPER1 with DEP_IN_UPPER2
965 r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2;
966 }
967 }
968 else if (is_M(i1)) {
969 if (is_N(i2)) {
970 // b > 0, x <= b, c <= y <= d <= 0 --> b*c <= x*y (uses the fact that d is not positive)
971 // a < 0, a <= x, c <= y <= d <= 0 --> x*y <= a*c (uses the fact that d is not positive)
972 r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
973 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
974 }
975 else if (is_M(i2)) {
976 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
977 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
978 }
979 else {
980 // a < 0, a <= x, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that c is not negative)
981 // b > 0, x <= b, 0 <= c <= y <= d --> x*y <= b*d (uses the fact that c is not negative)
982 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2 | DEP_IN_LOWER2;
983 r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER2;
984 }
985 }
986 else {
987 SASSERT(is_P(i1));
988 if (is_N(i2)) {
989 // 0 <= a <= x <= b, c <= y <= d <= 0 --> x*y <= b*c (uses the fact that x is pos (a is not neg) or y is neg (d is not pos))
990 // 0 <= a <= x, y <= d <= 0 --> a*d <= x*y
991 r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_UPPER2
992 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2;
993 }
994 else if (is_M(i2)) {
995 // 0 <= a <= x <= b, c <= y --> b*c <= x*y (uses the fact that a is not negative)
996 // 0 <= a <= x <= b, y <= d --> x*y <= b*d (uses the fact that a is not negative)
997 r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_LOWER1;
998 r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1;
999 }
1000 else {
1001 SASSERT(is_P(i2));
1002 // 0 <= a <= x, 0 <= c <= y --> a*c <= x*y
1003 // x <= b, y <= d --> x*y <= b*d (uses the fact that x is pos (a is not negative) or y is pos (c is not negative))
1004 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2;
1005 r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2 | DEP_IN_LOWER1; // we can replace DEP_IN_LOWER1 with DEP_IN_LOWER2
1006 }
1007 }
1008 }
1009
1010 template<typename C>
mul(interval const & i1,interval const & i2,interval & r,interval_deps_combine_rule & r_deps)1011 void interval_manager<C>::mul(interval const & i1, interval const & i2, interval & r, interval_deps_combine_rule & r_deps) {
1012 mul_jst(i1, i2, r_deps);
1013 mul(i1, i2, r);
1014 }
1015
1016 template<typename C>
mul(interval const & i1,interval const & i2,interval & r)1017 void interval_manager<C>::mul(interval const & i1, interval const & i2, interval & r) {
1018 #ifdef _TRACE
1019 static unsigned call_id = 0;
1020 #endif
1021 #if Z3DEBUG
1022 bool i1_contains_zero = contains_zero(i1);
1023 bool i2_contains_zero = contains_zero(i2);
1024 #endif
1025 if (is_zero(i1)) {
1026 set(r, i1);
1027 return;
1028 }
1029 if (is_zero(i2)) {
1030 set(r, i2);
1031 return;
1032 }
1033
1034 numeral const & a = lower(i1); ext_numeral_kind a_k = lower_kind(i1);
1035 numeral const & b = upper(i1); ext_numeral_kind b_k = upper_kind(i1);
1036 numeral const & c = lower(i2); ext_numeral_kind c_k = lower_kind(i2);
1037 numeral const & d = upper(i2); ext_numeral_kind d_k = upper_kind(i2);
1038
1039 bool a_o = lower_is_open(i1);
1040 bool b_o = upper_is_open(i1);
1041 bool c_o = lower_is_open(i2);
1042 bool d_o = upper_is_open(i2);
1043
1044 numeral & new_l_val = m_result_lower;
1045 numeral & new_u_val = m_result_upper;
1046 ext_numeral_kind new_l_kind, new_u_kind;
1047
1048 if (is_N(i1)) {
1049 if (is_N(i2)) {
1050 // x <= b <= 0, y <= d <= 0 --> b*d <= x*y
1051 // a <= x <= b <= 0, c <= y <= d <= 0 --> x*y <= a*c (we can use the fact that x or y is always negative (i.e., b is neg or d is neg))
1052 TRACE("interval_bug", tout << "(N, N) #" << call_id << "\n"; display(tout, i1); tout << "\n"; display(tout, i2); tout << "\n";
1053 tout << "a: "; m().display(tout, a); tout << "\n";
1054 tout << "b: "; m().display(tout, b); tout << "\n";
1055 tout << "c: "; m().display(tout, c); tout << "\n";
1056 tout << "d: "; m().display(tout, d); tout << "\n";
1057 tout << "is_N0(i1): " << is_N0(i1) << "\n";
1058 tout << "is_N0(i2): " << is_N0(i2) << "\n";
1059 );
1060 set_lower_is_open(r, (is_N0(i1) || is_N0(i2)) ? false : (b_o || d_o));
1061 set_upper_is_open(r, a_o || c_o);
1062 // if b = 0 (and the interval is closed), then the lower bound is closed
1063
1064 round_to_minus_inf();
1065 ::mul(m(), b, b_k, d, d_k, new_l_val, new_l_kind);
1066 round_to_plus_inf();
1067 ::mul(m(), a, a_k, c, c_k, new_u_val, new_u_kind);
1068 }
1069 else if (is_M(i2)) {
1070 // a <= x <= b <= 0, y <= d, d > 0 --> a*d <= x*y (uses the fact that b is not positive)
1071 // a <= x <= b <= 0, c <= y, c < 0 --> x*y <= a*c (uses the fact that b is not positive)
1072 TRACE("interval_bug", tout << "(N, M) #" << call_id << "\n";);
1073
1074 set_lower_is_open(r, a_o || d_o);
1075 set_upper_is_open(r, a_o || c_o);
1076
1077 round_to_minus_inf();
1078 ::mul(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
1079 round_to_plus_inf();
1080 ::mul(m(), a, a_k, c, c_k, new_u_val, new_u_kind);
1081 }
1082 else {
1083 // a <= x <= b <= 0, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that x is neg (b is not positive) or y is pos (c is not negative))
1084 // x <= b <= 0, 0 <= c <= y --> x*y <= b*c
1085 TRACE("interval_bug", tout << "(N, P) #" << call_id << "\n";);
1086 SASSERT(is_P(i2));
1087
1088 // must update upper_is_open first, since value of is_N0(i1) and is_P0(i2) may be affected by update
1089 set_upper_is_open(r, (is_N0(i1) || is_P0(i2)) ? false : (b_o || c_o));
1090 set_lower_is_open(r, a_o || d_o);
1091
1092 round_to_minus_inf();
1093 ::mul(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
1094 round_to_plus_inf();
1095 ::mul(m(), b, b_k, c, c_k, new_u_val, new_u_kind);
1096 }
1097 }
1098 else if (is_M(i1)) {
1099 if (is_N(i2)) {
1100 // b > 0, x <= b, c <= y <= d <= 0 --> b*c <= x*y (uses the fact that d is not positive)
1101 // a < 0, a <= x, c <= y <= d <= 0 --> x*y <= a*c (uses the fact that d is not positive)
1102 TRACE("interval_bug", tout << "(M, N) #" << call_id << "\n";);
1103
1104 set_lower_is_open(r, b_o || c_o);
1105 set_upper_is_open(r, a_o || c_o);
1106
1107 round_to_minus_inf();
1108 ::mul(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
1109 round_to_plus_inf();
1110 ::mul(m(), a, a_k, c, c_k, new_u_val, new_u_kind);
1111 }
1112 else if (is_M(i2)) {
1113 numeral & ad = m_mul_ad; ext_numeral_kind ad_k;
1114 numeral & bc = m_mul_bc; ext_numeral_kind bc_k;
1115 numeral & ac = m_mul_ac; ext_numeral_kind ac_k;
1116 numeral & bd = m_mul_bd; ext_numeral_kind bd_k;
1117
1118 bool ad_o = a_o || d_o;
1119 bool bc_o = b_o || c_o;
1120 bool ac_o = a_o || c_o;
1121 bool bd_o = b_o || d_o;
1122
1123 round_to_minus_inf();
1124 ::mul(m(), a, a_k, d, d_k, ad, ad_k);
1125 ::mul(m(), b, b_k, c, c_k, bc, bc_k);
1126 round_to_plus_inf();
1127 ::mul(m(), a, a_k, c, c_k, ac, ac_k);
1128 ::mul(m(), b, b_k, d, d_k, bd, bd_k);
1129
1130 if (::lt(m(), ad, ad_k, bc, bc_k) || (::eq(m(), ad, ad_k, bc, bc_k) && !ad_o && bc_o)) {
1131 m().swap(new_l_val, ad);
1132 new_l_kind = ad_k;
1133 set_lower_is_open(r, ad_o);
1134 }
1135 else {
1136 m().swap(new_l_val, bc);
1137 new_l_kind = bc_k;
1138 set_lower_is_open(r, bc_o);
1139 }
1140
1141
1142 if (::gt(m(), ac, ac_k, bd, bd_k) || (::eq(m(), ac, ac_k, bd, bd_k) && !ac_o && bd_o)) {
1143 m().swap(new_u_val, ac);
1144 new_u_kind = ac_k;
1145 set_upper_is_open(r, ac_o);
1146 }
1147 else {
1148 m().swap(new_u_val, bd);
1149 new_u_kind = bd_k;
1150 set_upper_is_open(r, bd_o);
1151 }
1152 }
1153 else {
1154 // a < 0, a <= x, 0 <= c <= y <= d --> a*d <= x*y (uses the fact that c is not negative)
1155 // b > 0, x <= b, 0 <= c <= y <= d --> x*y <= b*d (uses the fact that c is not negative)
1156 TRACE("interval_bug", tout << "(M, P) #" << call_id << "\n";);
1157 SASSERT(is_P(i2));
1158
1159 set_lower_is_open(r, a_o || d_o);
1160 set_upper_is_open(r, b_o || d_o);
1161
1162 round_to_minus_inf();
1163 ::mul(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
1164 round_to_plus_inf();
1165 ::mul(m(), b, b_k, d, d_k, new_u_val, new_u_kind);
1166 }
1167 }
1168 else {
1169 SASSERT(is_P(i1));
1170 if (is_N(i2)) {
1171 // 0 <= a <= x <= b, c <= y <= d <= 0 --> x*y <= b*c (uses the fact that x is pos (a is not neg) or y is neg (d is not pos))
1172 // 0 <= a <= x, y <= d <= 0 --> a*d <= x*y
1173 TRACE("interval_bug", tout << "(P, N) #" << call_id << "\n";);
1174
1175 // must update upper_is_open first, since value of is_P0(i1) and is_N0(i2) may be affected by update
1176 set_upper_is_open(r, (is_P0(i1) || is_N0(i2)) ? false : a_o || d_o);
1177 set_lower_is_open(r, b_o || c_o);
1178
1179 round_to_minus_inf();
1180 ::mul(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
1181 round_to_plus_inf();
1182 ::mul(m(), a, a_k, d, d_k, new_u_val, new_u_kind);
1183 }
1184 else if (is_M(i2)) {
1185 // 0 <= a <= x <= b, c <= y --> b*c <= x*y (uses the fact that a is not negative)
1186 // 0 <= a <= x <= b, y <= d --> x*y <= b*d (uses the fact that a is not negative)
1187 TRACE("interval_bug", tout << "(P, M) #" << call_id << "\n";);
1188
1189 set_lower_is_open(r, b_o || c_o);
1190 set_upper_is_open(r, b_o || d_o);
1191
1192 round_to_minus_inf();
1193 ::mul(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
1194 round_to_plus_inf();
1195 ::mul(m(), b, b_k, d, d_k, new_u_val, new_u_kind);
1196 }
1197 else {
1198 SASSERT(is_P(i2));
1199 // 0 <= a <= x, 0 <= c <= y --> a*c <= x*y
1200 // x <= b, y <= d --> x*y <= b*d (uses the fact that x is pos (a is not negative) or y is pos (c is not negative))
1201 TRACE("interval_bug", tout << "(P, P) #" << call_id << "\n";);
1202
1203 set_lower_is_open(r, (is_P0(i1) || is_P0(i2)) ? false : a_o || c_o);
1204 set_upper_is_open(r, b_o || d_o);
1205
1206 round_to_minus_inf();
1207 ::mul(m(), a, a_k, c, c_k, new_l_val, new_l_kind);
1208 round_to_plus_inf();
1209 ::mul(m(), b, b_k, d, d_k, new_u_val, new_u_kind);
1210 }
1211 }
1212
1213 m().swap(lower(r), new_l_val);
1214 m().swap(upper(r), new_u_val);
1215 set_lower_is_inf(r, new_l_kind == EN_MINUS_INFINITY);
1216 set_upper_is_inf(r, new_u_kind == EN_PLUS_INFINITY);
1217 SASSERT(!(i1_contains_zero || i2_contains_zero) || contains_zero(r));
1218 TRACE("interval_bug", tout << "result: "; display(tout, r); tout << "\n";);
1219 #ifdef _TRACE
1220 call_id++;
1221 #endif
1222 }
1223
1224 template<typename C>
power_jst(interval const & a,unsigned n,interval_deps_combine_rule & b_deps)1225 void interval_manager<C>::power_jst(interval const & a, unsigned n, interval_deps_combine_rule & b_deps) {
1226 if (n == 1) {
1227 b_deps.m_lower_combine = DEP_IN_LOWER1;
1228 b_deps.m_upper_combine = DEP_IN_UPPER1;
1229 }
1230 else if (n % 2 == 0) {
1231 if (lower_is_pos(a)) {
1232 // [l, u]^n = [l^n, u^n] if l > 0
1233 // 0 < l <= x --> l^n <= x^n (lower bound guarantees that is positive)
1234 // 0 < l <= x <= u --> x^n <= u^n (use lower and upper bound -- need the fact that x is positive)
1235 b_deps.m_lower_combine = DEP_IN_LOWER1;
1236 if (upper_is_inf(a))
1237 b_deps.m_upper_combine = 0;
1238 else
1239 b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
1240 }
1241 else if (upper_is_neg(a)) {
1242 // [l, u]^n = [u^n, l^n] if u < 0
1243 // l <= x <= u < 0 --> x^n <= l^n (use lower and upper bound -- need the fact that x is negative)
1244 // x <= u < 0 --> u^n <= x^n
1245 b_deps.m_lower_combine = DEP_IN_UPPER1;
1246 if (lower_is_inf(a))
1247 b_deps.m_upper_combine = 0;
1248 else
1249 b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
1250 }
1251 else {
1252 // [l, u]^n = [0, max{l^n, u^n}] otherwise
1253 // we need both bounds to justify upper bound
1254 b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
1255 b_deps.m_lower_combine = 0;
1256 }
1257 }
1258 else {
1259 // Remark: when n is odd x^n is monotonic.
1260 if (lower_is_inf(a))
1261 b_deps.m_lower_combine = 0;
1262 else
1263 b_deps.m_lower_combine = DEP_IN_LOWER1;
1264
1265 if (upper_is_inf(a))
1266 b_deps.m_upper_combine = 0;
1267 else
1268 b_deps.m_upper_combine = DEP_IN_UPPER1;
1269 }
1270 }
1271
1272 template<typename C>
power(interval const & a,unsigned n,interval & b,interval_deps_combine_rule & b_deps)1273 void interval_manager<C>::power(interval const & a, unsigned n, interval & b, interval_deps_combine_rule & b_deps) {
1274 power_jst(a, n, b_deps);
1275 power(a, n, b);
1276 }
1277
1278 template<typename C>
power(interval const & a,unsigned n,interval & b)1279 void interval_manager<C>::power(interval const & a, unsigned n, interval & b) {
1280 #ifdef _TRACE
1281 static unsigned call_id = 0;
1282 #endif
1283 if (n == 1) {
1284 set(b, a);
1285 }
1286 else if (n % 2 == 0) {
1287 if (lower_is_pos(a)) {
1288 // [l, u]^n = [l^n, u^n] if l > 0
1289 // 0 < l <= x --> l^n <= x^n (lower bound guarantees that is positive)
1290 // 0 < l <= x <= u --> x^n <= u^n (use lower and upper bound -- need the fact that x is positive)
1291 SASSERT(!lower_is_inf(a));
1292 round_to_minus_inf();
1293 m().power(lower(a), n, lower(b));
1294 set_lower_is_inf(b, false);
1295 set_lower_is_open(b, lower_is_open(a));
1296
1297 if (upper_is_inf(a)) {
1298 reset_upper(b);
1299 }
1300 else {
1301 round_to_plus_inf();
1302 m().power(upper(a), n, upper(b));
1303 set_upper_is_inf(b, false);
1304 set_upper_is_open(b, upper_is_open(a));
1305 }
1306 }
1307 else if (upper_is_neg(a)) {
1308 // [l, u]^n = [u^n, l^n] if u < 0
1309 // l <= x <= u < 0 --> x^n <= l^n (use lower and upper bound -- need the fact that x is negative)
1310 // x <= u < 0 --> u^n <= x^n
1311 SASSERT(!upper_is_inf(a));
1312 bool lower_a_open = lower_is_open(a), upper_a_open = upper_is_open(a);
1313 bool lower_a_inf = lower_is_inf(a);
1314
1315 m().set(lower(b), lower(a));
1316 m().set(upper(b), upper(a));
1317 m().swap(lower(b), upper(b)); // we use a swap because a and b can be aliased
1318
1319
1320 round_to_minus_inf();
1321 m().power(lower(b), n, lower(b));
1322
1323 set_lower_is_open(b, upper_a_open);
1324 set_lower_is_inf(b, false);
1325
1326 if (lower_a_inf) {
1327 reset_upper(b);
1328 }
1329 else {
1330 round_to_plus_inf();
1331 m().power(upper(b), n, upper(b));
1332 set_upper_is_inf(b, false);
1333 set_upper_is_open(b, lower_a_open);
1334 }
1335 }
1336 else {
1337 // [l, u]^n = [0, max{l^n, u^n}] otherwise
1338 // we need both bounds to justify upper bound
1339 TRACE("interval_bug", tout << "(M) #" << call_id << "\n"; display(tout, a); tout << "\nn:" << n << "\n";);
1340
1341 ext_numeral_kind un1_kind = lower_kind(a), un2_kind = upper_kind(a);
1342 numeral & un1 = m_result_lower;
1343 numeral & un2 = m_result_upper;
1344 m().set(un1, lower(a));
1345 m().set(un2, upper(a));
1346 round_to_plus_inf();
1347 ::power(m(), un1, un1_kind, n);
1348 ::power(m(), un2, un2_kind, n);
1349
1350 if (::gt(m(), un1, un1_kind, un2, un2_kind) || (::eq(m(), un1, un1_kind, un2, un2_kind) && !lower_is_open(a) && upper_is_open(a))) {
1351 m().swap(upper(b), un1);
1352 set_upper_is_inf(b, un1_kind == EN_PLUS_INFINITY);
1353 set_upper_is_open(b, lower_is_open(a));
1354 }
1355 else {
1356 m().swap(upper(b), un2);
1357 set_upper_is_inf(b, un2_kind == EN_PLUS_INFINITY);
1358 set_upper_is_open(b, upper_is_open(a));
1359 }
1360
1361 m().reset(lower(b));
1362 set_lower_is_inf(b, false);
1363 set_lower_is_open(b, false);
1364 }
1365 }
1366 else {
1367 // Remark: when n is odd x^n is monotonic.
1368 if (lower_is_inf(a)) {
1369 reset_lower(b);
1370 }
1371 else {
1372 m().power(lower(a), n, lower(b));
1373 set_lower_is_inf(b, false);
1374 set_lower_is_open(b, lower_is_open(a));
1375 }
1376
1377 if (upper_is_inf(a)) {
1378 reset_upper(b);
1379 }
1380 else {
1381 m().power(upper(a), n, upper(b));
1382 set_upper_is_inf(b, false);
1383 set_upper_is_open(b, upper_is_open(a));
1384 }
1385 }
1386 TRACE("interval_bug", tout << "result: "; display(tout, b); tout << "\n";);
1387 #ifdef _TRACE
1388 call_id++;
1389 #endif
1390 }
1391
1392
1393 template<typename C>
nth_root(interval const & a,unsigned n,numeral const & p,interval & b,interval_deps_combine_rule & b_deps)1394 void interval_manager<C>::nth_root(interval const & a, unsigned n, numeral const & p, interval & b, interval_deps_combine_rule & b_deps) {
1395 nth_root_jst(a, n, p, b_deps);
1396 nth_root(a, n, p, b);
1397 }
1398
1399 template<typename C>
nth_root(interval const & a,unsigned n,numeral const & p,interval & b)1400 void interval_manager<C>::nth_root(interval const & a, unsigned n, numeral const & p, interval & b) {
1401 SASSERT(n % 2 != 0 || !lower_is_neg(a));
1402 if (n == 1) {
1403 set(b, a);
1404 return;
1405 }
1406
1407 if (lower_is_inf(a)) {
1408 SASSERT(n % 2 != 0); // n must not be even.
1409 m().reset(lower(b));
1410 set_lower_is_inf(b, true);
1411 set_lower_is_open(b, true);
1412 }
1413 else {
1414 numeral & lo = m_result_lower;
1415 numeral & hi = m_result_upper;
1416 nth_root(lower(a), n, p, lo, hi);
1417 set_lower_is_inf(b, false);
1418 set_lower_is_open(b, lower_is_open(a) && m().eq(lo, hi));
1419 m().set(lower(b), lo);
1420 }
1421
1422 if (upper_is_inf(a)) {
1423 m().reset(upper(b));
1424 set_upper_is_inf(b, true);
1425 set_upper_is_open(b, true);
1426 }
1427 else {
1428 numeral & lo = m_result_lower;
1429 numeral & hi = m_result_upper;
1430 nth_root(upper(a), n, p, lo, hi);
1431 set_upper_is_inf(b, false);
1432 set_upper_is_open(b, upper_is_open(a) && m().eq(lo, hi));
1433 m().set(upper(b), hi);
1434 }
1435 TRACE("interval_nth_root", display(tout, a); tout << " --> "; display(tout, b); tout << "\n";);
1436 }
1437
1438 template<typename C>
nth_root_jst(interval const & a,unsigned n,numeral const & p,interval_deps_combine_rule & b_deps)1439 void interval_manager<C>::nth_root_jst(interval const & a, unsigned n, numeral const & p, interval_deps_combine_rule & b_deps) {
1440 b_deps.m_lower_combine = DEP_IN_LOWER1;
1441 if (n % 2 == 0)
1442 b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
1443 else
1444 b_deps.m_upper_combine = DEP_IN_UPPER1;
1445 }
1446
1447 template<typename C>
xn_eq_y(interval const & y,unsigned n,numeral const & p,interval & x,interval_deps_combine_rule & x_deps)1448 void interval_manager<C>::xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x, interval_deps_combine_rule & x_deps) {
1449 xn_eq_y_jst(y, n, p, x_deps);
1450 xn_eq_y(y, n, p, x);
1451 }
1452
1453 template<typename C>
xn_eq_y(interval const & y,unsigned n,numeral const & p,interval & x)1454 void interval_manager<C>::xn_eq_y(interval const & y, unsigned n, numeral const & p, interval & x) {
1455 SASSERT(n % 2 != 0 || !lower_is_neg(y));
1456 if (n % 2 == 0) {
1457 SASSERT(!lower_is_inf(y));
1458 if (upper_is_inf(y)) {
1459 reset(x);
1460 }
1461 else {
1462 numeral & lo = m_result_lower;
1463 numeral & hi = m_result_upper;
1464 nth_root(upper(y), n, p, lo, hi);
1465 // result is [-hi, hi]
1466 // result is open if upper(y) is open and lo == hi
1467 TRACE("interval_xn_eq_y", tout << "x^n = "; display(tout, y); tout << "\n";
1468 tout << "sqrt(y) in "; m().display(tout, lo); tout << " "; m().display(tout, hi); tout << "\n";);
1469 bool open = upper_is_open(y) && m().eq(lo, hi);
1470 set_lower_is_inf(x, false);
1471 set_upper_is_inf(x, false);
1472 set_lower_is_open(x, open);
1473 set_upper_is_open(x, open);
1474 m().set(upper(x), hi);
1475 round_to_minus_inf();
1476 m().set(lower(x), hi);
1477 m().neg(lower(x));
1478 TRACE("interval_xn_eq_y", tout << "interval for x: "; display(tout, x); tout << "\n";);
1479 }
1480 }
1481 else {
1482 SASSERT(n % 2 == 1); // n is odd
1483 nth_root(y, n, p, x);
1484 }
1485 }
1486
1487 template<typename C>
xn_eq_y_jst(interval const & y,unsigned n,numeral const & p,interval_deps_combine_rule & x_deps)1488 void interval_manager<C>::xn_eq_y_jst(interval const & y, unsigned n, numeral const & p, interval_deps_combine_rule & x_deps) {
1489 if (n % 2 == 0) {
1490 x_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
1491 x_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
1492 }
1493 else {
1494 x_deps.m_lower_combine = DEP_IN_LOWER1;
1495 x_deps.m_upper_combine = DEP_IN_UPPER1;
1496 }
1497 }
1498
1499 template<typename C>
inv_jst(interval const & a,interval_deps_combine_rule & b_deps)1500 void interval_manager<C>::inv_jst(interval const & a, interval_deps_combine_rule & b_deps) {
1501 SASSERT(!contains_zero(a));
1502 if (is_P1(a)) {
1503 b_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
1504 b_deps.m_upper_combine = DEP_IN_LOWER1;
1505 }
1506 else if (is_N1(a)) {
1507 // x <= u < 0 --> 1/u <= 1/x
1508 // l <= x <= u < 0 --> 1/l <= 1/x (use lower and upper bounds)
1509 b_deps.m_lower_combine = DEP_IN_UPPER1;
1510 b_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER1;
1511 }
1512 else {
1513 UNREACHABLE();
1514 }
1515 }
1516
1517 template<typename C>
inv(interval const & a,interval & b,interval_deps_combine_rule & b_deps)1518 void interval_manager<C>::inv(interval const & a, interval & b, interval_deps_combine_rule & b_deps) {
1519 inv_jst(a, b_deps);
1520 inv(a, b);
1521 }
1522
1523 template<typename C>
inv(interval const & a,interval & b)1524 void interval_manager<C>::inv(interval const & a, interval & b) {
1525 #ifdef _TRACE
1526 static unsigned call_id = 0;
1527 #endif
1528 // If the interval [l,u] does not contain 0, then 1/[l,u] = [1/u, 1/l]
1529 SASSERT(!contains_zero(a));
1530 TRACE("interval_bug", tout << "(inv) #" << call_id << "\n"; display(tout, a); tout << "\n";);
1531
1532 numeral & new_l_val = m_result_lower;
1533 numeral & new_u_val = m_result_upper;
1534 ext_numeral_kind new_l_kind, new_u_kind;
1535
1536 if (is_P1(a)) {
1537 // 0 < l <= x --> 1/x <= 1/l
1538 // 0 < l <= x <= u --> 1/u <= 1/x (use lower and upper bounds)
1539
1540 round_to_minus_inf();
1541 m().set(new_l_val, upper(a)); new_l_kind = upper_kind(a);
1542 ::inv(m(), new_l_val, new_l_kind);
1543 SASSERT(new_l_kind == EN_NUMERAL);
1544 bool new_l_open = upper_is_open(a);
1545
1546 if (lower_is_zero(a)) {
1547 SASSERT(lower_is_open(a));
1548 m().reset(upper(b));
1549 set_upper_is_inf(b, true);
1550 set_upper_is_open(b, true);
1551 }
1552 else {
1553 round_to_plus_inf();
1554 m().set(new_u_val, lower(a));
1555 m().inv(new_u_val);
1556 m().swap(upper(b), new_u_val);
1557 set_upper_is_inf(b, false);
1558 set_upper_is_open(b, lower_is_open(a));
1559 }
1560
1561 m().swap(lower(b), new_l_val);
1562 set_lower_is_inf(b, false);
1563 set_lower_is_open(b, new_l_open);
1564 }
1565 else if (is_N1(a)) {
1566 // x <= u < 0 --> 1/u <= 1/x
1567 // l <= x <= u < 0 --> 1/l <= 1/x (use lower and upper bounds)
1568
1569 round_to_plus_inf();
1570 m().set(new_u_val, lower(a)); new_u_kind = lower_kind(a);
1571 ::inv(m(), new_u_val, new_u_kind);
1572 SASSERT(new_u_kind == EN_NUMERAL);
1573 bool new_u_open = lower_is_open(a);
1574
1575 if (upper_is_zero(a)) {
1576 SASSERT(upper_is_open(a));
1577 m().reset(lower(b));
1578 set_lower_is_open(b, true);
1579 set_lower_is_inf(b, true);
1580 }
1581 else {
1582 round_to_minus_inf();
1583 m().set(new_l_val, upper(a));
1584 m().inv(new_l_val);
1585 m().swap(lower(b), new_l_val);
1586 set_lower_is_inf(b, false);
1587 set_lower_is_open(b, upper_is_open(a));
1588 }
1589
1590 m().swap(upper(b), new_u_val);
1591 set_upper_is_inf(b, false);
1592 set_upper_is_open(b, new_u_open);
1593 }
1594 else {
1595 UNREACHABLE();
1596 }
1597 TRACE("interval_bug", tout << "result: "; display(tout, b); tout << "\n";);
1598 #ifdef _TRACE
1599 call_id++;
1600 #endif
1601 }
1602
1603 template<typename C>
div_jst(interval const & i1,interval const & i2,interval_deps_combine_rule & r_deps)1604 void interval_manager<C>::div_jst(interval const & i1, interval const & i2, interval_deps_combine_rule & r_deps) {
1605 SASSERT(!contains_zero(i2));
1606 if (is_zero(i1)) {
1607 if (is_P1(i2)) {
1608 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2;
1609 r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2;
1610 }
1611 else {
1612 r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2;
1613 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2;
1614 }
1615 }
1616 else {
1617 if (is_N(i1)) {
1618 if (is_N1(i2)) {
1619 // x <= b <= 0, c <= y <= d < 0 --> b/c <= x/y
1620 // a <= x <= b <= 0, y <= d < 0 --> x/y <= a/d
1621 r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
1622 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2;
1623 }
1624 else {
1625 // a <= x, a < 0, 0 < c <= y --> a/c <= x/y
1626 // x <= b <= 0, 0 < c <= y <= d --> x/y <= b/d
1627 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2;
1628 r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
1629 }
1630 }
1631 else if (is_M(i1)) {
1632 if (is_N1(i2)) {
1633 // 0 < a <= x <= b < 0, y <= d < 0 --> b/d <= x/y
1634 // 0 < a <= x <= b < 0, y <= d < 0 --> x/y <= a/d
1635 r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2;
1636 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_UPPER2;
1637 }
1638 else {
1639 // 0 < a <= x <= b < 0, 0 < c <= y --> a/c <= x/y
1640 // 0 < a <= x <= b < 0, 0 < c <= y --> x/y <= b/c
1641 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2;
1642 r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2;
1643 }
1644 }
1645 else {
1646 SASSERT(is_P(i1));
1647 if (is_N1(i2)) {
1648 // b > 0, x <= b, c <= y <= d < 0 --> b/d <= x/y
1649 // 0 <= a <= x, c <= y <= d < 0 --> x/y <= a/c
1650 r_deps.m_lower_combine = DEP_IN_UPPER1 | DEP_IN_UPPER2;
1651 r_deps.m_upper_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
1652 }
1653 else {
1654 SASSERT(is_P1(i2));
1655 // 0 <= a <= x, 0 < c <= y <= d --> a/d <= x/y
1656 // b > 0 x <= b, 0 < c <= y --> x/y <= b/c
1657 r_deps.m_lower_combine = DEP_IN_LOWER1 | DEP_IN_LOWER2 | DEP_IN_UPPER2;
1658 r_deps.m_upper_combine = DEP_IN_UPPER1 | DEP_IN_LOWER2;
1659 }
1660 }
1661 }
1662 }
1663
1664 template<typename C>
div(interval const & i1,interval const & i2,interval & r,interval_deps_combine_rule & r_deps)1665 void interval_manager<C>::div(interval const & i1, interval const & i2, interval & r, interval_deps_combine_rule & r_deps) {
1666 div_jst(i1, i2, r_deps);
1667 div(i1, i2, r);
1668 }
1669
1670 template<typename C>
div(interval const & i1,interval const & i2,interval & r)1671 void interval_manager<C>::div(interval const & i1, interval const & i2, interval & r) {
1672 #ifdef _TRACE
1673 static unsigned call_id = 0;
1674 #endif
1675 SASSERT(!contains_zero(i2));
1676 SASSERT(&i1 != &r);
1677
1678 if (is_zero(i1)) {
1679 TRACE("interval_bug", tout << "div #" << call_id << "\n"; display(tout, i1); tout << "\n"; display(tout, i2); tout << "\n";);
1680
1681 // 0/other = 0 if other != 0
1682 m().reset(lower(r));
1683 m().reset(upper(r));
1684 set_lower_is_inf(r, false);
1685 set_upper_is_inf(r, false);
1686 set_lower_is_open(r, false);
1687 set_upper_is_open(r, false);
1688 }
1689 else {
1690 numeral const & a = lower(i1); ext_numeral_kind a_k = lower_kind(i1);
1691 numeral const & b = upper(i1); ext_numeral_kind b_k = upper_kind(i1);
1692 numeral const & c = lower(i2); ext_numeral_kind c_k = lower_kind(i2);
1693 numeral const & d = upper(i2); ext_numeral_kind d_k = upper_kind(i2);
1694
1695 bool a_o = lower_is_open(i1);
1696 bool b_o = upper_is_open(i1);
1697 bool c_o = lower_is_open(i2);
1698 bool d_o = upper_is_open(i2);
1699
1700 numeral & new_l_val = m_result_lower;
1701 numeral & new_u_val = m_result_upper;
1702 ext_numeral_kind new_l_kind, new_u_kind;
1703
1704 TRACE("interval_bug", tout << "div #" << call_id << "\n"; display(tout, i1); tout << "\n"; display(tout, i2); tout << "\n";
1705 tout << "a: "; m().display(tout, a); tout << "\n";
1706 tout << "b: "; m().display(tout, b); tout << "\n";
1707 tout << "c: "; m().display(tout, c); tout << "\n";
1708 tout << "d: "; m().display(tout, d); tout << "\n";
1709 );
1710
1711 if (is_N(i1)) {
1712 if (is_N1(i2)) {
1713 // x <= b <= 0, c <= y <= d < 0 --> b/c <= x/y
1714 // a <= x <= b <= 0, y <= d < 0 --> x/y <= a/d
1715 TRACE("interval_bug", tout << "(N, N) #" << call_id << "\n";);
1716
1717 set_lower_is_open(r, is_N0(i1) ? false : b_o || c_o);
1718 set_upper_is_open(r, a_o || d_o);
1719
1720 round_to_minus_inf();
1721 ::div(m(), b, b_k, c, c_k, new_l_val, new_l_kind);
1722 if (m().is_zero(d)) {
1723 SASSERT(d_o);
1724 m().reset(new_u_val);
1725 new_u_kind = EN_PLUS_INFINITY;
1726 }
1727 else {
1728 round_to_plus_inf();
1729 ::div(m(), a, a_k, d, d_k, new_u_val, new_u_kind);
1730 }
1731 }
1732 else {
1733 // a <= x, a < 0, 0 < c <= y --> a/c <= x/y
1734 // x <= b <= 0, 0 < c <= y <= d --> x/y <= b/d
1735 TRACE("interval_bug", tout << "(N, P) #" << call_id << "\n";);
1736 SASSERT(is_P1(i2));
1737
1738 set_upper_is_open(r, is_N0(i1) ? false : (b_o || d_o));
1739 set_lower_is_open(r, a_o || c_o);
1740
1741 if (m().is_zero(c)) {
1742 SASSERT(c_o);
1743 m().reset(new_l_val);
1744 new_l_kind = EN_MINUS_INFINITY;
1745 }
1746 else {
1747 round_to_minus_inf();
1748 ::div(m(), a, a_k, c, c_k, new_l_val, new_l_kind);
1749 }
1750 round_to_plus_inf();
1751 ::div(m(), b, b_k, d, d_k, new_u_val, new_u_kind);
1752 }
1753 }
1754 else if (is_M(i1)) {
1755 if (is_N1(i2)) {
1756 // 0 < a <= x <= b < 0, y <= d < 0 --> b/d <= x/y
1757 // 0 < a <= x <= b < 0, y <= d < 0 --> x/y <= a/d
1758 TRACE("interval_bug", tout << "(M, N) #" << call_id << "\n";);
1759
1760 set_lower_is_open(r, b_o || d_o);
1761 set_upper_is_open(r, a_o || d_o);
1762
1763 if (m().is_zero(d)) {
1764 SASSERT(d_o);
1765 m().reset(new_l_val); m().reset(new_u_val);
1766 new_l_kind = EN_MINUS_INFINITY;
1767 new_u_kind = EN_PLUS_INFINITY;
1768 }
1769 else {
1770 round_to_minus_inf();
1771 ::div(m(), b, b_k, d, d_k, new_l_val, new_l_kind);
1772 round_to_plus_inf();
1773 ::div(m(), a, a_k, d, d_k, new_u_val, new_u_kind);
1774 TRACE("interval_bug", tout << "new_l_kind: " << new_l_kind << ", new_u_kind: " << new_u_kind << "\n";);
1775 }
1776 }
1777 else {
1778 // 0 < a <= x <= b < 0, 0 < c <= y --> a/c <= x/y
1779 // 0 < a <= x <= b < 0, 0 < c <= y --> x/y <= b/c
1780
1781 TRACE("interval_bug", tout << "(M, P) #" << call_id << "\n";);
1782 SASSERT(is_P1(i2));
1783
1784 set_lower_is_open(r, a_o || c_o);
1785 set_upper_is_open(r, b_o || c_o);
1786
1787 if (m().is_zero(c)) {
1788 SASSERT(c_o);
1789 m().reset(new_l_val); m().reset(new_u_val);
1790 new_l_kind = EN_MINUS_INFINITY;
1791 new_u_kind = EN_PLUS_INFINITY;
1792 }
1793 else {
1794 round_to_minus_inf();
1795 ::div(m(), a, a_k, c, c_k, new_l_val, new_l_kind);
1796 round_to_plus_inf();
1797 ::div(m(), b, b_k, c, c_k, new_u_val, new_u_kind);
1798 }
1799 }
1800 }
1801 else {
1802 SASSERT(is_P(i1));
1803 if (is_N1(i2)) {
1804 // b > 0, x <= b, c <= y <= d < 0 --> b/d <= x/y
1805 // 0 <= a <= x, c <= y <= d < 0 --> x/y <= a/c
1806 TRACE("interval_bug", tout << "(P, N) #" << call_id << "\n";);
1807
1808 set_upper_is_open(r, is_P0(i1) ? false : a_o || c_o);
1809 set_lower_is_open(r, b_o || d_o);
1810
1811 if (m().is_zero(d)) {
1812 SASSERT(d_o);
1813 m().reset(new_l_val);
1814 new_l_kind = EN_MINUS_INFINITY;
1815 }
1816 else {
1817 round_to_minus_inf();
1818 ::div(m(), b, b_k, d, d_k, new_l_val, new_l_kind);
1819 }
1820 round_to_plus_inf();
1821 ::div(m(), a, a_k, c, c_k, new_u_val, new_u_kind);
1822 }
1823 else {
1824 SASSERT(is_P1(i2));
1825 // 0 <= a <= x, 0 < c <= y <= d --> a/d <= x/y
1826 // b > 0 x <= b, 0 < c <= y --> x/y <= b/c
1827 TRACE("interval_bug", tout << "(P, P) #" << call_id << "\n";);
1828
1829 set_lower_is_open(r, is_P0(i1) ? false : a_o || d_o);
1830 set_upper_is_open(r, b_o || c_o);
1831
1832 round_to_minus_inf();
1833 ::div(m(), a, a_k, d, d_k, new_l_val, new_l_kind);
1834 if (m().is_zero(c)) {
1835 SASSERT(c_o);
1836 m().reset(new_u_val);
1837 new_u_kind = EN_PLUS_INFINITY;
1838 }
1839 else {
1840 round_to_plus_inf();
1841 ::div(m(), b, b_k, c, c_k, new_u_val, new_u_kind);
1842 }
1843 }
1844 }
1845
1846 m().swap(lower(r), new_l_val);
1847 m().swap(upper(r), new_u_val);
1848 set_lower_is_inf(r, new_l_kind == EN_MINUS_INFINITY);
1849 set_upper_is_inf(r, new_u_kind == EN_PLUS_INFINITY);
1850 }
1851 TRACE("interval_bug", tout << "result: "; display(tout, r); tout << "\n";);
1852 #ifdef _TRACE
1853 call_id++;
1854 #endif
1855 }
1856
1857 template<typename C>
pi_series(int x,numeral & r,bool up)1858 void interval_manager<C>::pi_series(int x, numeral & r, bool up) {
1859 // Store in r the value: 1/16^x (4/(8x + 1) - 2/(8x + 4) - 1/(8x + 5) - 1/(8x + 6))
1860 _scoped_numeral<numeral_manager> f(m());
1861 set_rounding(up);
1862 m().set(r, 4, 8*x + 1);
1863 set_rounding(!up);
1864 m().set(f, 2, 8*x + 4);
1865 set_rounding(up);
1866 m().sub(r, f, r);
1867 set_rounding(!up);
1868 m().set(f, 1, 8*x + 5);
1869 set_rounding(up);
1870 m().sub(r, f, r);
1871 set_rounding(!up);
1872 m().set(f, 1, 8*x + 6);
1873 set_rounding(up);
1874 m().sub(r, f, r);
1875 m().set(f, 1, 16);
1876 m().power(f, x, f);
1877 m().mul(r, f, r);
1878 }
1879
1880 template<typename C>
pi(unsigned n,interval & r)1881 void interval_manager<C>::pi(unsigned n, interval & r) {
1882 // Compute an interval that contains pi using the series
1883 // P[0] + P[1] + ... + P[n]
1884 // where
1885 // P[n] := 1/16^x (4/(8x + 1) - 2/(8x + 4) - 1/(8x + 5) - 1/(8x + 6))
1886 //
1887 // The size of the interval is 1/15 * 1/(16^n)
1888 //
1889 // Lower is P[0] + P[1] + ... + P[n]
1890 // Upper is Lower + 1/15 * 1/(16^n)
1891
1892 // compute size of the resulting interval
1893 round_to_plus_inf(); // overestimate size of the interval
1894 _scoped_numeral<numeral_manager> len(m());
1895 _scoped_numeral<numeral_manager> p(m());
1896 m().set(len, 1, 16);
1897 m().power(len, n, len);
1898 m().set(p, 1, 15);
1899 m().mul(p, len, len);
1900
1901 // compute lower bound
1902 numeral & l_val = m_result_lower;
1903 m().reset(l_val);
1904 for (unsigned i = 0; i <= n; i++) {
1905 pi_series(i, p, false);
1906 round_to_minus_inf();
1907 m().add(l_val, p, l_val);
1908 }
1909
1910 // computer upper bound
1911 numeral & u_val = m_result_upper;
1912 if (m().precise()) {
1913 // the numeral manager is precise, so we do not need to recompute the series
1914 m().add(l_val, len, u_val);
1915 }
1916 else {
1917 // recompute the sum rounding to plus infinite
1918 m().reset(u_val);
1919 for (unsigned i = 0; i <= n; i++) {
1920 pi_series(i, p, true);
1921 round_to_plus_inf();
1922 m().add(u_val, p, u_val);
1923 }
1924 round_to_plus_inf();
1925 m().add(u_val, len, u_val);
1926 }
1927
1928 set_lower_is_open(r, false);
1929 set_upper_is_open(r, false);
1930 set_lower_is_inf(r, false);
1931 set_upper_is_inf(r, false);
1932 m().set(lower(r), l_val);
1933 m().set(upper(r), u_val);
1934 }
1935
1936 template<typename C>
set_pi_prec(unsigned n)1937 void interval_manager<C>::set_pi_prec(unsigned n) {
1938 SASSERT(n > 0);
1939 m_pi_n = n;
1940 pi(n, m_pi);
1941 mul(1, 2, m_pi, m_pi_div_2);
1942 mul(3, 2, m_pi, m_3_pi_div_2);
1943 mul(2, 1, m_pi, m_2_pi);
1944 }
1945
1946 template<typename C>
set_pi_at_least_prec(unsigned n)1947 void interval_manager<C>::set_pi_at_least_prec(unsigned n) {
1948 if (n > m_pi_n)
1949 set_pi_prec(n);
1950 }
1951
1952 template<typename C>
e_series(unsigned k,bool upper,numeral & o)1953 void interval_manager<C>::e_series(unsigned k, bool upper, numeral & o) {
1954 _scoped_numeral<numeral_manager> d(m()), a(m());
1955 m().set(o, 2);
1956 m().set(d, 1);
1957 for (unsigned i = 2; i <= k; i++) {
1958 set_rounding(!upper);
1959 m().set(a, static_cast<int>(i));
1960 m().mul(d, a, d); // d == i!
1961 m().set(a, d);
1962 set_rounding(upper);
1963 m().inv(a); // a == 1/i!
1964 m().add(o, a, o);
1965 }
1966 }
1967
1968 template<typename C>
e(unsigned k,interval & r)1969 void interval_manager<C>::e(unsigned k, interval & r) {
1970 // Store in r lower and upper bounds for Euler's constant.
1971 //
1972 // The procedure uses the series
1973 //
1974 // V = 1 + 1/1 + 1/2! + 1/3! + ... + 1/k!
1975 //
1976 // The error in the approximation above is <= E = 4/(k+1)!
1977 // Thus, e must be in the interval [V, V+E]
1978 numeral & lo = m_result_lower;
1979 numeral & hi = m_result_upper;
1980
1981 e_series(k, false, lo);
1982
1983 _scoped_numeral<numeral_manager> error(m()), aux(m());
1984 round_to_minus_inf();
1985 fact(k+1, error);
1986 round_to_plus_inf();
1987 m().inv(error); // error == 1/(k+1)!
1988 m().set(aux, 4);
1989 m().mul(aux, error, error); // error == 4/(k+1)!
1990
1991 if (m().precise()) {
1992 m().set(hi, lo);
1993 m().add(hi, error, hi);
1994 }
1995 else {
1996 e_series(k, true, hi);
1997 round_to_plus_inf();
1998 m().add(hi, error, hi);
1999 }
2000
2001 set_lower_is_open(r, false);
2002 set_upper_is_open(r, false);
2003 set_lower_is_inf(r, false);
2004 set_upper_is_inf(r, false);
2005 m().set(lower(r), lo);
2006 m().set(upper(r), hi);
2007 }
2008
2009