mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-15 13:28:47 +00:00
Code Activity

Fix memory management bugs

Browse source 
Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>
This commit is contained in:
Leonardo de Moura 2013-01-07 17:31:53 -08:00
parent 5873a59769
commit e01a7b6268
1 changed files with 329 additions and 225 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

554
src/math/realclosure/realclosure.cpp
View file

@ -1273,9 +1273,10 @@ namespace realclosure {
if (n == 2) { if (n == 2) {
// we don't need a field extension for linear polynomials. // we don't need a field extension for linear polynomials.
numeral r; numeral r;
value_ref neg_as_0(*this); value_ref v(*this);
neg_as_0 = neg(as[0]); neg(as[0], v);
set(r, div(neg_as_0, as[1])); div(v, as[1], v);
set(r, v);
roots.push_back(r); roots.push_back(r);
} }
else { else {
@ -1369,6 +1370,10 @@ namespace realclosure {
return is_real(a.m_value); return is_real(a.m_value);
} }
static void swap(mpbqi & a, mpbqi & b) {
realclosure::swap(a, b);
}
void mpq_to_mpbqi(mpq const & q, mpbqi & interval, unsigned k) { void mpq_to_mpbqi(mpq const & q, mpbqi & interval, unsigned k) {
interval.set_lower_is_inf(false); interval.set_lower_is_inf(false);
interval.set_upper_is_inf(false); interval.set_upper_is_inf(false);
@ -1489,7 +1494,9 @@ namespace realclosure {
// create the polynomial p of the form x^k - a // create the polynomial p of the form x^k - a
value_ref_buffer p(*this); value_ref_buffer p(*this);
p.push_back(neg(a.m_value)); value_ref neg_a(*this);
neg(a.m_value, neg_a);
p.push_back(neg_a);
for (unsigned i = 0; i < k - 1; i++) for (unsigned i = 0; i < k - 1; i++)
p.push_back(0); p.push_back(0);
p.push_back(one()); p.push_back(one());
@ -1500,15 +1507,17 @@ namespace realclosure {
void power(numeral const & a, unsigned k, numeral & b) { void power(numeral const & a, unsigned k, numeral & b) {
unsigned mask = 1; unsigned mask = 1;
value_ref power(*this); value_ref power(*this);
value_ref _b(*this);
power = a.m_value; power = a.m_value;
set(b, one()); _b = one();
while (mask <= k) { while (mask <= k) {
checkpoint(); checkpoint();
if (mask & k) if (mask & k)
set(b, mul(b.m_value, power)); mul(_b, power, _b);
power = mul(power, power); mul(power, power, power);
mask = mask << 1; mask = mask << 1;
} }
set(b, _b);
} }
/** /**
@ -1525,10 +1534,13 @@ namespace realclosure {
*/ */
void add(unsigned sz1, value * const * p1, unsigned sz2, value * const * p2, value_ref_buffer & r) { void add(unsigned sz1, value * const * p1, unsigned sz2, value * const * p2, value_ref_buffer & r) {
r.reset(); r.reset();
value_ref a_i(*this);
unsigned min = std::min(sz1, sz2); unsigned min = std::min(sz1, sz2);
unsigned i = 0; unsigned i = 0;
for (; i < min; i++) for (; i < min; i++) {
r.push_back(add(p1[i], p2[i])); add(p1[i], p2[i], a_i);
r.push_back(a_i);
}
for (; i < sz1; i++) for (; i < sz1; i++)
r.push_back(p1[i]); r.push_back(p1[i]);
for (; i < sz2; i++) for (; i < sz2; i++)
@ -1543,7 +1555,9 @@ namespace realclosure {
void add(unsigned sz, value * const * p, value * a, value_ref_buffer & r) { void add(unsigned sz, value * const * p, value * a, value_ref_buffer & r) {
SASSERT(sz > 0); SASSERT(sz > 0);
r.reset(); r.reset();
r.push_back(add(p[0], a)); value_ref a_0(*this);
add(p[0], a, a_0);
r.push_back(a_0);
r.append(sz - 1, p + 1); r.append(sz - 1, p + 1);
adjust_size(r); adjust_size(r);
} }
@ -1553,14 +1567,19 @@ namespace realclosure {
*/ */
void sub(unsigned sz1, value * const * p1, unsigned sz2, value * const * p2, value_ref_buffer & r) { void sub(unsigned sz1, value * const * p1, unsigned sz2, value * const * p2, value_ref_buffer & r) {
r.reset(); r.reset();
value_ref a_i(*this);
unsigned min = std::min(sz1, sz2); unsigned min = std::min(sz1, sz2);
unsigned i = 0; unsigned i = 0;
for (; i < min; i++) for (; i < min; i++) {
r.push_back(sub(p1[i], p2[i])); sub(p1[i], p2[i], a_i);
r.push_back(a_i);
}
for (; i < sz1; i++) for (; i < sz1; i++)
r.push_back(p1[i]); r.push_back(p1[i]);
for (; i < sz2; i++) for (; i < sz2; i++) {
r.push_back(neg(p2[i])); neg(p2[i], a_i);
r.push_back(a_i);
}
SASSERT(r.size() == std::max(sz1, sz2)); SASSERT(r.size() == std::max(sz1, sz2));
adjust_size(r); adjust_size(r);
} }
@ -1571,7 +1590,9 @@ namespace realclosure {
void sub(unsigned sz, value * const * p, value * a, value_ref_buffer & r) { void sub(unsigned sz, value * const * p, value * a, value_ref_buffer & r) {
SASSERT(sz > 0); SASSERT(sz > 0);
r.reset(); r.reset();
r.push_back(sub(p[0], a)); value_ref a_0(*this);
sub(p[0], a, a_0);
r.push_back(a_0);
r.append(sz - 1, p + 1); r.append(sz - 1, p + 1);
adjust_size(r); adjust_size(r);
} }
@ -1583,8 +1604,11 @@ namespace realclosure {
r.reset(); r.reset();
if (a == 0) if (a == 0)
return; return;
for (unsigned i = 0; i < sz; i++) value_ref a_i(*this);
r.push_back(mul(a, p[i])); for (unsigned i = 0; i < sz; i++) {
mul(a, p[i], a_i);
r.push_back(a_i);
}
} }
/** /**
@ -1605,8 +1629,9 @@ namespace realclosure {
continue; continue;
for (unsigned j = 0; j < sz2; j++) { for (unsigned j = 0; j < sz2; j++) {
// r[i+j] <- r[i+j] + p1[i]*p2[j] // r[i+j] <- r[i+j] + p1[i]*p2[j]
tmp = mul(p1[i], p2[j]); mul(p1[i], p2[j], tmp);
r.set(i+j, add(r[i+j], tmp)); add(r[i+j], tmp, tmp);
r.set(i+j, tmp);
} }
} }
adjust_size(r); adjust_size(r);
@ -1619,9 +1644,12 @@ namespace realclosure {
SASSERT(!is_zero(a)); SASSERT(!is_zero(a));
if (is_rational_one(a)) if (is_rational_one(a))
return; return;
value_ref a_i(*this);
unsigned sz = p.size(); unsigned sz = p.size();
for (unsigned i = 0; i < sz; i++) for (unsigned i = 0; i < sz; i++) {
p.set(i, div(p[i], a)); div(p[i], a, a_i);
p.set(i, a_i);
}
} }
/** /**
@ -1648,6 +1676,7 @@ namespace realclosure {
value * b_n = p2[sz2-1]; value * b_n = p2[sz2-1];
SASSERT(!is_zero(b_n)); SASSERT(!is_zero(b_n));
value_ref ratio(*this); value_ref ratio(*this);
value_ref aux(*this);
while (true) { while (true) {
checkpoint(); checkpoint();
sz1 = r.size(); sz1 = r.size();
@ -1656,13 +1685,15 @@ namespace realclosure {
break; break;
} }
unsigned m_n = sz1 - sz2; // m-n unsigned m_n = sz1 - sz2; // m-n
ratio = div(r[sz1 - 1], b_n); div(r[sz1 - 1], b_n, ratio);
// q[m_n] <- q[m_n] + r[sz1 - 1]/b_n // q[m_n] <- q[m_n] + r[sz1 - 1]/b_n
q.set(m_n, add(q[m_n], ratio)); add(q[m_n], ratio, aux);
q.set(m_n, aux);
for (unsigned i = 0; i < sz2 - 1; i++) { for (unsigned i = 0; i < sz2 - 1; i++) {
// r[i + m_n] <- r[i + m_n] - ratio * p2[i] // r[i + m_n] <- r[i + m_n] - ratio * p2[i]
ratio = mul(ratio, p2[i]); mul(ratio, p2[i], aux);
r.set(i + m_n, sub(r[i + m_n], ratio)); sub(r[i + m_n], aux, aux);
r.set(i + m_n, aux);
} }
r.shrink(sz1 - 1); r.shrink(sz1 - 1);
adjust_size(r); adjust_size(r);
@ -1683,8 +1714,10 @@ namespace realclosure {
\brief r <- p/a \brief r <- p/a
*/ */
void div(unsigned sz, value * const * p, value * a, value_ref_buffer & r) { void div(unsigned sz, value * const * p, value * a, value_ref_buffer & r) {
value_ref a_i(*this);
for (unsigned i = 0; i < sz; i++) { for (unsigned i = 0; i < sz; i++) {
r.push_back(div(p[i], a)); div(p[i], a, a_i);
r.push_back(a_i);
} }
} }
@ -1715,10 +1748,11 @@ namespace realclosure {
return; return;
} }
unsigned m_n = sz1 - sz2; unsigned m_n = sz1 - sz2;
ratio = div(r[sz1 - 1], b_n); div(r[sz1 - 1], b_n, ratio);
for (unsigned i = 0; i < sz2 - 1; i++) { for (unsigned i = 0; i < sz2 - 1; i++) {
new_a = mul(ratio, p2[i]); mul(ratio, p2[i], new_a);
r.set(i + m_n, sub(r[i + m_n], new_a)); sub(r[i + m_n], new_a, new_a);
r.set(i + m_n, new_a);
} }
r.shrink(sz1 - 1); r.shrink(sz1 - 1);
adjust_size(r); adjust_size(r);
@ -1730,29 +1764,36 @@ namespace realclosure {
*/ */
void neg(unsigned sz, value * const * p, value_ref_buffer & r) { void neg(unsigned sz, value * const * p, value_ref_buffer & r) {
r.reset(); r.reset();
for (unsigned i = 0; i < sz; i++) value_ref a_i(*this);
r.push_back(neg(p[i])); for (unsigned i = 0; i < sz; i++) {
neg(p[i], a_i);
r.push_back(a_i);
}
} }
/** /**
\brief r <- -r \brief r <- -r
*/ */
void neg(value_ref_buffer & r) { void neg(value_ref_buffer & r) {
value_ref a_i(*this);
unsigned sz = r.size(); unsigned sz = r.size();
for (unsigned i = 0; i < sz; i++) for (unsigned i = 0; i < sz; i++) {
r.set(i, neg(r[i])); neg(r[i], a_i);
r.set(i, a_i);
}
} }
/** /**
\brief p <- -p \brief p <- -p
*/ */
void neg(polynomial & p) { void neg(polynomial & p) {
value_ref a_i(*this);
unsigned sz = p.size(); unsigned sz = p.size();
for (unsigned i = 0; i < sz; i++) { for (unsigned i = 0; i < sz; i++) {
value * v = neg(p[i]); neg(p[i], a_i);
inc_ref(v); inc_ref(a_i.get());
dec_ref(p[i]); dec_ref(p[i]);
p[i] = v; p[i] = a_i.get();
} }
} }
@ -1771,10 +1812,12 @@ namespace realclosure {
void mk_monic(value_ref_buffer & p) { void mk_monic(value_ref_buffer & p) {
unsigned sz = p.size(); unsigned sz = p.size();
if (sz > 0) { if (sz > 0) {
value_ref a_i(*this);
SASSERT(p[sz-1] != 0); SASSERT(p[sz-1] != 0);
if (!is_rational_one(p[sz-1])) { if (!is_rational_one(p[sz-1])) {
for (unsigned i = 0; i < sz - 1; i++) { for (unsigned i = 0; i < sz - 1; i++) {
p.set(i, div(p[i], p[sz-1])); div(p[i], p[sz-1], a_i);
p.set(i, a_i);
} }
p.set(sz-1, one()); p.set(sz-1, one());
} }
@ -1825,9 +1868,10 @@ namespace realclosure {
if (sz > 1) { if (sz > 1) {
for (unsigned i = 1; i < sz; i++) { for (unsigned i = 1; i < sz; i++) {
mpq i_mpq(i); mpq i_mpq(i);
value_ref i_value(*this); value_ref a_i(*this);
i_value = mk_rational(i_mpq); a_i = mk_rational(i_mpq);
r.push_back(mul(i_value, p[i])); mul(a_i, p[i], a_i);
r.push_back(a_i);
} }
adjust_size(r); adjust_size(r);
} }
@ -2297,6 +2341,11 @@ namespace realclosure {
} }
} }
bool determine_sign(value_ref & r) {
SASSERT(is_rational_function(r.get()));
return determine_sign(to_rational_function(r.get()));
}
/** /**
\brief Set new_p1 and new_p2 using the following normalization rules: \brief Set new_p1 and new_p2 using the following normalization rules:
- new_p1 <- p1/p2[0]; new_p2 <- one IF sz2 == 1 - new_p1 <- p1/p2[0]; new_p2 <- one IF sz2 == 1
@ -2357,30 +2406,32 @@ namespace realclosure {
\brief Create a new value using the a->ext(), and the given numerator and denominator. \brief Create a new value using the a->ext(), and the given numerator and denominator.
Use interval(a) + interval(b) as an initial approximation for the interval of the result, and invoke determine_sign() Use interval(a) + interval(b) as an initial approximation for the interval of the result, and invoke determine_sign()
*/ */
value * mk_add_value(rational_function_value * a, value * b, unsigned num_sz, value * const * num, unsigned den_sz, value * const * den) { void mk_add_value(rational_function_value * a, value * b, unsigned num_sz, value * const * num, unsigned den_sz, value * const * den, value_ref & r) {
SASSERT(num_sz > 0 && den_sz > 0); SASSERT(num_sz > 0 && den_sz > 0);
if (num_sz == 1 && den_sz == 1) { if (num_sz == 1 && den_sz == 1) {
// In this case, the normalization rules guarantee that den is one. // In this case, the normalization rules guarantee that den is one.
SASSERT(is_rational_one(den[0])); SASSERT(is_rational_one(den[0]));
return num[0]; r = num[0];
}
rational_function_value * r = mk_rational_function_value_core(a->ext(), num_sz, num, den_sz, den);
bqim().add(interval(a), interval(b), r->interval());
if (determine_sign(r)) {
SASSERT(!contains_zero(r->interval()));
return r;
} }
else { else {
// The new value is 0 scoped_mpbqi ri(bqim());
del_rational_function(r); bqim().add(interval(a), interval(b), ri);
return 0; r = mk_rational_function_value_core(a->ext(), num_sz, num, den_sz, den);
swap(r->interval(), ri);
if (determine_sign(r)) {
SASSERT(!contains_zero(r->interval()));
}
else {
// The new value is 0
r = 0;
}
} }
} }
/** /**
\brief Add a value of 'a' the form n/1 with b where rank(a) > rank(b) \brief Add a value of 'a' the form n/1 with b where rank(a) > rank(b)
*/ */
value * add_p_v(rational_function_value * a, value * b) { void add_p_v(rational_function_value * a, value * b, value_ref & r) {
SASSERT(is_rational_one(a->den())); SASSERT(is_rational_one(a->den()));
SASSERT(compare_rank(a, b) > 0); SASSERT(compare_rank(a, b) > 0);
polynomial const & an = a->num(); polynomial const & an = a->num();
@ -2389,35 +2440,40 @@ namespace realclosure {
value_ref_buffer new_num(*this); value_ref_buffer new_num(*this);
add(an.size(), an.c_ptr(), b, new_num); add(an.size(), an.c_ptr(), b, new_num);
SASSERT(new_num.size() == an.size()); SASSERT(new_num.size() == an.size());
return mk_add_value(a, b, new_num.size(), new_num.c_ptr(), one.size(), one.c_ptr()); mk_add_value(a, b, new_num.size(), new_num.c_ptr(), one.size(), one.c_ptr(), r);
} }
/** /**
\brief Add a value 'a' of the form n/d with b where rank(a) > rank(b) \brief Add a value 'a' of the form n/d with b where rank(a) > rank(b)
*/ */
value * add_rf_v(rational_function_value * a, value * b) { void add_rf_v(rational_function_value * a, value * b, value_ref & r) {
value_ref_buffer b_ad(*this); value_ref_buffer b_ad(*this);
value_ref_buffer num(*this); value_ref_buffer num(*this);
polynomial const & an = a->num(); polynomial const & an = a->num();
polynomial const & ad = a->den(); polynomial const & ad = a->den();
if (is_rational_one(ad)) if (is_rational_one(ad)) {
return add_p_v(a, b); add_p_v(a, b, r);
// b_ad <- b * ad }
mul(b, ad.size(), ad.c_ptr(), b_ad); else {
// num <- a + b * ad // b_ad <- b * ad
add(an.size(), an.c_ptr(), b_ad.size(), b_ad.c_ptr(), num); mul(b, ad.size(), ad.c_ptr(), b_ad);
if (num.empty()) // num <- a + b * ad
return 0; add(an.size(), an.c_ptr(), b_ad.size(), b_ad.c_ptr(), num);
value_ref_buffer new_num(*this); if (num.empty())
value_ref_buffer new_den(*this); r = 0;
normalize(num.size(), num.c_ptr(), ad.size(), ad.c_ptr(), new_num, new_den); else {
return mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr()); value_ref_buffer new_num(*this);
value_ref_buffer new_den(*this);
normalize(num.size(), num.c_ptr(), ad.size(), ad.c_ptr(), new_num, new_den);
mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r);
}
}
} }
/** /**
\brief Add values 'a' and 'b' of the form n/1 and rank(a) == rank(b) \brief Add values 'a' and 'b' of the form n/1 and rank(a) == rank(b)
*/ */
value * add_p_p(rational_function_value * a, rational_function_value * b) { void add_p_p(rational_function_value * a, rational_function_value * b, value_ref & r) {
SASSERT(is_rational_one(a->den())); SASSERT(is_rational_one(a->den()));
SASSERT(is_rational_one(b->den())); SASSERT(is_rational_one(b->den()));
SASSERT(compare_rank(a, b) == 0); SASSERT(compare_rank(a, b) == 0);
@ -2427,109 +2483,120 @@ namespace realclosure {
value_ref_buffer new_num(*this); value_ref_buffer new_num(*this);
add(an.size(), an.c_ptr(), bn.size(), bn.c_ptr(), new_num); add(an.size(), an.c_ptr(), bn.size(), bn.c_ptr(), new_num);
if (new_num.empty()) if (new_num.empty())
return 0; r = 0;
return mk_add_value(a, b, new_num.size(), new_num.c_ptr(), one.size(), one.c_ptr()); else
mk_add_value(a, b, new_num.size(), new_num.c_ptr(), one.size(), one.c_ptr(), r);
} }
/** /**
\brief Add values 'a' and 'b' of the form n/d and rank(a) == rank(b) \brief Add values 'a' and 'b' of the form n/d and rank(a) == rank(b)
*/ */
value * add_rf_rf(rational_function_value * a, rational_function_value * b) { void add_rf_rf(rational_function_value * a, rational_function_value * b, value_ref & r) {
SASSERT(compare_rank(a, b) == 0); SASSERT(compare_rank(a, b) == 0);
polynomial const & an = a->num(); polynomial const & an = a->num();
polynomial const & ad = a->den(); polynomial const & ad = a->den();
polynomial const & bn = b->num(); polynomial const & bn = b->num();
polynomial const & bd = b->den(); polynomial const & bd = b->den();
if (is_rational_one(ad) && is_rational_one(bd)) if (is_rational_one(ad) && is_rational_one(bd)) {
return add_p_p(a, b); add_p_p(a, b, r);
value_ref_buffer an_bd(*this);
value_ref_buffer bn_ad(*this);
mul(an.size(), an.c_ptr(), bd.size(), bd.c_ptr(), an_bd);
mul(bn.size(), bn.c_ptr(), ad.size(), ad.c_ptr(), bn_ad);
value_ref_buffer num(*this);
add(an_bd.size(), an_bd.c_ptr(), bn_ad.size(), bn_ad.c_ptr(), num);
if (num.empty())
return 0;
value_ref_buffer den(*this);
mul(ad.size(), ad.c_ptr(), bd.size(), bd.c_ptr(), den);
value_ref_buffer new_num(*this);
value_ref_buffer new_den(*this);
normalize(num.size(), num.c_ptr(), den.size(), den.c_ptr(), new_num, new_den);
return mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr());
}
value * add(value * a, value * b) {
if (a == 0)
return b;
else if (b == 0)
return a;
else if (is_nz_rational(a) && is_nz_rational(b)) {
scoped_mpq r(qm());
qm().add(to_mpq(a), to_mpq(b), r);
if (qm().is_zero(r))
return 0;
else
return mk_rational(r);
} }
else { else {
switch (compare_rank(a, b)) { value_ref_buffer an_bd(*this);
case -1: return add_rf_v(to_rational_function(b), a); value_ref_buffer bn_ad(*this);
case 0: return add_rf_rf(to_rational_function(a), to_rational_function(b)); mul(an.size(), an.c_ptr(), bd.size(), bd.c_ptr(), an_bd);
case 1: return add_rf_v(to_rational_function(a), b); mul(bn.size(), bn.c_ptr(), ad.size(), ad.c_ptr(), bn_ad);
default: UNREACHABLE(); value_ref_buffer num(*this);
return 0; add(an_bd.size(), an_bd.c_ptr(), bn_ad.size(), bn_ad.c_ptr(), num);
if (num.empty()) {
r = 0;
}
else {
value_ref_buffer den(*this);
mul(ad.size(), ad.c_ptr(), bd.size(), bd.c_ptr(), den);
value_ref_buffer new_num(*this);
value_ref_buffer new_den(*this);
normalize(num.size(), num.c_ptr(), den.size(), den.c_ptr(), new_num, new_den);
mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r);
} }
} }
} }
value * sub(value * a, value * b) { void add(value * a, value * b, value_ref & r) {
if (a == 0) if (a == 0) {
return neg(b); r = b;
else if (b == 0) }
return a; else if (b == 0) {
r = a;
}
else if (is_nz_rational(a) && is_nz_rational(b)) { else if (is_nz_rational(a) && is_nz_rational(b)) {
scoped_mpq r(qm()); scoped_mpq v(qm());
qm().sub(to_mpq(a), to_mpq(b), r); qm().add(to_mpq(a), to_mpq(b), v);
if (qm().is_zero(r)) if (qm().is_zero(v))
return 0; r = 0;
else else
return mk_rational(r); r = mk_rational(v);
}
else {
switch (compare_rank(a, b)) {
case -1: add_rf_v(to_rational_function(b), a, r); break;
case 0: add_rf_rf(to_rational_function(a), to_rational_function(b), r); break;
case 1: add_rf_v(to_rational_function(a), b, r); break;
default: UNREACHABLE();
}
}
}
void sub(value * a, value * b, value_ref & r) {
if (a == 0) {
neg(b, r);
}
else if (b == 0) {
r = a;
}
else if (is_nz_rational(a) && is_nz_rational(b)) {
scoped_mpq v(qm());
qm().sub(to_mpq(a), to_mpq(b), v);
if (qm().is_zero(v))
r = 0;
else
r = mk_rational(v);
} }
else { else {
value_ref neg_b(*this); value_ref neg_b(*this);
neg_b = neg(b); neg(b, neg_b);
switch (compare_rank(a, neg_b)) { switch (compare_rank(a, neg_b)) {
case -1: return add_rf_v(to_rational_function(neg_b), a); case -1: add_rf_v(to_rational_function(neg_b), a, r); break;
case 0: return add_rf_rf(to_rational_function(a), to_rational_function(neg_b)); case 0: add_rf_rf(to_rational_function(a), to_rational_function(neg_b), r); break;
case 1: return add_rf_v(to_rational_function(a), b); case 1: add_rf_v(to_rational_function(a), neg_b, r); break;
default: UNREACHABLE(); default: UNREACHABLE();
return 0;
} }
} }
} }
value * neg_rf(rational_function_value * a) { void neg_rf(rational_function_value * a, value_ref & r) {
polynomial const & an = a->num(); polynomial const & an = a->num();
polynomial const & ad = a->den(); polynomial const & ad = a->den();
value_ref_buffer new_num(*this); value_ref_buffer new_num(*this);
neg(an.size(), an.c_ptr(), new_num); neg(an.size(), an.c_ptr(), new_num);
rational_function_value * r = mk_rational_function_value_core(a->ext(), new_num.size(), new_num.c_ptr(), ad.size(), ad.c_ptr()); scoped_mpbqi ri(bqim());
bqim().neg(interval(a), r->interval()); bqim().neg(interval(a), ri);
r = mk_rational_function_value_core(a->ext(), new_num.size(), new_num.c_ptr(), ad.size(), ad.c_ptr());
swap(r->interval(), ri);
SASSERT(!contains_zero(r->interval())); SASSERT(!contains_zero(r->interval()));
return r;
} }
value * neg(value * a) { void neg(value * a, value_ref & r) {
if (a == 0) if (a == 0) {
return 0; r = 0;
}
else if (is_nz_rational(a)) { else if (is_nz_rational(a)) {
scoped_mpq r(qm()); scoped_mpq v(qm());
qm().set(r, to_mpq(a)); qm().set(v, to_mpq(a));
qm().neg(r); qm().neg(v);
return mk_rational(r); r = mk_rational(v);
} }
else { else {
return neg_rf(to_rational_function(a)); neg_rf(to_rational_function(a), r);
} }
} }
@ -2537,30 +2604,32 @@ namespace realclosure {
\brief Create a new value using the a->ext(), and the given numerator and denominator. \brief Create a new value using the a->ext(), and the given numerator and denominator.
Use interval(a) * interval(b) as an initial approximation for the interval of the result, and invoke determine_sign() Use interval(a) * interval(b) as an initial approximation for the interval of the result, and invoke determine_sign()
*/ */
value * mk_mul_value(rational_function_value * a, value * b, unsigned num_sz, value * const * num, unsigned den_sz, value * const * den) { void mk_mul_value(rational_function_value * a, value * b, unsigned num_sz, value * const * num, unsigned den_sz, value * const * den, value_ref & r) {
SASSERT(num_sz > 0 && den_sz > 0); SASSERT(num_sz > 0 && den_sz > 0);
if (num_sz == 1 && den_sz == 1) { if (num_sz == 1 && den_sz == 1) {
// In this case, the normalization rules guarantee that den is one. // In this case, the normalization rules guarantee that den is one.
SASSERT(is_rational_one(den[0])); SASSERT(is_rational_one(den[0]));
return num[0]; r = num[0];
}
rational_function_value * r = mk_rational_function_value_core(a->ext(), num_sz, num, den_sz, den);
bqim().mul(interval(a), interval(b), r->interval());
if (determine_sign(r)) {
SASSERT(!contains_zero(r->interval()));
return r;
} }
else { else {
// The new value is 0 scoped_mpbqi ri(bqim());
del_rational_function(r); bqim().mul(interval(a), interval(b), ri);
return 0; r = mk_rational_function_value_core(a->ext(), num_sz, num, den_sz, den);
swap(ri, r->interval());
if (determine_sign(r)) {
SASSERT(!contains_zero(r->interval()));
}
else {
// The new value is 0
r = 0;
}
} }
} }
/** /**
\brief Multiply a value of 'a' the form n/1 with b where rank(a) > rank(b) \brief Multiply a value of 'a' the form n/1 with b where rank(a) > rank(b)
*/ */
value * mul_p_v(rational_function_value * a, value * b) { void mul_p_v(rational_function_value * a, value * b, value_ref & r) {
SASSERT(is_rational_one(a->den())); SASSERT(is_rational_one(a->den()));
SASSERT(b != 0); SASSERT(b != 0);
SASSERT(compare_rank(a, b) > 0); SASSERT(compare_rank(a, b) > 0);
@ -2570,31 +2639,34 @@ namespace realclosure {
value_ref_buffer new_num(*this); value_ref_buffer new_num(*this);
mul(b, an.size(), an.c_ptr(), new_num); mul(b, an.size(), an.c_ptr(), new_num);
SASSERT(new_num.size() == an.size()); SASSERT(new_num.size() == an.size());
return mk_mul_value(a, b, new_num.size(), new_num.c_ptr(), one.size(), one.c_ptr()); mk_mul_value(a, b, new_num.size(), new_num.c_ptr(), one.size(), one.c_ptr(), r);
} }
/** /**
\brief Multiply a value 'a' of the form n/d with b where rank(a) > rank(b) \brief Multiply a value 'a' of the form n/d with b where rank(a) > rank(b)
*/ */
value * mul_rf_v(rational_function_value * a, value * b) { void mul_rf_v(rational_function_value * a, value * b, value_ref & r) {
polynomial const & an = a->num(); polynomial const & an = a->num();
polynomial const & ad = a->den(); polynomial const & ad = a->den();
if (is_rational_one(ad)) if (is_rational_one(ad)) {
return mul_p_v(a, b); mul_p_v(a, b, r);
value_ref_buffer num(*this); }
// num <- b * an else {
mul(b, an.size(), an.c_ptr(), num); value_ref_buffer num(*this);
SASSERT(num.size() == an.size()); // num <- b * an
value_ref_buffer new_num(*this); mul(b, an.size(), an.c_ptr(), num);
value_ref_buffer new_den(*this); SASSERT(num.size() == an.size());
normalize(num.size(), num.c_ptr(), ad.size(), ad.c_ptr(), new_num, new_den); value_ref_buffer new_num(*this);
return mk_mul_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr()); value_ref_buffer new_den(*this);
normalize(num.size(), num.c_ptr(), ad.size(), ad.c_ptr(), new_num, new_den);
mk_mul_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r);
}
} }
/** /**
\brief Multiply values 'a' and 'b' of the form n/1 and rank(a) == rank(b) \brief Multiply values 'a' and 'b' of the form n/1 and rank(a) == rank(b)
*/ */
value * mul_p_p(rational_function_value * a, rational_function_value * b) { void mul_p_p(rational_function_value * a, rational_function_value * b, value_ref & r) {
SASSERT(is_rational_one(a->den())); SASSERT(is_rational_one(a->den()));
SASSERT(is_rational_one(b->den())); SASSERT(is_rational_one(b->den()));
SASSERT(compare_rank(a, b) == 0); SASSERT(compare_rank(a, b) == 0);
@ -2604,107 +2676,119 @@ namespace realclosure {
value_ref_buffer new_num(*this); value_ref_buffer new_num(*this);
mul(an.size(), an.c_ptr(), bn.size(), bn.c_ptr(), new_num); mul(an.size(), an.c_ptr(), bn.size(), bn.c_ptr(), new_num);
SASSERT(!new_num.empty()); SASSERT(!new_num.empty());
return mk_mul_value(a, b, new_num.size(), new_num.c_ptr(), one.size(), one.c_ptr()); mk_mul_value(a, b, new_num.size(), new_num.c_ptr(), one.size(), one.c_ptr(), r);
} }
/** /**
\brief Multiply values 'a' and 'b' of the form n/d and rank(a) == rank(b) \brief Multiply values 'a' and 'b' of the form n/d and rank(a) == rank(b)
*/ */
value * mul_rf_rf(rational_function_value * a, rational_function_value * b) { void mul_rf_rf(rational_function_value * a, rational_function_value * b, value_ref & r) {
SASSERT(compare_rank(a, b) == 0); SASSERT(compare_rank(a, b) == 0);
polynomial const & an = a->num(); polynomial const & an = a->num();
polynomial const & ad = a->den(); polynomial const & ad = a->den();
polynomial const & bn = b->num(); polynomial const & bn = b->num();
polynomial const & bd = b->den(); polynomial const & bd = b->den();
if (is_rational_one(ad) && is_rational_one(bd)) if (is_rational_one(ad) && is_rational_one(bd)) {
return mul_p_p(a, b); mul_p_p(a, b, r);
value_ref_buffer num(*this); }
value_ref_buffer den(*this); else {
mul(an.size(), an.c_ptr(), bn.size(), bn.c_ptr(), num); value_ref_buffer num(*this);
mul(ad.size(), ad.c_ptr(), bd.size(), bd.c_ptr(), den); value_ref_buffer den(*this);
SASSERT(!num.empty()); SASSERT(!den.empty()); mul(an.size(), an.c_ptr(), bn.size(), bn.c_ptr(), num);
value_ref_buffer new_num(*this); mul(ad.size(), ad.c_ptr(), bd.size(), bd.c_ptr(), den);
value_ref_buffer new_den(*this); SASSERT(!num.empty()); SASSERT(!den.empty());
normalize(num.size(), num.c_ptr(), den.size(), den.c_ptr(), new_num, new_den); value_ref_buffer new_num(*this);
return mk_mul_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr()); value_ref_buffer new_den(*this);
normalize(num.size(), num.c_ptr(), den.size(), den.c_ptr(), new_num, new_den);
mk_mul_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r);
}
} }
value * mul(value * a, value * b) { void mul(value * a, value * b, value_ref & r) {
if (a == 0 || b == 0) if (a == 0 || b == 0) {
return 0; r = 0;
else if (is_rational_one(a)) }
return b; else if (is_rational_one(a)) {
else if (is_rational_one(b)) r = b;
return a; }
else if (is_rational_minus_one(a)) else if (is_rational_one(b)) {
return neg(b); r = a;
else if (is_rational_minus_one(b)) }
return neg(a); else if (is_rational_minus_one(a)) {
neg(b, r);
}
else if (is_rational_minus_one(b)) {
neg(a, r);
}
else if (is_nz_rational(a) && is_nz_rational(b)) { else if (is_nz_rational(a) && is_nz_rational(b)) {
scoped_mpq r(qm()); scoped_mpq v(qm());
qm().mul(to_mpq(a), to_mpq(b), r); qm().mul(to_mpq(a), to_mpq(b), v);
return mk_rational(r); r = mk_rational(v);
} }
else { else {
switch (compare_rank(a, b)) { switch (compare_rank(a, b)) {
case -1: return mul_rf_v(to_rational_function(b), a); case -1: mul_rf_v(to_rational_function(b), a, r); break;
case 0: return mul_rf_rf(to_rational_function(a), to_rational_function(b)); case 0: mul_rf_rf(to_rational_function(a), to_rational_function(b), r); break;
case 1: return mul_rf_v(to_rational_function(a), b); case 1: mul_rf_v(to_rational_function(a), b, r); break;
default: UNREACHABLE(); default: UNREACHABLE();
return 0;
} }
} }
} }
value * div(value * a, value * b) { void div(value * a, value * b, value_ref & r) {
if (a == 0) if (a == 0) {
return 0; r = 0;
else if (b == 0) }
else if (b == 0) {
throw exception("division by zero"); throw exception("division by zero");
else if (is_rational_one(b)) }
return a; else if (is_rational_one(b)) {
else if (is_rational_one(a)) r = a;
return inv(b); }
else if (is_rational_minus_one(b)) else if (is_rational_one(a)) {
return neg(a); inv(b, r);
}
else if (is_rational_minus_one(b)) {
neg(a, r);
}
else if (is_nz_rational(a) && is_nz_rational(b)) { else if (is_nz_rational(a) && is_nz_rational(b)) {
scoped_mpq r(qm()); scoped_mpq v(qm());
qm().div(to_mpq(a), to_mpq(b), r); qm().div(to_mpq(a), to_mpq(b), v);
return mk_rational(r); r = mk_rational(v);
} }
else { else {
value_ref inv_b(*this); value_ref inv_b(*this);
inv_b = inv(b); inv(b, inv_b);
switch (compare_rank(a, inv_b)) { switch (compare_rank(a, inv_b)) {
case -1: return mul_rf_v(to_rational_function(inv_b), a); case -1: mul_rf_v(to_rational_function(inv_b), a, r); break;
case 0: return mul_rf_rf(to_rational_function(a), to_rational_function(inv_b)); case 0: mul_rf_rf(to_rational_function(a), to_rational_function(inv_b), r); break;
case 1: return mul_rf_v(to_rational_function(a), inv_b); case 1: mul_rf_v(to_rational_function(a), inv_b, r); break;
default: UNREACHABLE(); default: UNREACHABLE();
return 0;
} }
} }
} }
value * inv_rf(rational_function_value * a) { void inv_rf(rational_function_value * a, value_ref & r) {
polynomial const & an = a->num(); polynomial const & an = a->num();
polynomial const & ad = a->den(); polynomial const & ad = a->den();
rational_function_value * r = mk_rational_function_value_core(a->ext(), ad.size(), ad.c_ptr(), an.size(), an.c_ptr()); scoped_mpbqi ri(bqim());
bqim().inv(interval(a), r->interval()); bqim().inv(interval(a), ri);
r = mk_rational_function_value_core(a->ext(), ad.size(), ad.c_ptr(), an.size(), an.c_ptr());
swap(r->interval(), ri);
SASSERT(!contains_zero(r->interval())); SASSERT(!contains_zero(r->interval()));
return r;
} }
value * inv(value * a) { void inv(value * a, value_ref & r) {
if (a == 0) { if (a == 0) {
throw exception("division by zero"); throw exception("division by zero");
} }
if (is_nz_rational(a)) { if (is_nz_rational(a)) {
scoped_mpq r(qm()); scoped_mpq v(qm());
qm().inv(to_mpq(a), r); qm().inv(to_mpq(a), v);
return mk_rational(r); r = mk_rational(v);
} }
else { else {
return inv_rf(to_rational_function(a)); inv_rf(to_rational_function(a), r);
} }
} }
@ -2714,36 +2798,56 @@ namespace realclosure {
n.m_value = v; n.m_value = v;
} }
void set(numeral & n, value_ref const & v) {
set(n, v.get());
}
void neg(numeral & a) { void neg(numeral & a) {
set(a, neg(a.m_value)); value_ref r(*this);
neg(a.m_value, r);
set(a, r);
} }
void neg(numeral const & a, numeral & b) { void neg(numeral const & a, numeral & b) {
set(b, neg(a.m_value)); value_ref r(*this);
neg(a.m_value, r);
set(b, r);
} }
void inv(numeral & a) { void inv(numeral & a) {
set(a, inv(a.m_value)); value_ref r(*this);
inv(a.m_value, r);
set(a, r);
} }
void inv(numeral const & a, numeral & b) { void inv(numeral const & a, numeral & b) {
set(b, inv(a.m_value)); value_ref r(*this);
inv(a.m_value, r);
set(b, r);
} }
void add(numeral const & a, numeral const & b, numeral & c) { void add(numeral const & a, numeral const & b, numeral & c) {
set(c, add(a.m_value, b.m_value)); value_ref r(*this);
add(a.m_value, b.m_value, r);
set(c, r);
} }
void sub(numeral const & a, numeral const & b, numeral & c) { void sub(numeral const & a, numeral const & b, numeral & c) {
set(c, sub(a.m_value, b.m_value)); value_ref r(*this);
sub(a.m_value, b.m_value, r);
set(c, r);
} }
void mul(numeral const & a, numeral const & b, numeral & c) { void mul(numeral const & a, numeral const & b, numeral & c) {
set(c, mul(a.m_value, b.m_value)); value_ref r(*this);
mul(a.m_value, b.m_value, r);
set(c, r);
} }
void div(numeral const & a, numeral const & b, numeral & c) { void div(numeral const & a, numeral const & b, numeral & c) {
set(c, div(a.m_value, b.m_value)); value_ref r(*this);
div(a.m_value, b.m_value, r);
set(c, r);
} }
int compare(value * a, value * b) { int compare(value * a, value * b) {
@ -2761,7 +2865,7 @@ namespace realclosure {
return 1; return 1;
else { else {
value_ref diff(*this); value_ref diff(*this);
diff = sub(a, b); sub(a, b, diff);
return sign(diff); return sign(diff);
} }
} }