Head

GCC Code Coverage Report

Directory:	.		Exec	Total	Coverage
File:	src/theory/arith/nl/transcendental/sine_solver.cpp	Lines:	223	225	99.1 %
Date:	2021-05-22	Branches:	547	1146	47.7 %


/******************************************************************************
 * Top contributors (to current version):
 *   Gereon Kremer, Andrew Reynolds, Tim King
 *
 * 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/sine_solver.h"

#include <cmath>
#include <set>

#include "expr/node_algorithm.h"
#include "expr/node_builder.h"
#include "expr/proof.h"
#include "expr/skolem_manager.h"
#include "options/arith_options.h"
#include "theory/arith/arith_msum.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/transcendental_state.h"
#include "theory/rewriter.h"

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

/**
 * Ensure a is in the main phase:
 *   -pi <= a <= pi
 */
inline Node mkValidPhase(TNode a, TNode pi)
{
  return mkBounded(
      NodeManager::currentNM()->mkNode(Kind::MULT, mkRationalNode(-1), pi),
      a,
      pi);
}
}  // namespace

SineSolver::SineSolver(TranscendentalState* tstate) : d_data(tstate) {}

SineSolver::~SineSolver() {}

void SineSolver::doPhaseShift(TNode a, TNode new_a, TNode y)
{
  NodeManager* nm = NodeManager::currentNM();
  SkolemManager* sm = nm->getSkolemManager();
  Assert(a.getKind() == Kind::SINE);
  Trace("nl-ext-tf") << "Basis sine : " << new_a << " for " << a << std::endl;
  Assert(!d_data->d_pi.isNull());
  Node shift = sm->mkDummySkolem("s", nm->integerType(), "number of shifts");
  // TODO (cvc4-projects #47) : do not introduce shift here, instead needs model-based
  // refinement for constant shifts (cvc4-projects #1284)
  Node lem = nm->mkNode(
      Kind::AND,
      mkValidPhase(y, d_data->d_pi),
      nm->mkNode(Kind::ITE,
                 mkValidPhase(a[0], d_data->d_pi),
                 a[0].eqNode(y),
                 a[0].eqNode(nm->mkNode(Kind::PLUS,
                                        y,
                                        nm->mkNode(Kind::MULT,
                                                   nm->mkConst(Rational(2)),
                                                   shift,
                                                   d_data->d_pi)))),
      new_a.eqNode(a));
  CDProof* proof = nullptr;
  if (d_data->isProofEnabled())
  {
    proof = d_data->getProof();
    proof->addStep(lem, PfRule::ARITH_TRANS_SINE_SHIFT, {}, {a[0], y, shift});
  }
  // note we must do preprocess on this lemma
  Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : purify : " << lem
                        << std::endl;
  d_data->d_im.addPendingLemma(lem, InferenceId::ARITH_NL_T_PURIFY_ARG, proof);
}

