carl  24.04
Computer ARithmetic Library
UnivariatePolynomial.tpp
Go to the documentation of this file.
1 /**
2  * @file: UnivariatePolynomial.tpp
3  * @author: Sebastian Junges
4  *
5  * @since August 26, 2013
6  */
7 
8 #pragma once
9 
10 #include <carl-common/debug/debug.h>
11 #include <carl-common/meta/platform.h>
12 #include <carl-common/meta/SFINAE.h>
13 #include <carl-logging/carl-logging.h>
14 #include "MultivariatePolynomial.h"
15 #include <carl-arith/core/Sign.h>
16 
17 #include "functions/Derivative.h"
18 #include "functions/Division.h"
19 
20 #include <algorithm>
21 #include <iomanip>
22 
23 namespace carl
24 {
25 
26 template<typename Coeff>
27 UnivariatePolynomial<Coeff>::UnivariatePolynomial(const UnivariatePolynomial& p):
28  mMainVar(p.mMainVar), mCoefficients(p.mCoefficients)
29 {
30  assert(this->is_consistent());
31 }
32 
33 template<typename Coeff>
34 UnivariatePolynomial<Coeff>::UnivariatePolynomial(UnivariatePolynomial&& p) noexcept:
35  mMainVar(p.mMainVar), mCoefficients()
36 {
37  mCoefficients = std::move(p.mCoefficients);
38  assert(is_consistent());
39 }
40 
41 template<typename Coeff>
42 UnivariatePolynomial<Coeff>& UnivariatePolynomial<Coeff>::operator=(const UnivariatePolynomial& p) {
43  mMainVar = p.mMainVar;
44  mCoefficients = p.mCoefficients;
45  assert(is_consistent());
46  return *this;
47 }
48 
49 template<typename Coeff>
50 UnivariatePolynomial<Coeff>& UnivariatePolynomial<Coeff>::operator=(UnivariatePolynomial&& p) noexcept {
51  mMainVar = p.mMainVar;
52  mCoefficients = std::move(p.mCoefficients);
53  assert(is_consistent());
54  return *this;
55 }
56 
57 template<typename Coeff>
58 UnivariatePolynomial<Coeff>::UnivariatePolynomial(Variable mainVar)
59 : mMainVar(mainVar), mCoefficients()
60 {
61  assert(this->is_consistent());
62 }
63 template<typename Coeff>
64 UnivariatePolynomial<Coeff>::UnivariatePolynomial(Variable mainVar, const Coeff& coeff, std::size_t degree) :
65 mMainVar(mainVar),
66 mCoefficients(degree+1, Coeff(0)) // We would like to use 0 here, but Coeff(0) is not always constructable (some methods need more parameter)
67 {
68  if(coeff != Coeff(0))
69  {
70  mCoefficients[degree] = coeff;
71  }
72  else
73  {
74  mCoefficients.clear();
75  }
76  strip_leading_zeroes();
77  assert(is_consistent());
78 }
79 
80 template<typename Coeff>
81 UnivariatePolynomial<Coeff>::UnivariatePolynomial(Variable mainVar, std::initializer_list<Coeff> coefficients)
82 : mMainVar(mainVar), mCoefficients(coefficients)
83 {
84  this->strip_leading_zeroes();
85  assert(this->is_consistent());
86 }
87 
88 template<typename Coeff>
89 template<typename C, DisableIf<std::is_same<C, typename UnderlyingNumberType<C>::type>>>
90 UnivariatePolynomial<Coeff>::UnivariatePolynomial(Variable mainVar, std::initializer_list<typename UnderlyingNumberType<C>::type> coefficients)
91 : mMainVar(mainVar), mCoefficients()
92 {
93  for (auto c: coefficients) {
94  this->mCoefficients.push_back(Coeff(c));
95  }
96  this->strip_leading_zeroes();
97  assert(this->is_consistent());
98 }
99 
100 template<typename Coeff>
101 UnivariatePolynomial<Coeff>::UnivariatePolynomial(Variable mainVar, const std::vector<Coeff>& coefficients)
102 : mMainVar(mainVar), mCoefficients(coefficients)
103 {
104  this->strip_leading_zeroes();
105  assert(this->is_consistent());
106 }
107 
108 template<typename Coeff>
109 UnivariatePolynomial<Coeff>::UnivariatePolynomial(Variable mainVar, std::vector<Coeff>&& coefficients)
110 : mMainVar(mainVar), mCoefficients(coefficients)
111 {
112  this->strip_leading_zeroes();
113  assert(this->is_consistent());
114 }
115 
116 template<typename Coeff>
117 UnivariatePolynomial<Coeff>::UnivariatePolynomial(Variable mainVar, const std::map<uint, Coeff>& coefficients)
118 : mMainVar(mainVar)
119 {
120  mCoefficients.reserve(coefficients.rbegin()->first);
121  for (const auto& expAndCoeff : coefficients)
122  {
123  if(expAndCoeff.first != mCoefficients.size())
124  {
125  mCoefficients.resize(expAndCoeff.first, Coeff(0));
126  }
127  mCoefficients.push_back(expAndCoeff.second);
128  }
129  this->strip_leading_zeroes();
130  assert(this->is_consistent());
131 }
132 
133 template<typename Coeff>
134 Coeff UnivariatePolynomial<Coeff>::evaluate(const Coeff& value) const
135 {
136  Coeff result(0);
137  Coeff var(1);
138  for(const Coeff& coeff : mCoefficients)
139  {
140  result += (coeff * var);
141  var *= value;
142  }
143  return result;
144 }
145 
146 template<typename Coeff>
147 bool UnivariatePolynomial<Coeff>::is_normal() const
148 {
149  return unit_part() == Coeff(1);
150 }
151 
152 template<typename Coefficient>
153 UnivariatePolynomial<Coefficient>& UnivariatePolynomial<Coefficient>::mod(const Coefficient& modulus)
154 {
155  for(Coefficient& coeff : mCoefficients)
156  {
157  coeff = carl::mod(coeff, modulus);
158  }
159  return *this;
160 }
161 
162 template<typename Coefficient>
163 UnivariatePolynomial<Coefficient> UnivariatePolynomial<Coefficient>::mod(const Coefficient& modulus) const
164 {
165  UnivariatePolynomial<Coefficient> result;
166  result.mCoefficients.reserve(mCoefficients.size());
167  for(const Coefficient& coeff : mCoefficients)
168  {
169  result.mCoefficients.push_back(mod(coeff, modulus));
170  }
171  result.strip_leading_zeroes();
172  return result;
173 }
174 
175 template<typename Coeff>
176 UnivariatePolynomial<Coeff> UnivariatePolynomial<Coeff>::pow(std::size_t exp) const
177 {
178  if (exp == 0) return UnivariatePolynomial(mMainVar, constant_one<Coeff>::get());
179  if (carl::is_zero(*this)) return UnivariatePolynomial(mMainVar);
180  UnivariatePolynomial<Coeff> res(mMainVar, constant_one<Coeff>::get());
181  UnivariatePolynomial<Coeff> mult(*this);
182  while(exp > 0) {
183  if ((exp & 1) != 0) res *= mult;
184  exp /= 2;
185  if(exp > 0) mult = mult * mult;
186  }
187  return res;
188 }
189 
190 template<typename Coeff>
191 template<typename C, DisableIf<is_number_type<C>>>
192 UnivariatePolynomial<typename UnivariatePolynomial<Coeff>::NumberType> UnivariatePolynomial<Coeff>::toNumberCoefficients() const {
193  std::vector<NumberType> coeffs;
194  coeffs.reserve(this->mCoefficients.size());
195  for (auto c: this->mCoefficients) {
196  assert(c.is_constant());
197  coeffs.push_back(c.constant_part());
198  }
199  return UnivariatePolynomial<NumberType>(this->mMainVar, coeffs);
200 }
201 
202 template<typename Coeff>
203 template<typename NewCoeff>
204 UnivariatePolynomial<NewCoeff> UnivariatePolynomial<Coeff>::convert() const {
205  std::vector<NewCoeff> coeffs;
206  coeffs.resize(this->mCoefficients.size());
207  for (std::size_t i = 0; i < this->mCoefficients.size(); i++) {
208  coeffs[i] = NewCoeff(this->mCoefficients[i]);
209  }
210  return UnivariatePolynomial<NewCoeff>(this->mMainVar, coeffs);
211 }
212 
213 template<typename Coeff>
214 template<typename NewCoeff>
215 UnivariatePolynomial<NewCoeff> UnivariatePolynomial<Coeff>::convert(const std::function<NewCoeff(const Coeff&)>& f) const {
216  std::vector<NewCoeff> coeffs;
217  coeffs.resize(this->mCoefficients.size());
218  for (std::size_t i = 0; i < this->mCoefficients.size(); i++) {
219  coeffs[i] = f(this->mCoefficients[i]);
220  }
221  return UnivariatePolynomial<NewCoeff>(this->mMainVar, coeffs);
222 }
223 
224 template<typename Coeff>
225 UnivariatePolynomial<Coeff> UnivariatePolynomial<Coeff>::normalized() const
226 {
227  if(carl::is_zero(*this))
228  {
229  return *this;
230  }
231  return *this/unit_part();
232 }
233 
234 template<typename Coeff>
235 Coeff UnivariatePolynomial<Coeff>::unit_part() const {
236  if constexpr (is_field_type<Coeff>::value) {
237  // Coeffs from a field
238  return lcoeff();
239  } else if constexpr (is_number_type<Coeff>::value) {
240  // Coeffs from a number ring
241  if (carl::is_zero(*this) || lcoeff() > Coeff(0)) {
242  return Coeff(1);
243  } else {
244  return Coeff(-1);
245  }
246  } else {
247  // Coeffs from a polynomial ring
248  if (carl::is_zero(*this) || carl::is_zero(lcoeff()) || lcoeff().lcoeff() > NumberType(0)) {
249  return Coeff(1);
250  } else {
251  return Coeff(-1);
252  }
253  }
254 }
255 
256 template<typename Coeff>
257 template<typename C, EnableIf<is_subset_of_rationals_type<C>>>
258 Coeff UnivariatePolynomial<Coeff>::coprime_factor() const
259 {
260  assert(!carl::is_zero(*this));
261  auto it = mCoefficients.begin();
262  IntNumberType num = get_num(*it);
263  IntNumberType den = get_denom(*it);
264  for (++it; it != mCoefficients.end(); ++it) {
265  num = carl::gcd(num, get_num(*it));
266  den = carl::lcm(den, get_denom(*it));
267  }
268  return Coeff(den)/Coeff(num);
269 }
270 
271 template<typename Coeff>
272 template<typename C, DisableIf<is_subset_of_rationals_type<C>>>
273 typename UnderlyingNumberType<Coeff>::type UnivariatePolynomial<Coeff>::coprime_factor() const
274 {
275  assert(!carl::is_zero(*this));
276  auto it = mCoefficients.begin();
277  typename UnderlyingNumberType<Coeff>::type factor = it->coprime_factor();
278  for (++it; it != mCoefficients.end(); ++it) {
279  factor = carl::lcm(factor, it->coprime_factor());
280  }
281  return factor;
282 }
283 
284 template<typename Coeff>
285 template<typename C, EnableIf<is_subset_of_rationals_type<C>>>
286 UnivariatePolynomial<typename IntegralType<Coeff>::type> UnivariatePolynomial<Coeff>::coprime_coefficients() const
287 {
288  CARL_LOG_TRACE("carl.core", *this << " .coprime_coefficients()");
289  // Notice that even if factor is 1, we create a new polynomial
290  UnivariatePolynomial<typename IntegralType<Coeff>::type> result(mMainVar);
291  if (carl::is_zero(*this)) {
292  return result;
293  }
294  result.mCoefficients.reserve(mCoefficients.size());
295  Coeff factor = this->coprime_factor();
296  for (const Coeff& coeff: mCoefficients) {
297  assert(get_denom(coeff * factor) == 1);
298  result.mCoefficients.push_back(get_num(coeff * factor));
299  }
300  return result;
301 }
302 
303 template<typename Coeff>
304 template<typename C, DisableIf<is_subset_of_rationals_type<C>>>
305 UnivariatePolynomial<Coeff> UnivariatePolynomial<Coeff>::coprime_coefficients() const {
306  if (carl::is_zero(*this)) return *this;
307  UnivariatePolynomial<Coeff> result(mMainVar);
308  result.mCoefficients.reserve(mCoefficients.size());
309  auto factor = this->coprime_factor();
310  for (const Coeff& c: mCoefficients) {
311  result.mCoefficients.push_back(factor * c);
312  }
313  return result;
314 }
315 
316 template<typename Coeff>
317 template<typename C, EnableIf<is_subset_of_rationals_type<C>>>
318 UnivariatePolynomial<typename IntegralType<Coeff>::type> UnivariatePolynomial<Coeff>::coprime_coefficients_sign_preserving() const
319 {
320  CARL_LOG_TRACE("carl.core", *this << " .coprime_coefficients_sign_preserving()");
321  // Notice that even if factor is 1, we create a new polynomial
322  UnivariatePolynomial<typename IntegralType<Coeff>::type> result(mMainVar);
323  if (carl::is_zero(*this)) {
324  return result;
325  }
326  result.mCoefficients.reserve(mCoefficients.size());
327  Coeff factor = carl::abs(this->coprime_factor());
328  for (const Coeff& coeff: mCoefficients) {
329  assert(get_denom(coeff * factor) == 1);
330  result.mCoefficients.push_back(get_num(coeff * factor));
331  }
332  return result;
333 }
334 
335 template<typename Coeff>
336 template<typename C, DisableIf<is_subset_of_rationals_type<C>>>
337 UnivariatePolynomial<Coeff> UnivariatePolynomial<Coeff>::coprime_coefficients_sign_preserving() const {
338  if (this->is_zero()) return *this;
339  UnivariatePolynomial<Coeff> result(mMainVar);
340  result.mCoefficients.reserve(mCoefficients.size());
341  auto factor = carl::abs(this->coprime_factor());
342  for (const Coeff& c: mCoefficients) {
343  result.mCoefficients.push_back(factor * c);
344  }
345  return result;
346 }
347 
348 template<typename Coeff>
349 bool UnivariatePolynomial<Coeff>::divides(const UnivariatePolynomial& divisor) const
350 {
351  ///@todo Is this correct?
352  return carl::is_zero(divisor.divideBy(*this).remainder);
353 }
354 
355 template<typename Coeff>
356 template<typename C, EnableIf<is_instantiation_of<GFNumber, C>>>
357 UnivariatePolynomial<typename IntegralType<Coeff>::type> UnivariatePolynomial<Coeff>::to_integer_domain() const
358 {
359  UnivariatePolynomial<typename IntegralType<Coeff>::type> res(mMainVar);
360  res.mCoefficients.reserve(mCoefficients.size());
361  for(const Coeff& c : mCoefficients)
362  {
363  assert(carl::is_integer(c));
364  res.mCoefficients.push_back(c.representing_integer());
365  }
366  res.strip_leading_zeroes();
367  return res;
368 }
369 
370 template<typename Coeff>
371 template<typename C, DisableIf<is_instantiation_of<GFNumber, C>>>
372 UnivariatePolynomial<typename IntegralType<Coeff>::type> UnivariatePolynomial<Coeff>::to_integer_domain() const
373 {
374  UnivariatePolynomial<typename IntegralType<Coeff>::type> res(mMainVar);
375  res.mCoefficients.reserve(mCoefficients.size());
376  for(const Coeff& c : mCoefficients)
377  {
378  assert(carl::is_integer(c));
379  res.mCoefficients.push_back(get_num(c));
380  }
381  res.strip_leading_zeroes();
382  return res;
383 }
384 
385 template<typename Coeff>
386 //template<typename T = Coeff, EnableIf<!std::is_same<IntegralType<Coeff>, bool>::value>>
387 UnivariatePolynomial<GFNumber<typename IntegralType<Coeff>::type>> UnivariatePolynomial<Coeff>::toFiniteDomain(const GaloisField<typename IntegralType<Coeff>::type>* galoisField) const
388 {
389  UnivariatePolynomial<GFNumber<typename IntegralType<Coeff>::type>> res(mMainVar);
390  res.mCoefficients.reserve(mCoefficients.size());
391  for(const Coeff& c : mCoefficients)
392  {
393  assert(carl::is_integer(c));
394  res.mCoefficients.push_back(GFNumber<typename IntegralType<Coeff>::type>(c,galoisField));
395  }
396  res.strip_leading_zeroes();
397  return res;
398 
399 }
400 
401 template<typename Coeff>
402 template<typename N, EnableIf<is_subset_of_rationals_type<N>>>
403 typename UnivariatePolynomial<Coeff>::NumberType UnivariatePolynomial<Coeff>::numeric_content() const
404 {
405  if (carl::is_zero(*this)) return NumberType(0);
406  // Obtain main denominator for all coefficients.
407  IntNumberType mainDenom = this->main_denom();
408 
409  // now, some coefficient * mainDenom is always integral.
410  // we convert such a product to an integral data type by get_num()
411  UnivariatePolynomial<Coeff>::NumberType c = this->numeric_content(0) * mainDenom;
412  assert(get_denom(c) == 1);
413  IntNumberType res = get_num(c);
414  for (std::size_t i = 1; i < this->mCoefficients.size(); i++) {
415  c = this->numeric_content(i) * mainDenom;
416  assert(get_denom(c) == 1);
417  res = carl::gcd(get_num(c), res);
418  }
419  return NumberType(res) / mainDenom;
420 }
421 
422 template<typename Coeff>
423 template<typename C, EnableIf<is_number_type<C>>>
424 typename UnivariatePolynomial<Coeff>::IntNumberType UnivariatePolynomial<Coeff>::main_denom() const
425 {
426  IntNumberType denom = 1;
427  for (const auto& c: mCoefficients) {
428  denom = carl::lcm(denom, get_denom(c));
429  }
430  CARL_LOG_TRACE("carl.core", "mainDenom of " << *this << " is " << denom);
431  return denom;
432 }
433 template<typename Coeff>
434 template<typename C, DisableIf<is_number_type<C>>>
435 typename UnivariatePolynomial<Coeff>::IntNumberType UnivariatePolynomial<Coeff>::main_denom() const
436 {
437  IntNumberType denom = 1;
438  for (const auto& c: mCoefficients) {
439  denom = carl::lcm(denom, c.main_denom());
440  }
441  return denom;
442 }
443 
444 template<typename Coeff>
445 Coeff UnivariatePolynomial<Coeff>::synthetic_division(const Coeff& zeroOfDivisor)
446 {
447  if(coefficients().empty()) return Coeff(0);
448  if(coefficients().size() == 1) return coefficients().back();
449  std::vector<Coeff> secondRow;
450  secondRow.reserve(coefficients().size());
451  secondRow.push_back(Coeff(0));
452  std::vector<Coeff> thirdRow(coefficients().size(), Coeff(0));
453  size_t posThirdRow = coefficients().size()-1;
454  auto coeff = coefficients().rbegin();
455  thirdRow[posThirdRow] = (*coeff) + secondRow.front();
456  ++coeff;
457  while(coeff != coefficients().rend())
458  {
459  secondRow.push_back(zeroOfDivisor*thirdRow[posThirdRow]);
460  --posThirdRow;
461  thirdRow[posThirdRow] = (*coeff) + secondRow.back();
462  ++coeff;
463  }
464  assert(posThirdRow == 0);
465  CARL_LOG_TRACE("carl.core.upoly", "UnivSynDiv: (" << *this << ")[x -> " << zeroOfDivisor << "] = " << thirdRow.front());
466  if(thirdRow.front() == 0)
467  {
468  thirdRow.erase(thirdRow.begin());
469  this->mCoefficients.swap(thirdRow);
470  CARL_LOG_TRACE("carl.core.upoly", "UnivSynDiv: reduced by ((" << carl::abs(get_denom(thirdRow.front())) << ")*" << main_var() << " + (" << (thirdRow.front()<0 ? "-" : "") << carl::abs(get_num(thirdRow.front())) << ")) -> " << *this);
471  return Coeff(0);
472  }
473  return thirdRow.front();
474 }
475 
476 template<typename Coeff>
477 UnivariatePolynomial<Coeff> UnivariatePolynomial<Coeff>::operator -() const
478 {
479  UnivariatePolynomial result(mMainVar);
480  result.mCoefficients.reserve(mCoefficients.size());
481  for(const auto& c : mCoefficients) {
482  result.mCoefficients.push_back(-c);
483  }
484 
485  return result;
486 }
487 
488 template<typename Coefficient>
489 UnivariatePolynomial<Coefficient>& UnivariatePolynomial<Coefficient>::operator+=(const Coefficient& rhs)
490 {
491  if(rhs == Coefficient(0)) return *this;
492  if(mCoefficients.empty())
493  {
494  // Adding non-zero rhs to zero.
495  mCoefficients.resize(1, rhs);
496  }
497  else
498  {
499  mCoefficients.front() += rhs;
500  if(mCoefficients.size() == 1 && mCoefficients.front() == Coefficient(0))
501  {
502  // Result is zero.
503  mCoefficients.clear();
504  }
505  }
506  return *this;
507 }
508 
509 template<typename Coeff>
510 UnivariatePolynomial<Coeff>& UnivariatePolynomial<Coeff>::operator+=(const UnivariatePolynomial& rhs)
511 {
512  assert(mMainVar == rhs.mMainVar);
513 
514  if(carl::is_zero(rhs))
515  {
516  return *this;
517  }
518 
519  if(mCoefficients.size() < rhs.mCoefficients.size())
520  {
521  for(std::size_t i = 0; i < mCoefficients.size(); ++i)
522  {
523  mCoefficients[i] += rhs.mCoefficients[i];
524  }
525  mCoefficients.insert(mCoefficients.end(), rhs.mCoefficients.end() - sint(rhs.mCoefficients.size() - mCoefficients.size()), rhs.mCoefficients.end());
526  }
527  else
528  {
529  for(std::size_t i = 0; i < rhs.mCoefficients.size(); ++i)
530  {
531  mCoefficients[i] += rhs.mCoefficients[i];
532  }
533  }
534  strip_leading_zeroes();
535  return *this;
536 }
537 
538 template<typename C>
539 UnivariatePolynomial<C> operator+(const UnivariatePolynomial<C>& lhs, const UnivariatePolynomial<C>& rhs)
540 {
541  UnivariatePolynomial<C> res(lhs);
542  res += rhs;
543  return res;
544 }
545 
546 template<typename C>
547 UnivariatePolynomial<C> operator+(const UnivariatePolynomial<C>& lhs, const C& rhs)
548 {
549  UnivariatePolynomial<C> res(lhs);
550  res += rhs;
551  return res;
552 }
553 
554 template<typename C>
555 UnivariatePolynomial<C> operator+(const C& lhs, const UnivariatePolynomial<C>& rhs)
556 {
557  return rhs + lhs;
558 }
559 
560 
561 template<typename Coefficient>
562 UnivariatePolynomial<Coefficient>& UnivariatePolynomial<Coefficient>::operator-=(const Coefficient& rhs)
563 {
564 // CARL_LOG_INEFFICIENT();
565  return *this += -rhs;
566 }
567 
568 template<typename Coeff>
569 UnivariatePolynomial<Coeff>& UnivariatePolynomial<Coeff>::operator-=(const UnivariatePolynomial& rhs)
570 {
571 // CARL_LOG_INEFFICIENT();
572  return *this += -rhs;
573 }
574 
575 
576 template<typename C>
577 UnivariatePolynomial<C> operator-(const UnivariatePolynomial<C>& lhs, const UnivariatePolynomial<C>& rhs)
578 {
579  UnivariatePolynomial<C> res(lhs);
580  res -= rhs;
581  return res;
582 }
583 
584 template<typename C>
585 UnivariatePolynomial<C> operator-(const UnivariatePolynomial<C>& lhs, const C& rhs)
586 {
587  UnivariatePolynomial<C> res(lhs);
588  res -= rhs;
589  return res;
590 }
591 
592 template<typename C>
593 UnivariatePolynomial<C> operator-(const C& lhs, const UnivariatePolynomial<C>& rhs)
594 {
595  return rhs - lhs;
596 }
597 
598 template<typename Coefficient>
599 template<typename C, EnableIf<is_number_type<C>>>
600 UnivariatePolynomial<Coefficient>& UnivariatePolynomial<Coefficient>::operator*=(Variable rhs) {
601  if (rhs == this->mMainVar) {
602  this->mCoefficients.insert(this->mCoefficients.begin(), Coefficient(0));
603  return *this;
604  }
605  assert(!carl::is_number_type<Coefficient>::value);
606  return *this;
607 }
608 template<typename Coefficient>
609 template<typename C, DisableIf<is_number_type<C>>>
610 UnivariatePolynomial<Coefficient>& UnivariatePolynomial<Coefficient>::operator*=(Variable rhs) {
611  if (rhs == this->mMainVar) {
612  this->mCoefficients.insert(this->mCoefficients.begin(), Coefficient(0));
613  return *this;
614  }
615  assert(!carl::is_number_type<Coefficient>::value);
616  for (auto& c: this->mCoefficients) c *= rhs;
617  return *this;
618 }
619 
620 template<typename Coefficient>
621 UnivariatePolynomial<Coefficient>& UnivariatePolynomial<Coefficient>::operator*=(const Coefficient& rhs)
622 {
623  if(rhs == Coefficient(0))
624  {
625  mCoefficients.clear();
626  return *this;
627  }
628  for(Coefficient& c : mCoefficients)
629  {
630  c *= rhs;
631  }
632 
633  if(is_finite_type<Coefficient>::value)
634  {
635  strip_leading_zeroes();
636  }
637 
638  return *this;
639 }
640 
641 
642 template<typename Coeff>
643 template<typename I, DisableIf<std::is_same<Coeff, I>>...>
644 UnivariatePolynomial<Coeff>& UnivariatePolynomial<Coeff>::operator*=(const typename IntegralType<Coeff>::type& rhs)
645 {
646  static_assert(std::is_same<Coeff, I>::value, "Do not provide template parameters");
647  if(rhs == I(0))
648  {
649  mCoefficients.clear();
650  return *this;
651  }
652  for(Coeff& c : mCoefficients)
653  {
654  c *= rhs;
655  }
656  return *this;
657 }
658 
659 template<typename Coeff>
660 UnivariatePolynomial<Coeff>& UnivariatePolynomial<Coeff>::operator*=(const UnivariatePolynomial& rhs)
661 {
662  assert(mMainVar == rhs.mMainVar);
663  if(carl::is_zero(rhs))
664  {
665  mCoefficients.clear();
666  return *this;
667  }
668 
669  std::vector<Coeff> newCoeffs;
670  newCoeffs.reserve(mCoefficients.size() + rhs.mCoefficients.size());
671  for(std::size_t e = 0; e < mCoefficients.size() + rhs.degree(); ++e)
672  {
673  newCoeffs.push_back(Coeff(0));
674  for(std::size_t i = 0; i < mCoefficients.size() && i <= e; ++i)
675  {
676  if(e - i < rhs.mCoefficients.size())
677  {
678  newCoeffs.back() += mCoefficients[i] * rhs.mCoefficients[e-i];
679  }
680  }
681  }
682  mCoefficients.swap(newCoeffs);
683  strip_leading_zeroes();
684  return *this;
685 }
686 
687 
688 template<typename C>
689 UnivariatePolynomial<C> operator*(const UnivariatePolynomial<C>& lhs, const UnivariatePolynomial<C>& rhs)
690 {
691  UnivariatePolynomial<C> res(lhs);
692  res *= rhs;
693  return res;
694 }
695 
696 template<typename C>
697 UnivariatePolynomial<C> operator*(const UnivariatePolynomial<C>& lhs, Variable rhs) {
698  return std::move(UnivariatePolynomial<C>(lhs) *= rhs);
699 }
700 template<typename C>
701 UnivariatePolynomial<C> operator*(Variable lhs, const UnivariatePolynomial<C>& rhs) {
702  return std::move(UnivariatePolynomial<C>(rhs) *= lhs);
703 }
704 
705 template<typename C>
706 UnivariatePolynomial<C> operator*(const UnivariatePolynomial<C>& lhs, const C& rhs)
707 {
708  UnivariatePolynomial<C> res(lhs);
709  res *= rhs;
710  return res;
711 }
712 
713 template<typename C>
714 UnivariatePolynomial<C> operator*(const C& lhs, const UnivariatePolynomial<C>& rhs)
715 {
716  return rhs * lhs;
717 }
718 
719 template<typename C>
720 UnivariatePolynomial<C> operator*(const UnivariatePolynomial<C>& lhs, const IntegralTypeIfDifferent<C>& rhs)
721 {
722  UnivariatePolynomial<C> res(lhs);
723  res *= rhs;
724  return res;
725 }
726 
727 template<typename C>
728 UnivariatePolynomial<C> operator*(const IntegralTypeIfDifferent<C>& lhs, const UnivariatePolynomial<C>& rhs)
729 {
730  return rhs * lhs;
731 }
732 
733 
734 
735 template<typename Coeff>
736 template<typename C, EnableIf<is_field_type<C>>>
737 UnivariatePolynomial<Coeff>& UnivariatePolynomial<Coeff>::operator/=(const Coeff& rhs)
738 {
739  assert(rhs != Coeff(0));
740  for(Coeff& c : mCoefficients)
741  {
742  c /= rhs;
743  }
744  return *this;
745 }
746 
747 template<typename Coeff>
748 template<typename C, DisableIf<is_field_type<C>>>
749 UnivariatePolynomial<Coeff>& UnivariatePolynomial<Coeff>::operator/=(const Coeff& rhs)
750 {
751  assert(rhs != Coeff(0));
752  for(Coeff& c : mCoefficients)
753  {
754  c = quotient(c, rhs);
755  }
756  /// TODO not fully sure whether this is necessary
757  this->strip_leading_zeroes();
758  return *this;
759 }
760 
761 
762 
763 template<typename C>
764 UnivariatePolynomial<C> operator/(const UnivariatePolynomial<C>& lhs, const C& rhs)
765 {
766  assert(rhs != C(0));
767  if(carl::is_zero(lhs)) return lhs;
768  UnivariatePolynomial<C> res(lhs);
769  return res /= rhs;
770 }
771 
772 template<typename C>
773 bool operator==(const UnivariatePolynomial<C>& lhs, const UnivariatePolynomial<C>& rhs)
774 {
775  assert(lhs.is_consistent());
776  assert(rhs.is_consistent());
777  if(lhs.mMainVar == rhs.mMainVar)
778  {
779  return lhs.mCoefficients == rhs.mCoefficients;
780  }
781  else
782  {
783  // in different variables, polynomials can still be equal if constant.
784  if(carl::is_zero(lhs) && carl::is_zero(rhs)) return true;
785  if ((lhs.mCoefficients.size() == 1) && (rhs.mCoefficients.size() == 1) && (lhs.mCoefficients == rhs.mCoefficients)) return true;
786  // Convert to multivariate and compare that.
787  return MultivariatePolynomial<typename carl::UnderlyingNumberType<C>::type>(lhs) == MultivariatePolynomial<typename carl::UnderlyingNumberType<C>::type>(rhs);
788  }
789 }
790 template<typename C>
791 bool operator==(const UnivariatePolynomialPtr<C>& lhs, const UnivariatePolynomialPtr<C>& rhs)
792 {
793  if (lhs == nullptr && rhs == nullptr) return true;
794  if (lhs == nullptr || rhs == nullptr) return false;
795  return *lhs == *rhs;
796 }
797 
798 template<typename C>
799 bool operator==(const UnivariatePolynomial<C>& lhs, const C& rhs)
800 {
801  if (lhs.coefficients().size() == 0) {
802  return carl::is_zero(rhs);
803  }
804  return (lhs.coefficients().size() == 1) && lhs.lcoeff() == rhs;
805 }
806 
807 template<typename C>
808 bool operator==(const C& lhs, const UnivariatePolynomial<C>& rhs)
809 {
810  return rhs == lhs;
811 }
812 
813 
814 template<typename C>
815 bool operator!=(const UnivariatePolynomial<C>& lhs, const UnivariatePolynomial<C>& rhs)
816 {
817  return !(lhs == rhs);
818 }
819 template<typename C>
820 bool operator!=(const UnivariatePolynomialPtr<C>& lhs, const UnivariatePolynomialPtr<C>& rhs)
821 {
822  if (lhs == nullptr && rhs == nullptr) return false;
823  if (lhs == nullptr || rhs == nullptr) return true;
824  return *lhs != *rhs;
825 }
826 
827 template<typename C>
828 bool UnivariatePolynomial<C>::less(const UnivariatePolynomial<C>& rhs, const PolynomialComparisonOrder& order) const {
829  switch (order) {
830  case PolynomialComparisonOrder::CauchyBound: /*{
831  C a = this->cauchyBound();
832  C b = rhs.cauchyBound();
833  if (a < b) return true;
834  return *this < rhs;
835  }*/
836  case PolynomialComparisonOrder::LowDegree:
837  if (carl::is_zero(rhs)) return false;
838  if (carl::is_zero(*this) || this->degree() < rhs.degree()) return true;
839  return *this < rhs;
840  case PolynomialComparisonOrder::Memory:
841  return this < &rhs;
842  }
843  return false;
844 }
845 
846 template<typename C>
847 bool operator<(const UnivariatePolynomial<C>& lhs, const UnivariatePolynomial<C>& rhs)
848 {
849  if(lhs.mMainVar == rhs.mMainVar)
850  {
851  if(lhs.coefficients().size() == rhs.coefficients().size())
852  {
853  auto iterLhs = lhs.coefficients().rbegin();
854  auto iterRhs = rhs.coefficients().rbegin();
855  while(iterLhs != lhs.coefficients().rend())
856  {
857  assert(iterRhs != rhs.coefficients().rend());
858  if(*iterLhs == *iterRhs)
859  {
860  ++iterLhs;
861  ++iterRhs;
862  }
863  else
864  {
865  return *iterLhs < *iterRhs;
866  }
867  }
868  }
869  return lhs.coefficients().size() < rhs.coefficients().size();
870  }
871  return lhs.mMainVar < rhs.mMainVar;
872 }
873 
874 template<typename C>
875 std::ostream& operator<<(std::ostream& os, const UnivariatePolynomial<C>& rhs)
876 {
877  if(carl::is_zero(rhs)) return os << "0";
878  for(size_t i = 0; i < rhs.mCoefficients.size()-1; ++i )
879  {
880  const C& c = rhs.mCoefficients[rhs.mCoefficients.size()-i-1];
881  if(!(c == 0))
882  {
883  if (!(c == 1)) os << "(" << c << ")*";
884  os << rhs.mMainVar << "^" << rhs.mCoefficients.size()-i-1 << " + ";
885  }
886  }
887  os << rhs.mCoefficients[0];
888  return os;
889 }
890 
891 template<typename Coefficient>
892 template<typename C, EnableIf<is_number_type<C>>>
893 bool UnivariatePolynomial<Coefficient>::is_consistent() const {
894  if (!mCoefficients.empty()) {
895  assert(!carl::is_zero(lcoeff()));
896  }
897  return true;
898 }
899 
900 template<typename Coefficient>
901 template<typename C, DisableIf<is_number_type<C>>>
902 bool UnivariatePolynomial<Coefficient>::is_consistent() const {
903  if (!mCoefficients.empty()) {
904  assert(!carl::is_zero(lcoeff()));
905  }
906  for (const auto& c: mCoefficients) {
907  assert(!c.has(main_var()));
908  assert(c.is_consistent());
909  }
910  return true;
911 }
912 
913 }