Cleanup in transcendental solver, add ApproximationBounds struct. (#5945)

This PR merges some cleanup in the transcendental solver from proof-new.
It adds a new struct ApproximationBounds that replaces an opaque std::vector and does some general refactoring in the TaylorGenerator class, removing dead code and using fixed-width integers.

3 files changed, 121 insertions, 194 deletions

diff --git a/src/theory/arith/nl/transcendental/taylor_generator.cpp b/src/theory/arith/nl/transcendental/taylor_generator.cpp
index 745cba17a..375f1757e 100644
--- a/src/theory/arith/nl/transcendental/taylor_generator.cpp
+++ b/src/theory/arith/nl/transcendental/taylor_generator.cpp
@@ -25,214 +25,138 @@ namespace transcendental {
 
 TaylorGenerator::TaylorGenerator()
     : d_nm(NodeManager::currentNM()),
-      d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType())),
-      d_taylor_real_fv_base(d_nm->mkBoundVar("a", d_nm->realType())),
-      d_taylor_real_fv_base_rem(d_nm->mkBoundVar("b", d_nm->realType()))
+      d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType()))
 {
 }
 
 TNode TaylorGenerator::getTaylorVariable() { return d_taylor_real_fv; }
 
-std::pair<Node, Node> TaylorGenerator::getTaylor(TNode fa, std::uint64_t n)
+std::pair<Node, Node> TaylorGenerator::getTaylor(Kind k, std::uint64_t n)
 {
-  NodeManager* nm = NodeManager::currentNM();
-  Node d_zero = nm->mkConst(Rational(0));
-
   Assert(n > 0);
-  Node fac;  // what term we cache for fa
-  if (fa[0] == d_zero)
-  {
-    // optimization : simpler to compute (x-fa[0])^n if we are centered around 0
-    fac = fa;
-  }
-  else
+  // check if we have already computed this Taylor series
+  auto itt = d_taylor_terms[k].find(n);
+  if (itt != d_taylor_terms[k].end())
   {
-    // otherwise we use a standard factor a in (x-a)^n
-    fac = nm->mkNode(fa.getKind(), d_taylor_real_fv_base);
+    return itt->second;
   }
-  Node taylor_rem;
-  Node taylor_sum;
-  // check if we have already computed this Taylor series
-  auto itt = s_taylor_sum[fac].find(n);
-  if (itt == s_taylor_sum[fac].end())
+
+  NodeManager* nm = NodeManager::currentNM();
+
+  // the current factorial `counter!`
+  Integer factorial = 1;
+  // the current variable power `x^counter`
+  Node varpow = nm->mkConst(Rational(1));
+  std::vector<Node> sum;
+  for (std::uint64_t counter = 1; counter <= n; ++counter)
   {
-    Node i_exp_base;
-    if (fa[0] == d_zero)
+    if (k == Kind::EXPONENTIAL)
     {
-      i_exp_base = d_taylor_real_fv;
+      // Maclaurin series for exponential:
+      //   \sum_{n=0}^\infty x^n / n!
+      sum.push_back(
+          nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
     }
-    else
+    else if (k == Kind::SINE)
     {
-      i_exp_base = Rewriter::rewrite(
-          nm->mkNode(Kind::MINUS, d_taylor_real_fv, d_taylor_real_fv_base));
-    }
-    Node i_derv = fac;
-    Node i_fact = nm->mkConst(Rational(1));
-    Node i_exp = i_fact;
-    int i_derv_status = 0;
-    unsigned counter = 0;
-    std::vector<Node> sum;
-    do
-    {
-      counter++;
-      if (fa.getKind() == Kind::EXPONENTIAL)
-      {
-        // unchanged
-      }
-      else if (fa.getKind() == Kind::SINE)
+      // Maclaurin series for exponential:
+      //   \sum_{n=0}^\infty (-1)^n / (2n+1)! * x^(2n+1)
+      if (counter % 2 == 0)
       {
-        if (i_derv_status % 2 == 1)
-        {
-          Node pi = nm->mkNullaryOperator(nm->realType(), Kind::PI);
-          Node pi_2 = Rewriter::rewrite(nm->mkNode(
-              Kind::MULT, pi, nm->mkConst(Rational(1) / Rational(2))));
-
-          Node arg = nm->mkNode(Kind::PLUS, pi_2, d_taylor_real_fv_base);
-          i_derv = nm->mkNode(Kind::SINE, arg);
-        }
-        else
-        {
-          i_derv = fa;
-        }
-        if (i_derv_status >= 2)
-        {
-          i_derv = nm->mkNode(Kind::MINUS, d_zero, i_derv);
-        }
-        i_derv = Rewriter::rewrite(i_derv);
-        i_derv_status = i_derv_status == 3 ? 0 : i_derv_status + 1;
+        int sign = (counter % 4 == 0 ? -1 : 1);
+        sum.push_back(nm->mkNode(Kind::MULT,
+                                 nm->mkNode(Kind::DIVISION,
+                                            nm->mkConst<Rational>(sign),
+                                            nm->mkConst<Rational>(factorial)),
+                                 varpow));
       }
-      if (counter == (n + 1))
-      {
-        TNode x = d_taylor_real_fv_base;
-        i_derv = i_derv.substitute(x, d_taylor_real_fv_base_rem);
-      }
-      Node curr = nm->mkNode(
-          Kind::MULT, nm->mkNode(Kind::DIVISION, i_derv, i_fact), i_exp);
-      if (counter == (n + 1))
-      {
-        taylor_rem = curr;
-      }
-      else
-      {
-        sum.push_back(curr);
-        i_fact = Rewriter::rewrite(
-            nm->mkNode(Kind::MULT, nm->mkConst(Rational(counter)), i_fact));
-        i_exp = Rewriter::rewrite(nm->mkNode(Kind::MULT, i_exp_base, i_exp));
-      }
-    } while (counter <= n);
-    taylor_sum = sum.size() == 1 ? sum[0] : nm->mkNode(Kind::PLUS, sum);
-
-    if (fac[0] != d_taylor_real_fv_base)
-    {
-      TNode x = d_taylor_real_fv_base;
-      taylor_sum = taylor_sum.substitute(x, fac[0]);
     }
-
-    // cache
-    s_taylor_sum[fac][n] = taylor_sum;
-    d_taylor_rem[fac][n] = taylor_rem;
-  }
-  else
-  {
-    taylor_sum = itt->second;
-    Assert(d_taylor_rem[fac].find(n) != d_taylor_rem[fac].end());
-    taylor_rem = d_taylor_rem[fac][n];
+    factorial *= counter;
+    varpow =
+        Rewriter::rewrite(nm->mkNode(Kind::MULT, d_taylor_real_fv, varpow));
   }
+  Node taylor_sum =
+      Rewriter::rewrite(sum.size() == 1 ? sum[0] : nm->mkNode(Kind::PLUS, sum));
+  Node taylor_rem = Rewriter::rewrite(
+      nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
 
-  // must substitute for the argument if we were using a different lookup
-  if (fa[0] != fac[0])
-  {
-    TNode x = d_taylor_real_fv_base;
-    taylor_sum = taylor_sum.substitute(x, fa[0]);
-  }
-  return std::pair<Node, Node>(taylor_sum, taylor_rem);
+  auto res = std::make_pair(taylor_sum, taylor_rem);
+
+  // put result in cache
+  d_taylor_terms[k][n] = res;
+
+  return res;
 }
 
 void TaylorGenerator::getPolynomialApproximationBounds(
-    Kind k, unsigned d, std::vector<Node>& pbounds)
+    Kind k, std::uint64_t d, ApproximationBounds& pbounds)
 {
-  if (d_poly_bounds[k][d].empty())
+  auto it = d_poly_bounds[k].find(d);
+  if (it == d_poly_bounds[k].end())
   {
     NodeManager* nm = NodeManager::currentNM();
-    Node tft = nm->mkNode(k, nm->mkConst(Rational(0)));
     // n is the Taylor degree we are currently considering
-    unsigned n = 2 * d;
+    std::uint64_t n = 2 * d;
     // n must be even
-    std::pair<Node, Node> taylor = getTaylor(tft, n);
-    Trace("nl-ext-tftp-debug2")
-        << "Taylor for " << k << " is : " << taylor.first << std::endl;
-    Node taylor_sum = Rewriter::rewrite(taylor.first);
-    Trace("nl-ext-tftp-debug2")
-        << "Taylor for " << k << " is (post-rewrite) : " << taylor_sum
-        << std::endl;
-    Assert(taylor.second.getKind() == Kind::MULT);
-    Assert(taylor.second.getNumChildren() == 2);
-    Assert(taylor.second[0].getKind() == Kind::DIVISION);
-    Trace("nl-ext-tftp-debug2")
-        << "Taylor remainder for " << k << " is " << taylor.second << std::endl;
+    std::pair<Node, Node> taylor = getTaylor(k, n);
+    Node taylor_sum = taylor.first;
     // ru is x^{n+1}/(n+1)!
-    Node ru = nm->mkNode(Kind::DIVISION, taylor.second[1], taylor.second[0][1]);
-    ru = Rewriter::rewrite(ru);
-    Trace("nl-ext-tftp-debug2")
-        << "Taylor remainder factor is (post-rewrite) : " << ru << std::endl;
+    Node ru = taylor.second;
+    Trace("nl-trans") << "Taylor for " << k << " is : " << taylor.first
+                      << std::endl;
+    Trace("nl-trans") << "Taylor remainder for " << k << " is " << taylor.second
+                      << std::endl;
     if (k == Kind::EXPONENTIAL)
     {
-      pbounds.push_back(taylor_sum);
-      pbounds.push_back(taylor_sum);
-      pbounds.push_back(Rewriter::rewrite(
+      pbounds.d_lower = taylor_sum;
+      pbounds.d_upperNeg =
+          Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
+      pbounds.d_upperPos = Rewriter::rewrite(
           nm->mkNode(Kind::MULT,
                      taylor_sum,
-                     nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru))));
-      pbounds.push_back(
-          Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru)));
+                     nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru)));
     }
     else
     {
       Assert(k == Kind::SINE);
       Node l = Rewriter::rewrite(nm->mkNode(Kind::MINUS, taylor_sum, ru));
       Node u = Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
-      pbounds.push_back(l);
-      pbounds.push_back(l);
-      pbounds.push_back(u);
-      pbounds.push_back(u);
+      pbounds.d_lower = l;
+      pbounds.d_upperNeg = u;
+      pbounds.d_upperPos = u;
     }
