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());
<< "*** 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]);
}
bool isNeg = csign == -1;
- std::vector<Node> pbounds;
+ ApproximationBounds pbounds;
getPolynomialApproximationBoundForArg(k, c, d, pbounds);
std::vector<Node> bounds;
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]) }
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();
/**
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
* 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
* @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
*
* 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
projects / cvc5.git / commitdiff