135 lines · plain
1// (C) Copyright John Maddock 2006, 20152// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_RELATIVE_ERROR7#define BOOST_MATH_RELATIVE_ERROR8 9#include <boost/math/special_functions/fpclassify.hpp>10#include <boost/math/tools/promotion.hpp>11#include <boost/math/tools/precision.hpp>12 13namespace boost{14 namespace math{15 16 template <class T, class U>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_USING23#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS24 //25 // If math.h has no long double support we can't rely26 // on the math functions generating exponents outside27 // 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#else36 result_type min_val = tools::min_value<result_type>();37 result_type max_val = tools::max_value<result_type>();38#endif39 // 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 infinities: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 else48 return max_val; // one infinity and one finite value implies infinite difference49 }50 else if(fabs(a) > max_val)51 return max_val; // one infinity and one finite value implies infinite difference52 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 as62 // 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__)) && (LDBL_MAX_EXP <= DBL_MAX_EXP)73 template <>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_USING77 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 outside83 // 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 infinities:92 if(fabs(b) > max_val)93 {94 if(fabs(a) > max_val)95 return 0; // one infinity is as good as another!96 else97 return max_val; // one infinity and one finite value implies infinite difference98 }99 else if(fabs(a) > max_val)100 return max_val; // one infinity and one finite value implies infinite difference101 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 as111 // 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#endif121 122 template <class T, class U>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 math132} // namespace boost133 134#endif135