Skip to content

Commit 0e71a9d

Browse files
levnachlevnach
committed
comment and restore
Signed-off-by: Lev Nachmanson <[email protected]>
1 parent 84c8a93 commit 0e71a9d

File tree

1 file changed

+160
-3
lines changed

1 file changed

+160
-3
lines changed

‎src/nlsat/nlsat_explain.cpp‎

Lines changed: 160 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -617,7 +617,10 @@ namespace nlsat {
617617
m_todo.insert(p);
618618
}
619619
}
620-
620+
// It is the only check I found to test for well-orientedness, that first appears in
621+
// "An improved projection operation for cylindrical algebraic decomposition of three-dimensional space",
622+
// by McCallum, Scott
623+
621624
bool is_well_oriented(polynomial_ref_vector &ps, var x) {
622625
polynomial_ref p(m_pm);
623626
polynomial_ref lc_poly(m_pm);
@@ -631,15 +634,14 @@ namespace nlsat {
631634
// p depends on x
632635
disc_poly = discriminant(p, x); // Use global helper
633636
if (sign(disc_poly) == 0) { // Discriminant is zero
634-
635637
TRACE(nlsat_explain, tout << "Global !WO: Discriminant of poly is zero. Poly: "; display(tout, p); tout << " Disc: "; display(tout, disc_poly) << "\\n";);
636638
return false;
637639
}
638640
}
639641
return true;
640642
}
641643