void SineSolver::checkInitialRefine()
{
  NodeManager* nm = NodeManager::currentNM();
  for (std::pair<const Kind, std::vector<Node> >& tfl : d_data->d_funcMap)
  {
    if (tfl.first != Kind::SINE)
    {
      continue;
    }
    Trace("nl-ext") << "Get initial (sine) refinement lemmas for "
                       "transcendental functions..."
                    << std::endl;
    for (const Node& t : tfl.second)
    {
      // initial refinements
      if (d_tf_initial_refine.find(t) == d_tf_initial_refine.end())
      {
        d_tf_initial_refine[t] = true;
        Node symn = nm->mkNode(Kind::SINE,
                               nm->mkNode(Kind::MULT, d_data->d_neg_one, t[0]));
        symn = Rewriter::rewrite(symn);
        // Can assume it is its own master since phase is split over 0,
        // hence  -pi <= t[0] <= pi implies -pi <= -t[0] <= pi.
        d_data->d_trMaster[symn] = symn;
        d_data->d_trSlaves[symn].insert(symn);
        Assert(d_data->d_trSlaves.find(t) != d_data->d_trSlaves.end());

        {
          // sine bounds: -1 <= sin(t) <= 1
          Node lem = nm->mkNode(Kind::AND,
                                nm->mkNode(Kind::LEQ, t, d_data->d_one),
                                nm->mkNode(Kind::GEQ, t, d_data->d_neg_one));
          CDProof* proof = nullptr;
          if (d_data->isProofEnabled())
          {
            proof = d_data->getProof();
            proof->addStep(lem, PfRule::ARITH_TRANS_SINE_BOUNDS, {}, {t[0]});
          }
          d_data->d_im.addPendingLemma(
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
        }
        {
          // sine symmetry: sin(t) - sin(-t) = 0
          Node lem = nm->mkNode(Kind::PLUS, t, symn).eqNode(d_data->d_zero);
          CDProof* proof = nullptr;
          if (d_data->isProofEnabled())
          {
            proof = d_data->getProof();
            proof->addStep(lem, PfRule::ARITH_TRANS_SINE_SYMMETRY, {}, {t[0]});
          }
          d_data->d_im.addPendingLemma(
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
        }
        {
          // sine zero tangent:
          //   t > 0  =>  sin(t) < t
          //   t < 0  =>  sin(t) > t
          Node lem =
              nm->mkNode(Kind::AND,
                         nm->mkNode(Kind::IMPLIES,
                                    nm->mkNode(Kind::GT, t[0], d_data->d_zero),
                                    nm->mkNode(Kind::LT, t, t[0])),
                         nm->mkNode(Kind::IMPLIES,
                                    nm->mkNode(Kind::LT, t[0], d_data->d_zero),
                                    nm->mkNode(Kind::GT, t, t[0])));
          CDProof* proof = nullptr;
          if (d_data->isProofEnabled())
          {
            proof = d_data->getProof();
            proof->addStep(
                lem, PfRule::ARITH_TRANS_SINE_TANGENT_ZERO, {}, {t[0]});
          }
          d_data->d_im.addPendingLemma(
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
        }
        {
          // sine pi tangent:
          //   t > -pi  =>  sin(t) > -pi-t
          //   t <  pi  =>  sin(t) <  pi-t
          Node lem = nm->mkNode(
              Kind::AND,
              nm->mkNode(
                  Kind::IMPLIES,
                  nm->mkNode(Kind::GT, t[0], d_data->d_pi_neg),
                  nm->mkNode(Kind::GT,
                             t,
                             nm->mkNode(Kind::MINUS, d_data->d_pi_neg, t[0]))),
              nm->mkNode(
                  Kind::IMPLIES,
                  nm->mkNode(Kind::LT, t[0], d_data->d_pi),
                  nm->mkNode(Kind::LT,
                             t,
                             nm->mkNode(Kind::MINUS, d_data->d_pi, t[0]))));
          CDProof* proof = nullptr;
          if (d_data->isProofEnabled())
          {
            proof = d_data->getProof();
            proof->addStep(
                lem, PfRule::ARITH_TRANS_SINE_TANGENT_PI, {}, {t[0]});
          }
          d_data->d_im.addPendingLemma(
              lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
        }
        {
          Node lem =
              nm->mkNode(Kind::AND,
                         // sign
                         nm->mkNode(Kind::EQUAL,
                                    nm->mkNode(Kind::LT, t[0], d_data->d_zero),
                                    nm->mkNode(Kind::LT, t, d_data->d_zero)),
                         // zero val
                         nm->mkNode(Kind::EQUAL,
                                    nm->mkNode(Kind::GT, t[0], d_data->d_zero),
                                    nm->mkNode(Kind::GT, t, d_data->d_zero)));
          d_data->d_im.addPendingLemma(lem,
                                       InferenceId::ARITH_NL_T_INIT_REFINE);
        }
      }
    }
  }
}

void SineSolver::checkMonotonic()
{

  auto it = d_data->d_funcMap.find(Kind::SINE);
  if (it == d_data->d_funcMap.end())
  {
    Trace("nl-ext-exp") << "No sine terms" << std::endl;
    return;
  }
  Trace("nl-ext")
      << "Get monotonicity lemmas for (sine) transcendental functions..."
      << std::endl;

  // sort arguments of all transcendentals
  std::vector<Node> tf_args;
  std::map<Node, Node> tf_arg_to_term;

  for (const Node& tf : it->second)
  {
    Node a = tf[0];
    Node mvaa = d_data->d_model.computeAbstractModelValue(a);
    if (mvaa.isConst())
    {
      Trace("nl-ext-tf-mono-debug") << "...tf term : " << a << std::endl;
      tf_args.push_back(a);
      tf_arg_to_term[a] = tf;
    }
  }

  if (tf_args.empty())
  {
    return;
  }

  sortByNlModel(
      tf_args.begin(), tf_args.end(), &d_data->d_model, true, false, true);

  std::vector<Node> mpoints = {d_data->d_pi,
                               d_data->d_pi_2,
                               d_data->d_zero,
                               d_data->d_pi_neg_2,
                               d_data->d_pi_neg};
  std::vector<Node> mpoints_vals;

  // get model values for points
  for (const auto& point : mpoints)
  {
    mpoints_vals.emplace_back(d_data->d_model.computeAbstractModelValue(point));
    Assert(mpoints_vals.back().isConst());
  }

  unsigned mdir_index = 0;
  int monotonic_dir = -1;
  Node mono_bounds[2];
  Node targ, targval, t, tval;
  for (const auto& sarg : tf_args)
  {
    Node sargval = d_data->d_model.computeAbstractModelValue(sarg);
    Assert(sargval.isConst());
    Node s = tf_arg_to_term[sarg];
    Node sval = d_data->d_model.computeAbstractModelValue(s);
    Assert(sval.isConst());

    // increment to the proper monotonicity region
    bool increment = true;
    while (increment && mdir_index < mpoints.size())
    {
      increment = false;
      Node pval = mpoints_vals[mdir_index];
      Assert(pval.isConst());
      if (sargval.getConst<Rational>() < pval.getConst<Rational>())
      {
        increment = true;
        Trace("nl-ext-tf-mono")
            << "...increment at " << sarg << " since model value is less than "
            << mpoints[mdir_index] << std::endl;
      }
      if (increment)
      {
        tval = Node::null();
        mono_bounds[1] = mpoints[mdir_index];
        mdir_index++;
        monotonic_dir = regionToMonotonicityDir(mdir_index);
        if (mdir_index < mpoints.size())
        {
          mono_bounds[0] = mpoints[mdir_index];
        }
        else
        {
          mono_bounds[0] = Node::null();
        }
      }
    }
    // store the concavity region
    d_data->d_tf_region[s] = mdir_index;
    Trace("nl-ext-concavity")
        << "Transcendental function " << s << " is in region #" << mdir_index;
    Trace("nl-ext-concavity") << ", arg model value = " << sargval << std::endl;

    if (!tval.isNull())
    {
      NodeManager* nm = NodeManager::currentNM();
      Node mono_lem;
      if (monotonic_dir == 1
          && sval.getConst<Rational>() > tval.getConst<Rational>())
      {
        mono_lem = nm->mkNode(Kind::IMPLIES,
                              nm->mkNode(Kind::GEQ, targ, sarg),
                              nm->mkNode(Kind::GEQ, t, s));
      }
      else if (monotonic_dir == -1
               && sval.getConst<Rational>() < tval.getConst<Rational>())
      {
        mono_lem = nm->mkNode(Kind::IMPLIES,
                              nm->mkNode(Kind::LEQ, targ, sarg),
                              nm->mkNode(Kind::LEQ, t, s));
      }
      if (!mono_lem.isNull())
      {
        if (!mono_bounds[0].isNull())
        {
          Assert(!mono_bounds[1].isNull());
          mono_lem = nm->mkNode(
              Kind::IMPLIES,
              nm->mkNode(Kind::AND,
                         mkBounded(mono_bounds[0], targ, mono_bounds[1]),
                         mkBounded(mono_bounds[0], sarg, mono_bounds[1])),
              mono_lem);
        }
        Trace("nl-ext-tf-mono")
            << "Monotonicity lemma : " << mono_lem << std::endl;

        d_data->d_im.addPendingLemma(mono_lem,
                                     InferenceId::ARITH_NL_T_MONOTONICITY);
      }
    }
    // store the previous values
    targ = sarg;
    targval = sargval;
    t = s;
    tval = sval;
  }
}

void SineSolver::doTangentLemma(
    TNode e, TNode c, TNode poly_approx, int region, std::uint64_t d)
{
  NodeManager* nm = NodeManager::currentNM();

  // compute tangent plane
  // Figure 3: T( x )
  // We use zero slope tangent planes, since the concavity of the Taylor
  // approximation cannot be easily established.
  Convexity convexity = regionToConvexity(region);
  int mdir = regionToMonotonicityDir(region);
  bool usec = (mdir == 1) == (convexity == Convexity::CONCAVE);
  Node lem = nm->mkNode(
      Kind::IMPLIES,
      nm->mkNode(
          Kind::AND,
          nm->mkNode(
              Kind::GEQ, e[0], usec ? regionToLowerBound(region) : Node(c)),
          nm->mkNode(
              Kind::LEQ, e[0], usec ? Node(c) : regionToUpperBound(region))),
      nm->mkNode(convexity == Convexity::CONVEX ? Kind::GEQ : Kind::LEQ,
                 e,
                 poly_approx));

  Trace("nl-ext-sine") << "*** Tangent plane lemma (pre-rewrite): " << lem
                       << std::endl;
  lem = Rewriter::rewrite(lem);
  Trace("nl-ext-sine") << "*** Tangent plane lemma : " << lem << std::endl;
  Assert(d_data->d_model.computeAbstractModelValue(lem) == d_data->d_false);
  // Figure 3 : line 9
  CDProof* proof = nullptr;
  if (d_data->isProofEnabled())
  {
    proof = d_data->getProof();
    if (convexity == Convexity::CONVEX)
    {
      if (usec)
      {
        proof->addStep(lem,
                       PfRule::ARITH_TRANS_SINE_APPROX_BELOW_NEG,
                       {},
                       {nm->mkConst<Rational>(2 * d),
                        e[0],
                        c,
                        regionToLowerBound(region),
                        c});
      }
      else
      {
        proof->addStep(lem,
                       PfRule::ARITH_TRANS_SINE_APPROX_BELOW_NEG,
                       {},
                       {nm->mkConst<Rational>(2 * d),
                        e[0],
                        c,
                        c,
                        regionToUpperBound(region)});
      }
    }
    else
    {
      if (usec)
      {
        proof->addStep(lem,
                       PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_POS,
                       {},
                       {nm->mkConst<Rational>(2 * d),
                        e[0],
                        c,
                        regionToLowerBound(region),
                        c});
      }
      else
      {
        proof->addStep(lem,
                       PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_POS,
                       {},
                       {nm->mkConst<Rational>(2 * d),
                        e[0],
                        c,
                        c,
                        regionToUpperBound(region)});
      }
    }
  }
  d_data->d_im.addPendingLemma(
      lem, InferenceId::ARITH_NL_T_TANGENT, proof, true);
}

void SineSolver::doSecantLemmas(TNode e,
                                TNode poly_approx,
                                TNode c,
                                TNode poly_approx_c,
                                unsigned d,
                                unsigned actual_d,
                                int region)
{
  d_data->doSecantLemmas(getSecantBounds(e, c, d, region),
                         poly_approx,
                         c,
                         poly_approx_c,
                         e,
                         regionToConvexity(region),
                         d,
                         actual_d);
}

std::pair<Node, Node> SineSolver::getSecantBounds(TNode e,
                                                  TNode c,
                                                  unsigned d,
                                                  int region)
{
  std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, c, d);

  // Check if we already have neighboring secant points
  if (bounds.first.isNull())
  {
    // lower boundary point for this concavity region
    bounds.first = regionToLowerBound(region);
  }
  if (bounds.second.isNull())
  {
    // upper boundary point for this concavity region
    bounds.second = regionToUpperBound(region);
  }
  return bounds;
}

}  // 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, Tim King
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/sine_solver.h"
17
18		#include <cmath>
19		#include <set>
20
21		#include "expr/node_algorithm.h"
22		#include "expr/node_builder.h"
23		#include "expr/proof.h"
24		#include "expr/skolem_manager.h"
25		#include "options/arith_options.h"
26		#include "theory/arith/arith_msum.h"
27		#include "theory/arith/arith_utilities.h"
28		#include "theory/arith/inference_manager.h"
29		#include "theory/arith/nl/nl_model.h"
30		#include "theory/arith/nl/transcendental/transcendental_state.h"
31		#include "theory/rewriter.h"
32
33		namespace cvc5 {
34		namespace theory {
35		namespace arith {
36		namespace nl {
37		namespace transcendental {
38		namespace {
39
40		/**
41		* Ensure a is in the main phase:
42		* -pi <= a <= pi
43		*/
44	284	inline Node mkValidPhase(TNode a, TNode pi)
45		{
46		return mkBounded(
47	568	NodeManager::currentNM()->mkNode(Kind::MULT, mkRationalNode(-1), pi),
48		a,
49	852	pi);
50		}
51		} // namespace
52
53	4914	SineSolver::SineSolver(TranscendentalState* tstate) : d_data(tstate) {}
54
55	4914	SineSolver::~SineSolver() {}
56
57	142	void SineSolver::doPhaseShift(TNode a, TNode new_a, TNode y)
58		{
59	142	NodeManager* nm = NodeManager::currentNM();
60	142	SkolemManager* sm = nm->getSkolemManager();
61	142	Assert(a.getKind() == Kind::SINE);
62	142	Trace("nl-ext-tf") << "Basis sine : " << new_a << " for " << a << std::endl;
63	142	Assert(!d_data->d_pi.isNull());
64	284	Node shift = sm->mkDummySkolem("s", nm->integerType(), "number of shifts");
65		// TODO (cvc4-projects #47) : do not introduce shift here, instead needs model-based
66		// refinement for constant shifts (cvc4-projects #1284)
67		Node lem = nm->mkNode(
68		Kind::AND,
69	284	mkValidPhase(y, d_data->d_pi),
70	852	nm->mkNode(Kind::ITE,
71	284	mkValidPhase(a[0], d_data->d_pi),
72	284	a[0].eqNode(y),
73	568	a[0].eqNode(nm->mkNode(Kind::PLUS,
74		y,
75	568	nm->mkNode(Kind::MULT,
76	284	nm->mkConst(Rational(2)),
77		shift,
78	142	d_data->d_pi)))),
79	852	new_a.eqNode(a));
80	142	CDProof* proof = nullptr;
81	142	if (d_data->isProofEnabled())
82		{
83	24	proof = d_data->getProof();
84	24	proof->addStep(lem, PfRule::ARITH_TRANS_SINE_SHIFT, {}, {a[0], y, shift});
85		}
86		// note we must do preprocess on this lemma
87	284	Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : purify : " << lem
88	142	<< std::endl;
89	142	d_data->d_im.addPendingLemma(lem, InferenceId::ARITH_NL_T_PURIFY_ARG, proof);
90	142	}
91
92	2569	void SineSolver::checkInitialRefine()
93		{
94	2569	NodeManager* nm = NodeManager::currentNM();
95	3463	for (std::pair<const Kind, std::vector<Node> >& tfl : d_data->d_funcMap)
96		{
97	894	if (tfl.first != Kind::SINE)
98		{
99	524	continue;
100		}
101	740	Trace("nl-ext") << "Get initial (sine) refinement lemmas for "
102	370	"transcendental functions..."
103	370	<< std::endl;
104	3312	for (const Node& t : tfl.second)
105		{
106		// initial refinements
107	2942	if (d_tf_initial_refine.find(t) == d_tf_initial_refine.end())
108		{
109	252	d_tf_initial_refine[t] = true;
110		Node symn = nm->mkNode(Kind::SINE,
111	504	nm->mkNode(Kind::MULT, d_data->d_neg_one, t[0]));
112	252	symn = Rewriter::rewrite(symn);
113		// Can assume it is its own master since phase is split over 0,
114		// hence -pi <= t[0] <= pi implies -pi <= -t[0] <= pi.
115	252	d_data->d_trMaster[symn] = symn;
116	252	d_data->d_trSlaves[symn].insert(symn);
117	252	Assert(d_data->d_trSlaves.find(t) != d_data->d_trSlaves.end());
118
119		{
120		// sine bounds: -1 <= sin(t) <= 1
121		Node lem = nm->mkNode(Kind::AND,
122	504	nm->mkNode(Kind::LEQ, t, d_data->d_one),
123	1008	nm->mkNode(Kind::GEQ, t, d_data->d_neg_one));
124	252	CDProof* proof = nullptr;
125	252	if (d_data->isProofEnabled())
126		{
127	40	proof = d_data->getProof();
128	40	proof->addStep(lem, PfRule::ARITH_TRANS_SINE_BOUNDS, {}, {t[0]});
129		}
130	252	d_data->d_im.addPendingLemma(
131		lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
132		}
133		{
134		// sine symmetry: sin(t) - sin(-t) = 0
135	504	Node lem = nm->mkNode(Kind::PLUS, t, symn).eqNode(d_data->d_zero);
136	252	CDProof* proof = nullptr;
137	252	if (d_data->isProofEnabled())
138		{
139	40	proof = d_data->getProof();
140	40	proof->addStep(lem, PfRule::ARITH_TRANS_SINE_SYMMETRY, {}, {t[0]});
141		}
142	252	d_data->d_im.addPendingLemma(
143		lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
144		}
145		{
146		// sine zero tangent:
147		// t > 0 => sin(t) < t
148		// t < 0 => sin(t) > t
149		Node lem =
150		nm->mkNode(Kind::AND,
151	1008	nm->mkNode(Kind::IMPLIES,
152	504	nm->mkNode(Kind::GT, t[0], d_data->d_zero),
153	504	nm->mkNode(Kind::LT, t, t[0])),
154	1008	nm->mkNode(Kind::IMPLIES,
155	504	nm->mkNode(Kind::LT, t[0], d_data->d_zero),
156	1512	nm->mkNode(Kind::GT, t, t[0])));
157	252	CDProof* proof = nullptr;
158	252	if (d_data->isProofEnabled())
159		{
160	40	proof = d_data->getProof();
161	80	proof->addStep(
162	40	lem, PfRule::ARITH_TRANS_SINE_TANGENT_ZERO, {}, {t[0]});
163		}
164	252	d_data->d_im.addPendingLemma(
165		lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
166		}
167		{
168		// sine pi tangent:
169		// t > -pi => sin(t) > -pi-t
170		// t < pi => sin(t) < pi-t
171		Node lem = nm->mkNode(
172		Kind::AND,
173	1008	nm->mkNode(
174		Kind::IMPLIES,
175	504	nm->mkNode(Kind::GT, t[0], d_data->d_pi_neg),
176	504	nm->mkNode(Kind::GT,
177		t,
178	504	nm->mkNode(Kind::MINUS, d_data->d_pi_neg, t[0]))),
179	1008	nm->mkNode(
180		Kind::IMPLIES,
181	504	nm->mkNode(Kind::LT, t[0], d_data->d_pi),
182	504	nm->mkNode(Kind::LT,
183		t,
184	1512	nm->mkNode(Kind::MINUS, d_data->d_pi, t[0]))));
185	252	CDProof* proof = nullptr;
186	252	if (d_data->isProofEnabled())
187		{
188	40	proof = d_data->getProof();
189	80	proof->addStep(
190	40	lem, PfRule::ARITH_TRANS_SINE_TANGENT_PI, {}, {t[0]});
191		}
192	252	d_data->d_im.addPendingLemma(
193		lem, InferenceId::ARITH_NL_T_INIT_REFINE, proof);
194		}
195		{
196		Node lem =
197		nm->mkNode(Kind::AND,
198		// sign
199	1008	nm->mkNode(Kind::EQUAL,
200	504	nm->mkNode(Kind::LT, t[0], d_data->d_zero),
201	504	nm->mkNode(Kind::LT, t, d_data->d_zero)),
202		// zero val
203	1008	nm->mkNode(Kind::EQUAL,
204	504	nm->mkNode(Kind::GT, t[0], d_data->d_zero),
205	1512	nm->mkNode(Kind::GT, t, d_data->d_zero)));
206	252	d_data->d_im.addPendingLemma(lem,
207		InferenceId::ARITH_NL_T_INIT_REFINE);
208		}
209		}
210		}
211		}
212	2569	}
213
214	1106	void SineSolver::checkMonotonic()
215		{
216
217	1106	auto it = d_data->d_funcMap.find(Kind::SINE);
218	1106	if (it == d_data->d_funcMap.end())
219		{
220	966	Trace("nl-ext-exp") << "No sine terms" << std::endl;
221	1932	return;
222		}
223	280	Trace("nl-ext")
224	140	<< "Get monotonicity lemmas for (sine) transcendental functions..."
225	140	<< std::endl;
226
227		// sort arguments of all transcendentals
228	280	std::vector<Node> tf_args;
229	280	std::map<Node, Node> tf_arg_to_term;
230
231	1400	for (const Node& tf : it->second)
232		{
233	2520	Node a = tf[0];
234	2520	Node mvaa = d_data->d_model.computeAbstractModelValue(a);
235	1260	if (mvaa.isConst())
236		{
237	1260	Trace("nl-ext-tf-mono-debug") << "...tf term : " << a << std::endl;
238	1260	tf_args.push_back(a);
239	1260	tf_arg_to_term[a] = tf;
240		}
241		}
242
243	140	if (tf_args.empty())
244		{
245		return;
246		}
247
248	140	sortByNlModel(
249	140	tf_args.begin(), tf_args.end(), &d_data->d_model, true, false, true);
250
251	140	std::vector<Node> mpoints = {d_data->d_pi,
252	140	d_data->d_pi_2,
253	140	d_data->d_zero,
254	140	d_data->d_pi_neg_2,
255	700	d_data->d_pi_neg};
256	280	std::vector<Node> mpoints_vals;
257
258		// get model values for points
259	840	for (const auto& point : mpoints)
260		{
261	700	mpoints_vals.emplace_back(d_data->d_model.computeAbstractModelValue(point));
262	700	Assert(mpoints_vals.back().isConst());
263		}
264
265	140	unsigned mdir_index = 0;
266	140	int monotonic_dir = -1;
267	280	Node mono_bounds[2];
268	280	Node targ, targval, t, tval;
269	1400	for (const auto& sarg : tf_args)
270		{
271	2520	Node sargval = d_data->d_model.computeAbstractModelValue(sarg);
272	1260	Assert(sargval.isConst());
273	2520	Node s = tf_arg_to_term[sarg];
274	2520	Node sval = d_data->d_model.computeAbstractModelValue(s);
275	1260	Assert(sval.isConst());
276
277		// increment to the proper monotonicity region
278	1260	bool increment = true;
279	4740	while (increment && mdir_index < mpoints.size())
280		{
281	1740	increment = false;
282	3480	Node pval = mpoints_vals[mdir_index];
283	1740	Assert(pval.isConst());
284	1740	if (sargval.getConst<Rational>() < pval.getConst<Rational>())
285		{
286	480	increment = true;
287	960	Trace("nl-ext-tf-mono")
288	480	<< "...increment at " << sarg << " since model value is less than "
289	480	<< mpoints[mdir_index] << std::endl;
290		}
291	1740	if (increment)
292		{
293	480	tval = Node::null();
294	480	mono_bounds[1] = mpoints[mdir_index];
295	480	mdir_index++;
296	480	monotonic_dir = regionToMonotonicityDir(mdir_index);
297	480	if (mdir_index < mpoints.size())
298		{
299	480	mono_bounds[0] = mpoints[mdir_index];
300		}
301		else
302		{
303		mono_bounds[0] = Node::null();
304		}
305		}
306		}
307		// store the concavity region
308	1260	d_data->d_tf_region[s] = mdir_index;
309	2520	Trace("nl-ext-concavity")
310	1260	<< "Transcendental function " << s << " is in region #" << mdir_index;
311	1260	Trace("nl-ext-concavity") << ", arg model value = " << sargval << std::endl;
312
313	1260	if (!tval.isNull())
314		{
315	884	NodeManager* nm = NodeManager::currentNM();
316	1768	Node mono_lem;
317	884	if (monotonic_dir == 1
318	884	&& sval.getConst<Rational>() > tval.getConst<Rational>())
319		{
320	180	mono_lem = nm->mkNode(Kind::IMPLIES,
321	120	nm->mkNode(Kind::GEQ, targ, sarg),
322	120	nm->mkNode(Kind::GEQ, t, s));
323		}
324	824	else if (monotonic_dir == -1
325	824	&& sval.getConst<Rational>() < tval.getConst<Rational>())
326		{
327	84	mono_lem = nm->mkNode(Kind::IMPLIES,
328	56	nm->mkNode(Kind::LEQ, targ, sarg),
329	56	nm->mkNode(Kind::LEQ, t, s));
330		}
331	884	if (!mono_lem.isNull())
332		{
333	88	if (!mono_bounds[0].isNull())
334		{
335	88	Assert(!mono_bounds[1].isNull());
336	264	mono_lem = nm->mkNode(
337		Kind::IMPLIES,
338	352	nm->mkNode(Kind::AND,
339	176	mkBounded(mono_bounds[0], targ, mono_bounds[1]),
340	176	mkBounded(mono_bounds[0], sarg, mono_bounds[1])),
341		mono_lem);
342		}
343	176	Trace("nl-ext-tf-mono")
344	88	<< "Monotonicity lemma : " << mono_lem << std::endl;
345
346	88	d_data->d_im.addPendingLemma(mono_lem,
347		InferenceId::ARITH_NL_T_MONOTONICITY);
348		}
349		}
350		// store the previous values
351	1260	targ = sarg;
352	1260	targval = sargval;
353	1260	t = s;
354	1260	tval = sval;
355		}
356		}
357
358	32	void SineSolver::doTangentLemma(
359		TNode e, TNode c, TNode poly_approx, int region, std::uint64_t d)
360		{
361	32	NodeManager* nm = NodeManager::currentNM();
362
363		// compute tangent plane
364		// Figure 3: T( x )
365		// We use zero slope tangent planes, since the concavity of the Taylor
366		// approximation cannot be easily established.
367	32	Convexity convexity = regionToConvexity(region);
368	32	int mdir = regionToMonotonicityDir(region);
369	32	bool usec = (mdir == 1) == (convexity == Convexity::CONCAVE);
370		Node lem = nm->mkNode(
371		Kind::IMPLIES,
372	128	nm->mkNode(
373		Kind::AND,
374	64	nm->mkNode(
375	64	Kind::GEQ, e[0], usec ? regionToLowerBound(region) : Node(c)),
376	64	nm->mkNode(
377	64	Kind::LEQ, e[0], usec ? Node(c) : regionToUpperBound(region))),
378	64	nm->mkNode(convexity == Convexity::CONVEX ? Kind::GEQ : Kind::LEQ,
379		e,
380	128	poly_approx));
381
382	64	Trace("nl-ext-sine") << "*** Tangent plane lemma (pre-rewrite): " << lem
383	32	<< std::endl;
384	32	lem = Rewriter::rewrite(lem);
385	32	Trace("nl-ext-sine") << "*** Tangent plane lemma : " << lem << std::endl;
386	32	Assert(d_data->d_model.computeAbstractModelValue(lem) == d_data->d_false);
387		// Figure 3 : line 9
388	32	CDProof* proof = nullptr;
389	32	if (d_data->isProofEnabled())
390		{
391	8	proof = d_data->getProof();
392	8	if (convexity == Convexity::CONVEX)
393		{
394	4	if (usec)
395		{
396	6	proof->addStep(lem,
397		PfRule::ARITH_TRANS_SINE_APPROX_BELOW_NEG,
398		{},
399		{nm->mkConst<Rational>(2 * d),
400		e[0],
401		c,
402		regionToLowerBound(region),
403	5	c});
404		}
405		else
406		{
407	18	proof->addStep(lem,
408		PfRule::ARITH_TRANS_SINE_APPROX_BELOW_NEG,
409		{},
410		{nm->mkConst<Rational>(2 * d),
411		e[0],
412		c,
413		c,
414	15	regionToUpperBound(region)});
415		}
416		}
417		else
418		{
419	4	if (usec)
420		{
421	18	proof->addStep(lem,
422		PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_POS,
423		{},
424		{nm->mkConst<Rational>(2 * d),
425		e[0],
426		c,
427		regionToLowerBound(region),
428	15	c});
429		}
430		else
431		{
432	6	proof->addStep(lem,
433		PfRule::ARITH_TRANS_SINE_APPROX_ABOVE_POS,
434		{},
435		{nm->mkConst<Rational>(2 * d),
436		e[0],
437		c,
438		c,
439	5	regionToUpperBound(region)});
440		}
441		}
442		}
443	32	d_data->d_im.addPendingLemma(
444		lem, InferenceId::ARITH_NL_T_TANGENT, proof, true);
445	32	}
446
447	72	void SineSolver::doSecantLemmas(TNode e,
448		TNode poly_approx,
449		TNode c,
450		TNode poly_approx_c,
451		unsigned d,
452		unsigned actual_d,
453		int region)
454		{
455	72	d_data->doSecantLemmas(getSecantBounds(e, c, d, region),
456		poly_approx,
457		c,
458		poly_approx_c,
459		e,
460		regionToConvexity(region),
461		d,
462		actual_d);
463	72	}
464
465	72	std::pair<Node, Node> SineSolver::getSecantBounds(TNode e,
466		TNode c,
467		unsigned d,
468		int region)
469		{
470	72	std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, c, d);
471
472		// Check if we already have neighboring secant points
473	72	if (bounds.first.isNull())
474		{
475		// lower boundary point for this concavity region
476	72	bounds.first = regionToLowerBound(region);
477		}
478	72	if (bounds.second.isNull())
479		{
480		// upper boundary point for this concavity region
481	72	bounds.second = regionToUpperBound(region);
482		}
483	72	return bounds;
484		}
485
486		} // namespace transcendental
487		} // namespace nl
488		} // namespace arith
489		} // namespace theory
490	28191	} // namespace cvc5