1 /*++
2 Copyright (c) 2011 Microsoft Corporation
3
4 Module Name:
5
6 mpbq.cpp
7
8 Abstract:
9
10 Binary Rational Numbers
11
12 A binary rational is a number of the form a/2^k.
13 All integers are binary rationals.
14 Binary rational numbers can be implemented more efficiently than rationals.
15 Binary rationals form a Ring.
16 They are not closed under division.
17 In Z3, they are used to implement algebraic numbers.
18 The root isolation operations only use division by 2.
19
20 Author:
21
22 Leonardo de Moura (leonardo) 2011-11-24.
23
24 Revision History:
25
26 --*/
27 #include<sstream>
28 #include "util/mpbq.h"
29
30 #ifdef Z3DEBUG
31 #define MPBQ_DEBUG
32 #endif
33
to_rational(mpbq const & v)34 rational to_rational(mpbq const & v) {
35 rational r(v.numerator());
36 rational twok;
37 twok = power(rational(2), v.k());
38 return r/twok;
39 }
40
mpbq_manager(unsynch_mpz_manager & m)41 mpbq_manager::mpbq_manager(unsynch_mpz_manager & m):
42 m_manager(m) {
43 }
44
~mpbq_manager()45 mpbq_manager::~mpbq_manager() {
46 del(m_addmul_tmp);
47 m_manager.del(m_tmp);
48 m_manager.del(m_tmp2);
49 m_manager.del(m_select_int_tmp1);
50 m_manager.del(m_select_int_tmp2);
51 m_manager.del(m_select_small_tmp);
52 del(m_select_small_tmp1);
53 del(m_select_small_tmp2);
54 m_manager.del(m_div_tmp1);
55 m_manager.del(m_div_tmp2);
56 m_manager.del(m_div_tmp3);
57 }
58
reset(mpbq_vector & v)59 void mpbq_manager::reset(mpbq_vector & v) {
60 unsigned sz = v.size();
61 for (unsigned i = 0; i < sz; i++)
62 reset(v[i]);
63 v.reset();
64 }
65
normalize(mpbq & a)66 void mpbq_manager::normalize(mpbq & a) {
67 if (a.m_k == 0)
68 return;
69 if (m_manager.is_zero(a.m_num)) {
70 a.m_k = 0;
71 return;
72 }
73 #ifdef MPBQ_DEBUG
74 rational r = to_rational(a);
75 #endif
76 unsigned k = m_manager.power_of_two_multiple(a.m_num);
77 if (k > a.m_k)
78 k = a.m_k;
79 m_manager.machine_div2k(a.m_num, k);
80 a.m_k -= k;
81 #ifdef MPBQ_DEBUG
82 rational new_r = to_rational(a);
83 SASSERT(r == new_r);
84 #endif
85 }
86
magnitude_lb(mpbq const & a)87 int mpbq_manager::magnitude_lb(mpbq const & a) {
88 if (m_manager.is_zero(a.m_num))
89 return 0;
90 if (m_manager.is_pos(a.m_num))
91 return m_manager.log2(a.m_num) - a.m_k;
92 return m_manager.mlog2(a.m_num) - a.m_k + 1;
93 }
94
magnitude_ub(mpbq const & a)95 int mpbq_manager::magnitude_ub(mpbq const & a) {
96 if (m_manager.is_zero(a.m_num))
97 return 0;
98 if (m_manager.is_pos(a.m_num))
99 return m_manager.log2(a.m_num) - a.m_k + 1;
100 return m_manager.mlog2(a.m_num) - a.m_k;
101 }
102
mul2(mpbq & a)103 void mpbq_manager::mul2(mpbq & a) {
104 if (a.m_k == 0)
105 m_manager.mul2k(a.m_num, 1);
106 else
107 a.m_k--;
108 }
109
mul2k(mpbq & a,unsigned k)110 void mpbq_manager::mul2k(mpbq & a, unsigned k) {
111 if (k == 0)
112 return;
113 if (a.m_k < k) {
114 m_manager.mul2k(a.m_num, k - a.m_k);
115 a.m_k = 0;
116 }
117 else {
118 SASSERT(a.m_k >= k);
119 a.m_k -= k;
120 }
121 }
122
add(mpbq const & a,mpbq const & b,mpbq & r)123 void mpbq_manager::add(mpbq const & a, mpbq const & b, mpbq & r) {
124 #ifdef MPBQ_DEBUG
125 rational _a = to_rational(a);
126 rational _b = to_rational(b);
127 #endif
128 if (a.m_k == b.m_k) {
129 m_manager.add(a.m_num, b.m_num, r.m_num);
130 r.m_k = a.m_k;
131 }
132 else if (a.m_k < b.m_k) {
133 m_manager.mul2k(a.m_num, b.m_k - a.m_k, m_tmp);
134 m_manager.add(b.m_num, m_tmp, r.m_num);
135 r.m_k = b.m_k;
136 }
137 else {
138 SASSERT(a.m_k > b.m_k);
139 m_manager.mul2k(b.m_num, a.m_k - b.m_k, m_tmp);
140 m_manager.add(a.m_num, m_tmp, r.m_num);
141 r.m_k = a.m_k;
142 }
143 normalize(r);
144 #ifdef MPBQ_DEBUG
145 rational _r = to_rational(r);
146 SASSERT(_a + _b == _r);
147 #endif
148 }
149
add(mpbq const & a,mpz const & b,mpbq & r)150 void mpbq_manager::add(mpbq const & a, mpz const & b, mpbq & r) {
151 #ifdef MPBQ_DEBUG
152 rational _a = to_rational(a);
153 rational _b(b);
154 #endif
155 if (a.m_k == 0) {
156 m_manager.add(a.m_num, b, r.m_num);
157 r.m_k = a.m_k;
158 }
159 else {
160 m_manager.mul2k(b, a.m_k, m_tmp);
161 m_manager.add(a.m_num, m_tmp, r.m_num);
162 r.m_k = a.m_k;
163 }
164 normalize(r);
165 #ifdef MPBQ_DEBUG
166 rational _r = to_rational(r);
167 TRACE("mpbq_bug", tout << "add a: " << _a << ", b: " << _b << ", r: " << _r << ", expected: " << (_a + _b) << "\n";);
168 SASSERT(_a + _b == _r);
169 #endif
170 }
171
sub(mpbq const & a,mpbq const & b,mpbq & r)172 void mpbq_manager::sub(mpbq const & a, mpbq const & b, mpbq & r) {
173 #ifdef MPBQ_DEBUG
174 rational _a = to_rational(a);
175 rational _b = to_rational(b);
176 #endif
177 if (a.m_k == b.m_k) {
178 m_manager.sub(a.m_num, b.m_num, r.m_num);
179 r.m_k = a.m_k;
180 }
181 else if (a.m_k < b.m_k) {
182 m_manager.mul2k(a.m_num, b.m_k - a.m_k, m_tmp);
183 m_manager.sub(m_tmp, b.m_num, r.m_num);
184 r.m_k = b.m_k;
185 }
186 else {
187 SASSERT(a.m_k > b.m_k);
188 m_manager.mul2k(b.m_num, a.m_k - b.m_k, m_tmp);
189 m_manager.sub(a.m_num, m_tmp, r.m_num);
190 r.m_k = a.m_k;
191 }
192 normalize(r);
193 #ifdef MPBQ_DEBUG
194 rational _r = to_rational(r);
195 TRACE("mpbq_bug", tout << "sub a: " << _a << ", b: " << _b << ", r: " << _r << ", expected: " << (_a - _b) << "\n";);
196 SASSERT(_a - _b == _r);
197 #endif
198 }
199
sub(mpbq const & a,mpz const & b,mpbq & r)200 void mpbq_manager::sub(mpbq const & a, mpz const & b, mpbq & r) {
201 #ifdef MPBQ_DEBUG
202 rational _a = to_rational(a);
203 rational _b(b);
204 #endif
205 if (a.m_k == 0) {
206 m_manager.sub(a.m_num, b, r.m_num);
207 r.m_k = a.m_k;
208 }
209 else {
210 m_manager.mul2k(b, a.m_k, m_tmp);
211 m_manager.sub(a.m_num, m_tmp, r.m_num);
212 r.m_k = a.m_k;
213 }
214 normalize(r);
215 #ifdef MPBQ_DEBUG
216 rational _r = to_rational(r);
217 SASSERT(_a - _b == _r);
218 #endif
219 }
220
mul(mpbq const & a,mpbq const & b,mpbq & r)221 void mpbq_manager::mul(mpbq const & a, mpbq const & b, mpbq & r) {
222 #ifdef MPBQ_DEBUG
223 rational _a = to_rational(a);
224 rational _b = to_rational(b);
225 #endif
226 m_manager.mul(a.m_num, b.m_num, r.m_num);
227 r.m_k = a.m_k + b.m_k;
228 if (a.m_k == 0 || b.m_k == 0) {
229 // if a.m_k and b.m_k are greater than 0, then there is no point in normalizing r.
230 normalize(r);
231 }
232 #ifdef MPBQ_DEBUG
233 rational _r = to_rational(r);
234 SASSERT(_a * _b == _r);
235 #endif
236 }
237
mul(mpbq const & a,mpz const & b,mpbq & r)238 void mpbq_manager::mul(mpbq const & a, mpz const & b, mpbq & r) {
239 #ifdef MPBQ_DEBUG
240 rational _a = to_rational(a);
241 rational _b(b);
242 #endif
243 m_manager.mul(a.m_num, b, r.m_num);
244 r.m_k = a.m_k;
245 normalize(r);
246 #ifdef MPBQ_DEBUG
247 rational _r = to_rational(r);
248 SASSERT(_a * _b == _r);
249 #endif
250 }
251
power(mpbq & a,unsigned k)252 void mpbq_manager::power(mpbq & a, unsigned k) {
253 SASSERT(static_cast<uint64_t>(k) * static_cast<uint64_t>(a.k()) <= static_cast<uint64_t>(UINT_MAX));
254 // We don't need to normalize because:
255 // If a.m_k == 0, then a is an integer, and the result be an integer
256 // If a.m_k > 0, then a.m_num must be odd, and the (a.m_num)^k will also be odd
257 a.m_k *= k;
258 m_manager.power(a.m_num, k, a.m_num);
259 }
260
root_lower(mpbq & a,unsigned n)261 bool mpbq_manager::root_lower(mpbq & a, unsigned n) {
262 bool r = m_manager.root(a.m_num, n);
263 if (!r)
264 m_manager.dec(a.m_num);
265 if (a.m_k % n == 0) {
266 a.m_k /= n;
267 normalize(a);
268 return r;
269 }
270 else if (m_manager.is_neg(a.m_num)) {
271 a.m_k /= n;
272 normalize(a);
273 return false;
274 }
275 else {
276 a.m_k /= n;
277 a.m_k++;
278 normalize(a);
279 return false;
280 }
281 }
282
root_upper(mpbq & a,unsigned n)283 bool mpbq_manager::root_upper(mpbq & a, unsigned n) {
284 bool r = m_manager.root(a.m_num, n);
285 if (a.m_k % n == 0) {
286 a.m_k /= n;
287 normalize(a);
288 return r;
289 }
290 else if (m_manager.is_neg(a.m_num)) {
291 a.m_k /= n;
292 a.m_k++;
293 normalize(a);
294 return false;
295 }
296 else {
297 a.m_k /= n;
298 normalize(a);
299 return false;
300 }
301 }
302
lt(mpbq const & a,mpbq const & b)303 bool mpbq_manager::lt(mpbq const & a, mpbq const & b) {
304 // TODO: try the following trick when k1 != k2
305 // Given, a = n1/2^k1 b = n2/2^k2
306 // Suppose n1 > 0 and n2 > 0,
307 // Then, we have, n1 <= 2^{log2(n1) - k1} 2^{log2(n2) - 1 - k2} <= n2
308 // Thus, log2(n1) - k1 < log2(n2) - 1 - k2 implies a < b
309 // Similarly: log2(n2) - k2 < log2(n1) - 1 - k1 implies b < a
310 // That is we compare the "magnitude" of the numbers before performing mul2k
311 //
312 // If n1 < 0 and n2 < 0, a similar trick can be implemented using mlog2 instead log2.
313 //
314 // It seems the trick is not useful when n1 and n2 are small
315 // numbers, and k1 and k2 very small < 8. Since, no bignumber
316 // computation is needed for mul2k.
317
318 if (a.m_k == b.m_k) {
319 return m_manager.lt(a.m_num, b.m_num);
320 }
321 else if (a.m_k < b.m_k) {
322 m_manager.mul2k(a.m_num, b.m_k - a.m_k, m_tmp);
323 return m_manager.lt(m_tmp, b.m_num);
324 }
325 else {
326 SASSERT(a.m_k > b.m_k);
327 m_manager.mul2k(b.m_num, a.m_k - b.m_k, m_tmp);
328 return m_manager.lt(a.m_num, m_tmp);
329 }
330 }
331
lt_1div2k(mpbq const & a,unsigned k)332 bool mpbq_manager::lt_1div2k(mpbq const & a, unsigned k) {
333 if (m_manager.is_nonpos(a.m_num))
334 return true;
335 if (a.m_k <= k) {
336 // since a.m_num >= 1
337 return false;
338 }
339 else {
340 SASSERT(a.m_k > k);
341 m_manager.mul2k(mpz(1), a.m_k - k, m_tmp);
342 return m_manager.lt(a.m_num, m_tmp);
343 }
344 }
345
eq(mpbq const & a,mpq const & b)346 bool mpbq_manager::eq(mpbq const & a, mpq const & b) {
347 if (is_int(a) && m_manager.is_one(b.denominator()))
348 return m_manager.eq(a.m_num, b.numerator());
349 m_manager.mul2k(b.numerator(), a.m_k, m_tmp);
350 m_manager.mul(a.m_num, b.denominator(), m_tmp2);
351 return m_manager.eq(m_tmp, m_tmp2);
352 }
353
lt(mpbq const & a,mpq const & b)354 bool mpbq_manager::lt(mpbq const & a, mpq const & b) {
355 if (is_int(a) && m_manager.is_one(b.denominator()))
356 return m_manager.lt(a.m_num, b.numerator());
357 m_manager.mul(a.m_num, b.denominator(), m_tmp);
358 m_manager.mul2k(b.numerator(), a.m_k, m_tmp2);
359 return m_manager.lt(m_tmp, m_tmp2);
360 }
361
le(mpbq const & a,mpq const & b)362 bool mpbq_manager::le(mpbq const & a, mpq const & b) {
363 if (is_int(a) && m_manager.is_one(b.denominator()))
364 return m_manager.le(a.m_num, b.numerator());
365 m_manager.mul(a.m_num, b.denominator(), m_tmp);
366 m_manager.mul2k(b.numerator(), a.m_k, m_tmp2);
367 return m_manager.le(m_tmp, m_tmp2);
368 }
369
lt(mpbq const & a,mpz const & b)370 bool mpbq_manager::lt(mpbq const & a, mpz const & b) {
371 if (is_int(a))
372 return m_manager.lt(a.m_num, b);
373 m_manager.mul2k(b, a.m_k, m_tmp);
374 return m_manager.lt(a.m_num, m_tmp);
375 }
376
le(mpbq const & a,mpz const & b)377 bool mpbq_manager::le(mpbq const & a, mpz const & b) {
378 if (is_int(a))
379 return m_manager.le(a.m_num, b);
380 m_manager.mul2k(b, a.m_k, m_tmp);
381 return m_manager.le(a.m_num, m_tmp);
382 }
383
to_string(mpbq const & a)384 std::string mpbq_manager::to_string(mpbq const & a) {
385 std::ostringstream buffer;
386 buffer << m_manager.to_string(a.m_num);
387 if (a.m_k == 1)
388 buffer << "/2";
389 else if (a.m_k > 1)
390 buffer << "/2^" << a.m_k;
391 return buffer.str();
392 }
393
display(std::ostream & out,mpbq const & a)394 void mpbq_manager::display(std::ostream & out, mpbq const & a) {
395 out << m_manager.to_string(a.m_num);
396 if (a.m_k > 0)
397 out << "/2";
398 if (a.m_k > 1)
399 out << "^" << a.m_k;
400 }
401
display_pp(std::ostream & out,mpbq const & a)402 void mpbq_manager::display_pp(std::ostream & out, mpbq const & a) {
403 out << m_manager.to_string(a.m_num);
404 if (a.m_k > 0)
405 out << "/2";
406 if (a.m_k > 1)
407 out << "<sup>" << a.m_k << "</sup>";
408 }
409
display_smt2(std::ostream & out,mpbq const & a,bool decimal)410 void mpbq_manager::display_smt2(std::ostream & out, mpbq const & a, bool decimal) {
411 if (a.m_k == 0) {
412 m_manager.display_smt2(out, a.m_num, decimal);
413 }
414 else {
415 out << "(/ ";
416 m_manager.display_smt2(out, a.m_num, decimal);
417 out << " ";
418 out << "(^ 2";
419 if (decimal) out << ".0";
420 out << " " << a.m_k;
421 if (decimal) out << ".0";
422 out << "))";
423 }
424 }
425
display_decimal(std::ostream & out,mpbq const & a,unsigned prec)426 void mpbq_manager::display_decimal(std::ostream & out, mpbq const & a, unsigned prec) {
427 if (is_int(a)) {
428 out << m_manager.to_string(a.m_num);
429 return;
430 }
431 mpz two(2);
432 mpz ten(10);
433 mpz two_k;
434 mpz n1, v1;
435 if (m_manager.is_neg(a.m_num))
436 out << "-";
437 m_manager.set(v1, a.m_num);
438 m_manager.abs(v1);
439 m_manager.power(two, a.m_k, two_k);
440 m_manager.rem(v1, two_k, n1);
441 m_manager.div(v1, two_k, v1);
442 SASSERT(!m_manager.is_zero(n1));
443 out << m_manager.to_string(v1);
444 out << ".";
445 for (unsigned i = 0; i < prec; i++) {
446 m_manager.mul(n1, ten, n1);
447 m_manager.div(n1, two_k, v1);
448 m_manager.rem(n1, two_k, n1);
449 out << m_manager.to_string(v1);
450 if (m_manager.is_zero(n1))
451 goto end;
452 }
453 out << "?";
454 end:
455 m_manager.del(n1);
456 m_manager.del(v1);
457 m_manager.del(two_k);
458 }
459
display_decimal(std::ostream & out,mpbq const & a,mpbq const & b,unsigned prec)460 void mpbq_manager::display_decimal(std::ostream & out, mpbq const & a, mpbq const & b, unsigned prec) {
461 mpz two(2);
462 mpz ten(10);
463 mpz two_k1, two_k2;
464 mpz n1, v1, n2, v2;
465 if (m_manager.is_neg(a.m_num) != m_manager.is_neg(b.m_num)) {
466 out << "?";
467 return;
468 }
469 if (m_manager.is_neg(a.m_num))
470 out << "-";
471 m_manager.set(v1, a.m_num);
472 m_manager.abs(v1);
473 m_manager.set(v2, b.m_num);
474 m_manager.abs(v2);
475 m_manager.power(two, a.m_k, two_k1);
476 m_manager.power(two, b.m_k, two_k2);
477 m_manager.rem(v1, two_k1, n1);
478 m_manager.rem(v2, two_k2, n2);
479 m_manager.div(v1, two_k1, v1);
480 m_manager.div(v2, two_k2, v2);
481 if (!m_manager.eq(v1, v2)) {
482 out << "?";
483 goto end;
484 }
485
486 out << m_manager.to_string(v1);
487 if (m_manager.is_zero(n1) && m_manager.is_zero(n2))
488 goto end; // number is an integer
489 out << ".";
490 for (unsigned i = 0; i < prec; i++) {
491 m_manager.mul(n1, ten, n1);
492 m_manager.mul(n2, ten, n2);
493 m_manager.div(n1, two_k1, v1);
494 m_manager.div(n2, two_k2, v2);
495 if (m_manager.eq(v1, v2)) {
496 out << m_manager.to_string(v1);
497 }
498 else {
499 out << "?";
500 goto end;
501 }
502 m_manager.rem(n1, two_k1, n1);
503 m_manager.rem(n2, two_k2, n2);
504 if (m_manager.is_zero(n1) && m_manager.is_zero(n2))
505 goto end; // number is precise
506 }
507 out << "?";
508 end:
509 m_manager.del(n1);
510 m_manager.del(v1);
511 m_manager.del(n2);
512 m_manager.del(v2);
513 m_manager.del(two_k1);
514 m_manager.del(two_k2);
515 }
516
to_mpbq(mpq const & q,mpbq & bq)517 bool mpbq_manager::to_mpbq(mpq const & q, mpbq & bq) {
518 mpz const & n = q.numerator();
519 mpz const & d = q.denominator();
520 unsigned shift;
521 if (m_manager.is_one(d)) {
522 set(bq, n);
523 SASSERT(eq(bq, q));
524 return true;
525 }
526 else if (m_manager.is_power_of_two(d, shift)) {
527 SASSERT(shift>=1);
528 unsigned k = shift;
529 set(bq, n, k);
530 SASSERT(eq(bq, q));
531 return true;
532 }
533 else {
534 unsigned k = m_manager.log2(d);
535 set(bq, n, k+1);
536 return false;
537 }
538 }
539
refine_upper(mpq const & q,mpbq & l,mpbq & u)540 void mpbq_manager::refine_upper(mpq const & q, mpbq & l, mpbq & u) {
541 SASSERT(lt(l, q) && gt(u, q));
542 SASSERT(!m_manager.is_power_of_two(q.denominator()));
543 // l < q < u
544 mpbq mid;
545 while (true) {
546 add(l, u, mid);
547 div2(mid);
548 if (gt(mid, q)) {
549 swap(u, mid);
550 del(mid);
551 SASSERT(lt(l, q) && gt(u, q));
552 return;
553 }
554 swap(l, mid);
555 }
556 }
557
refine_lower(mpq const & q,mpbq & l,mpbq & u)558 void mpbq_manager::refine_lower(mpq const & q, mpbq & l, mpbq & u) {
559 SASSERT(lt(l, q) && gt(u, q));
560 SASSERT(!m_manager.is_power_of_two(q.denominator()));
561 // l < q < u
562 mpbq mid;
563 while (true) {
564 add(l, u, mid);
565 div2(mid);
566 if (lt(mid, q)) {
567 swap(l, mid);
568 del(mid);
569 SASSERT(lt(l, q) && gt(u, q));
570 return;
571 }
572 swap(u, mid);
573 }
574 }
575
576 // sectect integer in [lower, upper]
select_integer(mpbq const & lower,mpbq const & upper,mpz & r)577 bool mpbq_manager::select_integer(mpbq const & lower, mpbq const & upper, mpz & r) {
578 if (is_int(lower)) {
579 m_manager.set(r, lower.m_num);
580 return true;
581 }
582 if (is_int(upper)) {
583 m_manager.set(r, upper.m_num);
584 return true;
585 }
586 mpz & ceil_lower = m_select_int_tmp1;
587 mpz & floor_upper = m_select_int_tmp2;
588 ceil(m_manager, lower, ceil_lower);
589 floor(m_manager, upper, floor_upper);
590 if (m_manager.le(ceil_lower, floor_upper)) {
591 m_manager.set(r, ceil_lower);
592 return true;
593 }
594 return false;
595 }
596
597 // select integer in (lower, upper]
select_integer(unsynch_mpq_manager & qm,mpq const & lower,mpbq const & upper,mpz & r)598 bool mpbq_manager::select_integer(unsynch_mpq_manager & qm, mpq const & lower, mpbq const & upper, mpz & r) {
599 if (is_int(upper)) {
600 m_manager.set(r, upper.m_num);
601 return true;
602 }
603 mpz & ceil_lower = m_select_int_tmp1;
604 mpz & floor_upper = m_select_int_tmp2;
605 if (qm.is_int(lower)) {
606 m_manager.set(ceil_lower, lower.numerator());
607 m_manager.inc(ceil_lower);
608 }
609 else {
610 scoped_mpz tmp(qm);
611 qm.ceil(lower, tmp);
612 m_manager.set(ceil_lower, tmp);
613 }
614 floor(m_manager, upper, floor_upper);
615 if (m_manager.le(ceil_lower, floor_upper)) {
616 m_manager.set(r, ceil_lower);
617 return true;
618 }
619 return false;
620 }
621
622 // sectect integer in [lower, upper)
select_integer(unsynch_mpq_manager & qm,mpbq const & lower,mpq const & upper,mpz & r)623 bool mpbq_manager::select_integer(unsynch_mpq_manager & qm, mpbq const & lower, mpq const & upper, mpz & r) {
624 if (is_int(lower)) {
625 m_manager.set(r, lower.m_num);
626 return true;
627 }
628 mpz & ceil_lower = m_select_int_tmp1;
629 mpz & floor_upper = m_select_int_tmp2;
630 ceil(m_manager, lower, ceil_lower);
631 if (qm.is_int(upper)) {
632 m_manager.set(floor_upper, upper.numerator());
633 m_manager.dec(floor_upper);
634 }
635 else {
636 scoped_mpz tmp(qm);
637 qm.floor(upper, tmp);
638 m_manager.set(floor_upper, tmp);
639 }
640 if (m_manager.le(ceil_lower, floor_upper)) {
641 m_manager.set(r, ceil_lower);
642 return true;
643 }
644 return false;
645 }
646
647 // sectect integer in (lower, upper)
select_integer(unsynch_mpq_manager & qm,mpq const & lower,mpq const & upper,mpz & r)648 bool mpbq_manager::select_integer(unsynch_mpq_manager & qm, mpq const & lower, mpq const & upper, mpz & r) {
649 mpz & ceil_lower = m_select_int_tmp1;
650 mpz & floor_upper = m_select_int_tmp2;
651
652 if (qm.is_int(lower)) {
653 m_manager.set(ceil_lower, lower.numerator());
654 m_manager.inc(ceil_lower);
655 }
656 else {
657 scoped_mpz tmp(qm);
658 qm.ceil(lower, tmp);
659 m_manager.set(ceil_lower, tmp);
660 }
661
662 if (qm.is_int(upper)) {
663 m_manager.set(floor_upper, upper.numerator());
664 m_manager.dec(floor_upper);
665 }
666 else {
667 scoped_mpz tmp(qm);
668 qm.floor(upper, tmp);
669 m_manager.set(floor_upper, tmp);
670 }
671
672 if (m_manager.le(ceil_lower, floor_upper)) {
673 m_manager.set(r, ceil_lower);
674 return true;
675 }
676 return false;
677 }
678
679 #define LINEAR_SEARCH_THRESHOLD 8
680
select_small_core(mpbq const & lower,mpbq const & upper,mpbq & r)681 void mpbq_manager::select_small_core(mpbq const & lower, mpbq const & upper, mpbq & r) {
682 SASSERT(le(lower, upper));
683 mpz & aux = m_select_small_tmp;
684 if (select_integer(lower, upper, aux)) {
685 set(r, aux);
686 return;
687 }
688
689 // At this point we know that k=0 does not work, since there is no integer
690 // in the interval [lower, upper]
691 unsigned min_k = 0;
692 unsigned max_k = std::min(lower.m_k, upper.m_k);
693
694 if (max_k <= LINEAR_SEARCH_THRESHOLD) {
695 unsigned k = 0;
696 mpbq & l2k = m_select_small_tmp1;
697 mpbq & u2k = m_select_small_tmp2;
698 set(l2k, lower);
699 set(u2k, upper);
700 while (true) {
701 k++;
702 mul2(l2k);
703 mul2(u2k);
704 if (select_integer(l2k, u2k, aux)) {
705 set(r, aux, k);
706 break;
707 }
708 }
709 }
710 else {
711 mpbq & l2k = m_select_small_tmp1;
712 mpbq & u2k = m_select_small_tmp2;
713 while (true) {
714 unsigned mid_k = min_k + (max_k - min_k)/2;
715 set(l2k, lower);
716 set(u2k, upper);
717 mul2k(l2k, mid_k);
718 mul2k(u2k, mid_k);
719 if (select_integer(l2k, u2k, aux))
720 max_k = mid_k;
721 else
722 min_k = mid_k + 1;
723 if (min_k == max_k) {
724 if (max_k == mid_k) {
725 set(r, aux, max_k);
726 }
727 else {
728 set(l2k, lower);
729 set(u2k, upper);
730 mul2k(l2k, max_k);
731 mul2k(u2k, max_k);
732 VERIFY(select_integer(l2k, u2k, aux));
733 set(r, aux, max_k);
734 }
735 break;
736 }
737 }
738 }
739 SASSERT(le(lower, r));
740 SASSERT(le(r, upper));
741 }
742
select_small(mpbq const & lower,mpbq const & upper,mpbq & r)743 bool mpbq_manager::select_small(mpbq const & lower, mpbq const & upper, mpbq & r) {
744 if (gt(lower, upper))
745 return false;
746 select_small_core(lower, upper, r);
747 return true;
748 }
749
750
select_small_core(unsynch_mpq_manager & qm,mpq const & lower,mpbq const & upper,mpbq & r)751 void mpbq_manager::select_small_core(unsynch_mpq_manager & qm, mpq const & lower, mpbq const & upper, mpbq & r) {
752 TRACE("select_small", tout << "lower (q): " << qm.to_string(lower) << ", upper (bq): " << to_string(upper) << "\n";);
753 SASSERT(gt(upper, lower));
754 mpz & aux = m_select_small_tmp;
755 if (select_integer(qm, lower, upper, aux)) {
756 set(r, aux);
757 return;
758 }
759
760 // At this point we know that k=0 does not work, since there is no integer
761 // in the interval [lower, upper]
762 unsigned k = 0;
763 scoped_mpq l2k(qm);
764 mpq two(2);
765 mpbq & u2k = m_select_small_tmp2;
766 qm.set(l2k, lower);
767 set(u2k, upper);
768 while (true) {
769 k++;
770 qm.mul(l2k, two, l2k);
771 mul2(u2k);
772 if (select_integer(qm, l2k, u2k, aux)) {
773 set(r, aux, k);
774 break;
775 }
776 }
777 }
778
select_small_core(unsynch_mpq_manager & qm,mpbq const & lower,mpq const & upper,mpbq & r)779 void mpbq_manager::select_small_core(unsynch_mpq_manager & qm, mpbq const & lower, mpq const & upper, mpbq & r) {
780 SASSERT(lt(lower, upper));
781 mpz & aux = m_select_small_tmp;
782 if (select_integer(qm, lower, upper, aux)) {
783 set(r, aux);
784 return;
785 }
786
787 // At this point we know that k=0 does not work, since there is no integer
788 // in the interval [lower, upper]
789 unsigned k = 0;
790 mpbq & l2k = m_select_small_tmp2;
791 scoped_mpq u2k(qm);
792 mpq two(2);
793 set(l2k, lower);
794 qm.set(u2k, upper);
795 while (true) {
796 k++;
797 mul2(l2k);
798 qm.mul(u2k, two, u2k);
799 if (select_integer(qm, l2k, u2k, aux)) {
800 set(r, aux, k);
801 break;
802 }
803 }
804 }
805
select_small_core(unsynch_mpq_manager & qm,mpq const & lower,mpq const & upper,mpbq & r)806 void mpbq_manager::select_small_core(unsynch_mpq_manager & qm, mpq const & lower, mpq const & upper, mpbq & r) {
807 SASSERT(qm.lt(lower, upper));
808 mpz & aux = m_select_small_tmp;
809 if (select_integer(qm, lower, upper, aux)) {
810 set(r, aux);
811 return;
812 }
813
814 // At this point we know that k=0 does not work, since there is no integer
815 // in the interval [lower, upper]
816 unsigned k = 0;
817 scoped_mpq l2k(qm);
818 scoped_mpq u2k(qm);
819 mpq two(2);
820 qm.set(l2k, lower);
821 qm.set(u2k, upper);
822 while (true) {
823 k++;
824 qm.mul(l2k, two, l2k);
825 qm.mul(u2k, two, u2k);
826 if (select_integer(qm, l2k, u2k, aux)) {
827 set(r, aux, k);
828 break;
829 }
830 }
831 }
832
approx(mpbq & a,unsigned k,bool to_plus_inf)833 void mpbq_manager::approx(mpbq & a, unsigned k, bool to_plus_inf) {
834 if (a.m_k <= k)
835 return;
836 #ifdef MPBQ_DEBUG
837 scoped_mpbq old_a(*this);
838 old_a = a;
839 #endif
840 bool sgn = m_manager.is_neg(a.m_num);
841 bool _inc = (sgn != to_plus_inf);
842 unsigned shift = a.m_k - k;
843 m_manager.abs(a.m_num);
844 m_manager.machine_div2k(a.m_num, shift);
845 if (_inc)
846 m_manager.inc(a.m_num);
847 if (sgn)
848 m_manager.neg(a.m_num);
849 a.m_k = k;
850 normalize(a);
851 #ifdef MPBQ_DEBUG
852 if (to_plus_inf) {
853 SASSERT(lt(old_a, a));
854 }
855 else {
856 SASSERT(lt(a, old_a));
857 }
858 #endif
859 }
860
approx_div(mpbq const & a,mpbq const & b,mpbq & c,unsigned k,bool to_plus_inf)861 void mpbq_manager::approx_div(mpbq const & a, mpbq const & b, mpbq & c, unsigned k, bool to_plus_inf) {
862 SASSERT(!is_zero(b));
863 unsigned k_prime;
864 if (m_manager.is_power_of_two(b.m_num, k_prime)) {
865 // The division is precise, so we ignore k and to_plus_inf
866 SASSERT(b.m_k == 0 || k_prime == 0); // remark: b.m_num is odd when b.m_k > 0, since b.m_num is a power of two we have that b.m_k == 0 or b.m_num == 1.
867 m_manager.set(c.m_num, a.m_num);
868 if (b.m_k > 0) {
869 SASSERT(k_prime == 0);
870 mpz & pw2 = m_div_tmp1;
871 m_manager.power(mpz(2), b.m_k, pw2);
872 m_manager.mul(c.m_num, pw2, c.m_num);
873 }
874 c.m_k = a.m_k + k_prime;
875 normalize(c);
876 }
877 else if (m_manager.divides(b.m_num, a.m_num)) {
878 // result is also precise
879 m_manager.div(a.m_num, b.m_num, c.m_num);
880 if (a.m_k >= b.m_k) {
881 c.m_k = a.m_k - b.m_k;
882 }
883 else {
884 m_manager.mul2k(c.m_num, b.m_k - a.m_k);
885 c.m_k = 0;
886 }
887 normalize(c);
888 }
889 else {
890 bool sgn = is_neg(a) != is_neg(b);
891 mpz & abs_a = m_div_tmp1;
892 mpz & norm_a = m_div_tmp2;
893 mpz & abs_b = m_div_tmp3;
894 m_manager.set(abs_a, a.m_num);
895 m_manager.abs(abs_a);
896 m_manager.set(abs_b, b.m_num);
897 m_manager.abs(abs_b);
898 if (a.m_k > b.m_k) {
899 if (k >= a.m_k - b.m_k)
900 m_manager.mul2k(abs_a, k - (a.m_k - b.m_k), norm_a);
901 else
902 m_manager.machine_div2k(abs_a, (a.m_k - b.m_k) - k, norm_a);
903 }
904 else {
905 m_manager.mul2k(abs_a, k + b.m_k - a.m_k, norm_a);
906 }
907 c.m_k = k;
908 m_manager.div(norm_a, abs_b, c.m_num);
909 if (sgn != to_plus_inf)
910 m_manager.inc(c.m_num);
911 if (sgn)
912 m_manager.neg(c.m_num);
913 normalize(c);
914 }
915 }
916