#ifndef INTERVAL_HPP
#define INTERVAL_HPP

#include <iostream>
#include <cmath>
#include <stdexcept>
#include <fenv.h>

class interval {

	double inf;
	double sup;

	public:

	interval() {
		inf = 0;
		sup = 0;
	}

	interval(const double& x) {
		inf = x;
		sup = x;
	}

	interval(const double& x, const double& y) {
		inf = x;
		sup = y;
	}

	friend interval operator+(const interval& x, const interval& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x.inf + y.inf;
		fesetround(FE_UPWARD);
		r.sup = x.sup + y.sup;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator+(const interval& x, const double& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x.inf + y;
		fesetround(FE_UPWARD);
		r.sup = x.sup + y;
		fesetround(FE_TONEAREST);
		return r;
	}

	friend interval operator+(const double& x, const interval& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x + y.inf;
		fesetround(FE_UPWARD);
		r.sup = x + y.sup;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval& operator+=(interval& x, const interval& y) {
		x = x + y;
		return x;
	}

	friend interval& operator+=(interval& x, const double& y) {

		fesetround(FE_DOWNWARD);
		x.inf = x.inf + y;
		fesetround(FE_UPWARD);
		x.sup = x.sup + y;
		fesetround(FE_TONEAREST);

		return x;
	}

	friend interval operator-(const interval& x, const interval& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x.inf - y.sup;
		fesetround(FE_UPWARD);
		r.sup = x.sup - y.inf;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator-(const interval& x, const double& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x.inf - y;
		fesetround(FE_UPWARD);
		r.sup = x.sup - y;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator-(const double& x, const interval& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x - y.sup;
		fesetround(FE_UPWARD);
		r.sup = x - y.inf;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval& operator-=(interval& x, const interval& y) {
		x = x - y;
		return x;
	}

	friend interval& operator-=(interval& x, const double& y) {

		fesetround(FE_DOWNWARD);
		x.inf = x.inf - y;
		fesetround(FE_UPWARD);
		x.sup = x.sup - y;
		fesetround(FE_TONEAREST);

		return x;
	}

	friend interval operator-(const interval& x) {
		interval r;

		r.sup = - x.inf;
		r.inf = - x.sup;

		return r;
	}

	friend interval operator*(const interval& x, const interval& y) {
		interval r;
		double tmp;

		if (x.inf >= 0) {
			if (y.inf >= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.sup * y.sup;
			} else if (y.sup <= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup * y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.sup;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup * y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.sup * y.sup;
			}
		} else if (x.sup <= 0) {
			if (y.inf >= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.sup * y.inf;
			} else if (y.sup <= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup * y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.inf;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.inf;
			}
		} else {
			if (y.inf >= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.sup * y.sup;
			} else if (y.sup <= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup * y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.inf;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.sup;
				tmp = x.sup * y.inf;
				if (tmp < r.inf) r.inf = tmp;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.inf;
				tmp = x.sup * y.sup;
				if (tmp > r.sup) r.sup = tmp;
			}
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator*(const interval& x, const double& y) {
		interval r;

		if (y >= 0) {
			fesetround(FE_DOWNWARD);
			r.inf = x.inf * y;
			fesetround(FE_UPWARD);
			r.sup = x.sup * y;
		} else {
			fesetround(FE_DOWNWARD);
			r.inf = x.sup * y;
			fesetround(FE_UPWARD);
			r.sup = x.inf * y;
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator*(const double& x, const interval& y) {
		interval r;

		if (x >= 0) {
			fesetround(FE_DOWNWARD);
			r.inf = x * y.inf;
			fesetround(FE_UPWARD);
			r.sup = x * y.sup;
		} else {
			fesetround(FE_DOWNWARD);
			r.inf = x * y.sup;
			fesetround(FE_UPWARD);
			r.sup = x * y.inf;
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval& operator*=(interval& x, const interval& y) {
		x = x * y;
		return x;
	}

	friend interval& operator*=(interval& x, const double& y) {
		x = x * y;
		return x;
	}

	friend interval operator/(const interval& x, const interval& y) {
		interval r;

		if (y.inf > 0) {
			if (x.inf >= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf / y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.sup /  y.inf;
			} else if (x.sup <= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf / y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.sup / y.sup;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf / y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.sup / y.inf;
			}
		} else if (y.sup < 0) {
			if (x.inf >= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup / y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.inf / y.inf;
			} else if (x.sup <= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup / y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.inf / y.sup;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup / y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.inf / y.sup;
			}
		} else {
			fesetround(FE_TONEAREST);
			throw std::domain_error("interval: division by 0");
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator/(const interval& x, const double& y) {
		interval r;

		if (y > 0) {
			fesetround(FE_DOWNWARD);
			r.inf = x.inf / y;
			fesetround(FE_UPWARD);
			r.sup = x.sup / y;
		} else if (y < 0) {
			fesetround(FE_DOWNWARD);
			r.inf = x.sup / y;
			fesetround(FE_UPWARD);
			r.sup = x.inf / y;
		} else {
			fesetround(FE_TONEAREST);
			throw std::domain_error("interval: division by 0");
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator/(const double& x, const interval& y) {
		interval r;

		if (y.inf > 0 || y.sup < 0) {
			if (x >= 0) {
				fesetround(FE_DOWNWARD);
				r.inf = x / y.sup;
				fesetround(FE_UPWARD);
				r.sup = x / y.inf;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x / y.inf;
				fesetround(FE_UPWARD);
				r.sup = x / y.sup;
			}
		} else {
			fesetround(FE_TONEAREST);
			throw std::domain_error("interval: division by 0");
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval& operator/=(interval& x, const interval& y) {
		x = x / y;
		return x;
	}

	friend interval& operator/=(interval& x, const double& y) {
		x = x / y;
		return x;
	}

	friend interval sqrt(const interval& x) {
		interval r;

		if (x.inf < 0) {
			throw std::domain_error("interval: sqrt of negative value");
		}

		fesetround(FE_DOWNWARD);
		r.inf = sqrt(x.inf);
		fesetround(FE_UPWARD);
		r.sup = sqrt(x.sup);
		fesetround(FE_TONEAREST);

		return r;
	}

	friend std::ostream& operator<<(std::ostream& s, const interval& x) {
		s << '[' << x.inf << ',' << x.sup << ']';
		return s;
	}

	const double& lower() const {
		return inf;
	}
	const double& upper() const {
		return sup;
	}
	double& lower() {
		return inf;
	}
	double& upper() {
		return sup;
	}
};

#endif // INTERVAL_HPP