-    Trace("nl-ext-tf-tplanes")
-        << "Polynomial approximation for " << k << " is: " << std::endl;
-    Trace("nl-ext-tf-tplanes") << " Lower (pos): " << pbounds[0] << std::endl;
-    Trace("nl-ext-tf-tplanes") << " Upper (pos): " << pbounds[2] << std::endl;
-    Trace("nl-ext-tf-tplanes") << " Lower (neg): " << pbounds[1] << std::endl;
-    Trace("nl-ext-tf-tplanes") << " Upper (neg): " << pbounds[3] << std::endl;
-    d_poly_bounds[k][d].insert(
-        d_poly_bounds[k][d].end(), pbounds.begin(), pbounds.end());
+    Trace("nl-trans") << "Polynomial approximation for " << k
+                      << " is: " << std::endl;
+    Trace("nl-trans") << " Lower: " << pbounds.d_lower << std::endl;
+    Trace("nl-trans") << " Upper (neg): " << pbounds.d_upperNeg << std::endl;
+    Trace("nl-trans") << " Upper (pos): " << pbounds.d_upperPos << std::endl;
+    d_poly_bounds[k].emplace(d, pbounds);
   }
   else
   {
-    pbounds.insert(
-        pbounds.end(), d_poly_bounds[k][d].begin(), d_poly_bounds[k][d].end());
+    pbounds = it->second;
   }
 }
 
