(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.

diff --git a/src/theory/arith/nl/cad_solver.cpp b/src/theory/arith/nl/cad_solver.cpp
index bb1ef9911969f3a45f8cba5c0a320a818b701d98..b1f0894fb085b9ce8b2ba71cc97d70bb3dbdca22 100644 (file)
--- 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 e3bf4712d473430a689bb33f7a5175625a849d4b..d8dc2270d9dadd57729c45941a562d8c6a707ee8 100644 (file)
--- 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 82ae565d27862acf4242168a2ca8546b6ad11716..d83bc1b2e484632269a2126f58c4462b43142bf8 100644 (file)
--- 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