Interval Arithmetic

CMakeLists.txt

cmake_minimum_required(VERSION 2.8.4)
project(ean_logic)

set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++1y -g -DNDEBUG=1")

set(
	SOURCE_FILES example.cpp
	IntervalArithmetic.cpp
	Interval.cpp
	)
add_executable(ean_logic ${SOURCE_FILES})

target_link_libraries(ean_logic mpfr)

include_directories(
	./
	)

ErrorCode.h

/*
* ErrorCode.h
*
*  Created on: 03-04-2015
*  Author: Jakub Kuderski
*/

#pragma once

#include <stdexcept>
#include <string>
#include <cassert>

namespace ean
{
	template<typename T>
	class ErrorCode final
	{
	public:
		ErrorCode(T result, char error)
			: m_result(std::move(result))
			, m_erorrCode(error)
		{}
		ErrorCode(T result)
			: ErrorCode(std::move(result), 0)
		{}

		static ErrorCode createError(char errorCode)
		{
			ErrorCode error;
			error.m_erorrCode = errorCode;
			return error;
		}

		ErrorCode(const ErrorCode&) = default;
		ErrorCode(ErrorCode&&) noexcept = default;

		const T& get() const
		{
			if (isError())
			{
				throw std::runtime_error("Computation error: " + std::to_string(m_erorrCode));
			}
			return m_result;
		}

		T& get()
		{
			if (isError())
			{
				throw std::runtime_error("Computation error: " + std::to_string(m_erorrCode));
			}
			return m_result;
		}

		void setError(char errorCode) noexcept { m_erorrCode = errorCode; }
		char getError() const noexcept { return m_erorrCode; }
		bool isError() const noexcept { return m_erorrCode != 0; }

		operator bool() const noexcept { return !isError(); }

	private:
		ErrorCode() {}

		T m_result;
		char m_erorrCode;
	};

} // namespace ean

example.cpp

#include <iostream>
#include <string>

#include "Number.h"
#include "Interval.h"

using namespace std;

int main()
{
	auto wtw = ean::Interval("1.223");
	auto scnd = wtw;
	cout << (wtw.sin().get() * ean::Interval::pi()).getInverse() << endl;

	ean::Extended ex(5.0l);
	auto res = ex.cos();
	cout << res.isError() << endl;
	//res.setError(2);
	//cout << res.isError() << endl;
	cout << (res.get() * ean::Extended{2.0l}).getInverse();

	wtw = ean::Interval(-1, 1);
	auto div = scnd / wtw; // throws

	return 0;
}

Interval.cpp

/*
* Interval.cpp
*
*  Created on: 02-04-2015
*  Author: Jakub Kuderski
*/

#include "Interval.h"
#include <cmath>
#include <cfloat>

using namespace std;

namespace ean
{
	Interval::Interval(const std::string& value) noexcept
		: m_value(IntervalArithmetic::IntRead(value))
	{}

	Interval::Interval(IntervalArithmetic::interval interval) noexcept
		: m_value(interval)
	{}

	long double Interval::getWidth() const noexcept
	{
		return IntervalArithmetic::IntWidth(m_value);
	}

	Interval& Interval::operator+=(const Interval& rhs) noexcept
	{
		m_value = IntervalArithmetic::IAdd(m_value, rhs.m_value);
		return *this;
	}

	Interval Interval::operator+(const Interval& rhs) const noexcept
	{
		return {IntervalArithmetic::IAdd(m_value, rhs.m_value)};
	}

	Interval& Interval::operator-=(const Interval& rhs) noexcept
	{
		m_value = IntervalArithmetic::ISub(m_value, rhs.m_value);
		return *this;
	}

	Interval Interval::operator-(const Interval& rhs) const noexcept
	{
		return {IntervalArithmetic::ISub(m_value, rhs.m_value)};
	}

	Interval& Interval::operator*=(const Interval& rhs) noexcept
	{
		m_value = IntervalArithmetic::IMul(m_value, rhs.m_value);
		return *this;
	}

	Interval Interval::operator*(const Interval& rhs) const noexcept
	{
		return {IntervalArithmetic::IMul(m_value, rhs.m_value)};
	}

	Interval& Interval::operator/=(const Interval& rhs)
	{
		m_value = IntervalArithmetic::IDiv(m_value, rhs.m_value);
		return *this;
	}

	Interval Interval::operator/(const Interval& rhs) const
	{
		return {IntervalArithmetic::IDiv(m_value, rhs.m_value)};
	}

