Skip to content

Instantly share code, notes, and snippets.

@kuhar

kuhar/CMakeLists.txt

Last active August 29, 2015 14:18
Show Gist options
  • Save kuhar/154974b51e96e07d0b8e to your computer and use it in GitHub Desktop.
Save kuhar/154974b51e96e07d0b8e to your computer and use it in GitHub Desktop.
Download ZIP
Interval Arithmetic
Raw
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(
./
)
Raw
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
Raw
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;
}
Raw
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
Raw
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
Raw
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();
}
Raw
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_ */
Raw
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment