Head

GCC Code Coverage Report

Directory:	.		Exec	Total	Coverage
File:	src/theory/arith/nl/transcendental/transcendental_state.cpp	Lines:	182	194	93.8 %
Date:	2021-09-16	Branches:	498	1072	46.5 %


/******************************************************************************
 * Top contributors (to current version):
 *   Gereon Kremer, Andrew Reynolds
 *
 * This file is part of the cvc5 project.
 *
 * Copyright (c) 2009-2021 by the authors listed in the file AUTHORS
 * in the top-level source directory and their institutional affiliations.
 * All rights reserved.  See the file COPYING in the top-level source
 * directory for licensing information.
 * ****************************************************************************
 *
 * Implementation of solver for handling transcendental functions.
 */

#include "theory/arith/nl/transcendental/transcendental_state.h"

#include "proof/proof.h"
#include "theory/arith/arith_utilities.h"
#include "theory/arith/inference_manager.h"
#include "theory/arith/nl/nl_model.h"
#include "theory/arith/nl/transcendental/taylor_generator.h"
#include "theory/rewriter.h"

namespace cvc5 {
namespace theory {
namespace arith {
namespace nl {
namespace transcendental {

TranscendentalState::TranscendentalState(InferenceManager& im,
                                         NlModel& model,
                                         Env& env)
    : d_im(im), d_model(model), d_env(env)
{
  d_true = NodeManager::currentNM()->mkConst(true);
  d_false = NodeManager::currentNM()->mkConst(false);
  d_zero = NodeManager::currentNM()->mkConst(Rational(0));
  d_one = NodeManager::currentNM()->mkConst(Rational(1));
  d_neg_one = NodeManager::currentNM()->mkConst(Rational(-1));
  if (d_env.isTheoryProofProducing())
  {
    d_proof.reset(new CDProofSet<CDProof>(
        d_env.getProofNodeManager(), d_env.getUserContext(), "nl-trans"));
    d_proofChecker.reset(new TranscendentalProofRuleChecker());
    d_proofChecker->registerTo(d_env.getProofNodeManager()->getChecker());
  }
}

bool TranscendentalState::isProofEnabled() const
{
  return d_proof.get() != nullptr;
}

CDProof* TranscendentalState::getProof()
{
  Assert(isProofEnabled());
  return d_proof->allocateProof(d_env.getUserContext());
}

void TranscendentalState::init(const std::vector<Node>& xts,
                               std::vector<Node>& needsMaster)
{
  d_funcCongClass.clear();
  d_funcMap.clear();
  d_tf_region.clear();

  bool needPi = false;
  // for computing congruence
  std::map<Kind, ArgTrie> argTrie;
  for (std::size_t i = 0, xsize = xts.size(); i < xsize; ++i)
  {
    // Ignore if it is not a transcendental
    if (!isTranscendentalKind(xts[i].getKind()))
    {
      continue;
    }
    Node a = xts[i];
    Kind ak = a.getKind();
    bool consider = true;
    // if we've already computed master for a
    if (d_trMaster.find(a) != d_trMaster.end())
    {
      // a master has at least one slave
      consider = (d_trSlaves.find(a) != d_trSlaves.end());
    }
    else
    {
      if (ak == Kind::SINE)
      {
        // always not a master
        consider = false;
      }
      else
      {
        for (const Node& ac : a)
        {
          if (isTranscendentalKind(ac.getKind()))
          {
            consider = false;
            break;
          }
        }
      }
      if (!consider)
      {
        // wait to assign a master below
        needsMaster.push_back(a);
      }
      else
      {
        d_trMaster[a] = a;
        d_trSlaves[a].insert(a);
      }
    }
    if (ak == Kind::EXPONENTIAL || ak == Kind::SINE)
    {
      needPi = needPi || (ak == Kind::SINE);
      // if we didn't indicate that it should be purified above
      if (consider)
      {
        ensureCongruence(a, argTrie);
      }
    }
    else if (ak == Kind::PI)
    {
      Assert(consider);
      needPi = true;
      d_funcMap[ak].push_back(a);
      d_funcCongClass[a].push_back(a);
    }
  }
  // initialize pi if necessary
  if (needPi && d_pi.isNull())
  {
    mkPi();
    getCurrentPiBounds();
  }

  if (Trace.isOn("nl-ext-mv"))
  {
    Trace("nl-ext-mv") << "Arguments of trancendental functions : "
                       << std::endl;
    for (std::pair<const Kind, std::vector<Node> >& tfl : d_funcMap)
    {
      Kind k = tfl.first;
      if (k == Kind::SINE || k == Kind::EXPONENTIAL)
      {
        for (const Node& tf : tfl.second)
        {
          Node v = tf[0];
          d_model.computeConcreteModelValue(v);
          d_model.computeAbstractModelValue(v);
          d_model.printModelValue("nl-ext-mv", v);
        }
      }
    }
  }
}

void TranscendentalState::ensureCongruence(TNode a,
                                           std::map<Kind, ArgTrie>& argTrie)
{
  NodeManager* nm = NodeManager::currentNM();
  std::vector<Node> repList;
  for (const Node& ac : a)
  {
    Node r = d_model.computeConcreteModelValue(ac);
    repList.push_back(r);
  }
  Node aa = argTrie[a.getKind()].add(a, repList);
  if (aa != a)
  {
    // apply congruence to pairs of terms that are disequal and congruent
    Assert(aa.getNumChildren() == a.getNumChildren());
    Node mvaa = d_model.computeAbstractModelValue(a);
    Node mvaaa = d_model.computeAbstractModelValue(aa);
    if (mvaa != mvaaa)
    {
      std::vector<Node> exp;
      for (unsigned j = 0, size = a.getNumChildren(); j < size; j++)
      {
        exp.push_back(a[j].eqNode(aa[j]));
      }
      Node expn = exp.size() == 1 ? exp[0] : nm->mkNode(Kind::AND, exp);
      Node cong_lemma = expn.impNode(a.eqNode(aa));
      d_im.addPendingLemma(cong_lemma, InferenceId::ARITH_NL_CONGRUENCE);
    }
  }
  else
  {
    // new representative of congruence class
    d_funcMap[a.getKind()].push_back(a);
  }
  // add to congruence class
  d_funcCongClass[aa].push_back(a);
}

void TranscendentalState::mkPi()
{
  NodeManager* nm = NodeManager::currentNM();
  if (d_pi.isNull())
  {
    d_pi = nm->mkNullaryOperator(nm->realType(), Kind::PI);
    d_pi_2 = Rewriter::rewrite(
        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(1) / Rational(2))));
    d_pi_neg_2 = Rewriter::rewrite(
        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(-1) / Rational(2))));
    d_pi_neg = Rewriter::rewrite(
        nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(-1))));
    // initialize bounds
    d_pi_bound[0] = nm->mkConst(Rational(103993) / Rational(33102));
    d_pi_bound[1] = nm->mkConst(Rational(104348) / Rational(33215));
  }
}

