(cad-solver) Fixed excluding lemma generation. (#4958)

The lemma generation for partial cad checks had a number of issues that have been fixed in this PR.
The previous version had both soundness issues and a very naive approach to handling algebraic numbers.
This new version is sound (fingers crossed) and allows to construct more precise, but also more complex lemmas.
To avoid constructing very large lemmas, a (somewhat arbitrary) limit based on the size of coefficients was added.

3 files changed, 241 insertions, 19 deletions

diff --git a/src/theory/arith/nl/cad_solver.cpp b/src/theory/arith/nl/cad_solver.cpp
index bb1ef9911..b1f0894fb 100644
--- a/src/theory/arith/nl/cad_solver.cpp
+++ b/src/theory/arith/nl/cad_solver.cpp
@@ -134,10 +134,12 @@ std::vector<NlLemma> CadSolver::checkPartial()
         premise = nm->mkNode(Kind::AND, interval.d_origins);
       }
       Node conclusion =
-          excluding_interval_to_lemma(first_var, interval.d_interval);
-      Node lemma = nm->mkNode(Kind::IMPLIES, premise, conclusion);
-      Trace("nl-cad") << "Excluding " << first_var << " -> " << interval.d_interval << " using " << lemma << std::endl;
-      lems.emplace_back(lemma, Inference::CAD_EXCLUDED_INTERVAL);
+          excluding_interval_to_lemma(first_var, interval.d_interval, false);
+      if (!conclusion.isNull()) {
+        Node lemma = nm->mkNode(Kind::IMPLIES, premise, conclusion);
+        Trace("nl-cad") << "Excluding " << first_var << " -> " << interval.d_interval << " using " << lemma << std::endl;
+        lems.emplace_back(lemma, Inference::CAD_EXCLUDED_INTERVAL);
+       }
     }
   }
   return lems;
diff --git a/src/theory/arith/nl/poly_conversion.cpp b/src/theory/arith/nl/poly_conversion.cpp
index e3bf4712d..d8dc2270d 100644
--- a/src/theory/arith/nl/poly_conversion.cpp
+++ b/src/theory/arith/nl/poly_conversion.cpp
@@ -346,12 +346,126 @@ Node value_to_node(const poly::Value& v, const Node& ran_variable)
   return nm->mkConst(Rational(0));
 }
 
