Fix Taylor overapproximation for large exponentials (#2538)

5 files changed, 82 insertions, 4 deletions

diff --git a/src/theory/arith/nonlinear_extension.cpp b/src/theory/arith/nonlinear_extension.cpp
index 4cd7896d3..929c7808d 100644
--- a/src/theory/arith/nonlinear_extension.cpp
+++ b/src/theory/arith/nonlinear_extension.cpp
@@ -4488,6 +4488,8 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf,
         Assert(rcoeff_n.isConst());
         Rational rcoeff = rcoeff_n.getConst<Rational>();
         Assert(rcoeff.sgn() != 0);
+        poly_approx_b = Rewriter::rewrite(poly_approx_b);
+        poly_approx_c = Rewriter::rewrite(poly_approx_c);
         splane = nm->mkNode(
             PLUS,
             poly_approx_b,
@@ -4881,6 +4883,48 @@ void NonlinearExtension::getPolynomialApproximationBounds(
   }
 }
 
+void NonlinearExtension::getPolynomialApproximationBoundForArg(
+    Kind k, Node c, unsigned d, std::vector<Node>& pbounds)
+{
+  getPolynomialApproximationBounds(k, d, pbounds);
+  Assert(c.isConst());
+  if (k == EXPONENTIAL && c.getConst<Rational>().sgn() == 1)
+  {
+    NodeManager* nm = NodeManager::currentNM();
+    Node tft = nm->mkNode(k, d_zero);
+    bool success = false;
+    unsigned ds = d;
+    TNode ttrf = d_taylor_real_fv;
+    TNode tc = c;
+    do
+    {
+      success = true;
+      unsigned n = 2 * ds;
+      std::pair<Node, Node> taylor = getTaylor(tft, n);
+      // check that 1-c^{n+1}/(n+1)! > 0
+      Node ru = nm->mkNode(DIVISION, taylor.second[1], taylor.second[0][1]);
+      Node rus = ru.substitute(ttrf, tc);
+      rus = Rewriter::rewrite(rus);
+      Assert(rus.isConst());
+      if (rus.getConst<Rational>() > d_one.getConst<Rational>())
+      {
+        success = false;
+        ds = ds + 1;
+      }
+    } while (!success);
+    if (ds > d)
+    {
+      Trace("nl-ext-exp-taylor")
+          << "*** Increase Taylor bound to " << ds << " > " << d << " for ("
+          << k << " " << c << ")" << std::endl;
+      // must use sound upper bound
+      std::vector<Node> pboundss;
+      getPolynomialApproximationBounds(k, ds, pboundss);
+      pbounds[2] = pboundss[2];
+    }
+  }
+}
+
 std::pair<Node, Node> NonlinearExtension::getTfModelBounds(Node tf, unsigned d)
 {
   // compute the model value of the argument
@@ -4891,7 +4935,7 @@ std::pair<Node, Node> NonlinearExtension::getTfModelBounds(Node tf, unsigned d)
   bool isNeg = csign == -1;
 
   std::vector<Node> pbounds;
-  getPolynomialApproximationBounds(tf.getKind(), d, pbounds);
+  getPolynomialApproximationBoundForArg(tf.getKind(), c, d, pbounds);
 
   std::vector<Node> bounds;
   TNode tfv = d_taylor_real_fv;
diff --git a/src/theory/arith/nonlinear_extension.h b/src/theory/arith/nonlinear_extension.h
index f2cdea334..cb74502d6 100644
--- a/src/theory/arith/nonlinear_extension.h
+++ b/src/theory/arith/nonlinear_extension.h
@@ -649,15 +649,32 @@ class NonlinearExtension {
   unsigned d_taylor_degree;
   /** polynomial approximation bounds
    *
-   * This adds P_l+, P_l-, P_u+, P_u- to pbounds, where these are polynomial
-   * approximations of the Taylor series of <k>( 0 ) for degree 2*d where
-   * k is SINE or EXPONENTIAL.
+   * This adds P_l+[x], P_l-[x], P_u+[x], P_u-[x] to pbounds, where x is
+   * d_taylor_real_fv. These are polynomial approximations of the Taylor series
+   * of <k>( 0 ) for degree 2*d where k is SINE or EXPONENTIAL.
    * These correspond to P_l and P_u from Figure 3 of Cimatti et al., CADE 2017,
    * for positive/negative (+/-) values of the argument of <k>( 0 ).
+   *
+   * Notice that for certain bounds (e.g. upper bounds for exponential), the
+   * Taylor approximation for a fixed degree is only sound up to a given
+   * upper bound on the argument. To obtain sound lower/upper bounds for a
+   * given <k>( c ), use the function below.
    */
   void getPolynomialApproximationBounds(Kind k,
                                         unsigned d,
                                         std::vector<Node>& pbounds);
+  /** polynomial approximation bounds
+   *
+   * This computes polynomial approximations P_l+[x], P_l-[x], P_u+[x], P_u-[x]
+   * that are sound (lower, upper) bounds for <k>( c ). Notice that these
+   * polynomials may depend on c. In particular, for P_u+[x] for <k>( c ) where
+   * c>0, we return the P_u+[x] from the function above for the minimum degree
+   * d' >= d such that (1-c^{2*d'+1}/(2*d'+1)!) is positive.
+   */
+  void getPolynomialApproximationBoundForArg(Kind k,
+                                             Node c,
+                                             unsigned d,
+                                             std::vector<Node>& pbounds);
   /** cache of the above function */
   std::map<Kind, std::map<unsigned, std::vector<Node> > > d_poly_bounds;
   /** get transcendental function model bounds
diff --git a/test/regress/CMakeLists.txt b/test/regress/CMakeLists.txt
index cec33c615..4e8405d5d 100644
--- a/test/regress/CMakeLists.txt
+++ b/test/regress/CMakeLists.txt
@@ -1166,6 +1166,7 @@ set(regress_1_tests
   regress1/nl/dist-big.smt2
   regress1/nl/div-mod-partial.smt2
   regress1/nl/dumortier_llibre_artes_ex_5_13.transcendental.k2.smt2
+  regress1/nl/exp-approx.smt2
   regress1/nl/exp-4.5-lt.smt2
   regress1/nl/exp1-lb.smt2
   regress1/nl/exp_monotone.smt2
diff --git a/test/regress/Makefile.tests b/test/regress/Makefile.tests
index ce7bc8736..d6d8d1a68 100644
--- a/test/regress/Makefile.tests
+++ b/test/regress/Makefile.tests
@@ -1164,6 +1164,7 @@ REG1_TESTS = \
 	regress1/nl/div-mod-partial.smt2 \
 	regress1/nl/dumortier_llibre_artes_ex_5_13.transcendental.k2.smt2 \
 	regress1/nl/exp-4.5-lt.smt2 \
+	regress1/nl/exp-approx.smt2 \
 	regress1/nl/exp1-lb.smt2 \
 	regress1/nl/exp_monotone.smt2 \
 	regress1/nl/factor_agg_s.smt2 \
diff --git a/test/regress/regress1/nl/exp-approx.smt2 b/test/regress/regress1/nl/exp-approx.smt2
new file mode 100644
index 000000000..f3faceda5
--- /dev/null
+++ b/test/regress/regress1/nl/exp-approx.smt2
@@ -0,0 +1,15 @@
+; COMMAND-LINE: --no-check-models
+; EXPECT: sat
+(set-logic QF_NRAT)
+(set-info :status sat)
+
+(assert (and (<= 3.0 (exp 1.1)) (<= (exp 1.1) 3.1)))
+(assert (and (<= 8.1 (exp 2.1)) (<= (exp 2.1) 8.2)))
+(assert (and (<= 22.1 (exp 3.1)) (<= (exp 3.1) 22.2)))
+(assert (and (<= 60.3 (exp 4.1)) (<= (exp 4.1) 60.4)))
+(assert (and (<= 164.0 (exp 5.1)) (<= (exp 5.1) 164.1)))
+
+; this takes ~10 seconds in production
+;(assert (and (<= 536190464.4 (exp 20.1)) (<= (exp 20.1) 536190464.5)))
+
+(check-sat)