-unsigned TaylorGenerator::getPolynomialApproximationBoundForArg(
-    Kind k, Node c, unsigned d, std::vector<Node>& pbounds)
+std::uint64_t TaylorGenerator::getPolynomialApproximationBoundForArg(
+    Kind k, Node c, std::uint64_t d, ApproximationBounds& pbounds)
 {
   getPolynomialApproximationBounds(k, d, pbounds);
+  Trace("nl-trans") << "c = " << c << std::endl;
   Assert(c.isConst());
   if (k == Kind::EXPONENTIAL && c.getConst<Rational>().sgn() == 1)
   {
-    NodeManager* nm = NodeManager::currentNM();
-    Node tft = nm->mkNode(k, nm->mkConst(Rational(0)));
     bool success = false;
-    unsigned ds = d;
+    std::uint64_t ds = d;
     TNode ttrf = getTaylorVariable();
     TNode tc = c;
     do
     {
       success = true;
-      unsigned n = 2 * ds;
-      std::pair<Node, Node> taylor = getTaylor(tft, n);
+      std::uint64_t n = 2 * ds;
+      std::pair<Node, Node> taylor = getTaylor(k, n);
       // check that 1-c^{n+1}/(n+1)! > 0
-      Node ru =
-          nm->mkNode(Kind::DIVISION, taylor.second[1], taylor.second[0][1]);
+      Node ru = taylor.second;
       Node rus = ru.substitute(ttrf, tc);
       rus = Rewriter::rewrite(rus);
       Assert(rus.isConst());
@@ -248,16 +172,18 @@ unsigned TaylorGenerator::getPolynomialApproximationBoundForArg(
           << "*** Increase Taylor bound to " << ds << " > " << d << " for ("
           << k << " " << c << ")" << std::endl;
       // must use sound upper bound
-      std::vector<Node> pboundss;
+      ApproximationBounds pboundss;
       getPolynomialApproximationBounds(k, ds, pboundss);
-      pbounds[2] = pboundss[2];
+      pbounds.d_upperPos = pboundss.d_upperPos;
     }
     return ds;
   }
   return d;
 }
 
-std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlModel& model)
+std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf,
+                                                        std::uint64_t d,
+                                                        NlModel& model)
 {
   // compute the model value of the argument
   Node c = model.computeAbstractModelValue(tf[0]);
@@ -279,7 +205,7 @@ std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlM
   }
   bool isNeg = csign == -1;
 
