Fix transcendental secant plane lemmas (#5525)

The refactoring of the transcendental solver introduced a subtle issue that lead to incorrect secant plane lemmas.
This PR fixes this issue, so that we now produce the correct lemmas again.

7 files changed, 143 insertions, 123 deletions

diff --git a/src/theory/arith/nl/transcendental/exponential_solver.cpp b/src/theory/arith/nl/transcendental/exponential_solver.cpp
index 4275b2487..7d7d0c23c 100644
--- a/src/theory/arith/nl/transcendental/exponential_solver.cpp
+++ b/src/theory/arith/nl/transcendental/exponential_solver.cpp
@@ -170,12 +170,11 @@ void ExponentialSolver::doTangentLemma(TNode e, TNode c, TNode poly_approx)
   d_data->d_im.addPendingArithLemma(lem, InferenceId::NL_T_TANGENT, nullptr, true);
 }
 
-void ExponentialSolver::doSecantLemmas(TNode e,
-                                       TNode c,
-                                       TNode poly_approx,
-                                       unsigned d)
+void ExponentialSolver::doSecantLemmas(
+    TNode e, TNode poly_approx, TNode c, TNode poly_approx_c, unsigned d)
 {
-  d_data->doSecantLemmas(getSecantBounds(e, c, d), c, poly_approx, e, d, 1);
+  d_data->doSecantLemmas(
+      getSecantBounds(e, c, d), poly_approx, c, poly_approx_c, e, d, 1);
 }
 
 std::pair<Node, Node> ExponentialSolver::getSecantBounds(TNode e,
diff --git a/src/theory/arith/nl/transcendental/exponential_solver.h b/src/theory/arith/nl/transcendental/exponential_solver.h
index f757e5892..b20c23e56 100644
--- a/src/theory/arith/nl/transcendental/exponential_solver.h
+++ b/src/theory/arith/nl/transcendental/exponential_solver.h
@@ -87,7 +87,8 @@ class ExponentialSolver
   void doTangentLemma(TNode e, TNode c, TNode poly_approx);
 
   /** Sent secant lemmas around c for e */
-  void doSecantLemmas(TNode e, TNode c, TNode poly_approx, unsigned d);
+  void doSecantLemmas(
+      TNode e, TNode poly_approx, TNode c, TNode poly_approx_c, unsigned d);
 
  private:
   /** Generate bounds for secant lemmas */
diff --git a/src/theory/arith/nl/transcendental/sine_solver.cpp b/src/theory/arith/nl/transcendental/sine_solver.cpp
index 012b8e7ec..31fd47503 100644
--- a/src/theory/arith/nl/transcendental/sine_solver.cpp
+++ b/src/theory/arith/nl/transcendental/sine_solver.cpp
@@ -289,12 +289,17 @@ void SineSolver::doTangentLemma(TNode e, TNode c, TNode poly_approx, int region)
   d_data->d_im.addPendingArithLemma(lem, InferenceId::NL_T_TANGENT, nullptr, true);
 }
 
-void SineSolver::doSecantLemmas(
-    TNode e, TNode c, TNode poly_approx, unsigned d, int region)
+void SineSolver::doSecantLemmas(TNode e,
+                                TNode poly_approx,
+                                TNode c,
+                                TNode poly_approx_c,
+                                unsigned d,
+                                int region)
 {
   d_data->doSecantLemmas(getSecantBounds(e, c, d, region),
-                         c,
                          poly_approx,
+                         c,
+                         poly_approx_c,
                          e,
                          d,
                          regionToConcavity(region));
diff --git a/src/theory/arith/nl/transcendental/sine_solver.h b/src/theory/arith/nl/transcendental/sine_solver.h
index 394c2d927..9f9102a53 100644
--- a/src/theory/arith/nl/transcendental/sine_solver.h
+++ b/src/theory/arith/nl/transcendental/sine_solver.h
@@ -88,8 +88,12 @@ class SineSolver
   void doTangentLemma(TNode e, TNode c, TNode poly_approx, int region);
 
   /** Sent secant lemmas around c for e */
-  void doSecantLemmas(
-      TNode e, TNode c, TNode poly_approx, unsigned d, int region);
+  void doSecantLemmas(TNode e,
+                      TNode poly_approx,
+                      TNode c,
+                      TNode poly_approx_c,
+                      unsigned d,
+                      int region);
 
  private:
   std::pair<Node, Node> getSecantBounds(TNode e,
diff --git a/src/theory/arith/nl/transcendental/transcendental_solver.cpp b/src/theory/arith/nl/transcendental/transcendental_solver.cpp
index 1f76cd833..8eb69e50b 100644
--- a/src/theory/arith/nl/transcendental/transcendental_solver.cpp
+++ b/src/theory/arith/nl/transcendental/transcendental_solver.cpp
@@ -362,11 +362,11 @@ bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d)
   {
     if (k == EXPONENTIAL)
     {
-      d_expSlv.doSecantLemmas(tf, c, poly_approx_c, d);
+      d_expSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d);
     }
     else if (k == Kind::SINE)
     {
-      d_sineSlv.doSecantLemmas(tf, c, poly_approx_c, d, region);
+      d_sineSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d, region);
     }
   }
   return true;
