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