math/special_functions/ulp.hpp

0001 //  (C) Copyright John Maddock 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_SPECIAL_ULP_HPP
0007 #define BOOST_MATH_SPECIAL_ULP_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/special_functions/math_fwd.hpp>
0014 #include <boost/math/policies/error_handling.hpp>
0015 #include <boost/math/special_functions/fpclassify.hpp>
0016 #include <boost/math/special_functions/next.hpp>
0017 #include <boost/math/tools/precision.hpp>
0018
0019 namespace boost{ namespace math{ namespace detail{
0020
0021 template <class T, class Policy>
0022 T ulp_imp(const T& val, const std::true_type&, const Policy& pol)
0023 {
0024    BOOST_MATH_STD_USING
0025    int expon;
0026    static const char* function = "ulp<%1%>(%1%)";
0027
0028    int fpclass = (boost::math::fpclassify)(val);
0029
0030    if(fpclass == FP_NAN)
0031    {
0032       return policies::raise_domain_error<T>(
0033          function,
0034          "Argument must be finite, but got %1%", val, pol);
0035    }
0036    else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
0037    {
0038       return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, nullptr, pol);
0039    }
0040    else if(fpclass == FP_ZERO)
0041       return detail::get_smallest_value<T>();
0042    //
0043    // This code is almost the same as that for float_next, except for negative integers,
0044    // where we preserve the relation ulp(x) == ulp(-x) as does Java:
0045    //
0046    frexp(fabs(val), &expon);
0047    T diff = ldexp(T(1), expon - tools::digits<T>());
0048    if(diff == 0)
0049       diff = detail::get_smallest_value<T>();
0050    return diff;
0051 }
0052 // non-binary version:
0053 template <class T, class Policy>
0054 T ulp_imp(const T& val, const std::false_type&, const Policy& pol)
0055 {
0056    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0057    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0058    BOOST_MATH_STD_USING
0059    int expon;
0060    static const char* function = "ulp<%1%>(%1%)";
0061
0062    int fpclass = (boost::math::fpclassify)(val);
0063
0064    if(fpclass == FP_NAN)
0065    {
0066       return policies::raise_domain_error<T>(
0067          function,
0068          "Argument must be finite, but got %1%", val, pol);
0069    }
0070    else if((fpclass == FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
0071    {
0072       return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, nullptr, pol);
0073    }
0074    else if(fpclass == FP_ZERO)
0075       return detail::get_smallest_value<T>();
0076    //
0077    // This code is almost the same as that for float_next, except for negative integers,
0078    // where we preserve the relation ulp(x) == ulp(-x) as does Java:
0079    //
0080    expon = 1 + ilogb(fabs(val));
0081    T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
0082    if(diff == 0)
0083       diff = detail::get_smallest_value<T>();
0084    return diff;
0085 }
0086
0087 }
0088
0089 template <class T, class Policy>
0090 inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)
0091 {
0092    typedef typename tools::promote_args<T>::type result_type;
0093    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);
0094 }
0095
0096 template <class T>
0097 inline typename tools::promote_args<T>::type ulp(const T& val)
0098 {
0099    return ulp(val, policies::policy<>());
0100 }
0101
0102
0103 }} // namespaces
0104
0105 #endif // BOOST_MATH_SPECIAL_ULP_HPP
0106