#include "interval.h"

interval interval_init(double x)
{
	interval z;
	z.inf = x;
	z.sup = x;
	return z;
}

interval interval_init2(double x, double y)
{
	interval z;
	z.inf = x;
	z.sup = y;
	return z;
}

interval interval_add(interval x, interval y)
{
	interval z;

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

	return z;
}

interval interval_sub(interval x, interval y)
{
	interval z;

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

	return z;
}

interval interval_minus(interval x)
{
	interval z;

	z.inf = - x.sup;
	z.sup = - x.inf;

	return z;
}

interval interval_mul(interval x, interval y)
{
	interval z;
	double tmp;

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

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

	fesetround(FE_TONEAREST);

	return z;
}

interval interval_div(interval x, interval y)
{
	interval z;
	double tmp;

	if (y.inf <= 0 && y.sup >= 0) {
		printf("zero divide!\n");
		exit(1);
	}

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

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

	fesetround(FE_TONEAREST);

	return z;
}

interval interval_sqrt(interval x)
{
	interval z;

	if (x.inf < 0) {
		printf("sqrt of negative value!\n");
		exit(1);
	}

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

	return z;
}

void interval_print(interval x)
{
	printf("[%.17g,%.17g]\n", x.inf, x.sup);
}