ulp.hpp 3.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103
  1. // (C) Copyright John Maddock 2015.
  2. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. #ifndef BOOST_MATH_SPECIAL_ULP_HPP
  6. #define BOOST_MATH_SPECIAL_ULP_HPP
  7. #ifdef _MSC_VER
  8. #pragma once
  9. #endif
  10. #include <boost/math/special_functions/math_fwd.hpp>
  11. #include <boost/math/policies/error_handling.hpp>
  12. #include <boost/math/special_functions/fpclassify.hpp>
  13. #include <boost/math/special_functions/next.hpp>
  14. #include <boost/math/tools/precision.hpp>
  15. namespace boost{ namespace math{ namespace detail{
  16. template <class T, class Policy>
  17. T ulp_imp(const T& val, const std::true_type&, const Policy& pol)
  18. {
  19. BOOST_MATH_STD_USING
  20. int expon;
  21. static const char* function = "ulp<%1%>(%1%)";
  22. int fpclass = (boost::math::fpclassify)(val);
  23. if(fpclass == FP_NAN)
  24. {
  25. return policies::raise_domain_error<T>(function, "Argument must be finite, but got %1%", val, pol);
  26. }
  27. else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
  28. {
  29. return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, nullptr, pol);
  30. }
  31. else if(fpclass == FP_ZERO)
  32. return detail::get_smallest_value<T>();
  33. //
  34. // This code is almost the same as that for float_next, except for negative integers,
  35. // where we preserve the relation ulp(x) == ulp(-x) as does Java:
  36. //
  37. frexp(fabs(val), &expon);
  38. T diff = ldexp(T(1), expon - tools::digits<T>());
  39. if(diff == 0)
  40. diff = detail::get_smallest_value<T>();
  41. return diff;
  42. }
  43. // non-binary version:
  44. template <class T, class Policy>
  45. T ulp_imp(const T& val, const std::false_type&, const Policy& pol)
  46. {
  47. static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
  48. static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
  49. BOOST_MATH_STD_USING
  50. int expon;
  51. static const char* function = "ulp<%1%>(%1%)";
  52. int fpclass = (boost::math::fpclassify)(val);
  53. if(fpclass == FP_NAN)
  54. {
  55. return policies::raise_domain_error<T>(function,"Argument must be finite, but got %1%", val, pol);
  56. }
  57. else if((fpclass == FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
  58. {
  59. return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, nullptr, pol);
  60. }
  61. else if(fpclass == FP_ZERO)
  62. return detail::get_smallest_value<T>();
  63. //
  64. // This code is almost the same as that for float_next, except for negative integers,
  65. // where we preserve the relation ulp(x) == ulp(-x) as does Java:
  66. //
  67. expon = 1 + ilogb(fabs(val));
  68. T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
  69. if(diff == 0)
  70. diff = detail::get_smallest_value<T>();
  71. return diff; // LCOV_EXCL_LINE previous lines are covered so this one must be too.
  72. }
  73. }
  74. template <class T, class Policy>
  75. inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)
  76. {
  77. typedef typename tools::promote_args<T>::type result_type;
  78. 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);
  79. }
  80. template <class T>
  81. inline typename tools::promote_args<T>::type ulp(const T& val)
  82. {
  83. return ulp(val, policies::policy<>());
  84. }
  85. }} // namespaces
  86. #endif // BOOST_MATH_SPECIAL_ULP_HPP