	bool Interval::operator==(const Interval& rhs) const
	{
		return (std::fabs(m_value.a - rhs.m_value.a) < LDBL_EPSILON)
			&& (std::fabs(m_value.b - rhs.m_value.b) < LDBL_EPSILON);
	}

	Interval& Interval::opposite() noexcept
	{
		m_value = IntervalArithmetic::Opposite(m_value);
		return *this;
	}

	Interval Interval::getOpposite() const noexcept
	{
		return {IntervalArithmetic::Opposite(m_value)};
	}

	Interval& Interval::invert() noexcept
	{
		m_value = IntervalArithmetic::Inverse(m_value);
		return *this;
	}

	Interval Interval::getInverse() const noexcept
	{
		return {IntervalArithmetic::Inverse(m_value)};
	}

	ErrorCode<Interval> Interval::sin() const noexcept
	{
		try
		{
			char st;
			const auto result = IntervalArithmetic::ISin(m_value, st);
			return {{result}, st};
		}
		catch (...)
		{
			return  ErrorCode<Interval>::createError(DivisionByZeroCode);
		}
	}

	ErrorCode<Interval> Interval::cos() const noexcept
	{
		try
		{
			char st;
			const auto result = IntervalArithmetic::ICos(m_value, st);
			return {{result}, st};
		}
		catch (...)
		{
			return  ErrorCode<Interval>::createError(DivisionByZeroCode);
		}
	}

	ErrorCode<Interval> Interval::exp() const noexcept
	{
		try
		{
			char st;
			const auto result = IntervalArithmetic::IExp(m_value, st);
			return {{result}, st};
		}
		catch (...)
		{
			return  ErrorCode<Interval>::createError(DivisionByZeroCode);
		}
	}

	ErrorCode<Interval> Interval::sqrt() const noexcept
	{
		try
		{
			char st;
			const auto result = IntervalArithmetic::ISqr(m_value, st);
			return {{result}, st};
		}
		catch (...)
		{
			return  ErrorCode<Interval>::createError(DivisionByZeroCode);
		}
	}

	std::string Interval::to_string() const
	{
		string left, right;
		IntervalArithmetic::IEndsToStrings(m_value, left, right);
		return "[" + std::move(left) + ", " + std::move(right) + "]";
	}

	std::string Interval::getLeftAsString() const
	{
		string left, right;
		IntervalArithmetic::IEndsToStrings(m_value, left, right);
		return left;
	}

	std::string Interval::getRightAsString() const
	{
		string left, right;
		IntervalArithmetic::IEndsToStrings(m_value, left, right);
		return right;
	}

} // namespace ean

Interval.h

/*
* Interval.h
*
*  Created on: 02-04-2015
*  Author: Jakub Kuderski
*/

#pragma once

#include <ostream>

#include "IntervalArithmetic.h"
#include "ErrorCode.h"

namespace ean
{
	class Interval final
	{
	public:
		static const char DivisionByZeroCode = 3;

		Interval() noexcept = default;
		Interval(long double left, long double right) noexcept
			: m_value{left, right}
		{}
		Interval(const std::string& value) noexcept;

		long double getLeft() const noexcept { return m_value.a;  }
		std::string getLeftAsString() const;
		long double getRight() const noexcept { return m_value.b; }
		std::string getRightAsString() const;
		long double getWidth() const noexcept;

		Interval& operator+=(const Interval& rhs) noexcept;
		Interval operator+(const Interval& rhs) const noexcept;

		Interval& operator-=(const Interval& rhs) noexcept;
		Interval operator-(const Interval& rhs) const noexcept;

		Interval& operator*=(const Interval& rhs) noexcept;
		Interval operator*(const Interval& rhs) const noexcept;

		Interval& operator/=(const Interval& rhs);
		Interval operator/(const Interval& rhs) const;

		bool operator==(const Interval& rhs) const;

		Interval& opposite() noexcept;
		Interval getOpposite() const noexcept;
		Interval& invert() noexcept;
		Interval getInverse() const noexcept;

		ErrorCode<Interval> sin() const noexcept;
		ErrorCode<Interval> cos() const noexcept;
		ErrorCode<Interval> exp() const noexcept;
		ErrorCode<Interval> sqrt() const noexcept;

		static Interval sqrt2() noexcept { return {IntervalArithmetic::ISqr2()}; }
		static Interval sqrt3() noexcept { return {IntervalArithmetic::ISqr3()}; }
		static Interval pi() noexcept { return {IntervalArithmetic::IPi()}; }

		std::string to_string() const;

	private:
		Interval(IntervalArithmetic::interval interval) noexcept;