+Node lower_bound_as_node(const Node& var,
+                         const poly::Value& lower,
+                         bool open,
+                         bool allowNonlinearLemma)
+{
+  auto* nm = NodeManager::currentNM();
+  if (!poly::is_algebraic_number(lower))
+  {
+    return nm->mkNode(open ? Kind::LEQ : Kind::LT,
+                      var,
+                      nm->mkConst(poly_utils::toRationalAbove(lower)));
+  }
+  if (poly::represents_rational(lower))
+  {
+    return nm->mkNode(
+        open ? Kind::LEQ : Kind::LT,
+        var,
+        nm->mkConst(poly_utils::toRationalAbove(poly::get_rational(lower))));
+  }
+  if (!allowNonlinearLemma)
+  {
+    return Node();
+  }
+
+  const poly::AlgebraicNumber& alg = as_algebraic_number(lower);
+
+  Node poly = as_cvc_upolynomial(get_defining_polynomial(alg), var);
+  Rational l = poly_utils::toRational(
+      poly::get_lower(poly::get_isolating_interval(alg)));
+  Rational u = poly_utils::toRational(
+      poly::get_upper(poly::get_isolating_interval(alg)));
+  int sl = poly::sign_at(get_defining_polynomial(alg),
+                         poly::get_lower(poly::get_isolating_interval(alg)));
+  int su = poly::sign_at(get_defining_polynomial(alg),
+                         poly::get_upper(poly::get_isolating_interval(alg)));
+  Assert(sl != 0 && su != 0 && sl != su);
+
+  // open:  var <= l  or  (var < u  and  sgn(poly(var)) == sl)
+  // !open: var <= l  or  (var < u  and  sgn(poly(var)) == sl/0)
+  Kind relation;
+  if (open)
+  {
+    relation = (sl < 0) ? Kind::LEQ : Kind::GEQ;
+  }
+  else
+  {
+    relation = (sl < 0) ? Kind::LT : Kind::GT;
+  }
+  return nm->mkNode(
+      Kind::OR,
+      nm->mkNode(Kind::LEQ, var, nm->mkConst(l)),
+      nm->mkNode(Kind::AND,
+                 nm->mkNode(Kind::LT, var, nm->mkConst(u)),
+                 nm->mkNode(relation, poly, nm->mkConst(Rational(0)))));
+}
+
+Node upper_bound_as_node(const Node& var,
+                         const poly::Value& upper,
+                         bool open,
+                         bool allowNonlinearLemma)
+{
+  auto* nm = NodeManager::currentNM();
+  if (!poly::is_algebraic_number(upper))
+  {
+    return nm->mkNode(open ? Kind::GEQ : Kind::GT,
+                      var,
+                      nm->mkConst(poly_utils::toRationalAbove(upper)));
+  }
+  if (poly::represents_rational(upper))
+  {
+    return nm->mkNode(
+        open ? Kind::GEQ : Kind::GT,
+        var,
+        nm->mkConst(poly_utils::toRationalAbove(poly::get_rational(upper))));
+  }
+  if (!allowNonlinearLemma)
+  {
+    return Node();
+  }
+
+  const poly::AlgebraicNumber& alg = as_algebraic_number(upper);
+
+  Node poly = as_cvc_upolynomial(get_defining_polynomial(alg), var);
+  Rational l = poly_utils::toRational(
+      poly::get_lower(poly::get_isolating_interval(alg)));
+  Rational u = poly_utils::toRational(
+      poly::get_upper(poly::get_isolating_interval(alg)));
+  int sl = poly::sign_at(get_defining_polynomial(alg),
+                         poly::get_lower(poly::get_isolating_interval(alg)));
+  int su = poly::sign_at(get_defining_polynomial(alg),
+                         poly::get_upper(poly::get_isolating_interval(alg)));
+  Assert(sl != 0 && su != 0 && sl != su);
+
+  // open:  var >= u  or  (var > l  and  sgn(poly(var)) == su)
+  // !open: var >= u  or  (var > l  and  sgn(poly(var)) == su/0)
+  Kind relation;
+  if (open)
+  {
+    relation = (su < 0) ? Kind::LEQ : Kind::GEQ;
+  }
+  else
+  {
+    relation = (su < 0) ? Kind::LT : Kind::GT;
+  }
+  return nm->mkNode(
+      Kind::OR,
+      nm->mkNode(Kind::GEQ, var, nm->mkConst(u)),
+      nm->mkNode(Kind::AND,
+                 nm->mkNode(Kind::GT, var, nm->mkConst(l)),
+                 nm->mkNode(relation, poly, nm->mkConst(Rational(0)))));
+}
+
 Node excluding_interval_to_lemma(const Node& variable,
-                                 const poly::Interval& interval)
+                                 const poly::Interval& interval,
+                                 bool allowNonlinearLemma)
 {
   auto* nm = NodeManager::currentNM();
   const auto& lv = poly::get_lower(interval);
   const auto& uv = poly::get_upper(interval);
+  if (bitsize(lv) > 100 || bitsize(uv) > 100) return Node();
   bool li = poly::is_minus_infinity(lv);
   bool ui = poly::is_plus_infinity(uv);
   if (li && ui) return nm->mkConst(true);
@@ -359,16 +473,35 @@ Node excluding_interval_to_lemma(const Node& variable,
   {
     if (is_algebraic_number(lv))
     {
+      const poly::AlgebraicNumber& alg = as_algebraic_number(lv);
+      if (poly::is_rational(alg))
+      {
+        Trace("nl-cad") << "Rational point interval: " << interval << std::endl;
+        return nm->mkNode(Kind::DISTINCT,
+                          variable,
+                          nm->mkConst(poly_utils::toRational(
+                              poly::to_rational_approximation(alg))));
+      }
+      Trace("nl-cad") << "Algebraic point interval: " << interval << std::endl;
+      // p(x) != 0 or x <= lb or ub <= x
+      if (allowNonlinearLemma)
+      {
+        Node poly = as_cvc_upolynomial(get_defining_polynomial(alg), variable);
       return nm->mkNode(
           Kind::OR,
-          nm->mkNode(
-              Kind::LT, variable, nm->mkConst(poly_utils::toRationalBelow(lv))),
+            nm->mkNode(Kind::DISTINCT, poly, nm->mkConst(Rational(0))),
+            nm->mkNode(Kind::LT,
+                       variable,
+                       nm->mkConst(poly_utils::toRationalBelow(lv))),
           nm->mkNode(Kind::GT,
                      variable,
                      nm->mkConst(poly_utils::toRationalAbove(lv))));
+      }
+      return Node();
     }
     else
     {
+      Trace("nl-cad") << "Rational point interval: " << interval << std::endl;
       return nm->mkNode(Kind::DISTINCT,
                         variable,
                         nm->mkConst(poly_utils::toRationalBelow(lv)));
@@ -376,20 +509,23 @@ Node excluding_interval_to_lemma(const Node& variable,
   }
   if (li)
   {
-    return nm->mkNode(
-        Kind::GT, variable, nm->mkConst(poly_utils::toRationalAbove(uv)));
+    Trace("nl-cad") << "Only upper bound: " << interval << std::endl;
+    return upper_bound_as_node(
+        variable, uv, poly::get_upper_open(interval), allowNonlinearLemma);
   }
   if (ui)
   {
-    return nm->mkNode(
-        Kind::LT, variable, nm->mkConst(poly_utils::toRationalBelow(lv)));
+    Trace("nl-cad") << "Only lower bound: " << interval << std::endl;
+    return lower_bound_as_node(
+        variable, lv, poly::get_lower_open(interval), allowNonlinearLemma);
   }
-  return nm->mkNode(
-      Kind::OR,
-      nm->mkNode(
-          Kind::GT, variable, nm->mkConst(poly_utils::toRationalAbove(uv))),
-      nm->mkNode(
-          Kind::LT, variable, nm->mkConst(poly_utils::toRationalBelow(lv))));
+  Trace("nl-cad") << "Proper interval: " << interval << std::endl;
+  Node lb = lower_bound_as_node(
+      variable, lv, poly::get_lower_open(interval), allowNonlinearLemma);
+  Node ub = upper_bound_as_node(
+      variable, uv, poly::get_upper_open(interval), allowNonlinearLemma);
+  if (lb.isNull() || ub.isNull()) return Node();
+  return nm->mkNode(Kind::OR, lb, ub);
 }
 
 Maybe<Rational> get_lower_bound(const Node& n)
@@ -503,6 +639,73 @@ poly::Value node_to_value(const Node& n, const Node& ran_variable)
   return node_to_poly_ran(n, ran_variable);
 }
 
+/** Bitsize of a dyadic rational. */
+std::size_t bitsize(const poly::DyadicRational& v)
+{
+  return bit_size(numerator(v)) + bit_size(denominator(v));
+}
+/** Bitsize of an integer. */
+std::size_t bitsize(const poly::Integer& v) { return bit_size(v); }
+/** Bitsize of a rational. */
+std::size_t bitsize(const poly::Rational& v)
+{
+  return bit_size(numerator(v)) + bit_size(denominator(v));
+}
+/** Bitsize of a univariate polynomial. */
+std::size_t bitsize(const poly::UPolynomial& v)
+{
+  std::size_t sum = 0;
+  for (const auto& c : coefficients(v))
+  {
+    sum += bitsize(c);
+  }
+  return sum;
+}
+/** Bitsize of an algebraic number. */
+std::size_t bitsize(const poly::AlgebraicNumber& v)
+{
+  if (is_rational(v))
+  {
+    return bitsize(to_rational_approximation(v));
+  }
+  return bitsize(get_lower_bound(v)) + bitsize(get_upper_bound(v))
+         + bitsize(get_defining_polynomial(v));
+}
+
+std::size_t bitsize(const poly::Value& v)
+{
+  if (is_algebraic_number(v))
+  {
+    return bitsize(as_algebraic_number(v));
+  }
+  else if (is_dyadic_rational(v))
+  {
+    return bitsize(as_dyadic_rational(v));
+  }
+  else if (is_integer(v))
+  {
+    return bitsize(as_integer(v));
+  }
+  else if (is_minus_infinity(v))
+  {
+    return 1;
+  }
+  else if (is_none(v))
+  {
+    return 0;
+  }
+  else if (is_plus_infinity(v))
+  {
+    return 1;
+  }
+  else if (is_rational(v))
+  {
+    return bitsize(as_rational(v));
+  }
+  Assert(false);
+  return 0;
+}
+
 }  // namespace nl
 }  // namespace arith
 }  // namespace theory