void TranscendentalState::getCurrentPiBounds()
{
  NodeManager* nm = NodeManager::currentNM();
  Node pi_lem = nm->mkNode(Kind::AND,
                           nm->mkNode(Kind::GEQ, d_pi, d_pi_bound[0]),
                           nm->mkNode(Kind::LEQ, d_pi, d_pi_bound[1]));
  CDProof* proof = nullptr;
  if (isProofEnabled())
  {
    proof = getProof();
    proof->addStep(
        pi_lem, PfRule::ARITH_TRANS_PI, {}, {d_pi_bound[0], d_pi_bound[1]});
  }
  d_im.addPendingLemma(pi_lem, InferenceId::ARITH_NL_T_PI_BOUND, proof);
}

std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e,
                                                                  TNode center,
                                                                  unsigned d)
{
  // bounds are the minimum and maximum previous secant points
  // should not repeat secant points: secant lemmas should suffice to
  // rule out previous assignment
  Assert(std::find(
             d_secant_points[e][d].begin(), d_secant_points[e][d].end(), center)
         == d_secant_points[e][d].end());
  // Insert into the (temporary) vector. We do not update this vector
  // until we are sure this secant plane lemma has been processed. We do
  // this by mapping the lemma to a side effect below.
  std::vector<Node> spoints = d_secant_points[e][d];
  spoints.push_back(center);

  // sort
  sortByNlModel(spoints.begin(), spoints.end(), &d_model);
  // get the resulting index of c
  unsigned index =
      std::find(spoints.begin(), spoints.end(), center) - spoints.begin();

  // bounds are the next closest upper/lower bound values
  return {index > 0 ? spoints[index - 1] : Node(),
          index < spoints.size() - 1 ? spoints[index + 1] : Node()};
}