		IntervalArithmetic::interval m_value = {0.0, 0.0};
	};

	inline std::ostream& operator<<(std::ostream& os, const Interval& value)
	{
		return os << value.to_string();
	}

} // namespace ean

IntervalArithmetic.cpp

/*
* IntervalArithmetic.cpp
*
*  Created on: 20-11-2012
*  Author: Tomasz Hoffmann and Andrzej Marciniak
*  Some utter crap fixed by Jakub Kuderski
*/

#include "IntervalArithmetic.h"

#include <sstream>
#include <iomanip>
#include <stdexcept>
#include <cfenv>
#include <cstdlib>
#include <cstdint>
#include <climits>
#include <cmath>

#include "mpfr.h"

using namespace std;
using namespace IntervalArithmetic;


// store the original rounding mode
const int originalRounding = fegetround();

long double IntervalArithmetic::IntWidth(const interval& x) noexcept
{
	fesetround(FE_UPWARD);
	long double w = x.b - x.a;
	fesetround(FE_TONEAREST);
	return w;
}

IntervalArithmetic::interval IntervalArithmetic::IAdd(const interval& x, const interval& y) noexcept
{
	interval r;
	fesetround(FE_DOWNWARD);
	r.a = x.a + y.a;
	fesetround(FE_UPWARD);
	r.b = x.b + y.b;
	fesetround(FE_TONEAREST);
	return r;
}

IntervalArithmetic::interval IntervalArithmetic::ISub(const interval& x, const interval& y) noexcept
{
	interval r;
	fesetround(FE_DOWNWARD);
	r.a = x.a - y.b;
	fesetround(FE_UPWARD);
	r.b = x.b - y.a;
	fesetround(FE_TONEAREST);
	return r;
}

IntervalArithmetic::interval IntervalArithmetic::IMul(const interval& x, const interval& y) noexcept
{
	interval r;
	long double x1y1, x1y2, x2y1;

	fesetround(FE_DOWNWARD);
	x1y1 = x.a * y.a;
	x1y2 = x.a * y.b;
	x2y1 = x.b * y.a;
	r.a = x.b * y.b;
	if (x2y1 < r.a)
		r.a = x2y1;
	if (x1y2 < r.a)
		r.a = x1y2;
	if (x1y1 < r.a)
		r.a = x1y1;

	fesetround(FE_UPWARD);
	x1y1 = x.a * y.a;
	x1y2 = x.a * y.b;
	x2y1 = x.b * y.a;

	r.b = x.b * y.b;
	if (x2y1 > r.b)
		r.b = x2y1;
	if (x1y2 > r.b)
		r.b = x1y2;
	if (x1y1 > r.b)
		r.b = x1y1;
	fesetround(FE_TONEAREST);
	return r;
}

interval IntervalArithmetic::IDiv(const interval& x, const interval& y)
{
	interval r;
	long double x1y1, x1y2, x2y1, t;

	if ((y.a <= 0) && (y.b >= 0))
	{
		throw runtime_error("Division by an interval containing 0");
	}
	else
	{
		fesetround(FE_DOWNWARD);
		x1y1 = x.a / y.a;
		x1y2 = x.a / y.b;
		x2y1 = x.b / y.a;
		r.a = x.b / y.b;
		t = r.a;
		if (x2y1 < t)
			r.a = x2y1;
		if (x1y2 < t)
			r.a = x1y2;
		if (x1y1 < t)
			r.a = x1y1;

		fesetround(FE_UPWARD);
		x1y1 = x.a / y.a;
		x1y2 = x.a / y.b;
		x2y1 = x.b / y.a;

		r.b = x.b / y.b;
		t = r.b;
		if (x2y1 > t)
			r.b = x2y1;
		if (x1y2 > t)
			r.b = x1y2;
		if (x1y1 > t)
			r.b = x1y1;
	}
	fesetround(FE_TONEAREST);
	return r;
}

long double IntervalArithmetic::DIntWidth(const interval& x) noexcept
{
	long double w1, w2;

	fesetround(FE_UPWARD);
	w1 = x.b - x.a;
	if (w1 < 0)
		w1 = -w1;
	fesetround(FE_DOWNWARD);
	w2 = x.b - x.a;
	if (w2 < 0)
		w2 = -w2;
	fesetround(FE_TONEAREST);
	if (w1 > w2)
		return w1;
	else
		return w2;
}

IntervalArithmetic::interval IntervalArithmetic::Projection(const interval& x) noexcept
{
	interval r;
	r = x;
	if (x.a > x.b)
	{
		r.a = x.b;
		r.b = x.a;
	}
	return r;
}

