relative_difference.hpp 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134
  1. // (C) Copyright John Maddock 2006, 2015
  2. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. #ifndef BOOST_MATH_RELATIVE_ERROR
  6. #define BOOST_MATH_RELATIVE_ERROR
  7. #include <boost/math/special_functions/fpclassify.hpp>
  8. #include <boost/math/tools/promotion.hpp>
  9. #include <boost/math/tools/precision.hpp>
  10. namespace boost{
  11. namespace math{
  12. template <class T, class U>
  13. typename boost::math::tools::promote_args<T,U>::type relative_difference(const T& arg_a, const U& arg_b)
  14. {
  15. typedef typename boost::math::tools::promote_args<T, U>::type result_type;
  16. result_type a = arg_a;
  17. result_type b = arg_b;
  18. BOOST_MATH_STD_USING
  19. #ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  20. //
  21. // If math.h has no long double support we can't rely
  22. // on the math functions generating exponents outside
  23. // the range of a double:
  24. //
  25. result_type min_val = (std::max)(
  26. tools::min_value<result_type>(),
  27. static_cast<result_type>((std::numeric_limits<double>::min)()));
  28. result_type max_val = (std::min)(
  29. tools::max_value<result_type>(),
  30. static_cast<result_type>((std::numeric_limits<double>::max)()));
  31. #else
  32. result_type min_val = tools::min_value<result_type>();
  33. result_type max_val = tools::max_value<result_type>();
  34. #endif
  35. // Screen out NaN's first, if either value is a NaN then the distance is "infinite":
  36. if((boost::math::isnan)(a) || (boost::math::isnan)(b))
  37. return max_val;
  38. // Screen out infinites:
  39. if(fabs(b) > max_val)
  40. {
  41. if(fabs(a) > max_val)
  42. return (a < 0) == (b < 0) ? 0 : max_val; // one infinity is as good as another!
  43. else
  44. return max_val; // one infinity and one finite value implies infinite difference
  45. }
  46. else if(fabs(a) > max_val)
  47. return max_val; // one infinity and one finite value implies infinite difference
  48. //
  49. // If the values have different signs, treat as infinite difference:
  50. //
  51. if(((a < 0) != (b < 0)) && (a != 0) && (b != 0))
  52. return max_val;
  53. a = fabs(a);
  54. b = fabs(b);
  55. //
  56. // Now deal with zero's, if one value is zero (or denorm) then treat it the same as
  57. // min_val for the purposes of the calculation that follows:
  58. //
  59. if(a < min_val)
  60. a = min_val;
  61. if(b < min_val)
  62. b = min_val;
  63. return (std::max)(fabs((a - b) / a), fabs((a - b) / b));
  64. }
  65. #if (defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)) && (LDBL_MAX_EXP <= DBL_MAX_EXP)
  66. template <>
  67. inline boost::math::tools::promote_args<double, double>::type relative_difference(const double& arg_a, const double& arg_b)
  68. {
  69. BOOST_MATH_STD_USING
  70. double a = arg_a;
  71. double b = arg_b;
  72. //
  73. // On Mac OS X we evaluate "double" functions at "long double" precision,
  74. // but "long double" actually has a very slightly narrower range than "double"!
  75. // Therefore use the range of "long double" as our limits since results outside
  76. // that range may have been truncated to 0 or INF:
  77. //
  78. double min_val = (std::max)((double)tools::min_value<long double>(), tools::min_value<double>());
  79. double max_val = (std::min)((double)tools::max_value<long double>(), tools::max_value<double>());
  80. // Screen out NaN's first, if either value is a NaN then the distance is "infinite":
  81. if((boost::math::isnan)(a) || (boost::math::isnan)(b))
  82. return max_val;
  83. // Screen out infinites:
  84. if(fabs(b) > max_val)
  85. {
  86. if(fabs(a) > max_val)
  87. return 0; // one infinity is as good as another!
  88. else
  89. return max_val; // one infinity and one finite value implies infinite difference
  90. }
  91. else if(fabs(a) > max_val)
  92. return max_val; // one infinity and one finite value implies infinite difference
  93. //
  94. // If the values have different signs, treat as infinite difference:
  95. //
  96. if(((a < 0) != (b < 0)) && (a != 0) && (b != 0))
  97. return max_val;
  98. a = fabs(a);
  99. b = fabs(b);
  100. //
  101. // Now deal with zero's, if one value is zero (or denorm) then treat it the same as
  102. // min_val for the purposes of the calculation that follows:
  103. //
  104. if(a < min_val)
  105. a = min_val;
  106. if(b < min_val)
  107. b = min_val;
  108. return (std::max)(fabs((a - b) / a), fabs((a - b) / b));
  109. }
  110. #endif
  111. template <class T, class U>
  112. inline typename boost::math::tools::promote_args<T, U>::type epsilon_difference(const T& arg_a, const U& arg_b)
  113. {
  114. typedef typename boost::math::tools::promote_args<T, U>::type result_type;
  115. result_type r = relative_difference(arg_a, arg_b);
  116. if(tools::max_value<result_type>() * boost::math::tools::epsilon<result_type>() < r)
  117. return tools::max_value<result_type>();
  118. return r / boost::math::tools::epsilon<result_type>();
  119. }
  120. } // namespace math
  121. } // namespace boost
  122. #endif