#ifndef __INTERVAL__
#define __INTERVAL__

// Depend on OS we will use \verb|_control87| (\verb|float.h|, Windows) 
// or \verb|fesetround| (\verb|fenv.h|, Linux)
#ifdef _WIN32
  #include <float.h>
#else
  #include <fenv.h>
#endif

// Only for \verb|std::min| and \verb|std::max|
#include <algorithm>
// Output into standard stream
#include <iostream>

namespace MATH
{
// Set rounding direction. To nearest number by default
#ifdef _WIN32
  inline
    void __set_round_down() { _control87(_RC_DOWN, _MCW_RC); }
  inline 
    void __set_round_up() { _control87(_RC_UP, _MCW_RC); }
  inline 
    void __set_round_default() { _control87(_CW_DEFAULT, _MCW_RC); }
#else
  inline 
    void __set_round_down() { fesetround(FE_DOWNWARD); }
  inline 
    void __set_round_up() { fesetround(FE_UPWARD); }
  inline 
    void __set_round_default() { fesetround(FE_TONEAREST); }
#endif

  template <class T = double>
  class interval_t
  {
  public:
    // by default "--- exact zero: $[0,0]$
    interval_t<T>() : _lb(0), _ub(0) {}
    // exact number $[x,x]$
    interval_t<T>(long number) : _lb(number), _ub(number) {}
    // number $[\underline{x},\overline{x}]$
    interval_t<T>(long numer, long denom);
    // interval $[\underline{x},\overline{y}]$
    interval_t<T>(long numer_lb, long denom_lb,
                  long numer_ub, long denom_ub);

    T lo() const { return _lb; } // low bound
    T hi() const { return _ub; } // high bound

    // unary minus (negative number)
    interval_t<T> operator-()
    {
      interval_t<T> neg;
      neg._lb = -_ub;
      neg._ub = -_lb;
      return neg;
    }

    interval_t<T>& operator+=(const interval_t<T>& x);
    interval_t<T>& operator-=(const interval_t<T>& x);
    interval_t<T>& operator*=(const interval_t<T>& x);
    interval_t<T>& operator/=(const interval_t<T>& x);

    interval_t<T> operator + (const interval_t<T>& x)
    {
      return interval_t<T>(*this) += x;
    }
    interval_t<T> operator - (const interval_t<T>& x)
    {
      return interval_t<T>(*this) -= x;
    }
    interval_t<T> operator * (const interval_t<T>& x)
    {
      return interval_t<T>(*this) *= x;
    }
    interval_t<T> operator / (const interval_t<T>& x)
    {
      return interval_t<T>(*this) /= x;
    }

  private:
    T _lb, _ub;
  };

  template <class T>
  interval_t<T>::interval_t(long numer, long denom)
  {
    __set_round_down();
    _lb = T(numer) / T(denom);
    __set_round_up();
    _ub = T(numer) / T(denom);
    __set_round_default();
  }
  template <class T>
  interval_t<T>::interval_t(long numer_lb, long denom_lb,
                            long numer_ub, long denom_ub)
  {
    __set_round_down();
    _lb = T(numer_lb) / T(denom_lb);
    __set_round_up();
    _ub = T(numer_ub) / T(denom_ub);
    __set_round_default();
  }

  template <class T>
  interval_t<T>& interval_t<T>::operator += (const interval_t<T>& x)
  {
    __set_round_down();
    _lb += x._lb;
    __set_round_up();
    _ub += x._ub;
    __set_round_default();

    return *this;
  }
  template <class T>
  interval_t<T>& interval_t<T>::operator -= (const interval_t<T>& x)
  {
    if (this == &x)
    {
      _lb = _ub = 0.0;
    }
    else
    {
      __set_round_down();
      long double l = _lb - x._ub;
      __set_round_up();
      long double u = _ub - x._lb;
      __set_round_default();

      _lb = l;
      _ub = u;
    }

    return *this;
  }
  template <class T>
  interval_t<T>& interval_t<T>::operator *= (const interval_t<T>& x)
  {
    if (this == &x) // $x$ is same number
    {
      if (_lb > T(0) || _ub < T(0))  // \verb|_lb| and \verb|_ub| have same signs
      {
        __set_round_down();
        T l = std::min(_lb*_lb, _ub*_ub);
        __set_round_up();
        T u = std::max(_lb*_lb, _ub*_ub);
        __set_round_default();
        _lb = l;
        _ub = u;
      }
      else // \verb|_lb|$\,\leqslant0\leqslant\,$\verb|_ub|
      {
        _lb = T(0);
        __set_round_up();
        _ub = std::max(_lb*_lb, _ub*_ub);
        __set_round_default();
      }
    }
    else // $x$ is different number
    {
      __set_round_down();
      T l = std::min(std::min(_lb*x._lb, _lb*x._ub),
                               std::min(_ub*x._lb, _ub*x._ub));
      __set_round_up();
      T u = std::max(std::max(_lb*x._lb, _lb*x._ub),
                               std::max(_ub*x._lb, _ub*x._ub));
      __set_round_default();

      _lb = l;
      _ub = u;
    }

    return *this;
  }
  template <class T>
  interval_t<T>& interval_t<T>::operator /= (const interval_t<T>& x)
  {
    if ((x._lb <= 0) && (0 <= x._ub)) throw "Division by zero";

    if (this == &x)
    {
      _lb = _ub = 1.0;
      return *this;
    }
    else
    {
      interval_t<T> x_inv;

      __set_round_down();
      x_inv._lb = T(1) / x._ub;
      __set_round_up();
      x_inv._ub = T(1) / x._lb;
      __set_round_default();

      return (*this) *= x_inv;
    }
  }

  template <class T>
  std::ostream& operator<< (std::ostream& os, interval_t<T>& x)
  {
    if (x.lo() == x.hi())
      os << x.lo();
    else
      os << "[" << x.lo() << ";" << x.hi() << "]";
    return os;
  }
};

#endif