IntervalArithmetic::interval IntervalArithmetic::Opposite(const interval& x) noexcept
{
	interval r;
	r.a = -x.a;
	r.b = -x.b;
	return r;
}

IntervalArithmetic::interval IntervalArithmetic::Inverse(const interval& x) noexcept
{
	interval z1, z2;

	fesetround(FE_DOWNWARD);
	z1.a = 1 / x.a;
	z2.b = 1 / x.b;
	fesetround(FE_UPWARD);
	z1.b = 1 / x.b;
	z2.a = 1 / x.a;
	fesetround(FE_TONEAREST);
	if (DIntWidth(z1) >= DIntWidth(z2))
		return z1;
	else
		return z2;
}

IntervalArithmetic::interval IntervalArithmetic::DIAdd(const interval& x, const interval& y) noexcept
{
	interval z1, z2;
	if ((x.a <= x.b) && (y.a <= y.b))
	{
		return IAdd(x, y);
	}
	else
	{
		fesetround(FE_DOWNWARD);
		z1.a = x.a + y.a;
		z2.b = x.b + y.b;
		fesetround(FE_UPWARD);
		z1.b = x.b + y.b;
		z2.a = x.a + y.a;
		fesetround(FE_TONEAREST);
		if (DIntWidth(z1) >= DIntWidth(z2))
			return z1;
		else
			return z2;
	}
}

IntervalArithmetic::interval IntervalArithmetic::DISub(const interval& x, const interval& y) noexcept
{
	interval z1, z2;
	if ((x.a <= x.b) && (y.a <= y.b))
	{
		return ISub(x, y);
	}
	else
	{
		fesetround(FE_DOWNWARD);
		z1.a = x.a - y.b;
		z2.b = x.b - y.a;
		fesetround(FE_UPWARD);
		z1.b = x.b - y.a;
		z2.a = x.a - y.b;
		fesetround(FE_TONEAREST);
		if (DIntWidth(z1) >= DIntWidth(z2))
			return z1;
		else
			return z2;
	}
}

IntervalArithmetic::interval IntervalArithmetic::DIMul(const interval& x, const interval& y) noexcept
{
	interval z1, z2, r;
	long double z;
	bool xn, xp, yn, yp, zero;

	if ((x.a <= x.b) && (y.a <= y.b))
		r = IMul(x, y);
	else
	{
		xn = (x.a < 0) && (x.b < 0);
		xp = (x.a > 0) && (x.b > 0);
		yn = (y.a < 0) && (y.b < 0);
		yp = (y.a > 0) && (y.b > 0);
		zero = false;
		// A, B in H-T
		if ((xn || xp) && (yn || yp))
			if (xp && yp)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a * y.a;
				z2.b = x.b * y.b;
				fesetround(FE_UPWARD);
				z1.b = x.b * y.b;
				z2.a = x.a * y.a;
			}
			else if (xp && yn)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b * y.a;
				z2.b = x.a * y.b;
				fesetround(FE_UPWARD);
				z1.b = x.a * y.b;
				z2.a = x.b * y.a;
			}
			else if (xn && yp)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a * y.b;
				z2.b = x.b * y.a;
				fesetround(FE_UPWARD);
				z1.b = x.b * y.a;
				z2.a = x.a * y.b;
			}
			else
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b * y.b;
				z2.b = x.a * y.a;
				fesetround(FE_UPWARD);
				z1.b = x.a * y.a;
				z2.a = x.b * y.b;
			}
			// A in H-T, B in T
		else if ((xn || xp) && (((y.a <= 0) && (y.b >= 0)) || ((y.a >= 0) && (y.b <= 0))))
			if (xp && (y.a <= y.b))
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b * y.a;
				z2.b = x.b * y.b;
				fesetround(FE_UPWARD);
				z1.b = x.b * y.b;
				z2.a = x.b * y.a;
			}
			else if (xp && (y.a > y.b))
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a * y.a;
				z2.b = x.a * y.b;
				fesetround(FE_UPWARD);
				z1.b = x.a * y.b;
				z2.a = x.a * y.a;
			}
			else if (xn && (y.a <= y.b))
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a * y.b;
				z2.b = x.a * y.a;
				fesetround(FE_UPWARD);
				z1.b = x.a * y.a;
				z2.a = x.a * y.b;
			}
			else
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b * y.b;
				z2.b = x.b * y.a;
				fesetround(FE_UPWARD);
				z1.b = x.b * y.a;
				z2.a = x.b * y.b;
			}
			// A in T, B in H-T
		else if ((((x.a <= 0) && (x.b >= 0)) || ((x.a >= 0) && (x.b <= 0))) && (yn || yp))
			if ((x.a <= x.b) && yp)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a * y.b;
				z2.b = x.b * y.b;
				fesetround(FE_UPWARD);
				z1.b = x.b * y.b;
				z2.a = x.a * y.b;
			}
			else if ((x.a <= 0) && yn)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b * y.a;
				z2.b = x.a * y.a;
				fesetround(FE_UPWARD);
				z1.b = x.a * y.a;
				z2.a = x.b * y.a;
			}
			else if ((x.a > x.b) && yp)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a * y.a;
				z2.b = x.b * y.a;
				fesetround(FE_UPWARD);
				z1.b = x.b * y.a;
				z2.a = x.a * y.a;
			}
			else
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b * y.b;
				z2.b = x.a * y.b;
				fesetround(FE_UPWARD);
				z1.b = x.a * y.b;
				z2.a = x.b * y.b;
			}
			// A, B in Z-
		else if ((x.a >= 0) && (x.b <= 0) && (y.a >= 0) && (y.b <= 0))
		{
			fesetround(FE_DOWNWARD);
			z1.a = x.a * y.a;
			z = x.b * y.b;
			if (z1.a < z)
				z1.a = z;
			z2.b = x.a * y.b;
			z = x.b * y.a;
			if (z < z2.b)
				z2.b = z;
			fesetround(FE_UPWARD);
			z1.b = x.a * y.b;
			z = x.b * y.a;
			if (z < z1.b)
				z1.b = z;
			z2.a = x.a * y.a;
			z = x.b * y.b;
			if (z2.a < z)
				z2.a = z;
		}
		// A in Z and B in Z- or A in Z- and B in Z
		else
			zero = true;
		if (zero)
		{
			r.a = 0;
			r.b = 0;
		}
		else if (DIntWidth(z1) >= DIntWidth(z2))
			r = z1;
		else
			r = z2;
	}

	fesetround(FE_TONEAREST);
	return r;
}

