Head

GCC Code Coverage Report

Directory:	.		Exec	Total	Coverage
File:	src/theory/arith/nl/transcendental/sine_solver.cpp	Lines:	225	227	99.1 %
Date:	2021-11-07	Branches:	547	1144	47.8 %


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