Node TranscendentalState::mkSecantPlane(
    TNode arg, TNode lower, TNode upper, TNode lval, TNode uval)
{
  NodeManager* nm = NodeManager::currentNM();
  // Figure 3: S_l( x ), S_u( x ) for s = 0,1
  Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, lower, upper));
  Assert(rcoeff_n.isConst());
  Rational rcoeff = rcoeff_n.getConst<Rational>();
  Assert(rcoeff.sgn() != 0);
  Node res =
      nm->mkNode(Kind::PLUS,
                 lval,
                 nm->mkNode(Kind::MULT,
                            nm->mkNode(Kind::DIVISION,
                                       nm->mkNode(Kind::MINUS, lval, uval),
                                       nm->mkNode(Kind::MINUS, lower, upper)),
                            nm->mkNode(Kind::MINUS, arg, lower)));
  Trace("nl-trans") << "Creating secant plane for transcendental function of "
                    << arg << std::endl;
  Trace("nl-trans") << "\tfrom ( " << lower << " ; " << lval << " ) to ( "
                    << upper << " ; " << uval << " )" << std::endl;
  Trace("nl-trans") << "\t" << res << std::endl;
  Trace("nl-trans") << "\trewritten: " << Rewriter::rewrite(res) << std::endl;
  return res;
}

NlLemma TranscendentalState::mkSecantLemma(TNode lower,
                                           TNode upper,
                                           TNode lapprox,
                                           TNode uapprox,
                                           int csign,
                                           Convexity convexity,
                                           TNode tf,
                                           TNode splane,
                                           unsigned actual_d)
{
  NodeManager* nm = NodeManager::currentNM();
  // With respect to Figure 3, this is slightly different.
  // In particular, we chose b to be the model value of bounds[s],
  // which is a constant although bounds[s] may not be (e.g. if it
  // contains PI).
  // To ensure that c...b does not cross an inflection point,
  // we guard with the symbolic version of bounds[s].
  // This leads to lemmas e.g. of this form:
  //   ( c <= x <= PI/2 ) => ( sin(x) < ( P( b ) - P( c ) )*( x -
  //   b ) + P( b ) )
  // where b = (PI/2)^M, the current value of PI/2 in the model.
  // This is sound since we are guarded by the symbolic
  // representation of PI/2.
  Node antec_n = nm->mkNode(Kind::AND,
                            nm->mkNode(Kind::GEQ, tf[0], lower),
                            nm->mkNode(Kind::LEQ, tf[0], upper));
  Trace("nl-trans") << "Bound for secant plane: " << lower << " <= " << tf[0]
                    << " <= " << upper << std::endl;
  Trace("nl-trans") << "\t" << antec_n << std::endl;
  // Convex: actual value is below the secant
  // Concave: actual value is above the secant
  Node lem = nm->mkNode(
      Kind::IMPLIES,
      antec_n,
      nm->mkNode(
          convexity == Convexity::CONVEX ? Kind::LEQ : Kind::GEQ, tf, splane));
  Trace("nl-trans-lemma") << "*** Secant plane lemma (pre-rewrite) : " << lem
                          << std::endl;
  lem = Rewriter::rewrite(lem);
  Trace("nl-trans-lemma") << "*** Secant plane lemma : " << lem << std::endl;
  Assert(d_model.computeAbstractModelValue(lem) == d_false);
  CDProof* proof = nullptr;
  if (isProofEnabled())
  {
    proof = getProof();
    if (tf.getKind() == Kind::EXPONENTIAL)
    {
      if (csign == 1)
      {
        proof->addStep(
            lem,
            PfRule::ARITH_TRANS_EXP_APPROX_ABOVE_POS,
            {},
            {nm->mkConst<Rational>(2 * actual_d), tf[0], lower, upper});
      }
      else
      {
        proof->addStep(
            lem,
            PfRule::ARITH_TRANS_EXP_APPROX_ABOVE_NEG,
            {},
            {nm->mkConst<Rational>(2 * actual_d), tf[0], lower, upper});
      }
    }
    else if (tf.getKind() == Kind::SINE)
    {
      if (convexity == Convexity::CONCAVE)
      {
        proof->addStep(lem,
                       PfRule::ARITH_TRANS_SINE_APPROX_BELOW_POS,
                       {},
                       {nm->mkConst<Rational>(2 * actual_d),
                        tf[0],
                        lower,
                        upper,
                        lapprox,
                        uapprox

                       });
      }
      else
      {
        proof->addStep(lem,
                       PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_NEG,
                       {},
                       {nm->mkConst<Rational>(2 * actual_d),
                        tf[0],
                        lower,
                        upper,
                        lapprox,
                        uapprox});
      }
    }
  }
  return NlLemma(
      InferenceId::ARITH_NL_T_SECANT, lem, LemmaProperty::NONE, proof);
}

