ulp.hpp 3.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105
  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. namespace boost{ namespace math{ namespace detail{
  15. template <class T, class Policy>
  16. T ulp_imp(const T& val, const std::true_type&, const Policy& pol)
  17. {
  18. BOOST_MATH_STD_USING
  19. int expon;
  20. static const char* function = "ulp<%1%>(%1%)";
  21. int fpclass = (boost::math::fpclassify)(val);
  22. if(fpclass == (int)FP_NAN)
  23. {
  24. return policies::raise_domain_error<T>(
  25. function,
  26. "Argument must be finite, but got %1%", val, pol);
  27. }
  28. else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
  29. {
  30. return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
  31. }
  32. else if(fpclass == FP_ZERO)
  33. return detail::get_smallest_value<T>();
  34. //
  35. // This code is almost the same as that for float_next, except for negative integers,
  36. // where we preserve the relation ulp(x) == ulp(-x) as does Java:
  37. //
  38. frexp(fabs(val), &expon);
  39. T diff = ldexp(T(1), expon - tools::digits<T>());
  40. if(diff == 0)
  41. diff = detail::get_smallest_value<T>();
  42. return diff;
  43. }
  44. // non-binary version:
  45. template <class T, class Policy>
  46. T ulp_imp(const T& val, const std::false_type&, const Policy& pol)
  47. {
  48. BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
  49. BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
  50. BOOST_MATH_STD_USING
  51. int expon;
  52. static const char* function = "ulp<%1%>(%1%)";
  53. int fpclass = (boost::math::fpclassify)(val);
  54. if(fpclass == (int)FP_NAN)
  55. {
  56. return policies::raise_domain_error<T>(
  57. function,
  58. "Argument must be finite, but got %1%", val, pol);
  59. }
  60. else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
  61. {
  62. return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
  63. }
  64. else if(fpclass == FP_ZERO)
  65. return detail::get_smallest_value<T>();
  66. //
  67. // This code is almost the same as that for float_next, except for negative integers,
  68. // where we preserve the relation ulp(x) == ulp(-x) as does Java:
  69. //
  70. expon = 1 + ilogb(fabs(val));
  71. T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
  72. if(diff == 0)
  73. diff = detail::get_smallest_value<T>();
  74. return diff;
  75. }
  76. }
  77. template <class T, class Policy>
  78. inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)
  79. {
  80. typedef typename tools::promote_args<T>::type result_type;
  81. 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);
  82. }
  83. template <class T>
  84. inline typename tools::promote_args<T>::type ulp(const T& val)
  85. {
  86. return ulp(val, policies::policy<>());
  87. }
  88. }} // namespaces
  89. #endif // BOOST_MATH_SPECIAL_ULP_HPP