interval IntervalArithmetic::DIDiv(const interval& x, const interval& y)
{

	interval z1, z2, r;
	bool xn, xp, yn, yp, zero;

	if ((x.a <= x.b) && (y.a <= y.b))
		r = IDiv(x, y);
	else
	{
		xn = (x.a < 0) && (x.b < 0);
		xp = (x.a > 0) && (x.b > 0);
		yn = (y.a < 0) && (y.b < 0);
		yp = (y.a > 0) && (y.b > 0);
		zero = false;
		// A, B in H-T
		if ((xn || xp) && (yn || yp))
			if (xp && yp)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a / y.b;
				z2.b = x.b / y.a;
				fesetround(FE_UPWARD);
				z1.b = x.b / y.a;
				z2.a = x.a / y.b;
			}
			else if (xp && yn)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b / y.b;
				z2.b = x.a / y.a;
				fesetround(FE_UPWARD);
				z1.b = x.a / y.a;
				z2.a = x.b / y.b;
			}
			else if (xn && yp)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a / y.a;
				z2.b = x.b / y.b;
				fesetround(FE_UPWARD);
				z1.b = x.b / y.b;
				z2.a = x.a / y.a;
			}
			else
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b / y.a;
				z2.b = x.a / y.b;
				fesetround(FE_UPWARD);
				z1.b = x.a / y.b;
				z2.a = x.b / y.a;
			}
			// A in T, B in H-T
		else if (((x.a <= 0) && (x.b >= 0)) || (((x.a >= 0) && (x.b <= 0)) && (yn || yp)))
			if ((x.a <= x.b) && yp)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a / y.a;
				z2.b = x.b / y.a;
				fesetround(FE_UPWARD);
				z1.b = x.b / y.a;
				z2.a = x.a / y.a;
			}
			else if ((x.a <= x.b) && yn)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b / y.b;
				z2.b = x.a / y.b;
				fesetround(FE_UPWARD);
				z1.b = x.a / y.b;
				z2.a = x.b / y.b;
			}
			else if ((x.a > x.b) && yp)
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.a / y.b;
				z2.b = x.b / y.b;
				fesetround(FE_UPWARD);
				z1.b = x.b / y.b;
				z2.a = x.a / y.b;
			}
			else
			{
				fesetround(FE_DOWNWARD);
				z1.a = x.b / y.a;
				z2.b = x.a / y.a;
				fesetround(FE_UPWARD);
				z1.b = x.a / y.a;
				z2.a = x.b / y.a;
			}
		else
			zero = true;
		if (zero)
			throw runtime_error("Division by an interval containing 0.");
		else if (DIntWidth(z1) >= DIntWidth(z2))
			r = z1;
		else
			r = z2;
		fesetround(FE_TONEAREST);
	}
	return r;
}

