carl  24.04
Computer ARithmetic Library
Interval.tpp
Go to the documentation of this file.
1 /**
2  * The implementation for the templated interval class.
3  *
4  * @file Interval.tpp
5  * @author Stefan Schupp <stefan.schupp@cs.rwth-aachen.de>
6  *
7  * @since 2013-12-20
8  * @version 2014-11-11
9  */
10 #pragma once
11 #include "Interval.h"
12 
13 #include <carl-arith/numbers/numbers.h>
14 
15 #include <iostream>
16 
17 namespace carl
18 {
19  namespace bn = boost::numeric;
20 
21 /*******************************************************************************
22  * Transformations and advanced getters/setters
23  ******************************************************************************/
24 
25 template<typename Number>
26 Sign Interval<Number>::sgn() const
27 {
28  assert(this->is_consistent());
29  if (this->is_infinite()) return Sign::ZERO;
30  if ((mLowerBoundType == BoundType::STRICT && mContent.lower() >= carl::constant_zero<Number>().get()) || (mLowerBoundType == BoundType::WEAK && mContent.lower() > carl::constant_zero<Number>().get())) return Sign::POSITIVE;
31  if ((mUpperBoundType == BoundType::STRICT && mContent.upper() <= carl::constant_zero<Number>().get()) || (mUpperBoundType == BoundType::WEAK && mContent.upper() < carl::constant_zero<Number>().get())) return Sign::NEGATIVE;
32  return Sign::ZERO;
33 }
34 
35 template<typename Number>
36 Interval<Number> Interval<Number>::integral_part() const
37 {
38  if(this->is_empty())
39  return *this;
40 
41  Number newLowerBound = 0;
42  Number newUpperBound = 0;
43  BoundType newLowerBoundType = mLowerBoundType;
44  BoundType newUpperBoundType = mUpperBoundType;
45 
46  switch(mLowerBoundType) {
47  case BoundType::WEAK:
48  newLowerBound = ceil(mContent.lower());
49  newLowerBoundType = BoundType::WEAK;
50  break;
51  case BoundType::STRICT:
52  newLowerBound = ceil(mContent.lower());
53  newLowerBoundType = BoundType::WEAK;
54  if(newLowerBound == mContent.lower())
55  newLowerBound += carl::constant_one<Number>::get();
56  break;
57  default:
58  break;
59  }
60  switch(mUpperBoundType) {
61  case BoundType::WEAK:
62  newUpperBound = floor(mContent.upper());
63  newUpperBoundType = BoundType::WEAK;
64  if(newLowerBoundType == BoundType::INFTY)
65  newLowerBound = newUpperBound;
66  break;
67  case BoundType::STRICT:
68  newUpperBound = floor(mContent.upper());
69  newUpperBoundType = BoundType::WEAK;
70  if(newUpperBound == mContent.upper())
71  newUpperBound -= carl::constant_one<Number>::get();
72  if(newLowerBoundType == BoundType::INFTY)
73  newLowerBound = newUpperBound;
74  break;
75  default:
76  if(newLowerBoundType != BoundType::INFTY)
77  newUpperBound = newLowerBound;
78  break;
79  }
80  return Interval<Number>(newLowerBound, newLowerBoundType, newUpperBound, newUpperBoundType);
81 }
82 
83 template<typename Number>
84  void Interval<Number>::integralPart_assign()
85 {
86  *this = integral_part();
87 }
88 
89 template<typename Number>
90 bool Interval<Number>::contains_integer() const
91  {
92  assert(this->is_consistent());
93  switch (mLowerBoundType) {
94  case BoundType::INFTY:
95  return true;
96  case BoundType::STRICT:
97  break;
98  case BoundType::WEAK:
99  if (carl::is_integer(mContent.lower())) return true;
100  }
101  switch (mUpperBoundType) {
102  case BoundType::INFTY:
103  return true;
104  case BoundType::STRICT:
105  break;
106  case BoundType::WEAK:
107  if (carl::is_integer(mContent.upper())) return true;
108  }
109  return carl::floor(mContent.lower()) + carl::constant_one<Number>::get() < mContent.upper();
110  }
111 
112 
113 template<typename Number>
114 Number Interval<Number>::diameter() const
115  {
116  assert(this->is_consistent());
117  return boost::numeric::width(mContent);
118  }
119 
120 template<typename Number>
121 void Interval<Number>::diameter_assign()
122  {
123  this->set(BoostInterval(this->diameter()));
124  }
125 
126 template<typename Number>
127  Number Interval<Number>::diameter_ratio(const Interval<Number>& rhs) const
128  {
129  assert(rhs.diameter() != carl::constant_zero<Number>().get());
130  return this->diameter()/rhs.diameter();
131  }
132 
133 template<typename Number>
134  void Interval<Number>::diameter_ratio_assign(const Interval<Number>& rhs)
135  {
136  this->set(BoostInterval(this->diameter_ratio(rhs)));
137  }
138 
139 template<typename Number>
140  Number Interval<Number>::magnitude() const
141  {
142  assert(this->is_consistent());
143  return boost::numeric::norm(mContent);
144  }
145 
146 template<typename Number>
147  void Interval<Number>::magnitude_assign()
148  {
149  this->set(BoostInterval(this->magnitude()));
150  }
151 
152 template<typename Number>
153  void Interval<Number>::center_assign()
154  {
155  this->set(BoostInterval(carl::center(*this)));
156  }
157 
158  template<typename Number>
159  bool Interval<Number>::contains(const Number& val) const
160  {
161  assert(this->is_consistent());
162  switch( mLowerBoundType )
163  {
164  case BoundType::INFTY:
165  break;
166  case BoundType::STRICT:
167  if( mContent.lower() >= val )
168  return false;
169  break;
170  case BoundType::WEAK:
171  if( mContent.lower() > val )
172  return false;
173  }
174  // Invariant: n is not conflicting with lower bound
175  switch( mUpperBoundType )
176  {
177  case BoundType::INFTY:
178  break;
179  case BoundType::STRICT:
180  if( mContent.upper() <= val )
181  return false;
182  break;
183  case BoundType::WEAK:
184  if( mContent.upper() < val )
185  return false;
186  break;
187  }
188  return true; // for open intervals: (lower() < n && upper() > n) || (n == carl::constant_zero<Number>().get() && lower() == cln::cl_RA( 0 ) && upper() == cln::cl_RA( 0 )) }
189  }
190 
191  template<typename Number>
192  bool Interval<Number>::contains(const Interval<Number>& rhs) const
193  {
194  assert(this->is_consistent());
195  if( rhs.is_empty() )
196  return true;
197  // if one bound is totally wrong, we can just return false
198  if((mContent.lower() > rhs.lower() && mLowerBoundType != BoundType::INFTY) || (mContent.upper() < rhs.upper() && mUpperBoundType != BoundType::INFTY))
199  {
200  return false;
201  }
202  // check the bounds
203  bool lowerOk = mContent.lower() < rhs.lower() && rhs.lower_bound_type() != BoundType::INFTY;
204  bool upperOk = mContent.upper() > rhs.upper() && rhs.upper_bound_type() != BoundType::INFTY;
205  // if both are ok, return true
206  if( lowerOk && upperOk )
207  {
208  return true;
209  }
210  // Note that from this point on at least one bound is equal
211  // to our bounds but no bound is outside of our bounds
212  // check the bound types
213  bool lowerBoundTypesOk = get_weakest_bound_type(mLowerBoundType, rhs.lower_bound_type()) == mLowerBoundType;
214  bool upperBoundTypesOk = get_weakest_bound_type(mUpperBoundType, rhs.upper_bound_type()) == mUpperBoundType;
215  // if upper bounds are ok and lower bound types are ok, return true
216  if (upperOk && lowerBoundTypesOk)
217  {
218  return true;
219  }
220  // if lower bounds are ok and upper bound types are ok, return true
221  if (lowerOk && upperBoundTypesOk)
222  {
223  return true;
224  }
225  // if both bound types are ok, return true
226  if (lowerBoundTypesOk && upperBoundTypesOk)
227  {
228  return true;
229  }
230  // otherwise return false
231  return false; // not less and not equal
232  }
233 
234  template<typename Number>
235  bool Interval<Number>::meets(const Number& n) const
236  {
237  assert(this->is_consistent());
238  return (mContent.lower() <= n || mLowerBoundType == BoundType::INFTY) && (mContent.upper() >= n || mUpperBoundType == BoundType::INFTY);
239  }
240 
241  template<typename Number>
242  void Interval<Number>::bloat_by(const Number& width) {
243  if(this->is_empty()) {
244  return;
245  }
246  if (width == carl::constant_zero<Number>().get()) {
247  mLowerBoundType = BoundType::WEAK;
248  mUpperBoundType = BoundType::WEAK;
249  mContent.assign(0, 0);
250  return;
251  }
252  if (!is_infinite()) {
253  BoundType lowerTmp = mLowerBoundType;
254  BoundType upperTmp = mUpperBoundType;
255  this->set(boost::numeric::widen(mContent, width));
256  mLowerBoundType = lowerTmp;
257  mUpperBoundType = upperTmp;
258  } else if (mLowerBoundType != BoundType::INFTY) {
259  this->set(boost::numeric::widen(mContent, width));
260  mUpperBoundType = BoundType::INFTY;
261  } else if (mUpperBoundType != BoundType::INFTY) {
262  this->set(boost::numeric::widen(mContent, width));
263  mLowerBoundType = BoundType::INFTY;
264  }
265  }
266 
267  template<typename Number>
268  void Interval<Number>::shrink_by(const Number& width)
269  {
270  this->bloat_by(Number(-1)*width);
271  }
272 
273  template<typename Number>
274  std::pair<Interval<Number>, Interval<Number>> Interval<Number>::split() const
275  {
276  std::pair<BoostInterval, BoostInterval> bisection = boost::numeric::bisect(mContent);
277  if( this->is_empty() || this->is_point_interval() )
278  {
279  return std::pair<Interval<Number>, Interval<Number> >(Interval<Number>::empty_interval(), Interval<Number>::empty_interval());
280  }
281  return std::pair<Interval<Number>, Interval<Number> >(Interval(bisection.first, mLowerBoundType, BoundType::STRICT), Interval(bisection.second, BoundType::WEAK, mUpperBoundType));
282  }
283 
284  template<typename Number>
285  std::list<Interval<Number>> Interval<Number>::split(unsigned n) const
286  {
287  std::list<Interval<Number> > result;
288  if(n == 0)
289  {
290  assert(false);
291  result.emplace_back(empty_interval());
292  return result;
293  } else if ( n == 1) {
294  result.push_back(*this);
295  return result;
296  }
297  Number diameter = this->diameter();
298  diameter /= Number(n);
299 
300  Interval<Number> tmp;
301  tmp.set(mContent.lower(), mContent.lower()+diameter);
302  tmp.set_lower_bound_type(mLowerBoundType);
303  tmp.set_upper_bound_type(BoundType::STRICT);
304  result.push_back(tmp);
305 
306  tmp.set_lower_bound_type(BoundType::WEAK);
307  for( unsigned i = 1; i < (n-1); ++i )
308  {
309  tmp += diameter;
310  result.push_back(tmp);
311  }
312 
313  tmp += diameter;
314  tmp.set_upper_bound_type(mUpperBoundType);
315  result.push_back(tmp);
316  return result;
317  }
318 
319  template<typename Number>
320  std::string Interval<Number>::toString() const
321  {
322  std::ostringstream oss;
323  switch (mLowerBoundType) {
324  case BoundType::INFTY:
325  oss << std::string("]-INF, ");
326  break;
327  case BoundType::STRICT:
328  oss << std::string("]") << mContent.lower() << ", ";
329  break;
330  case BoundType::WEAK:
331  oss << std::string("[") << mContent.lower() << ", ";
332  }
333  switch (mUpperBoundType) {
334  case BoundType::INFTY:
335  oss << std::string("INF[");
336  break;
337  case BoundType::STRICT:
338  oss << mContent.upper() << std::string("[");
339  break;
340  case BoundType::WEAK:
341  oss << mContent.upper() << std::string("]");
342  }
343  return oss.str();
344  }
345 
346 
347 /**
348 * Calculates the distance between two Intervals.
349 * @param intervalA Interval to which we want to know the distance.
350 * @return distance to intervalA
351 */
352 template<typename Number>
353 Number Interval<Number>::distance(const Interval<Number>& intervalA)
354 {
355  if (intervalA.upper_bound() < lower_bound()) {
356  return lower() - intervalA.upper();
357  }
358  if (upper_bound() < intervalA.lower_bound()) {
359  return intervalA.lower() - upper();
360  }
361  return carl::constant_zero<Number>::get();
362 }
363 
364 template<typename Number>
365 Interval<Number> Interval<Number>::convex_hull(const Interval<Number>& interval) const {
366  if(this->is_empty())
367  return interval;
368 
369  if(interval.is_empty())
370  return *this;
371 
372  BoundType newLowerBound = get_strictest_bound_type(this->lower_bound_type(), interval.lower_bound_type());
373  BoundType newUpperBound = get_strictest_bound_type(this->upper_bound_type(), interval.upper_bound_type());
374  Number newLower = interval.lower() < this->lower() ? interval.lower() : this->lower();
375  Number newUpper = interval.upper() > this->upper() ? interval.upper() : this->upper();
376 
377  return Interval(newLower, newLowerBound, newUpper, newUpperBound);
378 }
379 
380 /*******************************************************************************
381  * Arithmetic functions
382  ******************************************************************************/
383 
384 template<typename Number>
385 Interval<Number> Interval<Number>::add(const Interval<Number>& rhs) const
386 {
387  assert(this->is_consistent());
388  assert(rhs.is_consistent());
389  if (rhs.is_empty()) return empty_interval();
390  return Interval<Number>( mContent + rhs.content(),
391  get_strictest_bound_type( mLowerBoundType, rhs.lower_bound_type() ),
392  get_strictest_bound_type( mUpperBoundType, rhs.upper_bound_type() ) );
393 }
394 
395 template<typename Number>
396 void Interval<Number>::add_assign(const Interval<Number>& rhs)
397 {
398  *this = this->add(rhs);
399 }
400 
401 template<typename Number>
402 Interval<Number> Interval<Number>::sub(const Interval<Number>& rhs) const
403 {
404  assert(this->is_consistent());
405  assert(rhs.is_consistent());
406  return this->add(rhs.inverse());
407 }
408 
409 template<typename Number>
410 void Interval<Number>::sub_assign(const Interval<Number>& rhs)
411 {
412  *this = this->sub(rhs);
413 }
414 
415 template<typename Number>
416 Interval<Number> Interval<Number>::mul(const Interval<Number>& rhs) const
417 {
418  assert(this->is_consistent());
419  assert(rhs.is_consistent());
420  if( this->is_empty() || rhs.is_empty() ) return empty_interval();
421  BoundType lowerBoundType = BoundType::WEAK;
422  BoundType upperBoundType = BoundType::WEAK;
423  BoostInterval resultInterval;
424  // The following unfortunately copies the content of boost::interval::mul, but we need to get into the case distinction to find the right bound types
425  typename BoostIntervalPolicies::rounding rnd;
426  const Number& xl = this->lower();
427  const Number& xu = this->upper();
428  const Number& yl = rhs.lower();
429  const Number& yu = rhs.upper();
430  const BoundType& xlt = this->lower_bound_type();
431  const BoundType& xut = this->upper_bound_type();
432  const BoundType& ylt = rhs.lower_bound_type();
433  const BoundType& yut = rhs.upper_bound_type();
434 
435  if (xl < carl::constant_zero<Number>().get())
436  {
437  if (xu > carl::constant_zero<Number>().get())
438  {
439  if (yl < carl::constant_zero<Number>().get())
440  {
441  if (yu > carl::constant_zero<Number>().get()) // M * M
442  {
443  Number lowerA = rnd.mul_down(xl, yu);
444  Number lowerB = rnd.mul_down(xu, yl);
445  if( lowerA > lowerB )
446  {
447  lowerA = lowerB;
448  lowerBoundType = get_strictest_bound_type(xut, ylt);
449  }
450  else
451  {
452  lowerBoundType = get_strictest_bound_type(xlt, yut);
453  }
454  Number upperA = rnd.mul_up(xl, yl);
455  Number upperB = rnd.mul_up(xu, yu);
456  if( upperA < upperB )
457  {
458  upperA = upperB;
459  upperBoundType = get_strictest_bound_type(xut, yut);
460  }
461  else
462  {
463  upperBoundType = get_strictest_bound_type(xlt, ylt);
464  }
465  resultInterval = BoostInterval(lowerA, upperA, true);
466  }
467  else // M * N
468  {
469  lowerBoundType = get_strictest_bound_type(xut, ylt);
470  upperBoundType = get_strictest_bound_type(xlt, ylt);
471  resultInterval = BoostInterval(rnd.mul_down(xu, yl), rnd.mul_up(xl, yl), true);
472  }
473  }
474  else
475  {
476  if (yu > carl::constant_zero<Number>().get()) // M * P
477  {
478  lowerBoundType = get_strictest_bound_type(xlt, yut);
479  upperBoundType = get_strictest_bound_type(xut, yut);
480  resultInterval = BoostInterval(rnd.mul_down(xl, yu), rnd.mul_up(xu, yu), true);
481  }
482  else // M * Z
483  {
484  if( ylt != BoundType::INFTY )
485  lowerBoundType = ylt;
486  if( yut != BoundType::INFTY )
487  upperBoundType = yut;
488  resultInterval = BoostInterval(static_cast<Number>(0), static_cast<Number>(0), true);
489  }
490  }
491  }
492  else
493  {
494  if (yl < carl::constant_zero<Number>().get())
495  {
496  if (yu > carl::constant_zero<Number>().get()) // N * M
497  {
498  lowerBoundType = get_strictest_bound_type(xlt, yut);
499  upperBoundType = get_strictest_bound_type(xlt, ylt);
500  resultInterval = BoostInterval(rnd.mul_down(xl, yu), rnd.mul_up(xl, yl), true);
501  }
502  else // N * N
503  {
504  lowerBoundType = get_strictest_bound_type(xut, yut);
505  upperBoundType = get_strictest_bound_type(xlt, ylt);
506  resultInterval = BoostInterval(rnd.mul_down(xu, yu), rnd.mul_up(xl, yl), true);
507  }
508  }
509  else
510  {
511  if (yu > carl::constant_zero<Number>().get()) // N * P
512  {
513  lowerBoundType = get_strictest_bound_type(xlt, yut);
514  upperBoundType = get_strictest_bound_type(xut, ylt);
515  resultInterval = BoostInterval(rnd.mul_down(xl, yu), rnd.mul_up(xu, yl), true);
516  }
517  else // N * Z
518  {
519  lowerBoundType = get_strictest_bound_type(xut, yut);
520  upperBoundType = get_strictest_bound_type(xut, ylt);
521  resultInterval = BoostInterval(static_cast<Number>(0), static_cast<Number>(0), true);
522  }
523  }
524  }
525  }
526  else
527  {
528  if (xu > carl::constant_zero<Number>().get())
529  {
530  if (yl < carl::constant_zero<Number>().get())
531  {
532  if (yu > carl::constant_zero<Number>().get()) // P * M
533  {
534  lowerBoundType = get_strictest_bound_type(xut, ylt);
535  upperBoundType = get_strictest_bound_type(xut, yut);
536  resultInterval = BoostInterval(rnd.mul_down(xu, yl), rnd.mul_up(xu, yu), true);
537  }
538  else // P * N
539  {
540  lowerBoundType = get_strictest_bound_type(xut, ylt);
541  upperBoundType = get_strictest_bound_type(xlt, yut);
542  resultInterval = BoostInterval(rnd.mul_down(xu, yl), rnd.mul_up(xl, yu), true);
543  }
544  }
545  else
546  {
547  if (yu > carl::constant_zero<Number>().get()) // P * P
548  {
549  lowerBoundType = get_strictest_bound_type(xlt, ylt);
550  upperBoundType = get_strictest_bound_type(xut, yut);
551  resultInterval = BoostInterval(rnd.mul_down(xl, yl), rnd.mul_up(xu, yu), true);
552  }
553  else // P * Z
554  {
555  if( ylt != BoundType::INFTY )
556  lowerBoundType = ylt;
557  if( yut != BoundType::INFTY )
558  upperBoundType = yut;
559  resultInterval = BoostInterval(static_cast<Number>(0), static_cast<Number>(0), true);
560  }
561  }
562  }
563  else // Z * ?
564  {
565  if( xlt != BoundType::INFTY )
566  lowerBoundType = xlt;
567  if( xut != BoundType::INFTY )
568  upperBoundType = xut;
569  resultInterval = BoostInterval(static_cast<Number>(0), static_cast<Number>(0), true);
570  }
571  }
572  unsigned zeroBoundInvolved = 2;
573  if( (xlt == BoundType::INFTY && (yu > carl::constant_zero<Number>().get() || yut == BoundType::INFTY))
574  || (xut == BoundType::INFTY && (yl < carl::constant_zero<Number>().get() || ylt == BoundType::INFTY))
575  || (ylt == BoundType::INFTY && (xu > carl::constant_zero<Number>().get() || xut == BoundType::INFTY))
576  || (yut == BoundType::INFTY && (xu < carl::constant_zero<Number>().get() || (xl < carl::constant_zero<Number>().get() || xlt == BoundType::INFTY))) )
577  {
578  lowerBoundType = BoundType::INFTY;
579  }
580  else if( resultInterval.lower() == carl::constant_zero<Number>().get() )
581  {
582  if( zeroBoundInvolved == 2 )
583  zeroBoundInvolved = (this->contains( carl::constant_zero<Number>().get() ) || rhs.contains( carl::constant_zero<Number>().get() )) ? 1 : 0;
584  if( zeroBoundInvolved == 1 )
585  lowerBoundType = BoundType::WEAK;
586  }
587  if( (xlt == BoundType::INFTY && (yu < carl::constant_zero<Number>().get() || (yl < carl::constant_zero<Number>().get() || ylt == BoundType::INFTY)))
588  || (xut == BoundType::INFTY && (yl > carl::constant_zero<Number>().get() || (yu > carl::constant_zero<Number>().get() || yut == BoundType::INFTY)))
589  || (ylt == BoundType::INFTY && (xu < carl::constant_zero<Number>().get() || (xl < carl::constant_zero<Number>().get() || xlt == BoundType::INFTY)))
590  || (yut == BoundType::INFTY && (xl > carl::constant_zero<Number>().get() || (xu > carl::constant_zero<Number>().get() || xut == BoundType::INFTY))) )
591  {
592  upperBoundType = BoundType::INFTY;
593  }
594  else if( resultInterval.upper() == carl::constant_zero<Number>().get() )
595  {
596  if( zeroBoundInvolved == 2 )
597  zeroBoundInvolved = (this->contains( carl::constant_zero<Number>().get() ) || rhs.contains( carl::constant_zero<Number>().get() )) ? 1 : 0;
598  if( zeroBoundInvolved == 1 )
599  upperBoundType = BoundType::WEAK;
600  }
601  return Interval<Number>(std::move(resultInterval), lowerBoundType, upperBoundType );
602 }
603 
604 template<typename Number>
605 void Interval<Number>::mul_assign(const Interval<Number>& rhs)
606  {
607  *this = this->mul(rhs);
608  }
609 
610 template<typename Number>
611 Interval<Number> Interval<Number>::div(const Interval<Number>& rhs) const
612  {
613  assert(this->is_consistent());
614  assert(rhs.is_consistent());
615  assert(!rhs.contains(carl::constant_zero<Number>().get()));
616  BoundType lowerBoundType = BoundType::WEAK;
617  BoundType upperBoundType = BoundType::WEAK;
618  ///@todo Correctly determine if bounds are strict or weak.
619  if (this->is_open_interval() || rhs.is_open_interval()) {
620  // just a quick heuristic, by no means complete.
621  lowerBoundType = BoundType::STRICT;
622  upperBoundType = BoundType::STRICT;
623  }
624  const Number& xl = mContent.lower();
625  const Number& yl = rhs.lower();
626  const Number& yu = rhs.upper();
627  const BoundType& xlt = mLowerBoundType;
628  const BoundType& xut = mUpperBoundType;
629  const BoundType& ylt = rhs.lower_bound_type();
630  const BoundType& yut = rhs.upper_bound_type();
631  if( (xlt == BoundType::INFTY && (carl::is_positive(yu) || yut == BoundType::INFTY))
632  || (xut == BoundType::INFTY && (carl::is_negative(yl) || ylt == BoundType::INFTY))
633  || (ylt == BoundType::INFTY && (carl::is_positive(xl) || xut == BoundType::INFTY))
634  || (yut == BoundType::INFTY && (carl::is_negative(xl) || (carl::is_negative(xl) || xlt == BoundType::INFTY))) )
635  {
636  lowerBoundType = BoundType::INFTY;
637  }
638  if( (xlt == BoundType::INFTY && (carl::is_negative(yu) || (carl::is_negative(yl) || ylt == BoundType::INFTY)))
639  || (xut == BoundType::INFTY && (carl::is_positive(yl) || (carl::is_positive(yu) || yut == BoundType::INFTY)))
640  || (ylt == BoundType::INFTY && (carl::is_negative(xl) || (carl::is_negative(xl) || xlt == BoundType::INFTY)))
641  || (yut == BoundType::INFTY && (carl::is_positive(xl) || (carl::is_positive(xl) || xut == BoundType::INFTY))) )
642  {
643  upperBoundType = BoundType::INFTY;
644  }
645  return Interval<Number>(BoostInterval( mContent/rhs.content() ), lowerBoundType, upperBoundType );
646  }
647 
648 template<typename Number>
649 void Interval<Number>::div_assign(const Interval<Number>& rhs)
650  {
651  *this = this->div(rhs);
652  }
653 
654 template<typename Number>
655 bool Interval<Number>::div_ext(const Interval<Number>& rhs, Interval<Number>& a, Interval<Number>& b) const
656  {
657  // Special case: if both contain 0 we can directly skip and return the unbounded interval.
658  if(this->contains(carl::constant_zero<Number>().get()) && rhs.contains(carl::constant_zero<Number>().get()))
659  {
660  a = unbounded_interval();
661  return false;
662  }
663 
664  Interval<Number> reciprocalA, reciprocalB;
665  bool splitOccured;
666 
667  if( rhs.lower_bound_type() != BoundType::INFTY && rhs.lower() == carl::constant_zero<Number>().get() && rhs.upper_bound_type() != BoundType::INFTY && rhs.upper() == carl::constant_zero<Number>().get() ) // point interval 0
668  {
669  splitOccured = false;
670  if( this->contains( carl::constant_zero<Number>().get() ))
671  {
672  a = unbounded_interval();
673  }
674  else
675  {
676  a = empty_interval();
677  }
678  return false;
679  }
680  else
681  {
682  if( rhs.is_infinite() )
683  {
684  a = unbounded_interval();
685  return false;
686  } // rhs.unbounded
687  else
688  {
689  //default case
690  splitOccured = rhs.reciprocal( reciprocalA, reciprocalB );
691  if( !splitOccured )
692  {
693  a = this->mul( reciprocalA );
694  return false;
695  }
696  else
697  {
698  a = this->mul( reciprocalA );
699  b = this->mul( reciprocalB );
700 
701  if( a == b )
702  {
703  return false;
704  }
705  else
706  {
707  return true;
708  }
709 
710  }
711  } // !rhs.unbounded
712 
713  } // !pointinterval 0
714  }
715 
716 template<typename Number>
717 Interval<Number> Interval<Number>::inverse() const
718  {
719  assert(this->is_consistent());
720  return Interval<Number>( mContent.upper()*Number(-1), mUpperBoundType, mContent.lower()*Number(-1), mLowerBoundType );
721  }
722 
723 template<typename Number>
724 void Interval<Number>::inverse_assign()
725  {
726  *this = this->inverse();
727  }
728 
729  template<typename Number>
730  Interval<Number> Interval<Number>::abs() const
731  {
732  if(this->contains(carl::constant_zero<Number>().get()))
733  {
734  Number max = mContent.upper() > mContent.lower() ? mContent.upper() : mContent.lower();
735  BoundType bt = mContent.upper() > mContent.lower() ? mUpperBoundType : mLowerBoundType;
736  return Interval<Number>( carl::constant_zero<Number>().get(), BoundType::WEAK, max, bt );
737  }
738  else if( mContent.upper() < carl::constant_zero<Number>().get()) // interval is fully negative
739  {
740  return Interval<Number>(-mContent.upper(), mUpperBoundType, -mContent.lower(), mLowerBoundType);
741  }
742  // otherwise inteval is already fully positive
743  return *this;
744  }
745 
746  template<typename Number>
747  void Interval<Number>::abs_assign()
748  {
749  *this = this->abs();
750  }
751 
752 template<typename Number>
753 bool Interval<Number>::reciprocal(Interval<Number>& a, Interval<Number>& b) const {
754  if( this->is_infinite() )
755  {
756  a = empty_interval();
757  return false;
758  }
759  else if( this->contains( carl::constant_zero<Number>().get() ) && mContent.lower() != carl::constant_zero<Number>().get() && mContent.upper() != carl::constant_zero<Number>().get() )
760  {
761  if( mLowerBoundType == BoundType::INFTY )
762  {
763  a = Interval<Number>(carl::constant_zero<Number>().get(), BoundType::INFTY,carl::constant_zero<Number>().get(), BoundType::WEAK );
764  b = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) /BoostInterval( mContent.upper() ), BoundType::WEAK, BoundType::INFTY );
765  }
766  else if( mUpperBoundType == BoundType::INFTY )
767  {
768  a = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) /BoostInterval( mContent.lower() ), BoundType::INFTY, BoundType::WEAK );
769  b = Interval<Number>(carl::constant_zero<Number>().get(), BoundType::WEAK,carl::constant_zero<Number>().get(), BoundType::INFTY );
770  }
771  else if( mContent.lower() == carl::constant_zero<Number>().get() && mContent.upper() != carl::constant_zero<Number>().get() )
772  {
773  a = Interval<Number>(carl::constant_zero<Number>().get(), BoundType::INFTY, carl::constant_zero<Number>().get(), BoundType::INFTY );
774  b = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) /BoostInterval( mContent.upper() ), BoundType::WEAK, BoundType::INFTY );
775  }
776  else if( mContent.lower() != carl::constant_zero<Number>().get() && mContent.upper() == carl::constant_zero<Number>().get() )
777  {
778  a = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) /BoostInterval( mContent.lower() ), BoundType::INFTY, BoundType::WEAK );
779  b = unbounded_interval(); // todo: really the whole interval here?
780  }
781  else if( mContent.lower() == carl::constant_zero<Number>().get() && mContent.upper() == carl::constant_zero<Number>().get() )
782  {
783  a = unbounded_interval();
784  return false;
785  }
786  else
787  {
788  a = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) /BoostInterval( mContent.lower() ), BoundType::INFTY, BoundType::WEAK );
789  b = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) /BoostInterval( mContent.upper() ), BoundType::WEAK, BoundType::INFTY );
790  }
791  return true;
792  }
793  else
794  {
795  if( mLowerBoundType == BoundType::INFTY && mContent.upper() != carl::constant_zero<Number>().get() )
796  {
797  a = Interval<Number>( carl::constant_one<Number>().get() / mContent.upper() , mUpperBoundType,carl::constant_zero<Number>().get(), BoundType::WEAK );
798  }
799  else if( mLowerBoundType == BoundType::INFTY && mContent.upper() == carl::constant_zero<Number>().get() )
800  {
801  a = Interval<Number>(carl::constant_zero<Number>().get(), BoundType::INFTY, carl::constant_zero<Number>().get(), BoundType::WEAK );
802  }
803  else if( mUpperBoundType == BoundType::INFTY && mContent.lower() != carl::constant_zero<Number>().get() )
804  {
805  a = Interval<Number>( carl::constant_zero<Number>().get() , BoundType::WEAK, carl::constant_one<Number>().get() / mContent.lower(), mLowerBoundType );
806  }
807  else if( mUpperBoundType == BoundType::INFTY && mContent.lower() == carl::constant_zero<Number>().get() )
808  {
809  a = Interval<Number>(carl::constant_zero<Number>().get(), BoundType::WEAK,carl::constant_zero<Number>().get(), BoundType::INFTY );
810  }
811  else if( mContent.lower() != carl::constant_zero<Number>().get() && mContent.upper() != carl::constant_zero<Number>().get() )
812  {
813  a = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) / mContent, mUpperBoundType, mLowerBoundType );
814  }
815  else if( mContent.lower() == carl::constant_zero<Number>().get() && mContent.upper() != carl::constant_zero<Number>().get() )
816  {
817  a = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) /BoostInterval( mContent.upper() ), mUpperBoundType, BoundType::INFTY );
818  }
819  else if( mContent.lower() != carl::constant_zero<Number>().get() && mContent.upper() == carl::constant_zero<Number>().get() )
820  {
821  a = Interval<Number>(BoostInterval( carl::constant_one<Number>().get() ) /BoostInterval( mContent.lower() ), BoundType::INFTY, mLowerBoundType );
822  }
823 
824  return false;
825  }
826  }
827 
828 template<typename Number>
829 template<typename Num, EnableIf<std::is_floating_point<Num>>>
830 Interval<Number> Interval<Number>::root(int deg) const
831  {
832  assert(this->is_consistent());
833  if( deg % 2 == 0 )
834  {
835  if( mUpperBoundType != BoundType::INFTY && mContent.upper() < carl::constant_zero<Number>().get() )
836  return Interval<Number>::empty_interval();
837  if( mLowerBoundType == BoundType::INFTY || mContent.lower() < carl::constant_zero<Number>().get() )
838  {
839  if( mUpperBoundType == BoundType::INFTY )
840  {
841  return Interval<Number>(BoostInterval(carl::constant_zero<Number>().get()), BoundType::WEAK, mUpperBoundType);
842  }
843  else
844  {
845  return Interval<Number>(boost::numeric::nth_root(BoostInterval(carl::constant_zero<Number>().get(),mContent.upper()), deg), BoundType::WEAK, mUpperBoundType);
846  }
847  }
848  }
849  return Interval<Number>(boost::numeric::nth_root(mContent, deg), mLowerBoundType, mUpperBoundType);
850  }
851 
852 template<typename Number>
853 template<typename Num, EnableIf<std::is_floating_point<Num>>>
854 void Interval<Number>::root_assign(unsigned deg)
855  {
856  *this = this->root(deg);
857  }
858 
859 }