void TranscendentalState::doSecantLemmas(const std::pair<Node, Node>& bounds,
                                         TNode poly_approx,
                                         TNode center,
                                         TNode cval,
                                         TNode tf,
                                         Convexity convexity,
                                         unsigned d,
                                         unsigned actual_d)
{
  int csign = center.getConst<Rational>().sgn();
  Trace("nl-trans") << "...do secant lemma with center " << center << " val "
                    << cval << " sign " << csign << std::endl;
  Trace("nl-trans") << "...secant bounds are : " << bounds.first << " ... "
                    << bounds.second << std::endl;
  // take the model value of lower (since may contain PI)
  // Make secant from bounds.first to center
  Node lower = d_model.computeAbstractModelValue(bounds.first);
  Trace("nl-trans") << "...model value of bound " << bounds.first << " is "
                    << lower << std::endl;
  if (lower != center)
  {
    // Figure 3 : P(l), P(u), for s = 0
    Node lval = Rewriter::rewrite(
        poly_approx.substitute(d_taylor.getTaylorVariable(), lower));
    Node splane = mkSecantPlane(tf[0], lower, center, lval, cval);
    NlLemma nlem = mkSecantLemma(
        lower, center, lval, cval, csign, convexity, tf, splane, actual_d);
    // The side effect says that if lem is added, then we should add the
    // secant point c for (tf,d).
    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center));
    d_im.addPendingLemma(nlem, true);
  }

  // take the model value of upper (since may contain PI)
  // Make secant from center to bounds.second
  Node upper = d_model.computeAbstractModelValue(bounds.second);
  Trace("nl-trans") << "...model value of bound " << bounds.second << " is "
                    << upper << std::endl;
  if (center != upper)
  {
    // Figure 3 : P(l), P(u), for s = 1
    Node uval = Rewriter::rewrite(
        poly_approx.substitute(d_taylor.getTaylorVariable(), upper));
    Node splane = mkSecantPlane(tf[0], center, upper, cval, uval);
    NlLemma nlem = mkSecantLemma(
        center, upper, cval, uval, csign, convexity, tf, splane, actual_d);
    // The side effect says that if lem is added, then we should add the
    // secant point c for (tf,d).
    nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center));
    d_im.addPendingLemma(nlem, true);
  }
}

}  // namespace transcendental
}  // namespace nl
}  // namespace arith
}  // namespace theory
}  // namespace cvc5


Generated by: GCOVR (Version 3.2)