IntervalArithmetic::interval IntervalArithmetic::IntRead(const std::string& sa) noexcept
{
	interval r;
	mpfr_t rop;
	mpfr_init2(rop, 80);
	mpfr_set_str(rop, sa.c_str(), 10, MPFR_RNDD);

	long double le = mpfr_get_ld(rop, MPFR_RNDD);
	mpfr_set_str(rop, sa.c_str(), 10, MPFR_RNDU);
	long double re = mpfr_get_ld(rop, MPFR_RNDU);
	fesetround(FE_TONEAREST);

	r.a = le;
	r.b = re;
	return r;
}

long double IntervalArithmetic::LeftRead(const std::string& sa) noexcept
{
	interval int_number;
	int_number = IntRead(sa);
	return int_number.a;
}

long double IntervalArithmetic::RightRead(const std::string& sa) noexcept
{
	interval int_number;
	int_number = IntRead(sa);
	return int_number.b;
}

IntervalArithmetic::interval IntervalArithmetic::ISin(const interval& x, char& st)
{
	bool is_even, finished = false;
	int k;
	interval d, s, w, w1, x2;
	if (x.a > x.b)
		st = 1;
	else
	{
		s = x;
		w = x;
		x2 = IMul(x, x);
		k = 1;
		is_even = true;
		finished = false;
		st = 0;

		do
		{
			d.a = (k + 1) * (k + 2);
			d.b = d.a;
			s = IMul(s, IDiv(x2, d));
			if (is_even)
				w1 = ISub(w, s);
			else
				w1 = IAdd(w, s);
			if ((w.a != 0) && (w.b != 0))
			{
				if ((std::fabs(w.a - w1.a) / std::fabs(w.a) < 1e-18) && (std::fabs(w.b - w1.b) / std::fabs(w.b) < 1e-18))
					finished = true;
			}

			else if ((w.a == 0) && (w.b != 0))
			{
				if ((std::fabs(w.a - w1.a) < 1e-18) && (std::fabs(w.b - w1.b) / std::fabs(w.b) < 1e-18))
					finished = true;
			}

			else if (w.a != 0)
			{
				if ((std::fabs(w.a - w1.a) / std::fabs(w.a) < 1e-18) & (std::fabs(w.b - w1.b) < 1e-18))
					finished = true;
				else if ((std::fabs(w.a - w1.a) < 1e-18) & (std::fabs(w.b - w1.b) < 1e-18))
					finished = true;
			}

			if (finished)
			{
				if (w1.b > 1)
				{
					w1.b = 1;
					if (w1.a > 1)
						w1.a = 1;
				}
				if (w1.a < -1)
				{
					w1.a = -1;
					if (w1.b < -1)
						w1.b = -1;
				}
				return w1;
			}
			else
			{
				w = w1;
				k = k + 2;
				is_even = !is_even;
			}
		} while (!(finished || (k > INT_MAX / 2)));
	}
	if (!finished)
		st = 2;

	interval r;
	r.a = 0;
	r.b = 0;
	return r;
}
IntervalArithmetic::interval IntervalArithmetic::ICos(const interval& x, char& st)
{
	bool is_even, finished = false;
	int k;
	interval d, c, w, w1, x2;
	if (x.a > x.b)
		st = 1;
	else
	{
		c.a = 1;
		c.b = 1;
		w = c;
		x2 = IMul(x, x);
		k = 1;
		is_even = true;
		finished = false;
		st = 0;

		do
		{
			d.a = k * (k + 1);
			d.b = d.a;
			c = IMul(c, IDiv(x2, d));
			if (is_even)
				w1 = ISub(w, c);
			else
				w1 = IAdd(w, c);

			if ((w.a != 0) && (w.b != 0))
			{
				if ((std::fabs(w.a - w1.a) / std::fabs(w.a) < 1e-18) && (std::fabs(w.b - w1.b) / std::fabs(w.b) < 1e-18))
					finished = true;
			}

			else if ((w.a == 0) && (w.b != 0))
			{
				if ((std::fabs(w.a - w1.a) < 1e-18) && (std::fabs(w.b - w1.b) / std::fabs(w.b) < 1e-18))
					finished = true;
			}

			else if (w.a != 0)
			{
				if ((std::fabs(w.a - w1.a) / std::fabs(w.a) < 1e-18) & (std::fabs(w.b - w1.b) < 1e-18))
					finished = true;
				else if ((std::fabs(w.a - w1.a) < 1e-18) & (std::fabs(w.b - w1.b) < 1e-18))
					finished = true;
			}

			if (finished)
			{
				if (w1.b > 1)
				{
					w1.b = 1;
					if (w1.a > 1)
						w1.a = 1;
				}
				if (w1.a < -1)
				{
					w1.a = -1;
					if (w1.b < -1)
						w1.b = -1;
				}
				return w1;
			}
			else
			{
				w = w1;
				k = k + 2;
				is_even = !is_even;
			}
		} while (!(finished || (k > INT_MAX / 2)));
	}
	if (!finished)
		st = 2;

	interval r;
	r.a = 0;
	r.b = 0;
	return r;
}
IntervalArithmetic::interval IntervalArithmetic::IExp(const interval& x, char& st)
{
	bool finished;
	int k;
	interval d, e, w, w1;
	if (x.a > x.b)
		st = 1;
	else
	{
		e.a = 1;
		e.b = 1;
		w = e;
		k = 1;
		finished = false;
		st = 0;
		do
		{
			d.a = k;
			d.b = k;
			e = IMul(e, IDiv(x, d));
			w1 = IAdd(w, e);
			if ((std::fabs(w.a - w1.a) / std::fabs(w.a) < 1e-18) && (std::fabs(w.b - w1.b) / std::fabs(w.b) < 1e-18))
			{
				finished = true;
				return w1;
			}
			else
			{
				w = w1;
				k = k + 1;
			}
		} while (!(finished || (k > INT_MAX / 2)));
		if (!finished)
			st = 2;
	}
	interval r;
	r.a = 0;
	r.b = 0;
	return r;
}

IntervalArithmetic::interval IntervalArithmetic::ISqr(const interval& x, char& st) noexcept
{
	long double minx, maxx;
	interval r;
	r.a = 0;
	r.b = 0;
	if (x.a > x.b)
		st = 1;
	else
	{
		st = 0;
		if ((x.a <= 0) && (x.b >= 0))
			minx = 0;
		else if (x.a > 0)
			minx = x.a;
		else
			minx = x.b;
		if (std::fabs(x.a) > std::fabs(x.b))
			maxx = std::fabs(x.a);
		else
			maxx = std::fabs(x.b);
		fesetround(FE_DOWNWARD);
		r.a = minx * minx;
		fesetround(FE_UPWARD);
		r.b = maxx * maxx;
		fesetround(FE_TONEAREST);
	}
	return r;
}

interval IntervalArithmetic::ISqr2()
{
	string i2;
	interval r;
	i2 = "1.414213562373095048";
	r.a = LeftRead(i2);
	i2 = "1.414213562373095049";
	r.b = RightRead(i2);
	return r;
}

interval IntervalArithmetic::ISqr3()
{
	string i2;
	interval r;
	i2 = "1.732050807568877293";
	r.a = LeftRead(i2);
	i2 = "1.732050807568877294";
	r.b = RightRead(i2);
	return r;
}
interval IntervalArithmetic::IPi()
{
	string i2;
	interval r;
	i2 = "3.141592653589793238";
	r.a = LeftRead(i2);
	i2 = "3.141592653589793239";
	r.b = RightRead(i2);
	return r;
}

void IntervalArithmetic::IEndsToStrings(const interval& i, string& left, string& right)
{
	ostringstream oss;
	oss << setprecision(15) << i.a;
	left = oss.str();
	oss = ostringstream();
	oss << setprecision(15) << i.b;
	right = oss.str();
}

IntervalArithmetic.h

/*
 * IntervalArithmetic.h
 *
 *  Created on: 20-11-2012
 *  Author: Tomasz Hoffmann and Andrzej Marciniak
 *  Some utter crap fixed by Jakub Kuderski
 *
 *  Before you start use module, please install libraries:
 *  	- GNU MPFR ( http://www.mpfr.org/ )
 */

#ifndef INTERVALARITHMETIC_H_
#define INTERVALARITHMETIC_H_

#include <string>

namespace IntervalArithmetic
{
	struct interval
	{
		long double a;
		long double b;
	};