diff --git a/src/theory/arith/nl/poly_conversion.h b/src/theory/arith/nl/poly_conversion.h
index 82ae565d2..d83bc1b2e 100644
--- a/src/theory/arith/nl/poly_conversion.h
+++ b/src/theory/arith/nl/poly_conversion.h
@@ -91,14 +91,24 @@ Node value_to_node(const poly::Value& v, const Node& ran_variable);
 
 /**
  * Constructs a lemma that excludes a given interval from the feasible values of
- * a variable. The resulting lemma has the form
+ * a variable. Conceptually, the resulting lemma has the form
  * (OR
  *    (<= var interval.lower)
  *    (<= interval.upper var)
  * )
+ * This method honors the interval bound types (open or closed), but also deals
+ * with real algebraic endpoints. If allowNonlinearLemma is false, real
+ * algebraic endpoints are reflected by coarse (numeric) approximations and thus
+ * may lead to lemmas that are weaker than they could be. Also, lemma creation
+ * may fail altogether.
+ * If allowNonlinearLemma is true, it tries to construct better lemmas (based on
+ * the sign of the defining polynomial of the real algebraic number). These
+ * lemmas are nonlinear, though, and may thus be expensive to use in the
+ * subsequent solving process.
  */
 Node excluding_interval_to_lemma(const Node& variable,
-                                 const poly::Interval& interval);
+                                 const poly::Interval& interval,
+                                 bool allowNonlinearLemma);
 
 /**
  * Transforms a node to a poly::AlgebraicNumber.
@@ -119,6 +129,13 @@ RealAlgebraicNumber node_to_ran(const Node& n, const Node& ran_variable);
  */
 poly::Value node_to_value(const Node& n, const Node& ran_variable);
 
+/**
+ * Give a rough estimate of the bitsize of the representation of `v`.
+ * It can be used as a rough measure of the size of complexity of a value, for
+ * example to avoid divergence or disallow huge lemmas.
+ */
+std::size_t bitsize(const poly::Value& v);
+
 }  // namespace nl
 }  // namespace arith
 }  // namespace theory