Chris@16: //  Copyright John Maddock 2005-2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_TOOLS_PRECISION_INCLUDED
Chris@16: #define BOOST_MATH_TOOLS_PRECISION_INCLUDED
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/limits.hpp>
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <boost/static_assert.hpp>
Chris@16: #include <boost/mpl/int.hpp>
Chris@16: #include <boost/mpl/bool.hpp>
Chris@16: #include <boost/mpl/if.hpp>
Chris@16: #include <boost/math/policies/policy.hpp>
Chris@16: 
Chris@16: // These two are for LDBL_MAN_DIG:
Chris@16: #include <limits.h>
Chris@16: #include <math.h>
Chris@16: 
Chris@16: namespace boost{ namespace math
Chris@16: {
Chris@16: namespace tools
Chris@16: {
Chris@16: // If T is not specialized, the functions digits, max_value and min_value,
Chris@16: // all get synthesised automatically from std::numeric_limits.
Chris@16: // However, if numeric_limits is not specialised for type RealType,
Chris@16: // for example with NTL::RR type, then you will get a compiler error
Chris@16: // when code tries to use these functions, unless you explicitly specialise them.
Chris@16: 
Chris@16: // For example if the precision of RealType varies at runtime,
Chris@16: // then numeric_limits support may not be appropriate,
Chris@16: // see boost/math/tools/ntl.hpp  for examples like
Chris@16: // template <> NTL::RR max_value<NTL::RR> ...
Chris@16: // See  Conceptual Requirements for Real Number Types.
Chris@16: 
Chris@16: template <class T>
Chris@16: inline int digits(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
Chris@16: {
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:    BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
Chris@16:    BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::radix == 2 || ::std::numeric_limits<T>::radix == 10);
Chris@16: #else
Chris@16:    BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
Chris@16:    BOOST_ASSERT(::std::numeric_limits<T>::radix == 2 || ::std::numeric_limits<T>::radix == 10);
Chris@16: #endif
Chris@16:    return std::numeric_limits<T>::radix == 2 
Chris@16:       ? std::numeric_limits<T>::digits
Chris@16:       : ((std::numeric_limits<T>::digits + 1) * 1000L) / 301L;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:    BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
Chris@16: #else
Chris@16:    BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
Chris@16: #endif
Chris@16:    return (std::numeric_limits<T>::max)();
Chris@16: } // Also used as a finite 'infinite' value for - and +infinity, for example:
Chris@16: // -max_value<double> = -1.79769e+308, max_value<double> = 1.79769e+308.
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:    BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
Chris@16: #else
Chris@16:    BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
Chris@16: #endif
Chris@16:    return (std::numeric_limits<T>::min)();
Chris@16: }
Chris@16: 
Chris@16: namespace detail{
Chris@16: //
Chris@16: // Logarithmic limits come next, note that although
Chris@16: // we can compute these from the log of the max value
Chris@16: // that is not in general thread safe (if we cache the value)
Chris@16: // so it's better to specialise these:
Chris@16: //
Chris@16: // For type float first:
Chris@16: //
Chris@16: template <class T>
Chris@16: inline T log_max_value(const mpl::int_<128>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16:    return 88.0f;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T log_min_value(const mpl::int_<128>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16:    return -87.0f;
Chris@16: }
Chris@16: //
Chris@16: // Now double:
Chris@16: //
Chris@16: template <class T>
Chris@16: inline T log_max_value(const mpl::int_<1024>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16:    return 709.0;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T log_min_value(const mpl::int_<1024>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16:    return -708.0;
Chris@16: }
Chris@16: //
Chris@16: // 80 and 128-bit long doubles:
Chris@16: //
Chris@16: template <class T>
Chris@16: inline T log_max_value(const mpl::int_<16384>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16:    return 11356.0L;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T log_min_value(const mpl::int_<16384>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16:    return -11355.0L;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T log_max_value(const mpl::int_<0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:    BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
Chris@16: #else
Chris@16:    BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
Chris@16: #endif
Chris@16:    BOOST_MATH_STD_USING
Chris@101: #ifdef __SUNPRO_CC
Chris@101:    static const T m = (std::numeric_limits<T>::max)();
Chris@101:    static const T val = log(m);
Chris@101: #else
Chris@16:    static const T val = log((std::numeric_limits<T>::max)());
Chris@101: #endif
Chris@16:    return val;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T log_min_value(const mpl::int_<0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:    BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
Chris@16: #else
Chris@16:    BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
Chris@16: #endif
Chris@16:    BOOST_MATH_STD_USING
Chris@101: #ifdef __SUNPRO_CC
Chris@101:    static const T m = (std::numeric_limits<T>::min)();
Chris@101:    static const T val = log(m);
Chris@101: #else
Chris@16:    static const T val = log((std::numeric_limits<T>::min)());
Chris@101: #endif
Chris@16:    return val;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T epsilon(const mpl::true_& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16:    return std::numeric_limits<T>::epsilon();
Chris@16: }
Chris@16: 
Chris@101: #if defined(__GNUC__) && ((LDBL_MANT_DIG == 106) || (__LDBL_MANT_DIG__ == 106))
Chris@16: template <>
Chris@16: inline long double epsilon<long double>(const mpl::true_& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(long double))
Chris@16: {
Chris@101:    // numeric_limits on Darwin (and elsewhere) tells lies here:
Chris@101:    // the issue is that long double on a few platforms is
Chris@101:    // really a "double double" which has a non-contiguous
Chris@101:    // mantissa: 53 bits followed by an unspecified number of
Chris@101:    // zero bits, followed by 53 more bits.  Thus the apparent
Chris@101:    // precision of the type varies depending where it's been.
Chris@101:    // Set epsilon to the value that a 106 bit fixed mantissa
Chris@101:    // type would have, as that will give us sensible behaviour everywhere.
Chris@101:    //
Chris@16:    // This static assert fails for some unknown reason, so
Chris@16:    // disabled for now...
Chris@16:    // BOOST_STATIC_ASSERT(std::numeric_limits<long double>::digits == 106);
Chris@16:    return 2.4651903288156618919116517665087e-32L;
Chris@16: }
Chris@16: #endif
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T epsilon(const mpl::false_& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std names
Chris@16:    static const T eps = ldexp(static_cast<T>(1), 1-policies::digits<T, policies::policy<> >());
Chris@16:    return eps;
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(push)
Chris@16: #pragma warning(disable:4309)
Chris@16: #endif
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:    typedef typename mpl::if_c<
Chris@16:       (std::numeric_limits<T>::radix == 2) &&
Chris@16:       (std::numeric_limits<T>::max_exponent == 128
Chris@16:       || std::numeric_limits<T>::max_exponent == 1024
Chris@16:       || std::numeric_limits<T>::max_exponent == 16384),
Chris@16:       mpl::int_<(std::numeric_limits<T>::max_exponent > INT_MAX ? INT_MAX : static_cast<int>(std::numeric_limits<T>::max_exponent))>,
Chris@16:       mpl::int_<0>
Chris@16:    >::type tag_type;
Chris@16:    BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
Chris@16:    return detail::log_max_value<T>(tag_type());
Chris@16: #else
Chris@16:    BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    static const T val = log((std::numeric_limits<T>::max)());
Chris@16:    return val;
Chris@16: #endif
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))
Chris@16: {
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:    typedef typename mpl::if_c<
Chris@16:       (std::numeric_limits<T>::radix == 2) &&
Chris@16:       (std::numeric_limits<T>::max_exponent == 128
Chris@16:       || std::numeric_limits<T>::max_exponent == 1024
Chris@16:       || std::numeric_limits<T>::max_exponent == 16384),
Chris@16:       mpl::int_<(std::numeric_limits<T>::max_exponent > INT_MAX ? INT_MAX : static_cast<int>(std::numeric_limits<T>::max_exponent))>,
Chris@16:       mpl::int_<0>
Chris@16:    >::type tag_type;
Chris@16: 
Chris@16:    BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized);
Chris@16:    return detail::log_min_value<T>(tag_type());
Chris@16: #else
Chris@16:    BOOST_ASSERT(::std::numeric_limits<T>::is_specialized);
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    static const T val = log((std::numeric_limits<T>::min)());
Chris@16:    return val;
Chris@16: #endif
Chris@16: }
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(pop)
Chris@16: #endif
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
Chris@16: {
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:    return detail::epsilon<T>(mpl::bool_< ::std::numeric_limits<T>::is_specialized>());
Chris@16: #else
Chris@16:    return ::std::numeric_limits<T>::is_specialized ?
Chris@16:       detail::epsilon<T>(mpl::true_()) :
Chris@16:       detail::epsilon<T>(mpl::false_());
Chris@16: #endif
Chris@16: }
Chris@16: 
Chris@16: namespace detail{
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T root_epsilon_imp(const mpl::int_<24>&)
Chris@16: {
Chris@16:    return static_cast<T>(0.00034526698300124390839884978618400831996329879769945L);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T root_epsilon_imp(const T*, const mpl::int_<53>&)
Chris@16: {
Chris@16:    return static_cast<T>(0.1490116119384765625e-7L);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T root_epsilon_imp(const T*, const mpl::int_<64>&)
Chris@16: {
Chris@16:    return static_cast<T>(0.32927225399135962333569506281281311031656150598474e-9L);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T root_epsilon_imp(const T*, const mpl::int_<113>&)
Chris@16: {
Chris@16:    return static_cast<T>(0.1387778780781445675529539585113525390625e-16L);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Tag>
Chris@16: inline T root_epsilon_imp(const T*, const Tag&)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    static const T r_eps = sqrt(tools::epsilon<T>());
Chris@16:    return r_eps;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@101: inline T cbrt_epsilon_imp(const mpl::int_<24>&)
Chris@101: {
Chris@101:    return static_cast<T>(0.0049215666011518482998719164346805794944150447839903L);
Chris@101: }
Chris@101: 
Chris@101: template <class T>
Chris@101: inline T cbrt_epsilon_imp(const T*, const mpl::int_<53>&)
Chris@101: {
Chris@101:    return static_cast<T>(6.05545445239333906078989272793696693569753008995e-6L);
Chris@101: }
Chris@101: 
Chris@101: template <class T>
Chris@101: inline T cbrt_epsilon_imp(const T*, const mpl::int_<64>&)
Chris@101: {
Chris@101:    return static_cast<T>(4.76837158203125e-7L);
Chris@101: }
Chris@101: 
Chris@101: template <class T>
Chris@101: inline T cbrt_epsilon_imp(const T*, const mpl::int_<113>&)
Chris@101: {
Chris@101:    return static_cast<T>(5.7749313854154005630396773604745549542403508090496e-12L);
Chris@101: }
Chris@101: 
Chris@101: template <class T, class Tag>
Chris@101: inline T cbrt_epsilon_imp(const T*, const Tag&)
Chris@101: {
Chris@101:    BOOST_MATH_STD_USING;
Chris@101:    static const T cbrt_eps = pow(tools::epsilon<T>(), T(1) / 3);
Chris@101:    return cbrt_eps;
Chris@101: }
Chris@101: 
Chris@101: template <class T>
Chris@16: inline T forth_root_epsilon_imp(const T*, const mpl::int_<24>&)
Chris@16: {
Chris@16:    return static_cast<T>(0.018581361171917516667460937040007436176452688944747L);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T forth_root_epsilon_imp(const T*, const mpl::int_<53>&)
Chris@16: {
Chris@16:    return static_cast<T>(0.0001220703125L);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T forth_root_epsilon_imp(const T*, const mpl::int_<64>&)
Chris@16: {
Chris@16:    return static_cast<T>(0.18145860519450699870567321328132261891067079047605e-4L);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T forth_root_epsilon_imp(const T*, const mpl::int_<113>&)
Chris@16: {
Chris@16:    return static_cast<T>(0.37252902984619140625e-8L);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Tag>
Chris@16: inline T forth_root_epsilon_imp(const T*, const Tag&)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    static const T r_eps = sqrt(sqrt(tools::epsilon<T>()));
Chris@16:    return r_eps;
Chris@16: }
Chris@16: 
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T root_epsilon()
Chris@16: {
Chris@16:    typedef mpl::int_< (::std::numeric_limits<T>::radix == 2) ? std::numeric_limits<T>::digits : 0> tag_type;
Chris@16:    return detail::root_epsilon_imp(static_cast<T const*>(0), tag_type());
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@101: inline T cbrt_epsilon()
Chris@101: {
Chris@101:    typedef mpl::int_< (::std::numeric_limits<T>::radix == 2) ? std::numeric_limits<T>::digits : 0> tag_type;
Chris@101:    return detail::cbrt_epsilon_imp(static_cast<T const*>(0), tag_type());
Chris@101: }
Chris@101: 
Chris@101: template <class T>
Chris@16: inline T forth_root_epsilon()
Chris@16: {
Chris@16:    typedef mpl::int_< (::std::numeric_limits<T>::radix == 2) ? std::numeric_limits<T>::digits : 0> tag_type;
Chris@16:    return detail::forth_root_epsilon_imp(static_cast<T const*>(0), tag_type());
Chris@16: }
Chris@16: 
Chris@16: } // namespace tools
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_TOOLS_PRECISION_INCLUDED
Chris@16: