#ifndef __CS_MATH_ERROR_H
#define __CS_MATH_ERROR_H

// Standard Headers
#include<functional>

// CS Headers
#include<cs/core/types.h>
#include<cs/core/binary_traits.h>
#include<cs/math/functions.h>

//
// EPSILON FUNCTIONS
//

namespace cs { namespace math {

template<class T>
struct zero_test : public std::unary_function<T, bool> {
 result_type operator ()(argument_type x) { return equals_zero(x); }
};

template<class T>
struct zero_epsilon_test : public std::unary_function<T, bool> {
 T epsilon;
 zero_epsilon_test() : epsilon(binary_traits<T>::epsilon()) {}
 explicit zero_epsilon_test(T value) : epsilon(value) {}
 result_type operator ()(argument_type x) { return (abs(x) < epsilon); }
};

// simple absolute error equality testing
template<class T>
struct absolute_error_equality_test : public std::unary_function<T, bool> {
 result_type operator ()(argument_type x, argument_type y) {
  return (x == y);
 }
 result_type operator ()(argument_type x, argument_type y, argument_type ex, argument_type ey) {
  return (x == y);
 }
};

// simple absolute error equality testing (binary32)
template<>
struct absolute_error_equality_test<binary32> : public std::unary_function<binary32, bool> {
 result_type operator ()(argument_type x, argument_type y) {
  return (abs(x - y) < binary_traits<binary32>::epsilon());
 }
 result_type operator ()(argument_type x, argument_type y, argument_type ex, argument_type ey) {
  return (abs(x - y) < binary_traits<binary32>::epsilon());
 }
};

// simple absolute error equality testing (binary64)
template<>
struct absolute_error_equality_test<binary64> : public std::unary_function<binary64, bool> {
 result_type operator ()(argument_type x, argument_type y) {
  return (abs(x - y) < binary_traits<binary64>::epsilon());
 }
 result_type operator ()(argument_type x, argument_type y, argument_type ex, argument_type ey) {
  return (abs(x - y) < binary_traits<binary64>::epsilon());
 }
};

// simple absolute error equality testing
template<class T>
struct upper_lower_boundary_equality_test : public std::unary_function<T, bool> {
 result_type operator ()(argument_type x, argument_type y, argument_type x_upper, argument_type x_lower, argument_type y_upper, argument_type y_lower) {
  T diff = x - y;
  T diff_upper = x_upper - y_upper;
  T diff_lower = x_lower - y_lower;
 cout << "diff = " << diff << endl;
 cout << "diff_upper = " << diff_upper << endl;
 cout << "diff_lower = " << diff_lower << endl;
  if(diff_lower < diff_upper) return ((diff <= diff_upper) && (diff_lower <= diff));
  return ((diff <= diff_lower) && (diff_upper <= diff));
 }
};

}}

//
// SIMPLE ERROR FUNCTIONS
// ADDITION, MULTIPLICATION, DIVISION
//
namespace cs { namespace math {

template<class T>
inline T mul_absolute_error(T x, T y, T dx, T dy)
{
 // exact formula for (z + dz) = (x + dx)*(y + dy)
 return dx*(y + dy) + x*dy;
}

template<class T>
inline T mul_absolute_average_error(T x, T y, T dx, T dy)
{
 // assumes dx*dy is small
 return (x*y)*(dx/x + dy/y);
}

template<class T>
inline T mul_absolute_stddev_error(T x, T y, T dx, T dy)
{
 // statistical error propagation
 T a = dx/x;
 T b = dy/y;
 return x*y*sqrt(a*a + b*b);
}

template<class T>
inline T mul_relative_error(T x, T y, T dx, T dy)
{
 // exact formula for dz = dx/x + dy/y + (dx*dy)/(x*y)
 T a = dx/x;
 T b = dy/y;
 return (a + b + a*b);
}

template<class T>
inline T mul_relative_average_error(T x, T y, T dx, T dy)
{
 // assumes dx*dy is small
 return (dx/x + dy/y);
}

template<class T>
inline T mul_relative_stddev_error(T x, T y, T dx, T dy)
{
 // statistical error propagation
 T a = dx/x;
 T b = dy/y;
 return sqrt(a*a + b*b);
}

template<class T>
inline T div_absolute_error(T x, T y, T dx, T dy)
{
 // exact formula for dz = (x + dx)/(y + dy) - z
 return (dx*y - dy*x)/((y + dy)*y);
}

template<class T>
inline T div_relative_error(T x, T y, T dx, T dy)
{
 // exact formula for dz = (x + dx)/(y + dy) - z
 return (dx*y - dy*x)/((y + dy)*x);
}

}}

#endif