	long double IntWidth(const interval& x) noexcept;
	interval IAdd(const interval& x, const interval& y) noexcept;
	interval ISub(const interval& x, const interval& y) noexcept;
	interval IMul(const interval& x, const interval& y) noexcept;
	interval IDiv(const interval& x, const interval& y);
	long double DIntWidth(const interval& x) noexcept;
	interval Projection(const interval& x) noexcept;
	interval Opposite(const interval& x) noexcept;
	interval Inverse(const interval& x) noexcept;
	interval DIAdd(const interval& x, const interval& y) noexcept;
	interval DISub(const interval& x, const interval& y) noexcept;
	interval DIMul(const interval& x, const interval& y) noexcept;
	interval DIDiv(const interval& x, const interval& y);
	interval IntRead(const std::string& sa) noexcept;
	long double LeftRead(const std::string& sa) noexcept;
	long double RightRead(const std::string& sa) noexcept;
	interval ISin(const interval& x, char& st);
	interval ICos(const interval& x, char& st);
	interval IExp(const interval& x, char& st);
	interval ISqr(const interval& x, char& st) noexcept;
	interval ISqr2();
	interval ISqr3();
	interval IPi();
	void IEndsToStrings(const interval& i, std::string& left, std::string& right);
} // namespace IntervalArithmetic

#endif /* INTERVALARITHMETIC_H_ */

Number.h

/*
* Interval.h
*
*  Created on: 03-04-2015
*  Author: Jakub Kuderski
*/

#pragma once

#include <string>
#include <ostream>
#include <sstream>
#include <iomanip>
#include <cmath>
#include <cfloat>
#include <type_traits>
#include <stdexcept>

#include "ErrorCode.h"

namespace ean
{
	template<typename T>
	class Number final
	{
		static_assert(std::is_floating_point<T>::value, "T must be a floating point type");
	public:
		Number() noexcept = default;
		Number(const std::string& value) noexcept
			: m_value(std::stold(value))
		{}

		explicit Number(T value) noexcept
			: m_value(value)
		{}

		T getWidth() const noexcept { return T(0); }

		Number& operator+=(const Number& rhs) noexcept
		{
			m_value += rhs.m_value;
			return *this;
		}

		Number operator+(const Number& rhs) const noexcept
		{
			Number temp = *this;
			temp += rhs;
			return temp;
		}

		Number& operator-=(const Number& rhs) noexcept
		{
			m_value -= rhs.m_value;
			return *this;
		}
		Number operator-(const Number& rhs) const noexcept
		{
			Number temp = *this;
			temp -= rhs;
			return temp;
		}

		Number& operator*=(const Number& rhs) noexcept
		{
			m_value *= rhs.m_value;
			return *this;
		}

		Number operator*(const Number& rhs) const noexcept
		{
			Number temp = *this;
			temp *= rhs;
			return temp;
		}

		Number& operator/=(const Number& rhs)
		{
			if (rhs.m_value == T(0)) // abs(x-y) < eps?
			{
				throw std::runtime_error("Division by 0");
			}

			m_value /= rhs.m_value;
			return *this;
		}

		Number operator/(const Number& rhs) const
		{
			Number temp = *this;
			temp /= rhs;
			return temp;
		}

		bool operator==(const Number& rhs) const
		{
			return std::fabs(m_value - rhs.m_value) < T(LDBL_EPSILON);
		}

		Number& opposite() noexcept
		{
			m_value = -m_value;
			return *this;
		}

		Number getOpposite() const noexcept
		{
			return Number{-m_value};
		}

		Number& invert() noexcept
		{
			m_value = 1 / m_value;
			return *this;
		}

		Number getInverse() const noexcept
		{
			Number temp = *this;
			temp.invert();
			return temp;
		}

		ErrorCode<Number> sin() const noexcept
		{
			return {Number{std::sin(m_value)}};
		}
		
		ErrorCode<Number> cos() const noexcept
		{
			return {Number{std::cos(m_value)}};
		}
		
		ErrorCode<Number> exp() const noexcept
		{
			return {Number{std::exp(m_value)}};
		}
		
		ErrorCode<Number> sqrt() const noexcept
		{
			return {Number{std::sqrt(m_value)}};
		}

		static Number sqrt2() noexcept { return Number{T(M_SQRT2)}; }
		static Number sqrt3() noexcept { return Number{T(std::sqrt(3.0l))}; }
		static Number pi() noexcept { return Number{T(M_PI)}; }

		std::string to_string() const 
		{
			std::ostringstream oss;
			oss << std::setprecision(15) << m_value;
			return oss.str();
		}

	private:
		T m_value = 0;
	};
	
	template<typename T>
	inline std::ostream& operator<<(std::ostream& os, const Number<T>& value)
	{
		return os << value.to_string();
	}

	using Extended = Number<long double>;

} // namespace ean