103 lines · plain
1// (C) Copyright John Maddock 2015.2// 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_SPECIAL_ULP_HPP7#define BOOST_MATH_SPECIAL_ULP_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/special_functions/math_fwd.hpp>14#include <boost/math/policies/error_handling.hpp>15#include <boost/math/special_functions/fpclassify.hpp>16#include <boost/math/special_functions/next.hpp>17#include <boost/math/tools/precision.hpp>18 19namespace boost{ namespace math{ namespace detail{20 21template <class T, class Policy>22T ulp_imp(const T& val, const std::true_type&, const Policy& pol)23{24 BOOST_MATH_STD_USING25 int expon;26 static const char* function = "ulp<%1%>(%1%)";27 28 int fpclass = (boost::math::fpclassify)(val);29 30 if(fpclass == FP_NAN)31 {32 return policies::raise_domain_error<T>(function, "Argument must be finite, but got %1%", val, pol);33 }34 else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))35 {36 return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, nullptr, pol);37 }38 else if(fpclass == FP_ZERO)39 return detail::get_smallest_value<T>();40 //41 // This code is almost the same as that for float_next, except for negative integers,42 // where we preserve the relation ulp(x) == ulp(-x) as does Java:43 //44 frexp(fabs(val), &expon);45 T diff = ldexp(T(1), expon - tools::digits<T>());46 if(diff == 0)47 diff = detail::get_smallest_value<T>();48 return diff;49}50// non-binary version:51template <class T, class Policy>52T ulp_imp(const T& val, const std::false_type&, const Policy& pol)53{54 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");55 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");56 BOOST_MATH_STD_USING57 int expon;58 static const char* function = "ulp<%1%>(%1%)";59 60 int fpclass = (boost::math::fpclassify)(val);61 62 if(fpclass == FP_NAN)63 {64 return policies::raise_domain_error<T>(function,"Argument must be finite, but got %1%", val, pol);65 }66 else if((fpclass == FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))67 {68 return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, nullptr, pol);69 }70 else if(fpclass == FP_ZERO)71 return detail::get_smallest_value<T>();72 //73 // This code is almost the same as that for float_next, except for negative integers,74 // where we preserve the relation ulp(x) == ulp(-x) as does Java:75 //76 expon = 1 + ilogb(fabs(val));77 T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);78 if(diff == 0)79 diff = detail::get_smallest_value<T>();80 return diff; // LCOV_EXCL_LINE previous lines are covered so this one must be too.81}82 83}84 85template <class T, class Policy>86inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)87{88 typedef typename tools::promote_args<T>::type result_type;89 return detail::ulp_imp(static_cast<result_type>(val), std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);90}91 92template <class T>93inline typename tools::promote_args<T>::type ulp(const T& val)94{95 return ulp(val, policies::policy<>());96}97 98 99}} // namespaces100 101#endif // BOOST_MATH_SPECIAL_ULP_HPP102 103