642-
// For each p in ps add the leading or all the coefficients of p to the projection,
644+
// For each p in ps add the leading or all the coefficients of p to the projection,
643645
// depending on the well-orientedness of ps.
644646
void add_lcs(polynomial_ref_vector &ps, var x) {
645647
polynomial_ref p(m_pm);
@@ -1918,6 +1920,161 @@ namespace nlsat {
19181920
scoped_anum_vector & roots = m_roots_tmp;
19191921
roots.reset();
19201922
m_am.isolate_roots(p, undef_var_assignment(m_assignment, x), roots);
1923+
for (auto const& r : roots) {
1924+
int s = m_am.compare(x_val, r);
1925+
SASSERT(s != 0);
1926+
1927+
if (s < 0 && (!lub_valid || m_am.lt(r, lub))) {
1928+
lub_index = i;
1929+
m_am.set(lub, r);
1930+
lub_valid = true;
1931+
}
1932+
1933+
if (s > 0 && (!glb_valid || m_am.lt(glb, r))) {
1934+
glb_index = i;
1935+
m_am.set(glb, r);
1936+
glb_valid = true;
1937+
}
1938+
if (s < 0) ++num_lub;
1939+
if (s > 0) ++num_glb;
1940+
}
1941+
}
1942+
TRACE(nlsat_explain, tout << "glb: " << num_glb << " lub: " << num_lub << "\n" << lub_index << "\n" << glb_index << "\n" << ps << "\n";);
1943+
1944+
if (num_lub == 0) {
1945+
project_plus_infinity(x, ps);
1946+
return;
1947+
}
1948+
1949+
if (num_glb == 0) {
1950+
project_minus_infinity(x, ps);
1951+
return;
1952+
}
1953+
1954+
if (num_lub <= num_glb) {
1955+
glb_index = lub_index;
1956+
}
1957+
1958+
project_pairs(x, glb_index, ps);
1959+
}
1960+
1961+
void project_plus_infinity(var x, polynomial_ref_vector const& ps) {
1962+
polynomial_ref p(m_pm), lc(m_pm);
1963+
for (unsigned i = 0; i < ps.size(); ++i) {
1964+
p = ps.get(i);
1965+
unsigned d = degree(p, x);
1966+
lc = m_pm.coeff(p, x, d);
1967+
if (!is_const(lc)) {
1968+
int s = sign(p);
1969+
SASSERT(s != 0);
1970+
atom::kind k = (s > 0)?(atom::GT):(atom::LT);
1971+
add_simple_assumption(k, lc);
1972+
}
1973+
}
1974+
}
1975+
1976+
void project_minus_infinity(var x, polynomial_ref_vector const& ps) {
1977+
polynomial_ref p(m_pm), lc(m_pm);
1978+
for (unsigned i = 0; i < ps.size(); ++i) {
1979+
p = ps.get(i);
1980+
unsigned d = degree(p, x);
1981+
lc = m_pm.coeff(p, x, d);
1982+
if (!is_const(lc)) {
1983+
int s = sign(p);
1984+
TRACE(nlsat_explain, tout << "degree: " << d << " " << lc << " sign: " << s << "\n";);
1985+
SASSERT(s != 0);
1986+
atom::kind k;
1987+
if (s > 0) {
1988+
k = (d % 2 == 0)?(atom::GT):(atom::LT);
1989+
}
1990+
else {
1991+
k = (d % 2 == 0)?(atom::LT):(atom::GT);
1992+
}
1993+
add_simple_assumption(k, lc);
1994+
}
1995+
}
1996+
}
1997+
1998+
void project_pairs(var x, unsigned idx, polynomial_ref_vector const& ps) {
1999+
TRACE(nlsat_explain, tout << "project pairs\n";);
2000+
polynomial_ref p(m_pm);
2001+
p = ps.get(idx);
2002+
for (unsigned i = 0; i < ps.size(); ++i) {
2003+
if (i != idx) {
2004+
project_pair(x, ps.get(i), p);
2005+
}
2006+
}
2007+
}
2008+
2009+
void project_pair(var x, polynomial::polynomial* p1, polynomial::polynomial* p2) {
2010+
m_ps2.reset();
2011+
m_ps2.push_back(p1);
2012+
m_ps2.push_back(p2);
2013+
project(m_ps2, x);
2014+
}
2015+
2016+
void project_single(var x, polynomial::polynomial* p) {
2017+
m_ps2.reset();
2018+
m_ps2.push_back(p);
2019+
project(m_ps2, x);
2020+
}
2021+
2022+
void solve_eq(var x, unsigned idx, polynomial_ref_vector const& ps) {
2023+
polynomial_ref p(m_pm), A(m_pm), B(m_pm), C(m_pm), D(m_pm), E(m_pm), q(m_pm), r(m_pm);
2024+
polynomial_ref_vector As(m_pm), Bs(m_pm);
2025+
p = ps.get(idx);
2026+
SASSERT(degree(p, x) == 1);
2027+
A = m_pm.coeff(p, x, 1);
2028+
B = m_pm.coeff(p, x, 0);
2029+
As.push_back(m_pm.mk_const(rational(1)));
2030+
Bs.push_back(m_pm.mk_const(rational(1)));
2031+
B = neg(B);
2032+
TRACE(nlsat_explain, tout << "p: " << p << " A: " << A << " B: " << B << "\n";);
2033+
// x = B/A
2034+
for (unsigned i = 0; i < ps.size(); ++i) {
2035+
if (i != idx) {
2036+
q = ps.get(i);
2037+
unsigned d = degree(q, x);
2038+
D = m_pm.mk_const(rational(1));
2039+
E = D;
2040+
r = m_pm.mk_zero();
2041+
for (unsigned j = As.size(); j <= d; ++j) {
2042+
D = As.back(); As.push_back(A * D);
2043+
D = Bs.back(); Bs.push_back(B * D);
2044+
}
2045+
for (unsigned j = 0; j <= d; ++j) {
2046+
// A^d*p0 + A^{d-1}*B*p1 + ... + B^j*A^{d-j}*pj + ... + B^d*p_d
2047+
C = m_pm.coeff(q, x, j);
2048+
TRACE(nlsat_explain, tout << "coeff: q" << j << ": " << C << "\n";);
2049+
if (!is_zero(C)) {
2050+
D = As.get(d - j);
2051+
E = Bs.get(j);
2052+
r = r + D*E*C;
2053+
}
2054+
}
2055+
TRACE(nlsat_explain, tout << "p: " << p << " q: " << q << " r: " << r << "\n";);
2056+
ensure_sign(r);
2057+
}
2058+
else {
2059+
ensure_sign(A);
2060+
}
2061+
}
2062+
2063+
}
2064+
2065+
void maximize(var x, unsigned num, literal const * ls, scoped_anum& val, bool& unbounded) {
2066+
svector<literal> lits;
2067+
polynomial_ref p(m_pm);
2068+
split_literals(x, num, ls, lits);
2069+
collect_polys(lits.size(), lits.data(), m_ps);
2070+
unbounded = true;
2071+
scoped_anum x_val(m_am);
2072+
x_val = m_assignment.value(x);
2073+
for (unsigned i = 0; i < m_ps.size(); ++i) {
2074+
p = m_ps.get(i);
2075+
scoped_anum_vector & roots = m_roots_tmp;
2076+
roots.reset();
2077+
m_am.isolate_roots(p, undef_var_assignment(m_assignment, x), roots);
19212078
for (unsigned j = 0; j < roots.size(); ++j) {
19222079
int s = m_am.compare(x_val, roots[j]);
19232080
if (s <= 0 && (unbounded || m_am.compare(roots[j], val) <= 0)) {

0 commit comments

Comments
 (0)