-  std::vector<Node> pbounds;
+  ApproximationBounds pbounds;
   getPolynomialApproximationBoundForArg(k, c, d, pbounds);
 
   std::vector<Node> bounds;
@@ -287,8 +213,8 @@ std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlM
   TNode tfs = tf[0];
   for (unsigned d2 = 0; d2 < 2; d2++)
   {
-    int index = d2 == 0 ? (isNeg ? 1 : 0) : (isNeg ? 3 : 2);
-    Node pab = pbounds[index];
+    Node pab = (d2 == 0 ? pbounds.d_lower
+                        : (isNeg ? pbounds.d_upperNeg : pbounds.d_upperPos));
     if (!pab.isNull())
     {
       // { x -> M_A(tf[0]) }
diff --git a/src/theory/arith/nl/transcendental/taylor_generator.h b/src/theory/arith/nl/transcendental/taylor_generator.h
index 261583ca9..ccdc80e0f 100644
--- a/src/theory/arith/nl/transcendental/taylor_generator.h
+++ b/src/theory/arith/nl/transcendental/taylor_generator.h
@@ -27,6 +27,17 @@ namespace transcendental {
 class TaylorGenerator
 {
  public:
+  /** Stores the approximation bounds for transcendental functions */
+  struct ApproximationBounds
+  {
+    /** Lower bound */
+    Node d_lower;
+    /** Upper bound for negative values */
+    Node d_upperNeg;
+    /** Upper bound for positive values */
+    Node d_upperPos;
+  };
+
   TaylorGenerator();
 
   /**
@@ -35,23 +46,17 @@ class TaylorGenerator
   TNode getTaylorVariable();
 
   /**
-   * Get Taylor series of degree n for function fa centered around point fa[0].
+   * Get Taylor series of degree n for function fa centered around zero.
    *
-   * Return value is ( P_{n,f(a)}( x ), R_{n+1,f(a)}( x ) ) where
+   * Return value is ( P_{n,f(0)}( x ), R_{n+1,f(0)}( x ) ) where
    * the first part of the pair is the Taylor series expansion :
-   *    P_{n,f(a)}( x ) = sum_{i=0}^n (f^i( a )/i!)*(x-a)^i
+   *    P_{n,f(0)}( x ) = sum_{i=0}^n (f^i(0)/i!)*x^i
    * and the second part of the pair is the Taylor series remainder :
-   *    R_{n+1,f(a),b}( x ) = (f^{n+1}( b )/(n+1)!)*(x-a)^{n+1}
+   *    R_{n+1,f(0)}( x ) = x^{n+1}/(n+1)!
    *
-   * The above values are cached for each (f,n) for a fixed variable "a".
-   * To compute the Taylor series for fa, we compute the Taylor series
-   *   for ( fa.getKind(), n ) then substitute { a -> fa[0] } if fa[0]!=0.
-   * We compute P_{n,f(0)}( x )/R_{n+1,f(0),b}( x ) for ( fa.getKind(), n )
-   *   if fa[0]=0.
-   * In the latter case, note we compute the exponential x^{n+1}
-   * instead of (x-a)^{n+1}, which can be done faster.
+   * The above values are cached for each (f,n).
    */
-  std::pair<Node, Node> getTaylor(TNode fa, std::uint64_t n);
+  std::pair<Node, Node> getTaylor(Kind k, std::uint64_t n);
 
   /**
    * polynomial approximation bounds
@@ -68,8 +73,8 @@ class TaylorGenerator
    * given <k>( c ), use the function below.
    */
   void getPolynomialApproximationBounds(Kind k,
-                                        unsigned d,
-                                        std::vector<Node>& pbounds);
+                                        std::uint64_t d,
+                                        ApproximationBounds& pbounds);
 
   /**
    * polynomial approximation bounds
@@ -82,10 +87,8 @@ class TaylorGenerator
    * @return the actual degree of the polynomial approximations (which may be
    * larger than d).
    */
-  unsigned getPolynomialApproximationBoundForArg(Kind k,
-                                                 Node c,
-                                                 unsigned d,
-                                                 std::vector<Node>& pbounds);
+  std::uint64_t getPolynomialApproximationBoundForArg(
+      Kind k, Node c, std::uint64_t d, ApproximationBounds& pbounds);
 
   /** get transcendental function model bounds
    *
@@ -93,22 +96,20 @@ class TaylorGenerator
    * function application tf based on Taylor of degree 2*d, which is dependent
    * on the model value of its argument.
    */
-  std::pair<Node, Node> getTfModelBounds(Node tf, unsigned d, NlModel& model);
+  std::pair<Node, Node> getTfModelBounds(Node tf,
+                                         std::uint64_t d,
+                                         NlModel& model);
 
  private:
   NodeManager* d_nm;
   const Node d_taylor_real_fv;
-  const Node d_taylor_real_fv_base;
-  const Node d_taylor_real_fv_base_rem;
-  std::unordered_map<Node,
-                     std::unordered_map<std::uint64_t, Node>,
-                     NodeHashFunction>
-      s_taylor_sum;
-  std::unordered_map<Node,
-                     std::unordered_map<std::uint64_t, Node>,
-                     NodeHashFunction>
-      d_taylor_rem;
-  std::map<Kind, std::map<unsigned, std::vector<Node>>> d_poly_bounds;
+
+  /**
+   * For every kind (EXP or SINE) and every degree we store the taylor series up
+   * to this degree and the next term in the series.
+   */
+  std::map<Kind, std::map<std::uint64_t, std::pair<Node, Node>>> d_taylor_terms;
+  std::map<Kind, std::map<std::uint64_t, ApproximationBounds>> d_poly_bounds;
 };
 
 }  // namespace transcendental
diff --git a/src/theory/arith/nl/transcendental/transcendental_solver.cpp b/src/theory/arith/nl/transcendental/transcendental_solver.cpp
index 54f430372..ee422ebb6 100644
--- a/src/theory/arith/nl/transcendental/transcendental_solver.cpp
+++ b/src/theory/arith/nl/transcendental/transcendental_solver.cpp
@@ -267,13 +267,13 @@ bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d)
   // Figure 3: P_l, P_u
   // mapped to for signs of c
   std::map<int, Node> poly_approx_bounds[2];
-  std::vector<Node> pbounds;
-  unsigned actual_d =
+  TaylorGenerator::ApproximationBounds pbounds;
+  std::uint64_t actual_d =
       d_tstate.d_taylor.getPolynomialApproximationBoundForArg(k, c, d, pbounds);
-  poly_approx_bounds[0][1] = pbounds[0];
-  poly_approx_bounds[0][-1] = pbounds[1];
-  poly_approx_bounds[1][1] = pbounds[2];
-  poly_approx_bounds[1][-1] = pbounds[3];
+  poly_approx_bounds[0][1] = pbounds.d_lower;
+  poly_approx_bounds[0][-1] = pbounds.d_lower;
+  poly_approx_bounds[1][1] = pbounds.d_upperPos;
+  poly_approx_bounds[1][-1] = pbounds.d_upperNeg;
 
   // Figure 3 : v
   Node v = d_tstate.d_model.computeAbstractModelValue(tf);