1 files changed, 36 insertions, 2 deletions
diff --git a/boost/math/special_functions/ulp.hpp b/boost/math/special_functions/ulp.hpp
index 3d78a1c7c2..bb0332a7e0 100644
--- a/boost/math/special_functions/ulp.hpp
+++ b/boost/math/special_functions/ulp.hpp
@@ -18,7 +18,7 @@
 namespace boost{ namespace math{ namespace detail{
 
 template <class T, class Policy>
-T ulp_imp(const T& val, const Policy& pol)
+T ulp_imp(const T& val, const mpl::true_&, const Policy& pol)
 {
    BOOST_MATH_STD_USING
    int expon;
@@ -48,6 +48,40 @@ T ulp_imp(const T& val, const Policy& pol)
       diff = detail::get_smallest_value<T>();
    return diff;
 }
+// non-binary version:
+template <class T, class Policy>
+T ulp_imp(const T& val, const mpl::false_&, const Policy& pol)
+{
+   BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
+   BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
+   BOOST_MATH_STD_USING
+   int expon;
+   static const char* function = "ulp<%1%>(%1%)";
+
+   int fpclass = (boost::math::fpclassify)(val);
+
+   if(fpclass == (int)FP_NAN)
+   {
+      return policies::raise_domain_error<T>(
+         function,
+         "Argument must be finite, but got %1%", val, pol);
+   }
+   else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
+   {
+      return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
+   }
+   else if(fpclass == FP_ZERO)
+      return detail::get_smallest_value<T>();
+   //
+   // This code is almost the same as that for float_next, except for negative integers,
+   // where we preserve the relation ulp(x) == ulp(-x) as does Java:
+   //
+   expon = 1 + ilogb(fabs(val));
+   T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
+   if(diff == 0)
+      diff = detail::get_smallest_value<T>();
+   return diff;
+}
 
 }
 
@@ -55,7 +89,7 @@ template <class T, class Policy>
 inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)
 {
    typedef typename tools::promote_args<T>::type result_type;
-   return detail::ulp_imp(static_cast<result_type>(val), pol);
+   return detail::ulp_imp(static_cast<result_type>(val), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
 }
 
 template <class T>