math/special_functions/relative_difference.hpp

0001 //  (C) Copyright John Maddock 2006, 2015
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_RELATIVE_ERROR
0007 #define BOOST_MATH_RELATIVE_ERROR
0008
0009 #include <boost/math/special_functions/fpclassify.hpp>
0010 #include <boost/math/tools/promotion.hpp>
0011 #include <boost/math/tools/precision.hpp>
0012
0013 namespace boost{
0014    namespace math{
0015
0016       template <class T, class U>
0017       typename boost::math::tools::promote_args<T,U>::type relative_difference(const T& arg_a, const U& arg_b)
0018       {
0019          typedef typename boost::math::tools::promote_args<T, U>::type result_type;
0020          result_type a = arg_a;
0021          result_type b = arg_b;
0022          BOOST_MATH_STD_USING
0023 #ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0024          //
0025          // If math.h has no long double support we can't rely
0026          // on the math functions generating exponents outside
0027          // the range of a double:
0028          //
0029          result_type min_val = (std::max)(
0030          tools::min_value<result_type>(),
0031          static_cast<result_type>((std::numeric_limits<double>::min)()));
0032          result_type max_val = (std::min)(
0033             tools::max_value<result_type>(),
0034             static_cast<result_type>((std::numeric_limits<double>::max)()));
0035 #else
0036          result_type min_val = tools::min_value<result_type>();
0037          result_type max_val = tools::max_value<result_type>();
0038 #endif
0039          // Screen out NaN's first, if either value is a NaN then the distance is "infinite":
0040          if((boost::math::isnan)(a) || (boost::math::isnan)(b))
0041             return max_val;
0042          // Screen out infinities:
0043          if(fabs(b) > max_val)
0044          {
0045             if(fabs(a) > max_val)
0046                return (a < 0) == (b < 0) ? 0 : max_val;  // one infinity is as good as another!
0047             else
0048                return max_val;  // one infinity and one finite value implies infinite difference
0049          }
0050          else if(fabs(a) > max_val)
0051             return max_val;    // one infinity and one finite value implies infinite difference
0052
0053          //
0054          // If the values have different signs, treat as infinite difference:
0055          //
0056          if(((a < 0) != (b < 0)) && (a != 0) && (b != 0))
0057             return max_val;
0058          a = fabs(a);
0059          b = fabs(b);
0060          //
0061          // Now deal with zero's, if one value is zero (or denorm) then treat it the same as
0062          // min_val for the purposes of the calculation that follows:
0063          //
0064          if(a < min_val)
0065             a = min_val;
0066          if(b < min_val)
0067             b = min_val;
0068
0069          return (std::max)(fabs((a - b) / a), fabs((a - b) / b));
0070       }
0071
0072 #if (defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)) && (LDBL_MAX_EXP <= DBL_MAX_EXP)
0073       template <>
0074       inline boost::math::tools::promote_args<double, double>::type relative_difference(const double& arg_a, const double& arg_b)
0075       {
0076          BOOST_MATH_STD_USING
0077          double a = arg_a;
0078          double b = arg_b;
0079          //
0080          // On Mac OS X we evaluate "double" functions at "long double" precision,
0081          // but "long double" actually has a very slightly narrower range than "double"!
0082          // Therefore use the range of "long double" as our limits since results outside
0083          // that range may have been truncated to 0 or INF:
0084          //
0085          double min_val = (std::max)((double)tools::min_value<long double>(), tools::min_value<double>());
0086          double max_val = (std::min)((double)tools::max_value<long double>(), tools::max_value<double>());
0087
0088          // Screen out NaN's first, if either value is a NaN then the distance is "infinite":
0089          if((boost::math::isnan)(a) || (boost::math::isnan)(b))
0090             return max_val;
0091          // Screen out infinities:
0092          if(fabs(b) > max_val)
0093          {
0094             if(fabs(a) > max_val)
0095                return 0;  // one infinity is as good as another!
0096             else
0097                return max_val;  // one infinity and one finite value implies infinite difference
0098          }
0099          else if(fabs(a) > max_val)
0100             return max_val;    // one infinity and one finite value implies infinite difference
0101
0102          //
0103          // If the values have different signs, treat as infinite difference:
0104          //
0105          if(((a < 0) != (b < 0)) && (a != 0) && (b != 0))
0106             return max_val;
0107          a = fabs(a);
0108          b = fabs(b);
0109          //
0110          // Now deal with zero's, if one value is zero (or denorm) then treat it the same as
0111          // min_val for the purposes of the calculation that follows:
0112          //
0113          if(a < min_val)
0114             a = min_val;
0115          if(b < min_val)
0116             b = min_val;
0117
0118          return (std::max)(fabs((a - b) / a), fabs((a - b) / b));
0119       }
0120 #endif
0121
0122       template <class T, class U>
0123       inline typename boost::math::tools::promote_args<T, U>::type epsilon_difference(const T& arg_a, const U& arg_b)
0124       {
0125          typedef typename boost::math::tools::promote_args<T, U>::type result_type;
0126          result_type r = relative_difference(arg_a, arg_b);
0127          if(tools::max_value<result_type>() * boost::math::tools::epsilon<result_type>() < r)
0128             return tools::max_value<result_type>();
0129          return r / boost::math::tools::epsilon<result_type>();
0130       }
0131 } // namespace math
0132 } // namespace boost
0133
0134 #endif