diff --git a/src/theory/arith/nl/transcendental/transcendental_state.cpp b/src/theory/arith/nl/transcendental/transcendental_state.cpp
index b6882d606..41ed2c53d 100644
--- a/src/theory/arith/nl/transcendental/transcendental_state.cpp
+++ b/src/theory/arith/nl/transcendental/transcendental_state.cpp
@@ -294,12 +294,9 @@ std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e,
 }
 
 Node TranscendentalState::mkSecantPlane(
-    TNode arg, TNode b, TNode c, TNode approx, TNode approx_c)
+    TNode arg, TNode b, TNode c, TNode approx_b, TNode approx_c)
 {
   NodeManager* nm = NodeManager::currentNM();
-  // Figure 3 : P(l), P(u), for s = 0,1
-  Node approx_b =
-      Rewriter::rewrite(approx.substitute(d_taylor.getTaylorVariable(), b));
   // Figure 3: S_l( x ), S_u( x ) for s = 0,1
   Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, b, c));
   Assert(rcoeff_n.isConst());
@@ -345,18 +342,26 @@ NlLemma TranscendentalState::mkSecantLemma(
 }
 
 void TranscendentalState::doSecantLemmas(const std::pair<Node, Node>& bounds,
+                                         TNode poly_approx,
                                          TNode c,
                                          TNode approx_c,
                                          TNode tf,
                                          unsigned d,
                                          int concavity)
 {
+  Trace("nl-ext-tftp-debug2") << "...secant bounds are : " << bounds.first
+                              << " ... " << bounds.second << std::endl;
   // take the model value of l or u (since may contain PI)
   // Make secant from bounds.first to c
   Node lval = d_model.computeAbstractModelValue(bounds.first);
+  Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds.first
+                              << " is " << lval << std::endl;
   if (lval != c)
   {
-    Node splane = mkSecantPlane(tf[0], lval, c, bounds.first, approx_c);
+    // Figure 3 : P(l), P(u), for s = 0
+    Node approx_l = Rewriter::rewrite(
+        poly_approx.substitute(d_taylor.getTaylorVariable(), lval));
+    Node splane = mkSecantPlane(tf[0], lval, c, approx_l, approx_c);
     NlLemma nlem = mkSecantLemma(lval, c, concavity, tf, splane);
     // The side effect says that if lem is added, then we should add the
     // secant point c for (tf,d).
@@ -366,9 +371,14 @@ void TranscendentalState::doSecantLemmas(const std::pair<Node, Node>& bounds,
 
   // Make secant from c to bounds.second
   Node uval = d_model.computeAbstractModelValue(bounds.second);
+  Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds.second
+                              << " is " << uval << std::endl;
   if (c != uval)
   {
-    Node splane = mkSecantPlane(tf[0], c, uval, approx_c, bounds.second);
+    // Figure 3 : P(l), P(u), for s = 1
+    Node approx_u = Rewriter::rewrite(
+        poly_approx.substitute(d_taylor.getTaylorVariable(), uval));
+    Node splane = mkSecantPlane(tf[0], c, uval, approx_c, approx_u);
     NlLemma nlem = mkSecantLemma(c, uval, concavity, tf, splane);
     // The side effect says that if lem is added, then we should add the
     // secant point c for (tf,d).
diff --git a/src/theory/arith/nl/transcendental/transcendental_state.h b/src/theory/arith/nl/transcendental/transcendental_state.h
index 2cf6400a3..0a4e46eca 100644
--- a/src/theory/arith/nl/transcendental/transcendental_state.h
+++ b/src/theory/arith/nl/transcendental/transcendental_state.h
@@ -89,6 +89,7 @@ struct TranscendentalState
   /**
    * Construct and send secant lemmas (if appropriate)
    * @param bounds Secant bounds
+   * @param poly_approx Polynomial approximation
    * @param c Current point
    * @param approx_c Approximation for c
    * @param tf Current transcendental term
@@ -96,118 +97,118 @@ struct TranscendentalState
    * @param concavity Concavity in region of c
    */
   void doSecantLemmas(const std::pair<Node, Node>& bounds,
+                      TNode poly_approx,
                       TNode c,
                       TNode approx_c,
                       TNode tf,
                       unsigned d,
-                      int concavity)
-    ;
-
-    Node d_true;
-    Node d_false;
-    Node d_zero;
-    Node d_one;
-    Node d_neg_one;
-
-    /** The inference manager that we push conflicts and lemmas to. */
-    InferenceManager& d_im;
-    /** Reference to the non-linear model object */
-    NlModel& d_model;
-    /** Utility to compute taylor approximations */
-    TaylorGenerator d_taylor;
-
-    /**
-     * Some transcendental functions f(t) are "purified", e.g. we add
-     * t = y ^ f(t) = f(y) where y is a fresh variable. Those that are not
-     * purified we call "master terms".
-     *
-     * The maps below maintain a master/slave relationship over
-     * transcendental functions (SINE, EXPONENTIAL, PI), where above
-     * f(y) is the master of itself and of f(t).
-     *
-     * This is used for ensuring that the argument y of SINE we process is on
-     * the interval [-pi .. pi], and that exponentials are not applied to
-     * arguments that contain transcendental functions.
-     */
-    std::map<Node, Node> d_trMaster;
-    std::map<Node, std::unordered_set<Node, NodeHashFunction>> d_trSlaves;
-
-    /** concavity region for transcendental functions
-     *
-     * This stores an integer that identifies an interval in
-     * which the current model value for an argument of an
-     * application of a transcendental function resides.
-     *
-     * For exp( x ):
-     *   region #1 is -infty < x < infty
-     * For sin( x ):
-     *   region #0 is pi < x < infty (this is an invalid region)
-     *   region #1 is pi/2 < x <= pi
-     *   region #2 is 0 < x <= pi/2
-     *   region #3 is -pi/2 < x <= 0
-     *   region #4 is -pi < x <= -pi/2
-     *   region #5 is -infty < x <= -pi (this is an invalid region)
-     * All regions not listed above, as well as regions 0 and 5
-     * for SINE are "invalid". We only process applications
-     * of transcendental functions whose arguments have model
-     * values that reside in valid regions.
-     */
-    std::unordered_map<Node, int, NodeHashFunction> d_tf_region;
-    /**
-     * Maps representives of a congruence class to the members of that class.
-     *
-     * In detail, a congruence class is a set of terms of the form
-     *   { f(t1), ..., f(tn) }
-     * such that t1 = ... = tn in the current context. We choose an arbitrary
-     * term among these to be the repesentative of this congruence class.
-     *
-     * Moreover, notice we compute congruence classes only over terms that
-     * are transcendental function applications that are "master terms",
-     * see d_trMaster/d_trSlave.
-     */
-    std::map<Node, std::vector<Node>> d_funcCongClass;
-    /**
-     * A list of all functions for each kind in { EXPONENTIAL, SINE, POW, PI }
-     * that are representives of their congruence class.
-     */
-    std::map<Kind, std::vector<Node>> d_funcMap;
-
-    /** secant points (sorted list) for transcendental functions
-     *
-     * This is used for tangent plane refinements for
-     * transcendental functions. This is the set
-     * "get-previous-secant-points" in "Satisfiability
-     * Modulo Transcendental Functions via Incremental
-     * Linearization" by Cimatti et al., CADE 2017, for
-     * each transcendental function application. We store this set for each
-     * Taylor degree.
-     */
-    std::unordered_map<Node,
-                       std::map<unsigned, std::vector<Node>>,
-                       NodeHashFunction>
-        d_secant_points;
-
-    /** PI
-     *
-     * Note that PI is a (symbolic, non-constant) nullary operator. This is
-     * because its value cannot be computed exactly. We constraint PI to
-     * concrete lower and upper bounds stored in d_pi_bound below.
-     */
-    Node d_pi;
-    /** PI/2 */
-    Node d_pi_2;
-    /** -PI/2 */
-    Node d_pi_neg_2;
-    /** -PI */
-    Node d_pi_neg;
-    /** the concrete lower and upper bounds for PI */
-    Node d_pi_bound[2];
-  };
+                      int concavity);
+
+  Node d_true;
+  Node d_false;
+  Node d_zero;
+  Node d_one;
+  Node d_neg_one;
+
+  /** The inference manager that we push conflicts and lemmas to. */
+  InferenceManager& d_im;
+  /** Reference to the non-linear model object */
+  NlModel& d_model;
+  /** Utility to compute taylor approximations */
+  TaylorGenerator d_taylor;
+
+  /**
+   * Some transcendental functions f(t) are "purified", e.g. we add
+   * t = y ^ f(t) = f(y) where y is a fresh variable. Those that are not
+   * purified we call "master terms".
+   *
+   * The maps below maintain a master/slave relationship over
+   * transcendental functions (SINE, EXPONENTIAL, PI), where above
+   * f(y) is the master of itself and of f(t).
+   *
+   * This is used for ensuring that the argument y of SINE we process is on
+   * the interval [-pi .. pi], and that exponentials are not applied to
+   * arguments that contain transcendental functions.
+   */
+  std::map<Node, Node> d_trMaster;
+  std::map<Node, std::unordered_set<Node, NodeHashFunction>> d_trSlaves;
+
+  /** concavity region for transcendental functions
+   *
+   * This stores an integer that identifies an interval in
+   * which the current model value for an argument of an
+   * application of a transcendental function resides.
+   *
+   * For exp( x ):
+   *   region #1 is -infty < x < infty
+   * For sin( x ):
+   *   region #0 is pi < x < infty (this is an invalid region)
+   *   region #1 is pi/2 < x <= pi
+   *   region #2 is 0 < x <= pi/2
+   *   region #3 is -pi/2 < x <= 0
+   *   region #4 is -pi < x <= -pi/2
+   *   region #5 is -infty < x <= -pi (this is an invalid region)
+   * All regions not listed above, as well as regions 0 and 5
+   * for SINE are "invalid". We only process applications
+   * of transcendental functions whose arguments have model
+   * values that reside in valid regions.
+   */
+  std::unordered_map<Node, int, NodeHashFunction> d_tf_region;
+  /**
+   * Maps representives of a congruence class to the members of that class.
+   *
+   * In detail, a congruence class is a set of terms of the form
+   *   { f(t1), ..., f(tn) }
+   * such that t1 = ... = tn in the current context. We choose an arbitrary
+   * term among these to be the repesentative of this congruence class.
+   *
+   * Moreover, notice we compute congruence classes only over terms that
+   * are transcendental function applications that are "master terms",
+   * see d_trMaster/d_trSlave.
+   */
+  std::map<Node, std::vector<Node>> d_funcCongClass;
+  /**
+   * A list of all functions for each kind in { EXPONENTIAL, SINE, POW, PI }
+   * that are representives of their congruence class.
+   */
+  std::map<Kind, std::vector<Node>> d_funcMap;
+
+  /** secant points (sorted list) for transcendental functions
+   *
+   * This is used for tangent plane refinements for
+   * transcendental functions. This is the set
+   * "get-previous-secant-points" in "Satisfiability
+   * Modulo Transcendental Functions via Incremental
+   * Linearization" by Cimatti et al., CADE 2017, for
+   * each transcendental function application. We store this set for each
+   * Taylor degree.
+   */
+  std::unordered_map<Node,
+                     std::map<unsigned, std::vector<Node>>,
+                     NodeHashFunction>
+      d_secant_points;
+
+  /** PI
+   *
+   * Note that PI is a (symbolic, non-constant) nullary operator. This is
+   * because its value cannot be computed exactly. We constraint PI to
+   * concrete lower and upper bounds stored in d_pi_bound below.
+   */
+  Node d_pi;
+  /** PI/2 */
+  Node d_pi_2;
+  /** -PI/2 */
+  Node d_pi_neg_2;
+  /** -PI */
+  Node d_pi_neg;
+  /** the concrete lower and upper bounds for PI */
+  Node d_pi_bound[2];
+};
 
-}  // namespace transcendental
 }  // namespace transcendental
 }  // namespace nl
 }  // namespace arith
 }  // namespace theory
+}  // namespace CVC4
 
 #endif /* CVC4__THEORY__ARITH__NL__TRANSCENDENTAL__TRANSCENDENTAL_STATE_H */