/*
 * Copyright (c) 2014-2026 Masahide Kashiwagi (kashi@waseda.jp)
 *
 * simplified version of interval.hpp from kv-0.4.59
 *  - no template, only double, single file, same functions
 */

#ifndef INTERVAL_HPP
#define INTERVAL_HPP

#include <iostream>
#include <iomanip>
#include <limits>
#include <stdexcept>
#include <cmath>
#include <string>
#include <sstream>
#include <cstdio>
#include <cctype>
#include <vector>
#include <list>
#include <cstdlib>
#include <algorithm>

#if defined(_MSC_VER)
#include <float.h>
#else
#include <fenv.h>
#endif


struct conv_double {

	static int get_sign_double(double x) {
		if (x == 0.) {
			x = 1. / x;
		}

		if (x > 0.) return 1;
		else return -1;
	}

	static int get_exponent(double x) {
		int i;

		if (x >= std::ldexp(1., 1023)) return 1023;
		if (x < std::ldexp(1., -1074)) return -1075;

		std::frexp(x, &i);

		return i - 1;
	}

	// convert double number to string
	// mode == -1 : down
	// mode ==  0 : nearest
	// mode ==  1 : up
	// format == 'e' : like %e of printf
	// format == 'f' : like %f of printf
	// format == 'g' : like %g of printf
	// format == 'a' : print all digits with no rounding

	static std::string dtostring(double x, int precision = 17, char format = 'g', int mode = 0) {
		int i, j;
		int sign, ex;
		double absx;

		if (x != x) return "nan";

		sign = get_sign_double(x);
		absx = std::fabs(x);

		if (absx == 0.) {
			if (sign == -1) {
				return "-0";
			} else {
				return "0";
			}
		}

		if (absx == std::numeric_limits<double>::infinity()) {
			if (sign == -1) {
				return "-inf";
			} else {
				return "inf";
			}
		}

		ex = get_exponent(absx);

		bool buf[1023 - (-1074) + 1];
		int offset = 1074;
		int emax, emin;
		double dtmp, dtmp2;

		dtmp = absx;
		dtmp2 = std::ldexp(1., ex);

		for (i=0; i<=52; i++) {
			if (dtmp >= dtmp2) {
				buf[offset + ex - i] = 1;
				dtmp -= dtmp2;
			} else {
				buf[offset + ex - i] = 0;
			}
			if (dtmp == 0) {
				emax = ex;
				emin = ex - i;
				break;
			}
			dtmp2 /= 2.;
		}

		if (emin > 0) {
			for (i=0; i<emin; i++) {
				buf[offset + i] = 0;
			}
			emin = 0;
		}

		if (emax < 0) {
			for (i=emax + 1; i<=0; i++) {
				buf[offset + i] = 0;
			}
			emax = 0;
		}

		std::list<int> result1, result2;
		int result_max, result_min, m, pm, carry, tmp;

		result_max = -1;

		while (emax >= 0) {
			if (emax >= 17) m= 5;
			else if (emax >= 14) m = 4;
			else if (emax >= 10) m = 3;
			else if (emax >= 7) m = 2;
			else  m = 1;

			pm = 1;
			for (i=0; i<m; i++) pm *= 10;

			carry = 0;
			for (i=emax; i>=0; i--) {
				tmp = carry * 2 + buf[offset + i];
				buf[offset + i] = tmp / pm;
				carry = tmp % pm;
			}

			for (i=0; i<m; i++) {
				result_max++;
				result1.push_back(carry  % 10);
				carry /= 10;
			}

			while (emax >= 0 && buf[offset + emax] == 0) {
				emax--;
			}
		}

		result_min = 0;

		while (emin < 0) {
			m = std::min(8, -emin);
			pm = 1;
			for (i=0; i<m; i++) pm *= 10;

			carry = 0;
			for (i=emin; i<=-1; i++) {
				tmp = buf[offset + i] * pm + carry;
				buf[offset + i] = tmp % 2;
				carry = tmp / 2;
			}

			for (i=0; i<m; i++) {
				result_min--;
				pm /= 10;
				result2.push_back(carry / pm);
				carry %= pm;
			}

			while (emin < 0 && buf[offset + emin] == 0) {
				emin++;
			}
		}

		std::vector<int> result;
		int offset2;

		// add 1byte margin to both ends of array
		result.resize(result_max - result_min + 1 + 2);
		offset2 = - result_min + 1;
		for (i=0; i<=result_max; i++) {
			result[offset2 + i] = result1.front();
			result1.pop_front();
		}
		for (i=-1; i>=result_min; i--) {
			result[offset2 + i] = result2.front();
			result2.pop_front();
		}

		#if 0
		for (i=result_min; i<=result_max; i++) {
			std::cout << i << ':' << result[offset2 + i] << "\n";
		}
		#endif

		std::ostringstream result_str;

		if (sign == 1) {
			result_str << "";
		} else {
			result_str << "-";
		}

		if (format == 'f') {
			// round to precision after decimal point
			if (-(precision+1) >= result_min) {
				result_min = -precision;
				tmp = result[offset2 + result_min - 1];
				if ((mode == 1 && sign == 1) || (mode == -1 && sign == -1) || (mode == 0 && tmp >= 5)) {
					result[offset2 + result_max + 1] = 0;
					result_max++;
					for (i=result_min; i<=result_max; i++) {
						result[offset2 + i]++;
						if (result[offset2 + i] != 10) break;
						result[offset2 + i] = 0;
					}
					if (result[offset2 + result_max] == 0) {
						result_max--;
					}
				}
			}

			// delete zeros of tail
			while (result_min < 0 && result[offset2 + result_min] == 0) {
				result_min++;
			}

			// make result string
			for (i=result_max; i>=result_min; i--) {
				if (i == -1) result_str << ".";
				result_str << result[offset2 + i];
			}

		} else if (format == 'e') {
			// delete zeros of head
			while (result[offset2 + result_max] == 0) {
				result_max--;
			}

			// round to precision
			if (result_max-precision-1 >= result_min) {
				result_min = result_max - precision;
				tmp = result[offset2 + result_min - 1];
				if ((mode == 1 && sign == 1) || (mode == -1 && sign == -1) || (mode == 0 && tmp >= 5)) {
					result[offset2 + result_max + 1] = 0;
					result_max++;
					for (i=result_min; i<=result_max; i++) {
						result[offset2 + i]++;
						if (result[offset2 + i] != 10) break;
						result[offset2 + i] = 0;
					}
					if (result[offset2 + result_max] == 0) {
						result_max--;
					} else {
						result_min++;
					}
				}
			}

			// delete zeros of tail
			while (result[offset2 + result_min] == 0) {
				result_min++;
			}

			// make result string
			for (i=result_max; i>=result_min; i--) {
				if (i == result_max -1) result_str << ".";
				result_str << result[offset2 + i];
			}
			// sprintf(stmp, "e%+03d", result_max);
			result_str << "e" << std::internal << std::showpos << std::setfill('0') << std::setw(3) << result_max << std::noshowpos;

		} else if (format == 'g') {
			// delete zeros of head
			while (result[offset2 + result_max] == 0) {
				result_max--;
			}

			// round to precision
			if (result_max-precision >= result_min) {
				result_min = result_max - precision + 1;
				tmp = result[offset2 + result_min - 1];
				if ((mode == 1 && sign == 1) || (mode == -1 && sign == -1) || (mode == 0 && tmp >= 5)) {
					result[offset2 + result_max + 1] = 0;
					result_max++;
					for (i=result_min; i<=result_max; i++) {
						result[offset2 + i]++;
						if (result[offset2 + i] != 10) break;
						result[offset2 + i] = 0;
					}
					if (result[offset2 + result_max] == 0) {
						result_max--;
					} else {
						result_min++;
					}
				}
			}

			if (-4 <= result_max && result_max <= precision -1) {
				// use 'f' like format

				// delete zeros of tail
				while (result_min < 0 && result[offset2 + result_min] == 0) {
					result_min++;
				}

				if (result_max < 0) {
					result_max = 0;
				}

				// make result string
				for (i=result_max; i>=result_min; i--) {
					if (i == -1) result_str << ".";
					result_str << result[offset2 + i];
				}

			} else {
				// use 'e' like format

				// delete zeros of tail
				while (result[offset2 + result_min] == 0) {
					result_min++;
				}

				// make result string
				for (i=result_max; i>=result_min; i--) {
					if (i == result_max -1) result_str << ".";
					result_str << result[offset2 + i];
				}
				// sprintf(stmp, "e%+03d", result_max);
				result_str << "e" << std::internal << std::showpos << std::setfill('0') << std::setw(3) << result_max << std::noshowpos;
			}

		} else if (format == 'a') {
			// make result string
			for (i=result_max; i>=result_min; i--) {
				if (i == -1) result_str << ".";
				result_str << result[offset2 + i];
			}
		}

		return result_str.str();
	}


