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