carl  24.04
Computer ARithmetic Library
Substitution.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "Degree.h"
4 #include "Evaluation.h"
5 #include "Power.h"
6 
7 #include "../Monomial.h"
8 
9 namespace carl {
10 
11 /**
12  * Applies the given substitutions to a monomial.
13  * Every variable may be substituted by some value.
14  * @param m The monomial.
15  * @param substitutions Maps variables to numbers.
16  * @return \f$ this[<substitutions>] \f$
17  */
18 template<typename Coeff>
19 Coeff substitute(const Monomial& m, const std::map<Variable, Coeff>& substitutions) {
20  CARL_LOG_FUNC("carl.core.monomial", m << ", " << substitutions);
22  for (const auto& ve : m) {
23  auto it = substitutions.find(ve.first);
24  if (it == substitutions.end()) {
25  res *= carl::pow(ve.first, ve.second);
26  } else {
27  res *= carl::pow(it->second, ve.second);
28  }
29  }
30  CARL_LOG_TRACE("carl.core.monomial", "Result: " << res);
31  return res;
32 }
33 
34 template<typename Coeff>
35 Term<Coeff> substitute(const Term<Coeff>& t, const std::map<Variable, Coeff>& substitutions) {
36  if (t.monomial()) {
38  Coeff coeff = t.coeff();
39  for (const auto& c: *t.monomial()) {
40  auto it = substitutions.find(c.first);
41  if (it == substitutions.end()) {
42  content.push_back(c);
43  } else {
44  coeff *= carl::pow(it->second, c.second);
45  }
46  }
47  if (content.empty()) return Term<Coeff>(coeff);
48  return Term<Coeff>(coeff, createMonomial(std::move(content)));
49  } else {
50  return Term<Coeff>(t.coeff());
51  }
52 }
53 
54 template<typename Coeff>
55 Term<Coeff> substitute(const Term<Coeff>& t, const std::map<Variable, Term<Coeff>>& substitutions) {
56  if (t.monomial()) {
57  return Coeff(t.coeff()) * substitute(*t.monomial(), substitutions);
58  } else {
59  return Term<Coeff>(t.coeff());
60  }
61 }
62 
63 template<typename C, typename O, typename P>
65  assert(p.is_consistent());
66  assert(value.is_consistent());
67  if (!p.has(var)) {
68  return;
69  }
71  // If we replace a variable by zero, just eliminate all terms containing the variable.
72  if(carl::is_zero(value))
73  {
74  bool removedLast = false;
75  for (const auto& term: p) {
76  if (!term.has(var)) {
77  newTerms.push_back(term);
78  removedLast = false;
79  } else removedLast = true;
80  }
81  p.terms().swap(newTerms);
82  CARL_LOG_TRACE("carl.core", p << " [ " << var << " -> " << value << " ] = " << p);
83  if (removedLast) {
84  p.reset_ordered();
85  p.template makeMinimallyOrdered<false, true>();
86  }
87  assert(p.is_consistent());
88  return;
89  }
90  // Find all exponents occurring with the variable to substitute as basis.
91  // expResults will finally be a mapping from every exponent e for which var^e occurs to the value^e and the number of times var^e occurs.
92  // Meanwhile, we store an upper bound on the expected number of terms of the result in expectedResultSize.
93  std::map<exponent, std::pair<MultivariatePolynomial<C,O,P>, size_t>> expResults;
94  size_t expectedResultSize = 0;
95  std::pair<MultivariatePolynomial<C,O,P>, unsigned> def( MultivariatePolynomial<C,O,P>(constant_one<C>::get()), 1 );
96  for(const auto& term: p)
97  {
98  if(term.monomial())
99  { // This is not the constant part.
100  exponent e = term.monomial()->exponent_of_variable(var);
101  if(e > 1)
102  { // Variable occurs with exponent at least two. Insert into map and increase counter in map.
103  auto iterBoolPair = expResults.insert(std::pair<exponent, std::pair<MultivariatePolynomial<C,O,P>, size_t>>(e, def));
104  if(!iterBoolPair.second)
105  {
106  ++(iterBoolPair.first->second.second);
107  }
108  }
109  else if(e == 1)
110  { // Variable occurs with exponent one.
111  expectedResultSize += value.nr_terms();
112  }
113  else
114  { // Variable does not occur in this term.
115  ++expectedResultSize;
116  }
117  }
118  else
119  { // This is the constant part.
120  ++expectedResultSize;
121  }
122  }
123  // Calculate the exponentiation of the multivariate polynomial to substitute the
124  // variable for, reusing the already calculated exponentiations.
125  if( !expResults.empty() )
126  {
127  // Last var^e
128  auto expResultA = expResults.begin();
129  // Next var^e
130  auto expResultB = expResultA;
131  // Calculate first one
132  expResultB->second.first = carl::pow(value, expResultB->first);
133  expectedResultSize += expResultB->second.second * expResultB->second.first.nr_terms();
134  ++expResultB;
135  while(expResultB != expResults.end())
136  {
137  // Calculate next var^e based on the last one.
138  expResultB->second.first = expResultA->second.first * carl::pow(value, expResultB->first - expResultA->first);
139  expectedResultSize += expResultB->second.second * expResultB->second.first.nr_terms();
140  ++expResultA;
141  ++expResultB;
142  }
143  }
144  // Substitute the variable.
146  auto id = tam.getId(expectedResultSize);
147  for (const auto& term: p)
148  {
149  if (term.monomial() == nullptr) {
150  tam.template addTerm<false>(id, term);
151  } else {
152  exponent e = term.monomial()->exponent_of_variable(var);
153  Monomial::Arg mon;
154  if (e > 0) mon = term.monomial()->drop_variable(var);
155  if (e == 1) {
156  for(auto vterm : value)
157  {
158  if (mon == nullptr) tam.template addTerm<false>(id, Term<C>(vterm.coeff() * term.coeff(), vterm.monomial()));
159  else if (vterm.monomial() == nullptr) tam.template addTerm<false>(id, Term<C>(vterm.coeff() * term.coeff(), mon));
160  else tam.template addTerm<false>(id, Term<C>(vterm.coeff() * term.coeff(), vterm.monomial() * mon));
161  }
162  } else if(e > 1) {
163  auto iter = expResults.find(e);
164  assert(iter != expResults.end());
165  for(auto vterm : iter->second.first)
166  {
167  if (mon == nullptr) tam.template addTerm<false>(id, Term<C>(vterm.coeff() * term.coeff(), vterm.monomial()));
168  else if (vterm.monomial() == nullptr) tam.template addTerm<false>(id, Term<C>(vterm.coeff() * term.coeff(), mon));
169  else tam.template addTerm<false>(id, Term<C>(vterm.coeff() * term.coeff(), vterm.monomial() * mon));
170  }
171  }
172  else
173  {
174  tam.template addTerm<false>(id, term);
175  }
176  }
177  }
178  tam.readTerms(id, p.terms());
179  p.reset_ordered();
180  p.template makeMinimallyOrdered<false, true>();
181  assert(p.nr_terms() <= expectedResultSize);
182  assert(p.is_consistent());
183 }
184 
185 template<typename C, typename O, typename P>
188  substitute_inplace(result, var, value);
189  return result;
190 }
191 
192 template<typename C, typename O, typename P, typename S>
193 MultivariatePolynomial<C,O,P> substitute(const MultivariatePolynomial<C,O,P>& p, const std::map<Variable,S>& substitutions) {
194  static_assert(!std::is_same<S, Term<C>>::value, "Terms are handled by a separate method.");
197  auto id = tam.getId(p.nr_terms());
198  for (const auto& term: p) {
199  Term<C> resultTerm = substitute(term, substitutions);
200  if( !carl::is_zero(resultTerm) )
201  {
202  tam.template addTerm<false>(id, resultTerm );
203  }
204  }
205  tam.readTerms(id, result.terms());
206  result.reset_ordered();
207  result.template makeMinimallyOrdered<false, true>();
208  assert(result.is_consistent());
209  return result;
210 }
211 
212 template<typename C, typename O, typename P>
216  auto id = tam.getId(p.nr_terms());
217  for (const auto& term: p) {
218  tam.template addTerm<false>(id, substitute(term, substitutions));
219  }
220  tam.readTerms(id, result.terms());
221  result.reset_ordered();
222  result.template makeMinimallyOrdered<false, true>();
223  assert(result.is_consistent());
224  return result;
225 }
226 
227 template<typename C, typename O, typename P>
230  if (is_constant(p) || substitutions.empty())
231  {
232  return result;
233  }
234  // Substitute the variables, which have to be replaced by 0, beforehand,
235  // as this could significantly simplify this multivariate polynomial.
236  for (const auto& sub: substitutions) {
237  if(carl::is_zero(sub.second))
238  {
239  substitute_inplace(result, sub.first, sub.second);
240  if(is_constant(result))
241  {
242  assert(result.is_consistent());
243  return result;
244  }
245  }
246  }
247  // Find and sort all exponents occurring for all variables to substitute as basis.
248  std::map<std::pair<Variable, exponent>, MultivariatePolynomial<C,O,P>> expResults;
249  for(const auto& term: result)
250  {
251  if(term.monomial())
252  {
253  const Monomial& m = *(term.monomial());
254  CARL_LOG_TRACE("carl.core.monomial", "Iterating over " << m);
255  for(unsigned i = 0; i < m.num_variables(); ++i)
256  {
257  CARL_LOG_TRACE("carl.core.monomial", "Iterating: " << m[i].first);
258  if(m[i].second > 1 && substitutions.find(m[i].first) != substitutions.end())
259  {
261  }
262  }
263  }
264  }
265  // Calculate the exponentiation of the multivariate polynomial to substitute the
266  // for variables for, reusing the already calculated exponentiations.
267  if(!expResults.empty())
268  {
269  auto expResultA = expResults.begin();
270  auto expResultB = expResultA;
271  auto sub = substitutions.begin();
272  while (sub->first != expResultB->first.first)
273  {
274  assert(sub != substitutions.end());
275  ++sub;
276  }
277  expResultB->second = carl::pow(sub->second, expResultB->first.second);
278  ++expResultB;
279  while(expResultB != expResults.end())
280  {
281  if(expResultA->first.first != expResultB->first.first)
282  {
283  ++sub;
284  assert(sub != substitutions.end());
285  // Go to the next variable.
286  while (sub->first != expResultB->first.first)
287  {
288  assert(sub != substitutions.end());
289  ++sub;
290  }
291  assert(sub->first == expResultB->first.first);
292  expResultB->second = pow(sub->second, expResultB->first.second);
293  }
294  else
295  {
296  expResultB->second = expResultA->second * pow(sub->second, expResultB->first.second-expResultA->first.second);
297  }
298  ++expResultA;
299  ++expResultB;
300  }
301  }
303  // Substitute the variable for which all occurring exponentiations are calculated.
304  for(const auto& term: result)
305  {
306  MultivariatePolynomial<C,O,P> termResult(term.coeff());
307  if (term.monomial())
308  {
309  const Monomial& m = *(term.monomial());
310  CARL_LOG_TRACE("carl.core.monomial", "Iterating over " << m);
311  for(unsigned i = 0; i < m.num_variables(); ++i)
312  {
313  CARL_LOG_TRACE("carl.core.monomial", "Iterating: " << m[i].first);
314  if(m[i].second == 1)
315  {
316  auto iter = substitutions.find(m[i].first);
317  if(iter != substitutions.end())
318  {
319  termResult *= iter->second;
320  }
321  else
322  {
323  termResult *= m[i].first;
324  }
325  }
326  else
327  {
328  auto iter = expResults.find(m[i]);
329  if(iter != expResults.end())
330  {
331  termResult *= iter->second;
332  }
333  else
334  {
335  termResult *= Term<C>(constant_one<C>::get(), m[i].first, m[i].second);
336  }
337  }
338  }
339  }
340  resultB += termResult;
341  }
342  assert(resultB.is_consistent());
343  return resultB;
344 }
345 
346 
347 template<typename Coeff>
349  if (carl::is_zero(p)) return;
350  if (var == p.main_var()) {
352  } else if constexpr (!is_number_type<Coeff>::value) {
353  // Coefficients from a polynomial ring
354  if (value.has(var)) {
355  // Fall back to multivariate substitution.
357  substitute_inplace(tmp, var, value);
359  } else {
360  // Safely substitute into each coefficient separately
361  for (auto& c: p.coefficients()) {
362  substitute_inplace(c, var, value);
363  }
364  }
365  }
367  assert(p.is_consistent());
368 }
369 
370 template<typename Coeff>
372  if constexpr (is_number_type<Coeff>::value) {
373  if (var == p.main_var()) {
374  return UnivariatePolynomial<Coeff>(p.main_var(), p.evaluate(value));
375  }
376  return p;
377  } else {
378  if (var == p.main_var()) {
380  for (const auto& c: p.coefficients()) {
381  res += substitute(c, var, value);
382  }
383  CARL_LOG_TRACE("carl.core.uvpolynomial", p << " [ " << var << " -> " << value << " ] = " << res);
384  return res;
385  } else {
386  if (value.has(var)) {
387  // Fall back to multivariate substitution.
389  substitute_inplace(tmp, var, value);
390  return to_univariate_polynomial(tmp, p.main_var());
391  } else {
392  std::vector<Coeff> res(p.coefficients().size());
393  for (std::size_t i = 0; i < res.size(); i++) {
394  res[i] = substitute(p.coefficients()[i], var, value);
395  }
396  UnivariatePolynomial<Coeff> resp(p.main_var(), res);
397  resp.strip_leading_zeroes();
398  CARL_LOG_TRACE("carl.core.uvpolynomial", p << " [ " << var << " -> " << value << " ] = " << resp);
399  return resp;
400  }
401  }
402  }
403 }
404 
405 /**
406  * Substitutes a variable with a rational within a polynomial.
407  */
408 template<typename Rational>
411 }
412 template<typename Poly, typename Rational>
414  carl::substitute_inplace(p, var, Poly(r));
415 }
416 
417 }
Implements utility functions concerning the (total) degree of monomials, terms and polynomials.
#define CARL_LOG_FUNC(channel, args)
Definition: carl-logging.h:46
#define CARL_LOG_TRACE(channel, msg)
Definition: carl-logging.h:44
mpq_class Rational
Definition: HornerTest.cpp:12
carl is the main namespace for the library.
Coeff content(const UnivariatePolynomial< Coeff > &p)
The content of a polynomial is the gcd of the coefficients of the normal part of a polynomial.
Definition: Content.h:22
Monomial::Arg createMonomial(T &&... t)
Definition: MonomialPool.h:168
UnivariatePolynomial< C > to_univariate_polynomial(const MultivariatePolynomial< C, O, P > &p)
Convert a univariate polynomial that is currently (mis)represented by a 'MultivariatePolynomial' into...
bool is_constant(const ContextPolynomial< Coeff, Ordering, Policies > &p)
std::size_t exponent
Type of an exponent.
Definition: Monomial.h:29
bool is_zero(const Interval< Number > &i)
Check if this interval is a point-interval containing 0.
Definition: Interval.h:1453
bool evaluate(const BasicConstraint< Poly > &c, const Assignment< Number > &m)
Definition: Evaluation.h:10
typename UnderlyingNumberType< P >::type Coeff
Coeff substitute(const Monomial &m, const std::map< Variable, Coeff > &substitutions)
Applies the given substitutions to a monomial.
Definition: Substitution.h:19
void substitute_inplace(MultivariateRoot< Poly > &mr, Variable var, const Poly &poly)
Create a copy of the underlying polynomial with the given variable replaced by the given polynomial.
Interval< Number > pow(const Interval< Number > &i, Integer exp)
Definition: Power.h:11
A Variable represents an algebraic variable that can be used throughout carl.
Definition: Variable.h:85
static const T & get()
Definition: constants.h:51
This class represents a univariate polynomial with coefficients of an arbitrary type.
const std::vector< Coefficient > & coefficients() const &
Retrieves the coefficients defining this polynomial.
Variable main_var() const
Retrieves the main variable of this polynomial.
bool is_consistent() const
Asserts that this polynomial over numeric coefficients complies with the requirements and assumptions...
Coefficient evaluate(const Coefficient &value) const
The general-purpose multivariate polynomial class.
const TermsType & terms() const
bool is_consistent() const
Asserts that this polynomial complies with the requirements and assumptions for MultivariatePolynomia...
std::vector< Term< Coeff > > TermsType
Type our terms vector.f.
std::size_t nr_terms() const
Calculate the number of terms.
States if a type is a number type.
Definition: typetraits.h:237
The general-purpose monomials.
Definition: Monomial.h:59
std::shared_ptr< const Monomial > Arg
Definition: Monomial.h:62
exponents_it end()
Returns past-the-end iterator.
Definition: Monomial.h:154
std::vector< std::pair< Variable, std::size_t > > Content
Definition: Monomial.h:63
std::size_t num_variables() const
Returns the number of variables that occur in the monomial.
Definition: Monomial.h:245
Coefficient & coeff()
Get the coefficient.
Definition: Term.h:80
Monomial::Arg & monomial()
Get the monomial.
Definition: Term.h:91