	static void ignore_space(std::string& s) {
		int p = 0;
		while (p < s.size() && isspace(s[p])) p++;
		s = s.substr(p);
	}

	static int get_sign(std::string& s) {
		int r, p;

		p = 0;
		r = 1;
		if (p < s.size() && s[p] == '-') {
			r = -1;
			p++;
		} else if (p < s.size() && s[p] == '+') {
			r = 1;
			p++;
		}

		s = s.substr(p);
		return r;
	}

	static std::string get_number(std::string& s) {
		std::string r;
		int p;

		p = 0;
		r = "";
		while (p < s.size() && isdigit(s[p])) {
			r += s[p];
			p++;
		}

		s = s.substr(p);
		return r;
	}

	// convert string to double number
	// mode == -1 : down
	// mode ==  0 : nearest
	// mode ==  1 : up

	static double stringtod(std::string s, int mode = 0) {
		int i, j, tmp;
		bool flag;
		int sign, e10, esign;
		std::string num1_s, num2_s, nume_s;

		ignore_space(s);
		sign = get_sign(s);
		num1_s = get_number(s);

		if (0 < s.size() && s[0] == '.') {
			s = s.substr(1);
			num2_s = get_number(s);
		}

		if (0 < s.size() && (s[0] == 'e' || s[0] == 'E')) {
			s = s.substr(1);
			esign = get_sign(s);
			nume_s = get_number(s);
			e10 = esign * atoi(nume_s.c_str());
		} else {
			e10 = 0;
		}

		// delete 0s from the head of num1_s
		while (0 < num1_s.size() && num1_s[0] == '0') {
			num1_s = num1_s.substr(1);
		}

		// delete 0s from the tail of num2_s
		while (0 < num2_s.size() && num2_s[num2_s.size() - 1] == '0') {
			num2_s = num2_s.substr(0, num2_s.size() - 1);
		}

		// set table and offset
		// |x| = \sum_{table_min}^{table_max} table[offset + i] * 10^i
		int table_max, table_min, offset;
		std::vector<int> table;

		table_max = num1_s.size() - 1 + e10;
		table_min = e10 - num2_s.size();
		table.resize(table_max - table_min + 1);
		offset = - table_min;

		for (i=0; i<num1_s.size(); i++) {
			table[offset + num1_s.size() - 1 - i + e10] = num1_s[i] - '0';
		}

		for (i=0; i<num2_s.size(); i++) {
			table[offset - i - 1 + e10] = num2_s[i] - '0';
		}

		// extend table
		if (table_min > 0) {
			tmp = table.size();
			table.resize(tmp + table_min);
			for (i=tmp-1; i>=0; i--) {
				table[i + table_min] = table[i];
			}
			for (i=0; i<table_min; i++) {
				table[i] = 0;
			}
			offset += table_min;
			table_min = 0;
		}

		if (table_max < -1) {
			tmp = table.size();
			table.resize(tmp + (-1-table_max));
			for (i=0; i<(-1-table_max); i++) {
				table[tmp + i] = 0;
			}
			table_max = -1;
		}

		#if 0
		for (i=table_max; i>=table_min; i--) {
			std::cout << i << ':' << table[offset + i] << "\n";
		}
		#endif

		// convert decimal number to binary number
		// set result and offset2
		// |x| = \sum_{result_min}^{reuslt_max} result[offset2 + i] * 2^i

		int result_min, result_max, m, pm, carry, carry2;
		std::list<bool> result1, result2;

		// integer part

		result_max = -1;

		while (table_max >= 0) {
			if (table_max >= 5) m = 16;
			else if (table_max >= 4) m = 13;
			else if (table_max >= 3) m = 9;
			else if (table_max >= 2) m = 6;
			else if (table_max >= 1) m = 3;
			else m = 1;
			pm = 1 << m;
			
			carry = 0;
			for (i=table_max; i>=0; i--) {
				tmp = carry * 10 + table[offset + i];
				table[offset + i] = tmp / pm;
				carry = tmp % pm;
			}
			for (i=0; i<m; i++) {
				result_max++;
				result1.push_back(carry % 2);
				carry = carry / 2;
			}
			while (table_max >= 0 && table[offset + table_max] == 0) {
				table_max--;
			}
		}

		// fraction part

		//  flag means whether most significant bit already found or not
		if (result_max >= 0) flag = true;
		else flag = false;

		result_min = 0;

		while (table_min < 0) {
			tmp = 53 - (result_max - result_min);
			if (flag && tmp <= 0) break;
			if (!flag) {
				m = 16;
			} else {
				m = std::min(16, tmp);
			}
			pm = 1 << m;

			carry = 0;
			for (i=table_min; i<=-1; i++) {
				tmp = table[offset + i] * pm + carry;
				table[offset + i] = tmp % 10;
				carry = tmp / 10;
			}

			for (i=0; i<m; i++) {
				result_min--;
				pm /= 2;
				carry2 = carry / pm;
				carry = carry % pm;

				if (flag) {
					result2.push_back(carry2);
				} else {
					if (carry2 != 0) {
						result2.push_back(carry2);
						result_max = result_min;
						flag = true;
					}
				}
			}

			while (table_min < 0 && table[offset + table_min] == 0) {
				table_min++;
			}
		}

		// append integer and fraction part

		std::vector<bool> result;
		int offset2;

		result.resize(result_max - result_min + 1);
		offset2 = - result_min;
		for (i=0; i<=result_max; i++) {
			result[offset2 + i] = result1.front();
			result1.pop_front();
		}
		for (i=std::min(-1, result_max); i>=result_min; i--) {
			result[offset2 + i] = result2.front();
			result2.pop_front();
		}

		#if 0
		for (i=result_max; i>=result_min; i--) {
			printf("%d %d\n", i, result[offset2 + i]);
		}
		#endif

		// convert binary to double number

		if (result_max > 1023) {
			if ((sign == 1 && mode == -1) || (sign == -1 && mode == 1)) {
                        	return sign * (std::numeric_limits<double>::max)();
			}

			return sign * std::numeric_limits<double>::infinity();
		}

		if (result_max < -1075) {
			if ((sign == 1 && mode == 1) || (sign == -1 && mode == -1)) {
				return sign * std::numeric_limits<double>::denorm_min();
			}
			return sign * 0.;
		}

		double r;

		r = 0.;
		for (i=result_max; i >= result_min; i--) {
			if (result_max - i == 53 || i == -1075) {
				if (sign == 1) {
					if (mode == -1) {
					} else if (mode == 0) {
						if (result[offset2 + i] == 0) {
						} else {
							r += std::ldexp(1., i+1);
						}
					} else {
						r += std::ldexp(1., i+1);
					}
				} else {
					if (mode == -1) {
						r += std::ldexp(1., i+1);
					} else if (mode == 0) {
						if (result[offset2 + i] == 0) {
						} else {
							r += std::ldexp(1., i+1);
						}
					} else {
					}
				}
				break;
			}
			r += std::ldexp((double)result[offset2 + i], i);
		}

		return sign * r;
	}
};


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;
	}

	interval(const std::string& x) {
		inf = fromstring_down(x);
		sup = fromstring_up(x);
	}

	interval(const std::string& x, const std::string& y) {
		inf = fromstring_down(x);
		sup = fromstring_up(y);
	}

	interval& operator=(const double& x) {
		inf = x;
		sup = x;
		return *this;
	}

	interval& operator=(const std::string& x) {
		inf = fromstring_down(x);
		sup = fromstring_up(x);
		return *this;
	}

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

		r.inf = add_down(x.inf, y.inf);
		r.sup = add_up(x.sup, y.sup);

		return r;
	}

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

		r.inf = add_down(x.inf, y);
		r.sup = add_up(x.sup, y);
		return r;
	}

	friend interval operator+(const interval& x, const std::string& y) {
		return x + interval(y);
	}

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

		r.inf = add_down(x, y.inf);
		r.sup = add_up(x, y.sup);

		return r;
	}

	friend interval operator+(const std::string& x, const interval& y) {
		return interval(x) + y;
	}

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

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

		x.inf = add_down(x.inf, y);
		x.sup = add_up(x.sup, y);

		return x;
	}

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

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

		r.inf = sub_down(x.inf, y.sup);
		r.sup = sub_up(x.sup, y.inf);

		return r;
	}

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

		r.inf = sub_down(x.inf, y);
		r.sup = sub_up(x.sup, y);

		return r;
	}

	friend interval operator-(const interval& x, const std::string& y) {
		return x - interval(y);
	}

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

		r.inf = sub_down(x, y.sup);
		r.sup = sub_up(x, y.inf);

		return r;
	}

	friend interval operator-(const std::string& x, const interval& y) {
		return interval(x) - y;
	}

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

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

		x.inf = sub_down(x.inf, y);
		x.sup = sub_up(x.sup, y);

		return x;
	}

	friend interval& operator-=(interval& x, const std::string& y) {
		x = x - interval(y);
		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;

		using std::abs;

		if (x.inf >= 0.) {
			if (x.sup == 0.) {
				r = interval(0., 0.);
			} else {
				if (y.inf >= 0.) {
					if (y.sup == 0.) {
						r = interval(0., 0.);
					} else {
						r.inf = mul_down(x.inf, y.inf);
						r.sup = mul_up(x.sup, y.sup);
					}
				} else if (y.sup <= 0.) {
					r.inf = mul_down(x.sup, y.inf);
					r.sup = mul_up(x.inf, y.sup);
				} else {
					r.inf = mul_down(x.sup, y.inf);
					r.sup = mul_up(x.sup, y.sup);
				}
			}
		} else if (x.sup <= 0.) {
			if (y.inf >= 0.) {
				if (y.sup == 0.) {
					r = interval(0., 0.);
				} else {
					r.inf = mul_down(x.inf, y.sup);
					r.sup = mul_up(x.sup, y.inf);
				}
			} else if (y.sup <= 0.) {
				r.inf = mul_down(x.sup, y.sup);
				r.sup = mul_up(x.inf, y.inf);
			} else {
				r.inf = mul_down(x.inf, y.sup);
				r.sup = mul_up(x.inf, y.inf);
			}
		} else {
			if (y.inf >= 0.) {
				if (y.sup == 0.) {
					r = interval(0., 0.);
				} else {
					r.inf = mul_down(x.inf, y.sup);
					r.sup = mul_up(x.sup, y.sup);
				}
			} else if (y.sup <= 0.) {
				r.inf = mul_down(x.sup, y.inf);
				r.sup = mul_up(x.inf, y.inf);
			} else {
				r.inf = mul_down(x.inf, y.sup);
				tmp = mul_down(x.sup, y.inf);
				if (tmp < r.inf) r.inf = tmp;
				r.sup = mul_up(x.inf, y.inf);
				tmp = mul_up(x.sup, y.sup);
				if (tmp > r.sup) r.sup = tmp;
			}
		}

		return r;
	}

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

		using std::abs;

		if (y > 0.) {
			r.inf = mul_down(x.inf, y);
			r.sup = mul_up(x.sup, y);
		} else if (y < 0.) {
			r.inf = mul_down(x.sup, y);
			r.sup = mul_up(x.inf, y);
		} else {
			r = interval(0., 0.);
		}

		return r;
	}

	friend interval operator*(const interval& x, const std::string& y) {
		return x * interval(y);
	}

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

		using std::abs;

		if (x > 0.) {
			r.inf = mul_down(x, y.inf);
			r.sup = mul_up(x, y.sup);
		} else if (x < 0.) {
			r.inf = mul_down(x, y.sup);
			r.sup = mul_up(x, y.inf);
		} else {
			r = interval(0., 0.);
		}

		return r;
	}

	friend interval operator*(const std::string& x, const interval& y) {
		return interval(x) * y;
	}

	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*=(interval& x, const std::string& y) {
		x = x * interval(y);
		return x;
	}

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

		if (y.inf > 0.) {
			if (x.inf >= 0.) {
				r.inf = div_down(x.inf, y.sup);
				r.sup = div_up(x.sup, y.inf);
			} else if (x.sup <= 0.) {
				r.inf = div_down(x.inf, y.inf);
				r.sup = div_up(x.sup, y.sup);
			} else {
				r.inf = div_down(x.inf, y.inf);
				r.sup = div_up(x.sup, y.inf);
			}
		} else if (y.sup < 0.) {
			if (x.inf >= 0.) {
				r.inf = div_down(x.sup, y.sup);
				r.sup = div_up(x.inf, y.inf);
			} else if (x.sup <= 0.) {
				r.inf = div_down(x.sup, y.inf);
				r.sup = div_up(x.inf, y.sup);
			} else {
				r.inf = div_down(x.sup, y.sup);
				r.sup = div_up(x.inf, y.sup);
			}
		} else {
			throw std::domain_error("interval: division by 0");
		}

		return r;
	}

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

		if (y > 0.) {
			r.inf = div_down(x.inf, y);
			r.sup = div_up(x.sup, y);
		} else if (y < 0.) {
			r.inf = div_down(x.sup, y);
			r.sup = div_up(x.inf, y);
		} else {
			throw std::domain_error("interval: division by 0");
		}

		return r;
	}

	friend interval operator/(const interval& x, const std::string& y) {
		return x / interval(y);
	}

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

		if (y.inf > 0. || y.sup < 0.) {
			if (x >= 0.) {
				r.inf = div_down(x, y.sup);
				r.sup = div_up(x, y.inf);
			} else {
				r.inf = div_down(x, y.inf);
				r.sup = div_up(x, y.sup);
			}
		} else {
			throw std::domain_error("interval: division by 0");
		}

		return r;
	}

	friend interval operator/(const std::string& x, const interval& y) {
		return interval(x) / y;
	}

	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/=(interval& x, const std::string& y) {
		x = x / interval(y);
		return x;
	}

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

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

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

		r.inf = sqrt_down(x.inf);
		r.sup = sqrt_up(x.sup);

		return r;
	}

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

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

	static interval whole() {
		return interval(-std::numeric_limits<double>::infinity(),
		                 std::numeric_limits<double>::infinity() );
	}

	static interval hull(const double& x, const double& y) {
		if (x < y) return interval(x, y);
		else return interval(y, x);
	}

	static interval hull(const interval& x, const interval& y) {
		double tmp1, tmp2;
		tmp1 = x.inf;
		if (y.inf < tmp1) tmp1 = y.inf;
		tmp2 = x.sup;
		if (y.sup > tmp2) tmp2 = y.sup;
		return interval(tmp1, tmp2);
	}

	static interval hull(const interval& x, const double& y) {
		double tmp1, tmp2;
		tmp1 = x.inf;
		if (y < tmp1) tmp1 = y;
		tmp2 = x.sup;
		if (y > tmp2) tmp2 = y;
		return interval(tmp1, tmp2);
	}

	static interval hull(const double& x, const interval& y) {
		double tmp1, tmp2;
		tmp1 = x;
		if (y.inf < tmp1) tmp1 = y.inf;
		tmp2 = x;
		if (y.sup > tmp2) tmp2 = y.sup;
		return interval(tmp1, tmp2);
	}


	friend double width(const interval& x) {
		double tmp;

		tmp = sub_up(x.sup, x.inf);

		return tmp;
	}

	friend double rad(const interval& x) {
		double tmp;

		tmp = mul_up(sub_up(x.sup, x.inf), 0.5);

		return tmp;
	}

	friend double median(const interval& x) {
		return (x.inf + x.sup) * 0.5;
	}

	friend double mid(const interval& x) {
		return median(x);
	}

	friend void midrad(const interval& x, double& m, double& r) {
		double tmp;
		m = mid(x);
		r = sub_up(x.sup, m);
		tmp = sub_up(m, x.inf);
		if (tmp > r) r = tmp;
	}

	friend double norm(const interval& x) {
		if (x.inf >= 0.) return x.sup;
		if (x.sup <= 0.) return -x.inf;
		double tmp = -x.inf;
		if (x.sup > tmp) tmp = x.sup;
		return tmp;
	}

	friend double mag(const interval& x) {
		return norm(x);
	}

	friend double mig(const interval& x) {
		if (zero_in(x)) return 0.;
		if (x.inf > 0.) return x.inf;
		else return -x.sup;
	}

	friend interval abs(const interval& x) {
		if (x.inf >= 0.) return x;
		if (x.sup <= 0.) return -x;
		double tmp = -x.inf;
		if (x.sup > tmp) tmp = x.sup;
		return interval(0., tmp);
	}

	friend interval max(const interval& x, const interval& y) {
		double z1, z2;
		z1 = x.inf;
		if (y.inf > z1) z1 = y.inf;
		z2 = x.sup;
		if (y.sup > z2) z2 = y.sup;
		return interval(z1, z2);
	}

	friend interval min(const interval& x, const interval& y) {
		double z1, z2;
		z1 = x.inf;
		if (y.inf < z1) z1 = y.inf;
		z2 = x.sup;
		if (y.sup < z2) z2 = y.sup;
		return interval(z1, z2);
	}

	friend interval floor(const interval& x) {
		using std::floor;
		return interval(floor(x.inf), floor(x.sup));
	}

        friend interval ceil(const interval& x) {
		using std::ceil;
		return interval(ceil(x.inf), ceil(x.sup));
	}

	friend bool in(const double& a, const interval& x) {
		return x.inf <= a && a <= x.sup;
	}

	friend bool zero_in(const interval& x) {
		return x.inf <= 0. && 0. <= x.sup;
	}

	friend bool subset(const interval& x, const interval& y) {
		return y.inf <= x.inf && x.sup <= y.sup;
	}

	friend bool proper_subset(const interval& x, const interval& y) {
		return subset(x, y) && (y.inf < x.inf || x.sup < y.sup);
	}

	friend bool overlap(const interval& x, const interval& y) {
		double tmp1, tmp2;
		tmp1 = x.inf;
		if (y.inf > tmp1) tmp1 = y.inf;
		tmp2 = x.sup;
		if (y.sup < tmp2) tmp2 = y.sup;
		return tmp1 <= tmp2;
	}

	friend interval intersect(const interval& x, const interval& y) {
		double tmp1, tmp2;
		tmp1 = x.inf;
		if (y.inf > tmp1) tmp1 = y.inf;
		tmp2 = x.sup;
		if (y.sup < tmp2) tmp2 = y.sup;
		return interval(tmp1, tmp2);
	}


	friend bool operator<(const interval& x, const interval& y) {
		return x.sup < y.inf;
	}

	friend bool operator<(const interval& x, const double& y) {
		return x.sup < y;
	}

	friend bool operator<(const interval& x, const std::string& y) {
		return x.sup < fromstring_down(y);
	}

	friend bool operator<(const double& x, const interval& y) {
		return x < y.inf;
	}

	friend bool operator<(const std::string& x, const interval& y) {
		return fromstring_up(x) < y.inf;
	}


	friend bool operator<=(const interval& x, const interval& y) {
		return x.sup <= y.inf;
	}

	friend bool operator<=(const interval& x, const double& y) {
		return x.sup <= y;
	}

	friend bool operator<=(const interval& x, const std::string& y) {
		return x.sup <= fromstring_down(y);
	}

	friend bool operator<=(const double& x, const interval& y) {
		return x <= y.inf;
	}

	friend bool operator<=(const std::string& x, const interval& y) {
		return fromstring_up(x) <= y.inf;
	}


	friend bool operator>(const interval& x, const interval& y) {
		return x.inf > y.sup;
	}

	friend bool operator>(const interval& x, const double& y) {
		return x.inf > y;
	}

	friend bool operator>(const interval& x, const std::string& y) {
		return x.inf > fromstring_up(y);
	}

	friend bool operator>(const double& x, const interval& y) {
		return x > y.sup;
	}

	friend bool operator>(const std::string& x, const interval& y) {
		return fromstring_down(x) > y.sup;
	}


	friend bool operator>=(const interval& x, const interval& y) {
		return x.inf >= y.sup;
	}

	friend bool operator>=(const interval& x, const double& y) {
		return x.inf >= y;
	}

	friend bool operator>=(const interval& x, const std::string& y) {
		return x.inf >= fromstring_up(y);
	}

	friend bool operator>=(const double& x, const interval& y) {
		return x >= y.sup;
	}

	friend bool operator>=(const std::string& x, const interval& y) {
		return fromstring_down(x) >= y.sup;
	}


	friend bool operator==(const interval& x, const interval& y) {
		return x.inf == x.sup && x.sup == y.inf && y.inf == y.sup;
	}

	friend bool operator==(const interval& x, const double& y) {
		return x.inf == x.sup && x.sup == y;
	}

	friend bool operator==(const interval& x, const std::string& y) {
		interval ytmp(y);
		return x.inf == x.sup && x.sup == ytmp.inf && ytmp.inf == ytmp.sup;
	}

	friend bool operator==(const double& x, const interval& y) {
		return y.inf == y.sup && y.sup == x;
	}

	friend bool operator==(const std::string& x, const interval& y) {
		interval xtmp(x);
		return xtmp.inf == xtmp.sup && xtmp.sup == y.inf && y.inf == y.sup;
	}


	friend bool operator!=(const interval& x, const interval& y) {
		return !overlap(x, y);
	}

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

	friend bool operator!=(const interval& x, const std::string& y) {
		return !overlap(x, interval(y));
	}

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

	friend bool operator!=(const std::string& x, const interval& y) {
		return !overlap(interval(x), y);
	}


	friend interval division_part1(const interval& x, const interval& y, bool& parted) {
		interval r;

		parted = false;

		if (y.inf > 0. || y.sup < 0.) {
			return x / y;
		}

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

		if (y.inf < 0. && y.sup == 0.) {
			if (x.inf == 0. && x.sup == 0) {
				r.inf = 0.;
				r.sup = 0.;
			} else if (x.inf >= 0.) {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = div_up(x.inf, y.inf);
			} else if (x.sup <= 0.) {
				r.inf = div_down(x.sup, y.inf);
				r.sup = std::numeric_limits<double>::infinity();
			} else {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = std::numeric_limits<double>::infinity();
			}
		} else if (y.inf == 0. && y.sup > 0.) {
			if (x.inf == 0. && x.sup == 0) {
				r.inf = 0.;
				r.sup = 0.;
			} else if (x.inf >= 0.) {
				r.inf = div_down(x.inf, y.sup);
				r.sup = std::numeric_limits<double>::infinity();
			} else if (x.sup <= 0.) {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = div_up(x.sup, y.sup);
			} else {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = std::numeric_limits<double>::infinity();
			}
		} else {
			if (x.inf == 0. && x.sup == 0) {
				r.inf = 0.;
				r.sup = 0.;
			} else if (x.inf > 0.) {
				parted = true;
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = div_up(x.inf, y.inf);
			} else if (x.sup < 0.) {
				parted = true;
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = div_up(x.sup, y.sup);
			} else {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = std::numeric_limits<double>::infinity();
			}
		}

		return r;
	}

	friend interval division_part2(const interval& x, const interval& y, bool parted = true) {
		interval r;

		if (y.inf > 0. || y.sup < 0.) {
			return x / y;
		}

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

		if (y.inf < 0. && y.sup == 0.) {
			if (x.inf == 0. && x.sup == 0) {
				r.inf = 0.;
				r.sup = 0.;
			} else if (x.inf >= 0.) {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = div_up(x.inf, y.inf);
			} else if (x.sup <= 0.) {
				r.inf = div_down(x.sup, y.inf);
				r.sup = std::numeric_limits<double>::infinity();
			} else {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = std::numeric_limits<double>::infinity();
			}
		} else if (y.inf == 0. && y.sup > 0.) {
			if (x.inf == 0. && x.sup == 0) {
				r.inf = 0.;
				r.sup = 0.;
			} else if (x.inf >= 0.) {
				r.inf = div_down(x.inf, y.sup);
				r.sup = std::numeric_limits<double>::infinity();
			} else if (x.sup <= 0.) {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = div_up(x.sup, y.sup);
			} else {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = std::numeric_limits<double>::infinity();
			}
		} else {
			if (x.inf == 0. && x.sup == 0) {
				r.inf = 0.;
				r.sup = 0.;
			} else if (x.inf > 0.) {
				r.inf = div_down(x.inf, y.sup);
				r.sup = std::numeric_limits<double>::infinity();
			} else if (x.sup < 0.) {
				r.inf = div_down(x.sup, y.inf);
				r.sup = std::numeric_limits<double>::infinity();
			} else {
				r.inf = -std::numeric_limits<double>::infinity();
				r.sup = std::numeric_limits<double>::infinity();
			}
		}

		return r;
	}

	static interval pow_point(const double& x, int y) {
		interval r, xp;

		r = 1.;
		xp = (interval)x;

		while (y != 0) {
			if (y % 2 != 0) {
				r *= xp;
			}
			y /= 2;
			xp *= xp;
		}

		return r;
	}

	friend interval pow(const interval& x, int y) {
		interval r, xp;
		int a;

		if (y == 0) return interval(1.);

		a = (y >= 0) ? y : -y;

		if (a % 2 == 0 && in(0., x)) {
			r = hull(0., pow_point(mag(x), a));
		} else {
			r = hull(pow_point(x.lower(), a), pow_point(x.upper(), a));
		}

		if (y < 0) {
			r = 1. / r;
		}

		return r;
	}

	friend interval pow(const interval& x, const interval& y) {
		return exp(y * log(x));
	}

	friend interval pow(const interval& x, const double& y) {
		return pow(x, interval(y));
	}

	friend interval pow(const interval& x, const std::string& y) {
		return pow(x, interval(y));
	}

	friend interval pow(const double& x, const interval& y) {
		return pow(interval(x), y);
	}

	friend interval pow(const std::string& x, const interval& y) {
		return pow(interval(x), y);
	}

	// power by integer
	static interval ipower(const interval& x, double i) {
		double tmp;
		interval xp = x;
		interval r(1.);

		if (i != i) {
			throw std::domain_error("interval: cannot calculate power of nan");
		}

		while (i != 0.) {
			i *= 0.5;
			using std::floor;
			tmp = floor(i);
			if (tmp != i) {
				i = tmp;
				r *= xp;
			}
			xp = xp * xp;
		}

		return r;
	}

	static interval exp_point(const double& x) {
		double x_i, x_f, tmp;
		interval r, y, remainder;
		int i;

		if (x == std::numeric_limits<double>::infinity()) {
			return interval((std::numeric_limits<double>::max)(), std::numeric_limits<double>::infinity());
		}
		if (x == -std::numeric_limits<double>::infinity()) {
			return interval(0., std::numeric_limits<double>::denorm_min());
		}

		remainder = interval((1./sqrt(e())).lower(), sqrt(e()).upper());

		using std::floor;
		if (x >= 0.) {
			x_i = floor(x);
			x_f = x - x_i;
			if (x_f >= 0.5) {
				x_f -= 1.;
				x_i += 1.;
			}
		} else {
			x_i = -floor(-x);
			x_f = x - x_i;
			if (x_f <= -0.5) {
				x_f += 1.;
				x_i -= 1.;
			}
		}

		r = 1.;
		y = 1.;
		for (i=1;  ; i++) {
			y *= x_f;
			y /= (double)i;
			if (mag(y) * remainder.upper() < std::numeric_limits<double>::epsilon()) {
				r += y * remainder;
				break;
			} else {
				r += y;
			}
		}

		if (x_i >= 0.) {
			r *= ipower(e(), x_i);
		} else {
			r /= ipower(e(), -x_i);
		}

		return r;
	}

	friend interval exp(const interval& I) {
		return interval(exp_point(I.lower()).lower(), exp_point(I.upper()).upper());
	}

	static interval expm1_origin(const double& x) {
		interval r, y, remainder;
		int i;

		remainder = interval((1./sqrt(e())).lower(), sqrt(e()).upper());

		r = 0.;
		y = 1.;
		for (i=1;  ; i++) {
			y *= x;
			y /= (double)i;
			if (mag(y) * remainder.upper() < std::numeric_limits<double>::epsilon()) {
				r += y * remainder;
				break;
			} else {
				r += y;
			}
		}

		return r;
	}

	static interval expm1_point(const double& x) {
		if (x >= -0.5 && x <= 0.5) {
			return expm1_origin(x);
		} else {
			return exp_point(x) - 1.;
		}
	}

	friend interval expm1(const interval& I) {
		return interval(expm1_point(I.lower()).lower(), expm1_point(I.upper()).upper());
	}

	static double log_point(const double& x, int round) {
		interval tmp;
		double x2, x2m1;
		interval cinv;
		interval r;
		interval xn, xn2;
		interval sqrt2 = sqrt(interval((double)2.));
		int p_i;
		double p;
		int i;

		if (x == std::numeric_limits<double>::infinity()) {
			if (round == 1) {
				return std::numeric_limits<double>::infinity();
			} else {
				return (std::numeric_limits<double>::max)();
			}
		}
		if (x == 0.) {
			if (round == 1) {
				return -(std::numeric_limits<double>::max)();
			} else {
				return -std::numeric_limits<double>::infinity();
			}
		}

		using std::frexp;
		x2 = frexp(x, &p_i);
		p = (double)p_i;

		while (x2 > 4. * std::sqrt(2.) - 4.) {
			x2 *= 0.5;
			p += 1.;
		}
		while (x2 > 4. - 2. * std::sqrt(2.)) {
			tmp = x2 / sqrt2;
			if (round == -1) x2 = tmp.lower();
			else x2 = tmp.upper();
			p += 0.5;
		}
		while (x2 < 2. - std::sqrt(2.)) {
			x2 *= 2.;
			p -= 1.;
		}
		while (x2 < 2. * std::sqrt(2.) - 2.) {
			tmp = x2 * sqrt2;
			if (round == -1) x2 = tmp.lower();
			else x2 = tmp.upper();
			p -= 0.5;
		}

		x2m1 = x2 - 1.;
		cinv = 1. / hull(x2, 1.);
		r = 0.;
		xn = -1.;
		xn2 = -1.;
		for (i=1;  ; i++) {
			xn = -xn * x2m1; 
			xn2 = -xn2 * cinv * x2m1;
			tmp = xn2 / (double)i;
			if (mag(tmp) < std::numeric_limits<double>::epsilon()) {
				r += tmp;
				break;
			} else {
				r += xn / (double)i;
			}
		}

		r += ln2() * p;

		if (round == -1) return r.lower();
		else return r.upper();
	}

	friend interval log(const interval& I) {
		if (I.inf < 0.) {
			throw std::domain_error("interval: log of negative value");
		}
		return interval(log_point(I.lower(), -1), log_point(I.upper(), 1));
	}

	static interval log1p_origin(const double& x) {
		interval tmp;
		interval cinv;
		interval r;
		interval xn, xn2;
		int i;

		cinv = 1. / hull(x + interval(1.), 1.);
		r = 0.;
		xn = -1.;
		xn2 = -1.;
		for (i=1;  ; i++) {
			xn = -xn * x; 
			xn2 = -xn2 * cinv * x;
			tmp = xn2 / (double)(i);
			if (mag(tmp) < std::numeric_limits<double>::epsilon()) {
				r += tmp;
				break;
			} else {
				r += xn / (double)(i);
			}
		}

		return r;
	}

	static double log1p_point(const double& x, int round) {
		interval tmp;

		if (x >= -(3. - 2. * std::sqrt(2.)) && x <= 3. - 2. * std::sqrt(2.)) {
			tmp = log1p_origin(x);
			if (round == -1) return tmp.lower();
			else return tmp.upper();
		} else {
			tmp = x + interval(1.);
			if (round == -1) {
				return log_point(tmp.lower(), -1);
			} else {
				return log_point(tmp.upper(), 1);
			}
		}
	}

	friend interval log1p(const interval& I) {
		if (I.inf < -1.) {
			throw std::domain_error("interval: log of negative value");
		}
		return interval(log1p_point(I.lower(), -1), log1p_point(I.upper(), 1));
	}

	static interval sin_origin(const interval& I) {
		interval r, y;
		int i;

		r = 0.;
		y = 1.;
		for (i=1;  ; i++) {
			y *= I;
			y /= (double)i;
			if (mag(y) < std::numeric_limits<double>::epsilon()) {
				r += y * interval(-1., 1.);
				break;
			} else {
				if (i % 2 != 0) {
					if (i % 4 == 1) {
						r += y;
					} else {
						r -= y;
					}
				}
			}
		}

		return r;
	}

	static interval cos_origin(const interval& I) {
		interval r, y;
		int i;

		r = 1.;
		y = 1.;
		for (i=1;  ; i++) {
			y *= I;
			y /= (double)i;
			if (mag(y) < std::numeric_limits<double>::epsilon()) {
				r += y * interval(-1., 1.);
				break;
			} else {
				if (i % 2 == 0) {
					if (i % 4 == 0) {
						r += y;
					} else {
						r -= y;
					}
				}
			}
		}

		return r;
	}

	static interval sin_point(const interval& I) {
		const interval pi = interval::pi();
		const double mpi = mid(pi);

		if (I.lower() >= mpi) {
			return sin_point(I - pi * 2.);
		}

		if (I.upper() <= -mpi * 3. / 4.) {
			return -sin_origin(I + pi);
		}
		if (I.upper() <= -mpi * 0.5) {
			return -cos_origin(-pi * 0.5 - I);
		}
		if (I.upper() <= -mpi * 0.25) {
			return -cos_origin(I + pi * 0.5);
		}
		if (I.upper() <= 0.) {
			return -sin_origin(-I);
		}
		if (I.upper() <= mpi * 0.25) {
			return sin_origin(I);
		}
		if (I.upper() <= mpi * 0.5) {
			return cos_origin(pi * 0.5 - I);
		}
		if (I.upper() <= mpi * 3. / 4.) {
			return cos_origin(I - pi * 0.5);
		}
		return sin_origin(pi - I);
	}

	static interval cos_point(const interval& I) {
		const interval pi = interval::pi();
		const double mpi = mid(pi);

		if (I.lower() >= mpi) {
			return cos_point(I - pi * 2.);
		}

		if (I.upper() <= -mpi * 3. / 4.) {
			return -cos_origin(I + pi);
		}
		if (I.upper() <= -mpi * 0.5) {
			return -sin_origin(-pi * 0.5 - I);
		}
		if (I.upper() <= -mpi * 0.25) {
			return sin_origin(I + pi * 0.5);
		}
		if (I.upper() <= 0.) {
			return cos_origin(-I);
		}
		if (I.upper() <= mpi * 0.25) {
			return cos_origin(I);
		}
		if (I.upper() <= mpi * 0.5) {
			return sin_origin(pi * 0.5 - I);
		}
		if (I.upper() <= mpi * 3. / 4.) {
			return -sin_origin(I - pi * 0.5);
		}
		return -cos_origin(pi - I);
	}

	friend interval sin(const interval& I) {
		const interval pi = interval::pi();
		const interval pi2 = pi * 2.;

		double n;
		interval r, I2;

		using std::abs;
		if (abs(I.lower()) == std::numeric_limits<double>::infinity()) {
			return interval(-1. , 1.);
		}
		if (abs(I.upper()) == std::numeric_limits<double>::infinity()) {
			return interval(-1. , 1.);
		}

		I2 = I;
		while (I2.lower() <= -pi.upper() || I2.lower() >= pi.upper()) {
			using std::floor;
			n = floor((I2.lower() / pi2.lower()) + 0.5);
			I2 -= n * pi2;
		}

		if ((interval(I2.upper()) - I2.lower()).lower() >= pi2.upper()) {
			return interval(-1., 1.);
		}

		r = hull(sin_point(interval(I2.lower())), sin_point(interval(I2.upper())));

		if (overlap(pi * 0.5, I2)) {
			r = hull(r, 1.);
		}

		if (overlap(pi * 2.5, I2)) {
			r = hull(r, 1.);
		}

		if (overlap(-pi * 0.5, I2)) {
			r = hull(r, -1.);
		}

		if (overlap(pi * 1.5, I2)) {
			r = hull(r, -1.);
		}

		return intersect(r, interval(-1., 1.));
	}

	friend interval cos(const interval& I) {
		const interval pi = interval::pi();
		const interval pi2 = pi * 2.;

		double n;
		interval r, I2;

		using std::abs;
		if (abs(I.lower()) == std::numeric_limits<double>::infinity()) {
			return interval(-1. , 1.);
		}
		if (abs(I.upper()) == std::numeric_limits<double>::infinity()) {
			return interval(-1. , 1.);
		}

		I2 = I;
		while (I2.lower() <= -pi.upper() || I2.lower() >= pi.upper()) {
			using std::floor;
			n = floor((I2.lower() / pi2.lower()) + 0.5);
			I2 -= n * pi2;
		}

		if ((interval(I2.upper()) - I2.lower()).lower() >= pi2.upper()) {
			return interval(-1., 1.);
		}

		r = hull(cos_point(interval(I2.lower())), cos_point(interval(I2.upper())));

		if (in(0., I2)) {
			r = hull(r, 1.);
		}

		if (overlap(pi2, I2)) {
			r = hull(r, 1.);
		}

		if (overlap(-pi, I2)) {
			r = hull(r, -1.);
		}

		if (overlap(pi, I2)) {
			r = hull(r, -1.);
		}

		if (overlap(pi * 3., I2)) {
			r = hull(r, -1.);
		}

		return intersect(r, interval(-1., 1.));
	}

	static interval tan_point(const double& x) {
		return sin_point(interval(x)) / cos_point(interval(x));
	}

	friend interval tan(const interval& I) {
		const interval pi = interval::pi();
		const interval pih = pi * 0.5;

		double n;
		interval I2;

		using std::abs;
		if (abs(I.lower()) == std::numeric_limits<double>::infinity()) {
			return whole();
		}
		if (abs(I.upper()) == std::numeric_limits<double>::infinity()) {
			return whole();
		}

		I2 = I;
		while (I2.lower() <= -pih.upper() || I2.lower() >= pih.upper()) {
			using std::floor;
			n = floor((I2.lower() / pi.lower()) + 0.5);
			I2 -= n * pi;
		}

		if (I2.upper() >= pih.upper()) {
			return whole();
		}

		return interval(tan_point(I2.lower()).lower(), tan_point(I2.upper()).upper());
	}

	static interval atan_origin(const interval& I) {
		interval tmp;
		interval r, y;
		int i;

		r = 0.;
		y = 1.;
		for (i=1;  ; i++) {
			y *= I;
			tmp = y * interval(-1., 1.) / (double)i;
			if (mag(tmp) < std::numeric_limits<double>::epsilon()) {
				r += tmp;
				break;
			} else {
				if (i % 2 != 0) {
					if (i % 4 == 1) {
						r += y / (double)i;
					} else {
						r -= y / (double)i;
					}
				}
			}
		}

		return r;
	}

	static interval atan_point(const double& x) {
		const interval pi = interval::pi();

		interval I = interval(x);

		if (x < -(std::sqrt(2.) + 1.)) {
			return -pi * 0.5 - atan_origin(1. / I);
		}
		if (x < -(std::sqrt(2.) - 1.)) {
			return -pi * 0.25 + atan_origin((1. + I)/(1 - I));
		}
		if (x < (std::sqrt(2.) - 1.)) {
			return atan_origin(I);
		}
		if (x < (std::sqrt(2.) + 1.)) {
			return pi * 0.25 + atan_origin((I - 1.)/(I + 1.));
		}
		return pi * 0.5 - atan_origin(1. / I);
	}

	friend interval atan(const interval& I) {
		return interval(atan_point(I.lower()).lower(), atan_point(I.upper()).upper());
	}

	static interval asin_point(const double& x) {
		const interval pih = interval::pi() * 0.5;

		if (x < -1. || x > 1.) {
			throw std::domain_error("interval: asin domain error");
		}

		if (x == 1) return pih;
		if (x == -1) return -pih;
		using std::abs;
		if (abs(x) < std::sqrt(6.) / 3.) {
			return atan(x / sqrt(1. - interval(x) * x));
		} else {
			if (x > 0.) {
				return atan(x / sqrt((1. + interval(x)) * (1. - x)));
			} else {
				return atan(x / sqrt((1. + x) * (1. - interval(x))));
			}
		}
	}

	friend interval asin(const interval& I) {
		return interval(asin_point(I.lower()).lower(), asin_point(I.upper()).upper());
	}

	static interval pih_m_atan_point(const interval& I) {
		const interval pi = interval::pi();

		if (I.lower() < -(std::sqrt(2.) + 1.)) {
			return pi + atan_origin(1. / I);
		}
		if (I.lower() < -(std::sqrt(2.) - 1.)) {
			return pi * 0.75 - atan_origin((1. + I)/(1 - I));
		}
		if (I.lower() < (std::sqrt(2.) - 1.)) {
			return pi * 0.5 - atan_origin(I);
		}
		if (I.lower() < (std::sqrt(2.) + 1.)) {
			return pi * 0.25 - atan_origin((I - 1.)/(I + 1.));
		}
		return atan_origin(1. / I);
	}

	static interval acos_point(const double& x) {
		const interval pi = interval::pi();

		if (x < -1. || x > 1.) {
			throw std::domain_error("interval: acos domain error");
		}
		if (x == 1.) return interval(0.);
		if (x == -1.) return pi;
		using std::abs;
		if (abs(x) < std::sqrt(6.) / 3.) {
			return pih_m_atan_point(x / sqrt(1. - interval(x) * x));
		} else {
			if (x > 0.) {
				return pih_m_atan_point(x / sqrt((1. + interval(x)) * (1. - x)));
			} else {
				return pih_m_atan_point(x / sqrt((1. + x) * (1. - interval(x))));
			}
		}
	}

	friend interval acos(const interval& I) {
		return interval(acos_point(I.upper()).lower(), acos_point(I.lower()).upper());
	}

	static interval atan2_point(const double& y, const double& x) {
		const interval pi = interval::pi();

		interval Ix = interval(x);
		interval Iy = interval(y);

		if (y <= x && y > -x) {
			return atan(Iy / Ix);
		}
		if (y > x && y > -x) {
			return pi * 0.5 - atan(Ix / Iy);
		}
		if (y > x && y <= -x) {
			if (y >= 0.) {
				return pi + atan(Iy / Ix);
			} else {
				return -pi + atan(Iy / Ix);
			}
		}
		return -pi * 0.5 - atan(Ix / Iy);
	}

	friend interval atan2(const interval& Iy, const interval& Ix) {
		const interval pi = interval::pi();

		if (zero_in(Ix)) {
			if (zero_in(Iy)) {
				return interval(-pi.upper(), pi.upper());
			} else {
				if (Iy > 0.) {
					return interval(atan2_point(Iy.lower(), Ix.upper()).lower(), atan2_point(Iy.lower(), Ix.lower()).upper());
				} else {
					return interval(atan2_point(Iy.upper(), Ix.lower()).lower(), atan2_point(Iy.upper(), Ix.upper()).upper());
				}
			}
		} else {
			if (zero_in(Iy)) {
				if (Ix > 0.) {
					return interval(atan2_point(Iy.lower(), Ix.lower()).lower(), atan2_point(Iy.upper(), Ix.lower()).upper());
				} else {
					if (Iy.lower() < 0)
						return interval(atan2_point(Iy.upper(), Ix.upper()).lower(), (pi * 2. + atan2_point(Iy.lower(), Ix.upper())).upper());
					else {
						return interval(atan2_point(Iy.upper(), Ix.upper()).lower(), (atan2_point(Iy.lower(), Ix.upper())).upper());
					}
				}
			} else {
				if (Ix > 0.) {
					if (Iy > 0.) {
						return interval(atan2_point(Iy.lower(), Ix.upper()).lower(), atan2_point(Iy.upper(), Ix.lower()).upper());
					} else {
						return interval(atan2_point(Iy.lower(), Ix.lower()).lower(), atan2_point(Iy.upper(), Ix.upper()).upper());
					}
				} else {
					if (Iy > 0.) {
						return interval(atan2_point(Iy.upper(), Ix.upper()).lower(), atan2_point(Iy.lower(), Ix.lower()).upper());
					} else {
						return interval(atan2_point(Iy.upper(), Ix.lower()).lower(), atan2_point(Iy.lower(), Ix.upper()).upper());
					}
				}
			}
		}
	}

	static interval sinh_origin(const double& x) {
		// cosh(0.5) = (exp(0.5)+exp(-0.5))/2
		const interval exph = sqrt(e());
		const interval coshh = (exph + 1. / exph) * 0.5;

		interval tmp;
		interval r, y;
		int i;

		r = 0.;
		y = 1.;
		for (i=1;  ; i++) {
			y *= x;
			y /= (double)i;
			tmp = y * interval(-coshh.upper(), coshh.upper());
			if (mag(tmp) < std::numeric_limits<double>::epsilon()) {
				r += tmp;
				break;
			} else {
				if (i % 2 != 0) {
					r += y;
				}
			}
		}

		return r;
	}

	static interval sinh_point(const double& x) {
		interval tmp;
		if (x >= -0.5 && x <= 0.5) {
			return sinh_origin(x);
		} else if (x > 0.) {
			tmp = exp_point(x);
			return (tmp - 1. / tmp) * 0.5;
		} else {
			tmp = exp_point(-x);
			return (1. / tmp - tmp) * 0.5;
		}
	}

	friend interval sinh(const interval& I) {
		return interval(sinh_point(I.lower()).lower(), sinh_point(I.upper()).upper());
	}

	static interval cosh_point(const double& x) {
		interval tmp;
		if (x >= 0.) {
			tmp = exp_point(x);
			return (tmp + 1. / tmp) * 0.5;
		} else {
			tmp = exp_point(-x);
			return (1. / tmp + tmp) * 0.5;
		}
	}

	friend interval cosh(const interval& I) {
		interval r;
		r = hull(cosh_point(I.lower()), cosh_point(I.upper()));
		if (zero_in(I)) {
			r = hull(r, 1.);
		}
		return r;
	}

	static interval tanh_point(const double& x) {
		if (x > 0.5) {
			return 1. - 2. / (1. + exp_point(2. * x));
		} else if (x < -0.5) {
			return 2. / (1. + exp_point(-2. * x)) - 1.;
		} else {
			return sinh_origin(x) / cosh_point(x);
		}
	}

	friend interval tanh(const interval& I) {
		return interval(tanh_point(I.lower()).lower(), tanh_point(I.upper()).upper());
	}

	static interval asinh_point(const double& x) {
		if (x < -0.5) {
			return -log(-x + sqrt(1. + interval(x) * x));
		} else if (x <= 0.5) {
			return log1p((1. + x / (1. + sqrt(1. + interval(x) * x))) * x);
		} else {
			return log(x + sqrt(1. + interval(x) * x));
		}
	}

	friend interval asinh(const interval& I) {
		return interval(asinh_point(I.lower()).lower(), asinh_point(I.upper()).upper());
	}

	static interval acosh_point(const double& x) {
		if (x < 1.) {
			throw std::domain_error("interval: acosh domain error");
		}
		if (x == 1.) {
			return interval(0.);
		} else if (x <= 1.5) {
			interval y(x - 1.);
			return log1p(y + sqrt(y * (interval(x) + 1.)));
		} else {
			return log(x + sqrt(interval(x) * x - 1.));
		}
	}

	friend interval acosh(const interval& I) {
		return interval(acosh_point(I.lower()).lower(), acosh_point(I.upper()).upper());
	}

	static interval atanh_point(const double& x) {
		if (x < -1. || x > 1.) {
			throw std::domain_error("interval: atanh domain error");
		}
		if (x == -1.) return interval(-std::numeric_limits<double>::infinity(), -(std::numeric_limits<double>::max)());
		if (x == 1.) return interval((std::numeric_limits<double>::max)(), std::numeric_limits<double>::infinity());
		if (x < -0.5) {
			return 0.5 * log((1. + x) / (1. - interval(x)));
		} else if (x <= 0.5) {
			return 0.5 * log1p(2. * x / (1. - interval(x)));
		} else {
			return 0.5 * log((1. + interval(x)) / (1. - x));
		}
	}

	friend interval atanh(const interval& I) {
		return interval(atanh_point(I.lower()).lower(), atanh_point(I.upper()).upper());
	}

	static interval pi() {
		// static is used so that string is evaluated only "one time"
		static const interval tmp(
			"3.1415926535897932384626433832795028841971693993751",
			"3.1415926535897932384626433832795028841971693993752"
		);
		return tmp;
	}

	static interval e() {
		static const interval tmp(
			"2.7182818284590452353602874713526624977572470936999",
			"2.7182818284590452353602874713526624977572470937000"
		);
		return tmp;
	}

	static interval ln2() {
		static const interval tmp(
			"0.69314718055994530941723212145817656807550013436025",
			"0.69314718055994530941723212145817656807550013436026"
		);
		return tmp;
	}

	static double add_up(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_UP, _MCW_RC);
		#else
		fesetround(FE_UPWARD);
		#endif
		r = x1 + y1;
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double add_down(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_DOWN, _MCW_RC);
		#else
		fesetround(FE_DOWNWARD);
		#endif
		r = x1 + y1;
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double sub_up(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_UP, _MCW_RC);
		#else
		fesetround(FE_UPWARD);
		#endif
		r = x1 - y1;
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double sub_down(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_DOWN, _MCW_RC);
		#else
		fesetround(FE_DOWNWARD);
		#endif
		r = x1 - y1;
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double mul_up(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_UP, _MCW_RC);
		#else
		fesetround(FE_UPWARD);
		#endif
		r = x1 * y1;
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double mul_down(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_DOWN, _MCW_RC);
		#else
		fesetround(FE_DOWNWARD);
		#endif
		r = x1 * y1;
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double div_up(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_UP, _MCW_RC);
		#else
		fesetround(FE_UPWARD);
		#endif
		r = x1 / y1;
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double div_down(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_DOWN, _MCW_RC);
		#else
		fesetround(FE_DOWNWARD);
		#endif
		r = x1 / y1;
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double sqrt_up(const double& x) {
		volatile double r, x1 = x;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_UP, _MCW_RC);
		#else
		fesetround(FE_UPWARD);
		#endif
		r = sqrt(x1);
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static double sqrt_down(const double& x) {
		volatile double r, x1 = x;
		#if defined(_MSC_VER)
		unsigned int cw = 0;
		_controlfp_s(&cw, _RC_DOWN, _MCW_RC);
		#else
		fesetround(FE_DOWNWARD);
		#endif
		r = sqrt(x1);
		#if defined(_MSC_VER)
		_controlfp_s(&cw, _RC_NEAR, _MCW_RC);
		#else
		fesetround(FE_TONEAREST);
		#endif
		return r;
	}

	static void print_up(const double& x, std::ostream& s) {
		char format;
		if (s.flags() & s.scientific) {
			if (s.flags() & s.fixed) {
				format = 'g';
			} else {
				format = 'e';
			}
		} else {
			if (s.flags() & s.fixed) {
				format = 'f';
			} else {
				format = 'g';
			}
		}
		s << conv_double::dtostring(x, s.precision(), format, 1);
	}

	static void print_down(const double& x, std::ostream& s) {
		char format;
		if (s.flags() & s.scientific) {
			if (s.flags() & s.fixed) {
				format = 'g';
			} else {
				format = 'e';
			}
		} else {
			if (s.flags() & s.fixed) {
				format = 'f';
			} else {
				format = 'g';
			}
		}
		s << conv_double::dtostring(x, s.precision(), format, -1);
	}

	static double fromstring_up(const std::string& s) {
		return conv_double::stringtod(s, 1);
	}

	static double fromstring_down(const std::string& s) {
		return conv_double::stringtod(s, -1);
	}
};

#endif // INTERVAL_HPP