Line	Exec	Source
1		/******************************************************************************
2		* Top contributors (to current version):
3		* Gereon Kremer, Andrew Reynolds
4		*
5		* This file is part of the cvc5 project.
6		*
7		* Copyright (c) 2009-2021 by the authors listed in the file AUTHORS
8		* in the top-level source directory and their institutional affiliations.
9		* All rights reserved. See the file COPYING in the top-level source
10		* directory for licensing information.
11		* ****************************************************************************
12		*
13		* Implementation of solver for handling transcendental functions.
14		*/
15
16		#include "theory/arith/nl/transcendental/transcendental_state.h"
17
18		#include "proof/proof.h"
19		#include "theory/arith/arith_utilities.h"
20		#include "theory/arith/inference_manager.h"
21		#include "theory/arith/nl/nl_model.h"
22		#include "theory/arith/nl/transcendental/taylor_generator.h"
23		#include "theory/rewriter.h"
24
25		namespace cvc5 {
26		namespace theory {
27		namespace arith {
28		namespace nl {
29		namespace transcendental {
30
31	5203	TranscendentalState::TranscendentalState(InferenceManager& im,
32		NlModel& model,
33	5203	Env& env)
34	5203	: d_im(im), d_model(model), d_env(env)
35		{
36	5203	d_true = NodeManager::currentNM()->mkConst(true);
37	5203	d_false = NodeManager::currentNM()->mkConst(false);
38	5203	d_zero = NodeManager::currentNM()->mkConst(Rational(0));
39	5203	d_one = NodeManager::currentNM()->mkConst(Rational(1));
40	5203	d_neg_one = NodeManager::currentNM()->mkConst(Rational(-1));
41	5203	if (d_env.isTheoryProofProducing())
42		{
43	1779	d_proof.reset(new CDProofSet<CDProof>(
44	1186	d_env.getProofNodeManager(), d_env.getUserContext(), "nl-trans"));
45	593	d_proofChecker.reset(new TranscendentalProofRuleChecker());
46	593	d_proofChecker->registerTo(d_env.getProofNodeManager()->getChecker());
47		}
48	5203	}
49
50	2176	bool TranscendentalState::isProofEnabled() const
51		{
52	2176	return d_proof.get() != nullptr;
53		}
54
55	315	CDProof* TranscendentalState::getProof()
56		{
57	315	Assert(isProofEnabled());
58	315	return d_proof->allocateProof(d_env.getUserContext());
59		}
60
61	3461	void TranscendentalState::init(const std::vector<Node>& xts,
62		std::vector<Node>& needsMaster)
63		{
64	3461	d_funcCongClass.clear();
65	3461	d_funcMap.clear();
66	3461	d_tf_region.clear();
67
68	3461	bool needPi = false;
69		// for computing congruence
70	6922	std::map<Kind, ArgTrie> argTrie;
71	30204	for (std::size_t i = 0, xsize = xts.size(); i < xsize; ++i)
72		{
73		// Ignore if it is not a transcendental
74	26743	if (!isTranscendentalKind(xts[i].getKind()))
75		{
76	19710	continue;
77		}
78	14066	Node a = xts[i];
79	7033	Kind ak = a.getKind();
80	7033	bool consider = true;
81		// if we've already computed master for a
82	7033	if (d_trMaster.find(a) != d_trMaster.end())
83		{
84		// a master has at least one slave
85	6565	consider = (d_trSlaves.find(a) != d_trSlaves.end());
86		}
87		else
88		{
89	468	if (ak == Kind::SINE)
90		{
91		// always not a master
92	312	consider = false;
93		}
94		else
95		{
96	240	for (const Node& ac : a)
97		{
98	84	if (isTranscendentalKind(ac.getKind()))
99		{
100		consider = false;
101		break;
102		}
103		}
104		}
105	468	if (!consider)
106		{
107		// wait to assign a master below
108	312	needsMaster.push_back(a);
109		}
110		else
111		{
112	156	d_trMaster[a] = a;
113	156	d_trSlaves[a].insert(a);
114		}
115		}
116	7033	if (ak == Kind::EXPONENTIAL \|\| ak == Kind::SINE)
117		{
118	6548	needPi = needPi \|\| (ak == Kind::SINE);
119		// if we didn't indicate that it should be purified above
120	13096	if (consider)
121		{
122	4264	ensureCongruence(a, argTrie);
123		}
124		}
125	485	else if (ak == Kind::PI)
126		{
127	485	Assert(consider);
128	485	needPi = true;
129	485	d_funcMap[ak].push_back(a);
130	485	d_funcCongClass[a].push_back(a);
131		}
132		}
133		// initialize pi if necessary
134	3461	if (needPi && d_pi.isNull())
135		{
136	74	mkPi();
137	74	getCurrentPiBounds();
138		}
139
140	3461	if (Trace.isOn("nl-ext-mv"))
141		{
142		Trace("nl-ext-mv") << "Arguments of trancendental functions : "
143		<< std::endl;
144		for (std::pair<const Kind, std::vector<Node> >& tfl : d_funcMap)
145		{
146		Kind k = tfl.first;
147		if (k == Kind::SINE \|\| k == Kind::EXPONENTIAL)
148		{
149		for (const Node& tf : tfl.second)
150		{
151		Node v = tf[0];
152		d_model.computeConcreteModelValue(v);
153		d_model.computeAbstractModelValue(v);
154		d_model.printModelValue("nl-ext-mv", v);
155		}
156		}
157		}
158		}
159	3461	}
160
161	4264	void TranscendentalState::ensureCongruence(TNode a,
162		std::map<Kind, ArgTrie>& argTrie)
163		{
164	4264	NodeManager* nm = NodeManager::currentNM();
165	8528	std::vector<Node> repList;
166	8528	for (const Node& ac : a)
167		{
168	8528	Node r = d_model.computeConcreteModelValue(ac);
169	4264	repList.push_back(r);
170		}
171	8528	Node aa = argTrie[a.getKind()].add(a, repList);
172	4264	if (aa != a)
173		{
174		// apply congruence to pairs of terms that are disequal and congruent
175	352	Assert(aa.getNumChildren() == a.getNumChildren());
176	704	Node mvaa = d_model.computeAbstractModelValue(a);
177	704	Node mvaaa = d_model.computeAbstractModelValue(aa);
178	352	if (mvaa != mvaaa)
179		{
180	148	std::vector<Node> exp;
181	148	for (unsigned j = 0, size = a.getNumChildren(); j < size; j++)
182		{
183	74	exp.push_back(a[j].eqNode(aa[j]));
184		}
185	148	Node expn = exp.size() == 1 ? exp[0] : nm->mkNode(Kind::AND, exp);
186	148	Node cong_lemma = expn.impNode(a.eqNode(aa));
187	74	d_im.addPendingLemma(cong_lemma, InferenceId::ARITH_NL_CONGRUENCE);
188		}
189		}
190		else
191		{
192		// new representative of congruence class
193	3912	d_funcMap[a.getKind()].push_back(a);
194		}
195		// add to congruence class
196	4264	d_funcCongClass[aa].push_back(a);
197	4264	}
198
199	74	void TranscendentalState::mkPi()
200		{
201	74	NodeManager* nm = NodeManager::currentNM();
202	74	if (d_pi.isNull())
203		{
204	74	d_pi = nm->mkNullaryOperator(nm->realType(), Kind::PI);
205	148	d_pi_2 = Rewriter::rewrite(
206	222	nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(1) / Rational(2))));
207	148	d_pi_neg_2 = Rewriter::rewrite(
208	222	nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(-1) / Rational(2))));
209	148	d_pi_neg = Rewriter::rewrite(
210	222	nm->mkNode(Kind::MULT, d_pi, nm->mkConst(Rational(-1))));
211		// initialize bounds
212	74	d_pi_bound[0] = nm->mkConst(Rational(103993) / Rational(33102));
213	74	d_pi_bound[1] = nm->mkConst(Rational(104348) / Rational(33215));
214		}
215	74	}
216
217	74	void TranscendentalState::getCurrentPiBounds()
218		{
219	74	NodeManager* nm = NodeManager::currentNM();
220		Node pi_lem = nm->mkNode(Kind::AND,
221	148	nm->mkNode(Kind::GEQ, d_pi, d_pi_bound[0]),
222	296	nm->mkNode(Kind::LEQ, d_pi, d_pi_bound[1]));
223	74	CDProof* proof = nullptr;
224	74	if (isProofEnabled())
225		{
226	16	proof = getProof();
227	48	proof->addStep(
228	32	pi_lem, PfRule::ARITH_TRANS_PI, {}, {d_pi_bound[0], d_pi_bound[1]});
229		}
230	74	d_im.addPendingLemma(pi_lem, InferenceId::ARITH_NL_T_PI_BOUND, proof);
231	74	}
232
233	123	std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e,
234		TNode center,
235		unsigned d)
236		{
237		// bounds are the minimum and maximum previous secant points
238		// should not repeat secant points: secant lemmas should suffice to
239		// rule out previous assignment
240	123	Assert(std::find(
241		d_secant_points[e][d].begin(), d_secant_points[e][d].end(), center)
242		== d_secant_points[e][d].end());
243		// Insert into the (temporary) vector. We do not update this vector
244		// until we are sure this secant plane lemma has been processed. We do
245		// this by mapping the lemma to a side effect below.
246	246	std::vector<Node> spoints = d_secant_points[e][d];
247	123	spoints.push_back(center);
248
249		// sort
250	123	sortByNlModel(spoints.begin(), spoints.end(), &d_model);
251		// get the resulting index of c
252		unsigned index =
253	123	std::find(spoints.begin(), spoints.end(), center) - spoints.begin();
254
255		// bounds are the next closest upper/lower bound values
256	246	return {index > 0 ? spoints[index - 1] : Node(),
257	246	index < spoints.size() - 1 ? spoints[index + 1] : Node()};
258		}
259
260	246	Node TranscendentalState::mkSecantPlane(
261		TNode arg, TNode lower, TNode upper, TNode lval, TNode uval)
262		{
263	246	NodeManager* nm = NodeManager::currentNM();
264		// Figure 3: S_l( x ), S_u( x ) for s = 0,1
265	492	Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, lower, upper));
266	246	Assert(rcoeff_n.isConst());
267	492	Rational rcoeff = rcoeff_n.getConst<Rational>();
268	246	Assert(rcoeff.sgn() != 0);
269		Node res =
270		nm->mkNode(Kind::PLUS,
271		lval,
272	984	nm->mkNode(Kind::MULT,
273	984	nm->mkNode(Kind::DIVISION,
274	492	nm->mkNode(Kind::MINUS, lval, uval),
275	492	nm->mkNode(Kind::MINUS, lower, upper)),
276	738	nm->mkNode(Kind::MINUS, arg, lower)));
277	492	Trace("nl-trans") << "Creating secant plane for transcendental function of "
278	246	<< arg << std::endl;
279	492	Trace("nl-trans") << "\tfrom ( " << lower << " ; " << lval << " ) to ( "
280	246	<< upper << " ; " << uval << " )" << std::endl;
281	246	Trace("nl-trans") << "\t" << res << std::endl;
282	246	Trace("nl-trans") << "\trewritten: " << Rewriter::rewrite(res) << std::endl;
283	492	return res;
284		}
285
286	246	NlLemma TranscendentalState::mkSecantLemma(TNode lower,
287		TNode upper,
288		TNode lapprox,
289		TNode uapprox,
290		int csign,
291		Convexity convexity,
292		TNode tf,
293		TNode splane,
294		unsigned actual_d)
295		{
296	246	NodeManager* nm = NodeManager::currentNM();
297		// With respect to Figure 3, this is slightly different.
298		// In particular, we chose b to be the model value of bounds[s],
299		// which is a constant although bounds[s] may not be (e.g. if it
300		// contains PI).
301		// To ensure that c...b does not cross an inflection point,
302		// we guard with the symbolic version of bounds[s].
303		// This leads to lemmas e.g. of this form:
304		// ( c <= x <= PI/2 ) => ( sin(x) < ( P( b ) - P( c ) )*( x -
305		// b ) + P( b ) )
306		// where b = (PI/2)^M, the current value of PI/2 in the model.
307		// This is sound since we are guarded by the symbolic
308		// representation of PI/2.
309		Node antec_n = nm->mkNode(Kind::AND,
310	492	nm->mkNode(Kind::GEQ, tf[0], lower),
311	984	nm->mkNode(Kind::LEQ, tf[0], upper));
312	492	Trace("nl-trans") << "Bound for secant plane: " << lower << " <= " << tf[0]
313	246	<< " <= " << upper << std::endl;
314	246	Trace("nl-trans") << "\t" << antec_n << std::endl;
315		// Convex: actual value is below the secant
316		// Concave: actual value is above the secant
317		Node lem = nm->mkNode(
318		Kind::IMPLIES,
319		antec_n,
320	492	nm->mkNode(
321	492	convexity == Convexity::CONVEX ? Kind::LEQ : Kind::GEQ, tf, splane));
322	492	Trace("nl-trans-lemma") << "*** Secant plane lemma (pre-rewrite) : " << lem
323	246	<< std::endl;
324	246	lem = Rewriter::rewrite(lem);
325	246	Trace("nl-trans-lemma") << "*** Secant plane lemma : " << lem << std::endl;
326	246	Assert(d_model.computeAbstractModelValue(lem) == d_false);
327	246	CDProof* proof = nullptr;
328	246	if (isProofEnabled())
329		{
330	42	proof = getProof();
331	42	if (tf.getKind() == Kind::EXPONENTIAL)
332		{
333	18	if (csign == 1)
334		{
335	80	proof->addStep(
336		lem,
337		PfRule::ARITH_TRANS_EXP_APPROX_ABOVE_POS,
338		{},
339	64	{nm->mkConst<Rational>(2 * actual_d), tf[0], lower, upper});
340		}
341		else
342		{
343	10	proof->addStep(
344		lem,
345		PfRule::ARITH_TRANS_EXP_APPROX_ABOVE_NEG,
346		{},
347	8	{nm->mkConst<Rational>(2 * actual_d), tf[0], lower, upper});
348		}
349		}
350	24	else if (tf.getKind() == Kind::SINE)
351		{
352	24	if (convexity == Convexity::CONCAVE)
353		{
354	84	proof->addStep(lem,
355		PfRule::ARITH_TRANS_SINE_APPROX_BELOW_POS,
356		{},
357		{nm->mkConst<Rational>(2 * actual_d),
358		tf[0],
359		lower,
360		upper,
361		lapprox,
362		uapprox
363
364	72	});
365		}
366		else
367		{
368	84	proof->addStep(lem,
369		PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_NEG,
370		{},
371		{nm->mkConst<Rational>(2 * actual_d),
372		tf[0],
373		lower,
374		upper,
375		lapprox,
376	72	uapprox});
377		}
378		}
379		}
380		return NlLemma(
381	492	InferenceId::ARITH_NL_T_SECANT, lem, LemmaProperty::NONE, proof);
382		}
383
384	123	void TranscendentalState::doSecantLemmas(const std::pair<Node, Node>& bounds,
385		TNode poly_approx,
386		TNode center,
387		TNode cval,
388		TNode tf,
389		Convexity convexity,
390		unsigned d,
391		unsigned actual_d)
392		{
393	123	int csign = center.getConst<Rational>().sgn();
394	246	Trace("nl-trans") << "...do secant lemma with center " << center << " val "
395	123	<< cval << " sign " << csign << std::endl;
396	246	Trace("nl-trans") << "...secant bounds are : " << bounds.first << " ... "
397	123	<< bounds.second << std::endl;
398		// take the model value of lower (since may contain PI)
399		// Make secant from bounds.first to center
400	246	Node lower = d_model.computeAbstractModelValue(bounds.first);
401	246	Trace("nl-trans") << "...model value of bound " << bounds.first << " is "
402	123	<< lower << std::endl;
403	123	if (lower != center)
404		{
405		// Figure 3 : P(l), P(u), for s = 0
406		Node lval = Rewriter::rewrite(
407	246	poly_approx.substitute(d_taylor.getTaylorVariable(), lower));
408	246	Node splane = mkSecantPlane(tf[0], lower, center, lval, cval);
409		NlLemma nlem = mkSecantLemma(
410	246	lower, center, lval, cval, csign, convexity, tf, splane, actual_d);
411		// The side effect says that if lem is added, then we should add the
412		// secant point c for (tf,d).
413	123	nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center));
414	123	d_im.addPendingLemma(nlem, true);
415		}
416
417		// take the model value of upper (since may contain PI)
418		// Make secant from center to bounds.second
419	246	Node upper = d_model.computeAbstractModelValue(bounds.second);
420	246	Trace("nl-trans") << "...model value of bound " << bounds.second << " is "
421	123	<< upper << std::endl;
422	123	if (center != upper)
423		{
424		// Figure 3 : P(l), P(u), for s = 1
425		Node uval = Rewriter::rewrite(
426	246	poly_approx.substitute(d_taylor.getTaylorVariable(), upper));
427	246	Node splane = mkSecantPlane(tf[0], center, upper, cval, uval);
428		NlLemma nlem = mkSecantLemma(
429	246	center, upper, cval, uval, csign, convexity, tf, splane, actual_d);
430		// The side effect says that if lem is added, then we should add the
431		// secant point c for (tf,d).
432	123	nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center));
433	123	d_im.addPendingLemma(nlem, true);
434		}
435	123	}
436
437		} // namespace transcendental
438		} // namespace nl
439		} // namespace arith
440		} // namespace theory
441	29577	} // namespace cvc5