GCC Code Coverage Report
Directory: . Exec Total Coverage
File: src/theory/arith/nl/transcendental/taylor_generator.cpp Lines: 118 126 93.7 %
Date: 2021-03-22 Branches: 238 552 43.1 %

Line Exec Source
1
/*********************                                                        */
2
/*! \file taylor_generator.cpp
3
 ** \verbatim
4
 ** Top contributors (to current version):
5
 **   Gereon Kremer, Andrew Reynolds
6
 ** This file is part of the CVC4 project.
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.\endverbatim
11
 **
12
 ** \brief Generate taylor approximations transcendental lemmas.
13
 **/
14
15
#include "theory/arith/nl/transcendental/taylor_generator.h"
16
17
#include "theory/arith/arith_utilities.h"
18
#include "theory/arith/nl/nl_model.h"
19
#include "theory/rewriter.h"
20
21
namespace CVC4 {
22
namespace theory {
23
namespace arith {
24
namespace nl {
25
namespace transcendental {
26
27
4456
TaylorGenerator::TaylorGenerator()
28
4456
    : d_nm(NodeManager::currentNM()),
29
4456
      d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType()))
30
{
31
4456
}
32
33
3351
TNode TaylorGenerator::getTaylorVariable() { return d_taylor_real_fv; }
34
35
3496
std::pair<Node, Node> TaylorGenerator::getTaylor(Kind k, std::uint64_t n)
36
{
37
3496
  Assert(n > 0);
38
  // check if we have already computed this Taylor series
39
3496
  auto itt = d_taylor_terms[k].find(n);
40
3496
  if (itt != d_taylor_terms[k].end())
41
  {
42
3254
    return itt->second;
43
  }
44
45
242
  NodeManager* nm = NodeManager::currentNM();
46
47
  // the current factorial `counter!`
48
484
  Integer factorial = 1;
49
  // the current variable power `x^counter`
50
484
  Node varpow = nm->mkConst(Rational(1));
51
484
  std::vector<Node> sum;
52
1642
  for (std::uint64_t counter = 1; counter <= n; ++counter)
53
  {
54
1400
    if (k == Kind::EXPONENTIAL)
55
    {
56
      // Maclaurin series for exponential:
57
      //   \sum_{n=0}^\infty x^n / n!
58
788
      sum.push_back(
59
1576
          nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
60
    }
61
612
    else if (k == Kind::SINE)
62
    {
63
      // Maclaurin series for exponential:
64
      //   \sum_{n=0}^\infty (-1)^n / (2n+1)! * x^(2n+1)
65
612
      if (counter % 2 == 0)
66
      {
67
306
        int sign = (counter % 4 == 0 ? -1 : 1);
68
918
        sum.push_back(nm->mkNode(Kind::MULT,
69
1224
                                 nm->mkNode(Kind::DIVISION,
70
612
                                            nm->mkConst<Rational>(sign),
71
612
                                            nm->mkConst<Rational>(factorial)),
72
                                 varpow));
73
      }
74
    }
75
1400
    factorial *= counter;
76
1400
    varpow =
77
2800
        Rewriter::rewrite(nm->mkNode(Kind::MULT, d_taylor_real_fv, varpow));
78
  }
79
  Node taylor_sum =
80
484
      Rewriter::rewrite(sum.size() == 1 ? sum[0] : nm->mkNode(Kind::PLUS, sum));
81
  Node taylor_rem = Rewriter::rewrite(
82
484
      nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
83
84
484
  auto res = std::make_pair(taylor_sum, taylor_rem);
85
86
  // put result in cache
87
242
  d_taylor_terms[k][n] = res;
88
89
242
  return res;
90
}
91
92
3240
void TaylorGenerator::getPolynomialApproximationBounds(
93
    Kind k, std::uint64_t d, ApproximationBounds& pbounds)
94
{
95
3240
  auto it = d_poly_bounds[k].find(d);
96
3240
  if (it == d_poly_bounds[k].end())
97
  {
98
242
    NodeManager* nm = NodeManager::currentNM();
99
    // n is the Taylor degree we are currently considering
100
242
    std::uint64_t n = 2 * d;
101
    // n must be even
102
484
    std::pair<Node, Node> taylor = getTaylor(k, n);
103
484
    Node taylor_sum = taylor.first;
104
    // ru is x^{n+1}/(n+1)!
105
484
    Node ru = taylor.second;
106
484
    Trace("nl-trans") << "Taylor for " << k << " is : " << taylor.first
107
242
                      << std::endl;
108
484
    Trace("nl-trans") << "Taylor remainder for " << k << " is " << taylor.second
109
242
                      << std::endl;
110
242
    if (k == Kind::EXPONENTIAL)
111
    {
112
130
      pbounds.d_lower = taylor_sum;
113
130
      pbounds.d_upperNeg =
114
260
          Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
115
260
      pbounds.d_upperPos = Rewriter::rewrite(
116
260
          nm->mkNode(Kind::MULT,
117
                     taylor_sum,
118
390
                     nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru)));
119
    }
120
    else
121
    {
122
112
      Assert(k == Kind::SINE);
123
224
      Node l = Rewriter::rewrite(nm->mkNode(Kind::MINUS, taylor_sum, ru));
124
224
      Node u = Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
125
112
      pbounds.d_lower = l;
126
112
      pbounds.d_upperNeg = u;
127
112
      pbounds.d_upperPos = u;
128
    }
129
484
    Trace("nl-trans") << "Polynomial approximation for " << k
130
242
                      << " is: " << std::endl;
131
242
    Trace("nl-trans") << " Lower: " << pbounds.d_lower << std::endl;
132
242
    Trace("nl-trans") << " Upper (neg): " << pbounds.d_upperNeg << std::endl;
133
242
    Trace("nl-trans") << " Upper (pos): " << pbounds.d_upperPos << std::endl;
134
242
    d_poly_bounds[k].emplace(d, pbounds);
135
  }
136
  else
137
  {
138
2998
    pbounds = it->second;
139
  }
140
3240
}
141
142
2533
std::uint64_t TaylorGenerator::getPolynomialApproximationBoundForArg(
143
    Kind k, Node c, std::uint64_t d, ApproximationBounds& pbounds)
144
{
145
2533
  getPolynomialApproximationBounds(k, d, pbounds);
146
2533
  Trace("nl-trans") << "c = " << c << std::endl;
147
2533
  Assert(c.isConst());
148
2533
  if (k == Kind::EXPONENTIAL && c.getConst<Rational>().sgn() == 1)
149
  {
150
1604
    bool success = false;
151
1604
    std::uint64_t ds = d;
152
3208
    TNode ttrf = getTaylorVariable();
153
3208
    TNode tc = c;
154
1650
    do
155
    {
156
3254
      success = true;
157
3254
      std::uint64_t n = 2 * ds;
158
6508
      std::pair<Node, Node> taylor = getTaylor(k, n);
159
      // check that 1-c^{n+1}/(n+1)! > 0
160
6508
      Node ru = taylor.second;
161
6508
      Node rus = ru.substitute(ttrf, tc);
162
3254
      rus = Rewriter::rewrite(rus);
163
3254
      Assert(rus.isConst());
164
3254
      if (rus.getConst<Rational>() > 1)
165
      {
166
1650
        success = false;
167
1650
        ds = ds + 1;
168
      }
169
3254
    } while (!success);
170
1604
    if (ds > d)
171
    {
172
1280
      Trace("nl-ext-exp-taylor")
173
640
          << "*** Increase Taylor bound to " << ds << " > " << d << " for ("
174
640
          << k << " " << c << ")" << std::endl;
175
      // must use sound upper bound
176
1280
      ApproximationBounds pboundss;
177
640
      getPolynomialApproximationBounds(k, ds, pboundss);
178
640
      pbounds.d_upperPos = pboundss.d_upperPos;
179
    }
180
1604
    return ds;
181
  }
182
929
  return d;
183
}
184
185
1402
std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf,
186
                                                        std::uint64_t d,
187
                                                        NlModel& model)
188
{
189
  // compute the model value of the argument
190
2804
  Node c = model.computeAbstractModelValue(tf[0]);
191
1402
  Assert(c.isConst());
192
1402
  int csign = c.getConst<Rational>().sgn();
193
1402
  Kind k = tf.getKind();
194
1402
  if (csign == 0)
195
  {
196
    NodeManager* nm = NodeManager::currentNM();
197
    // at zero, its trivial
198
    if (k == Kind::SINE)
199
    {
200
      Node zero = nm->mkConst(Rational(0));
201
      return std::pair<Node, Node>(zero, zero);
202
    }
203
    Assert(k == Kind::EXPONENTIAL);
204
    Node one = nm->mkConst(Rational(1));
205
    return std::pair<Node, Node>(one, one);
206
  }
207
1402
  bool isNeg = csign == -1;
208
209
2804
  ApproximationBounds pbounds;
210
1402
  getPolynomialApproximationBoundForArg(k, c, d, pbounds);
211
212
2804
  std::vector<Node> bounds;
213
2804
  TNode tfv = getTaylorVariable();
214
2804
  TNode tfs = tf[0];
215
4206
  for (unsigned d2 = 0; d2 < 2; d2++)
216
  {
217
    Node pab = (d2 == 0 ? pbounds.d_lower
218
5608
                        : (isNeg ? pbounds.d_upperNeg : pbounds.d_upperPos));
219
2804
    if (!pab.isNull())
220
    {
221
      // { x -> M_A(tf[0]) }
222
      // Notice that we compute the model value of tfs first, so that
223
      // the call to rewrite below does not modify the term, where notice that
224
      // rewrite( x*x { x -> M_A(t) } ) = M_A(t)*M_A(t)
225
      // is not equal to
226
      // M_A( x*x { x -> t } ) = M_A( t*t )
227
      // where M_A denotes the abstract model.
228
5608
      Node mtfs = model.computeAbstractModelValue(tfs);
229
2804
      pab = pab.substitute(tfv, mtfs);
230
2804
      pab = Rewriter::rewrite(pab);
231
2804
      Assert(pab.isConst());
232
2804
      bounds.push_back(pab);
233
    }
234
    else
235
    {
236
      bounds.push_back(Node::null());
237
    }
238
  }
239
1402
  return std::pair<Node, Node>(bounds[0], bounds[1]);
240
}
241
242
}  // namespace transcendental
243
}  // namespace nl
244
}  // namespace arith
245
}  // namespace theory
246
26676
}  // namespace CVC4