#ifndef INTERVAL_HPP
#define INTERVAL_HPP

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

struct interval {

	double inf;
	double sup;

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

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

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

	friend interval operator+(interval x, 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-(interval x, 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-(interval x) {
		interval r;

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

		return r;
	}

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

		fesetround(FE_DOWNWARD);
		r.inf = x.inf * y.inf;
		tmp = x.inf * y.sup;
		if (tmp < r.inf) r.inf = tmp;
		tmp = x.sup * y.inf;
		if (tmp < r.inf) r.inf = tmp;
		tmp = x.sup * y.sup;
		if (tmp < r.inf) r.inf = tmp;

		fesetround(FE_UPWARD);
		r.sup = x.inf * y.inf;
		tmp = x.inf * y.sup;
		if (tmp > r.sup) r.sup = tmp;
		tmp = x.sup * y.inf;
		if (tmp > r.sup) r.sup = tmp;
		tmp = x.sup * y.sup;
		if (tmp > r.sup) r.sup = tmp;

		fesetround(FE_TONEAREST);

		return r;
	}

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

		if (y.inf <= 0. && y.sup >= 0.) {
			throw std::domain_error("interval: division by 0");
		}

		fesetround(FE_DOWNWARD);
		r.inf = x.inf / y.inf;
		tmp = x.inf / y.sup;
		if (tmp < r.inf) r.inf = tmp;
		tmp = x.sup / y.inf;
		if (tmp < r.inf) r.inf = tmp;
		tmp = x.sup / y.sup;
		if (tmp < r.inf) r.inf = tmp;

		fesetround(FE_UPWARD);
		r.sup = x.inf / y.inf;
		tmp = x.inf / y.sup;
		if (tmp > r.sup) r.sup = tmp;
		tmp = x.sup / y.inf;
		if (tmp > r.sup) r.sup = tmp;
		tmp = x.sup / y.sup;
		if (tmp > r.sup) r.sup = tmp;

		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval sqrt(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, interval x) {
		s << '[' << x.inf << ',' << x.sup << ']';
		return s;
	}

};

#endif // INTERVAL_HPP