math/tools/precision.hpp

0001 //  Copyright John Maddock 2005-2006.
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_TOOLS_PRECISION_INCLUDED
0007 #define BOOST_MATH_TOOLS_PRECISION_INCLUDED
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/config.hpp>
0014 #include <boost/math/tools/assert.hpp>
0015 #include <boost/math/tools/type_traits.hpp>
0016 #include <boost/math/tools/numeric_limits.hpp>
0017 #include <boost/math/policies/policy.hpp>
0018
0019 #ifndef BOOST_MATH_HAS_NVRTC
0020 #include <type_traits>
0021 #include <limits>
0022 #include <climits>
0023 #include <cmath>
0024 #include <cstdint>
0025 #include <cfloat> // LDBL_MANT_DIG
0026 #endif
0027
0028 namespace boost{ namespace math
0029 {
0030 namespace tools
0031 {
0032 // If T is not specialized, the functions digits, max_value and min_value,
0033 // all get synthesised automatically from std::numeric_limits.
0034 // However, if numeric_limits is not specialised for type RealType,
0035 // for example with NTL::RR type, then you will get a compiler error
0036 // when code tries to use these functions, unless you explicitly specialise them.
0037
0038 // For example if the precision of RealType varies at runtime,
0039 // then numeric_limits support may not be appropriate,
0040 // see boost/math/tools/ntl.hpp  for examples like
0041 // template <> NTL::RR max_value<NTL::RR> ...
0042 // See  Conceptual Requirements for Real Number Types.
0043
0044 template <class T>
0045 BOOST_MATH_GPU_ENABLED inline constexpr int digits(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) noexcept
0046 {
0047    static_assert( ::boost::math::numeric_limits<T>::is_specialized, "Type T must be specialized");
0048    static_assert( ::boost::math::numeric_limits<T>::radix == 2 || ::boost::math::numeric_limits<T>::radix == 10, "Type T must have a radix of 2 or 10");
0049
0050    return boost::math::numeric_limits<T>::radix == 2
0051       ? boost::math::numeric_limits<T>::digits
0052       : ((boost::math::numeric_limits<T>::digits + 1) * 1000L) / 301L;
0053 }
0054
0055 template <class T>
0056 BOOST_MATH_GPU_ENABLED inline constexpr T max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))  noexcept(boost::math::is_floating_point<T>::value)
0057 {
0058    static_assert( ::boost::math::numeric_limits<T>::is_specialized, "Type T must be specialized");
0059    return (boost::math::numeric_limits<T>::max)();
0060 } // Also used as a finite 'infinite' value for - and +infinity, for example:
0061 // -max_value<double> = -1.79769e+308, max_value<double> = 1.79769e+308.
0062
0063 template <class T>
0064 BOOST_MATH_GPU_ENABLED inline constexpr T min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)
0065 {
0066    static_assert( ::boost::math::numeric_limits<T>::is_specialized, "Type T must be specialized");
0067
0068    return (boost::math::numeric_limits<T>::min)();
0069 }
0070
0071 namespace detail{
0072 //
0073 // Logarithmic limits come next, note that although
0074 // we can compute these from the log of the max value
0075 // that is not in general thread safe (if we cache the value)
0076 // so it's better to specialise these:
0077 //
0078 // For type float first:
0079 //
0080 template <class T>
0081 BOOST_MATH_GPU_ENABLED constexpr T log_max_value(const boost::math::integral_constant<int, 128>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)
0082 {
0083    return 88.0f;
0084 }
0085
0086 template <class T>
0087 BOOST_MATH_GPU_ENABLED constexpr T log_min_value(const boost::math::integral_constant<int, 128>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)
0088 {
0089    return -87.0f;
0090 }
0091 //
0092 // Now double:
0093 //
0094 template <class T>
0095 BOOST_MATH_GPU_ENABLED constexpr T log_max_value(const boost::math::integral_constant<int, 1024>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)
0096 {
0097    return 709.0;
0098 }
0099
0100 template <class T>
0101 BOOST_MATH_GPU_ENABLED constexpr T log_min_value(const boost::math::integral_constant<int, 1024>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)
0102 {
0103    return -708.0;
0104 }
0105 //
0106 // 80 and 128-bit long doubles:
0107 //
0108 template <class T>
0109 BOOST_MATH_GPU_ENABLED inline constexpr T log_max_value(const boost::math::integral_constant<int, 16384>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)
0110 {
0111    return 11356.0L;
0112 }
0113
0114 template <class T>
0115 BOOST_MATH_GPU_ENABLED inline constexpr T log_min_value(const boost::math::integral_constant<int, 16384>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)
0116 {
0117    return -11355.0L;
0118 }
0119
0120 template <class T>
0121 BOOST_MATH_GPU_ENABLED inline T log_max_value(const boost::math::integral_constant<int, 0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
0122 {
0123    BOOST_MATH_STD_USING
0124 #ifdef __SUNPRO_CC
0125    static const T m = boost::math::tools::max_value<T>();
0126    static const T val = log(m);
0127 #else
0128    static const T val = log(boost::math::tools::max_value<T>());
0129 #endif
0130    return val;
0131 }
0132
0133 template <class T>
0134 BOOST_MATH_GPU_ENABLED inline T log_min_value(const boost::math::integral_constant<int, 0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
0135 {
0136    BOOST_MATH_STD_USING
0137 #ifdef __SUNPRO_CC
0138    static const T m = boost::math::tools::min_value<T>();
0139    static const T val = log(m);
0140 #else
0141    static const T val = log(boost::math::tools::min_value<T>());
0142 #endif
0143    return val;
0144 }
0145
0146 template <class T>
0147 BOOST_MATH_GPU_ENABLED constexpr T epsilon(const boost::math::true_type& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)
0148 {
0149    return boost::math::numeric_limits<T>::epsilon();
0150 }
0151
0152 #if defined(__GNUC__) && ((LDBL_MANT_DIG == 106) || (__LDBL_MANT_DIG__ == 106))
0153 template <>
0154 BOOST_MATH_GPU_ENABLED inline constexpr long double epsilon<long double>(const boost::math::true_type& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(long double)) noexcept(boost::math::is_floating_point<long double>::value)
0155 {
0156    // numeric_limits on Darwin (and elsewhere) tells lies here:
0157    // the issue is that long double on a few platforms is
0158    // really a "double double" which has a non-contiguous
0159    // mantissa: 53 bits followed by an unspecified number of
0160    // zero bits, followed by 53 more bits.  Thus the apparent
0161    // precision of the type varies depending where it's been.
0162    // Set epsilon to the value that a 106 bit fixed mantissa
0163    // type would have, as that will give us sensible behaviour everywhere.
0164    //
0165    // This static assert fails for some unknown reason, so
0166    // disabled for now...
0167    // static_assert(std::numeric_limits<long double>::digits == 106);
0168    return 2.4651903288156618919116517665087e-32L;
0169 }
0170 #endif
0171
0172 template <class T>
0173 BOOST_MATH_GPU_ENABLED inline T epsilon(const boost::math::false_type& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
0174 {
0175    // Note: don't cache result as precision may vary at runtime:
0176    BOOST_MATH_STD_USING  // for ADL of std names
0177    return ldexp(static_cast<T>(1), 1-policies::digits<T, policies::policy<> >());
0178 }
0179
0180 template <class T>
0181 struct log_limit_traits
0182 {
0183    typedef typename boost::math::conditional<
0184       (boost::math::numeric_limits<T>::radix == 2) &&
0185       (
0186          (     boost::math::numeric_limits<T>::max_exponent == 128
0187             || boost::math::numeric_limits<T>::max_exponent == 1024
0188             || boost::math::numeric_limits<T>::max_exponent == 16384
0189          )
0190          && (-boost::math::numeric_limits<T>::min_exponent10 + 1 == boost::math::numeric_limits<T>::max_exponent10)
0191       ),
0192       boost::math::integral_constant<int, (boost::math::numeric_limits<T>::max_exponent > (boost::math::numeric_limits<int>::max)() ? (boost::math::numeric_limits<int>::max)() : static_cast<int>(boost::math::numeric_limits<T>::max_exponent))>,
0193       boost::math::integral_constant<int, 0>
0194    >::type tag_type;
0195    static constexpr bool value = (tag_type::value != 0);
0196    static_assert(::boost::math::numeric_limits<T>::is_specialized || !value, "Type T must be specialized or equal to 0");
0197 };
0198
0199 template <class T, bool b> struct log_limit_noexcept_traits_imp : public log_limit_traits<T> {};
0200 template <class T> struct log_limit_noexcept_traits_imp<T, false> : public boost::math::integral_constant<bool, false> {};
0201
0202 template <class T>
0203 struct log_limit_noexcept_traits : public log_limit_noexcept_traits_imp<T, boost::math::is_floating_point<T>::value> {};
0204
0205 } // namespace detail
0206
0207 #ifdef _MSC_VER
0208 #pragma warning(push)
0209 #pragma warning(disable:4309)
0210 #endif
0211
0212 template <class T>
0213 BOOST_MATH_GPU_ENABLED inline T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
0214 {
0215 #ifndef BOOST_MATH_HAS_NVRTC
0216    #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
0217       return detail::log_max_value<T>(typename detail::log_limit_traits<T>::tag_type());
0218    #else
0219       BOOST_MATH_ASSERT(::boost::math::numeric_limits<T>::is_specialized);
0220       BOOST_MATH_STD_USING
0221       static const T val = log((boost::math::numeric_limits<T>::max)());
0222       return val;
0223    #endif
0224 #else
0225    return log((boost::math::numeric_limits<T>::max)());
0226 #endif
0227 }
0228
0229 template <class T>
0230 BOOST_MATH_GPU_ENABLED inline T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
0231 {
0232 #ifndef BOOST_MATH_HAS_NVRTC
0233    #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
0234       return detail::log_min_value<T>(typename detail::log_limit_traits<T>::tag_type());
0235    #else
0236       BOOST_MATH_ASSERT(::boost::math::numeric_limits<T>::is_specialized);
0237       BOOST_MATH_STD_USING
0238       static const T val = log((boost::math::numeric_limits<T>::min)());
0239       return val;
0240    #endif
0241 #else
0242    return log((boost::math::numeric_limits<T>::min)());
0243 #endif
0244 }
0245
0246 #ifdef _MSC_VER
0247 #pragma warning(pop)
0248 #endif
0249
0250 template <class T>
0251 BOOST_MATH_GPU_ENABLED constexpr T epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) noexcept(boost::math::is_floating_point<T>::value)
0252 {
0253    // NVRTC does not like this dispatching method so we just skip to where we want to go
0254 #ifndef BOOST_MATH_HAS_NVRTC
0255    #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
0256       return detail::epsilon<T>(boost::math::integral_constant<bool, ::boost::math::numeric_limits<T>::is_specialized>());
0257    #else
0258       return ::boost::math::numeric_limits<T>::is_specialized ?
0259          detail::epsilon<T>(boost::math::true_type()) :
0260          detail::epsilon<T>(boost::math::false_type());
0261    #endif
0262 #else
0263    return boost::math::numeric_limits<T>::epsilon();
0264 #endif
0265 }
0266
0267 namespace detail{
0268
0269 template <class T>
0270 BOOST_MATH_GPU_ENABLED inline constexpr T root_epsilon_imp(const boost::math::integral_constant<int, 24>&) noexcept(boost::math::is_floating_point<T>::value)
0271 {
0272    return static_cast<T>(0.00034526698300124390839884978618400831996329879769945L);
0273 }
0274
0275 template <class T>
0276 BOOST_MATH_GPU_ENABLED inline constexpr T root_epsilon_imp(const T*, const boost::math::integral_constant<int, 53>&) noexcept(boost::math::is_floating_point<T>::value)
0277 {
0278    return static_cast<T>(0.1490116119384765625e-7L);
0279 }
0280
0281 template <class T>
0282 BOOST_MATH_GPU_ENABLED inline constexpr T root_epsilon_imp(const T*, const boost::math::integral_constant<int, 64>&) noexcept(boost::math::is_floating_point<T>::value)
0283 {
0284    return static_cast<T>(0.32927225399135962333569506281281311031656150598474e-9L);
0285 }
0286
0287 template <class T>
0288 BOOST_MATH_GPU_ENABLED inline constexpr T root_epsilon_imp(const T*, const boost::math::integral_constant<int, 113>&) noexcept(boost::math::is_floating_point<T>::value)
0289 {
0290    return static_cast<T>(0.1387778780781445675529539585113525390625e-16L);
0291 }
0292
0293 template <class T, class Tag>
0294 BOOST_MATH_GPU_ENABLED inline T root_epsilon_imp(const T*, const Tag&)
0295 {
0296    BOOST_MATH_STD_USING
0297    BOOST_MATH_STATIC_LOCAL_VARIABLE const T r_eps = sqrt(tools::epsilon<T>());
0298    return r_eps;
0299 }
0300
0301 template <class T>
0302 BOOST_MATH_GPU_ENABLED inline T root_epsilon_imp(const T*, const boost::math::integral_constant<int, 0>&)
0303 {
0304    BOOST_MATH_STD_USING
0305    return sqrt(tools::epsilon<T>());
0306 }
0307
0308 template <class T>
0309 BOOST_MATH_GPU_ENABLED inline constexpr T cbrt_epsilon_imp(const boost::math::integral_constant<int, 24>&) noexcept(boost::math::is_floating_point<T>::value)
0310 {
0311    return static_cast<T>(0.0049215666011518482998719164346805794944150447839903L);
0312 }
0313
0314 template <class T>
0315 BOOST_MATH_GPU_ENABLED inline constexpr T cbrt_epsilon_imp(const T*, const boost::math::integral_constant<int, 53>&) noexcept(boost::math::is_floating_point<T>::value)
0316 {
0317    return static_cast<T>(6.05545445239333906078989272793696693569753008995e-6L);
0318 }
0319
0320 template <class T>
0321 BOOST_MATH_GPU_ENABLED inline constexpr T cbrt_epsilon_imp(const T*, const boost::math::integral_constant<int, 64>&) noexcept(boost::math::is_floating_point<T>::value)
0322 {
0323    return static_cast<T>(4.76837158203125e-7L);
0324 }
0325
0326 template <class T>
0327 BOOST_MATH_GPU_ENABLED inline constexpr T cbrt_epsilon_imp(const T*, const boost::math::integral_constant<int, 113>&) noexcept(boost::math::is_floating_point<T>::value)
0328 {
0329    return static_cast<T>(5.7749313854154005630396773604745549542403508090496e-12L);
0330 }
0331
0332 template <class T, class Tag>
0333 BOOST_MATH_GPU_ENABLED inline T cbrt_epsilon_imp(const T*, const Tag&)
0334 {
0335    BOOST_MATH_STD_USING;
0336    static const T cbrt_eps = pow(tools::epsilon<T>(), T(1) / 3);
0337    return cbrt_eps;
0338 }
0339
0340 template <class T>
0341 BOOST_MATH_GPU_ENABLED inline T cbrt_epsilon_imp(const T*, const boost::math::integral_constant<int, 0>&)
0342 {
0343    BOOST_MATH_STD_USING;
0344    return pow(tools::epsilon<T>(), T(1) / 3);
0345 }
0346
0347 template <class T>
0348 BOOST_MATH_GPU_ENABLED inline constexpr T forth_root_epsilon_imp(const T*, const boost::math::integral_constant<int, 24>&) noexcept(boost::math::is_floating_point<T>::value)
0349 {
0350    return static_cast<T>(0.018581361171917516667460937040007436176452688944747L);
0351 }
0352
0353 template <class T>
0354 BOOST_MATH_GPU_ENABLED inline constexpr T forth_root_epsilon_imp(const T*, const boost::math::integral_constant<int, 53>&) noexcept(boost::math::is_floating_point<T>::value)
0355 {
0356    return static_cast<T>(0.0001220703125L);
0357 }
0358
0359 template <class T>
0360 BOOST_MATH_GPU_ENABLED inline constexpr T forth_root_epsilon_imp(const T*, const boost::math::integral_constant<int, 64>&) noexcept(boost::math::is_floating_point<T>::value)
0361 {
0362    return static_cast<T>(0.18145860519450699870567321328132261891067079047605e-4L);
0363 }
0364
0365 template <class T>
0366 BOOST_MATH_GPU_ENABLED inline constexpr T forth_root_epsilon_imp(const T*, const boost::math::integral_constant<int, 113>&) noexcept(boost::math::is_floating_point<T>::value)
0367 {
0368    return static_cast<T>(0.37252902984619140625e-8L);
0369 }
0370
0371 template <class T, class Tag>
0372 BOOST_MATH_GPU_ENABLED inline T forth_root_epsilon_imp(const T*, const Tag&)
0373 {
0374    BOOST_MATH_STD_USING
0375    static const T r_eps = sqrt(sqrt(tools::epsilon<T>()));
0376    return r_eps;
0377 }
0378
0379 template <class T>
0380 BOOST_MATH_GPU_ENABLED inline T forth_root_epsilon_imp(const T*, const boost::math::integral_constant<int, 0>&)
0381 {
0382    BOOST_MATH_STD_USING
0383    return sqrt(sqrt(tools::epsilon<T>()));
0384 }
0385
0386 template <class T>
0387 struct root_epsilon_traits
0388 {
0389    typedef boost::math::integral_constant<int, (::boost::math::numeric_limits<T>::radix == 2) && (::boost::math::numeric_limits<T>::digits != (boost::math::numeric_limits<int>::max)()) ? boost::math::numeric_limits<T>::digits : 0> tag_type;
0390    static constexpr bool has_noexcept = (tag_type::value == 113) || (tag_type::value == 64) || (tag_type::value == 53) || (tag_type::value == 24);
0391 };
0392
0393 }
0394
0395 template <class T>
0396 BOOST_MATH_GPU_ENABLED inline constexpr T root_epsilon() noexcept(boost::math::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)
0397 {
0398    return detail::root_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());
0399 }
0400
0401 template <class T>
0402 BOOST_MATH_GPU_ENABLED inline constexpr T cbrt_epsilon() noexcept(boost::math::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)
0403 {
0404    return detail::cbrt_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());
0405 }
0406
0407 template <class T>
0408 BOOST_MATH_GPU_ENABLED inline constexpr T forth_root_epsilon() noexcept(boost::math::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)
0409 {
0410    return detail::forth_root_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());
0411 }
0412
0413 } // namespace tools
0414 } // namespace math
0415 } // namespace boost
0416
0417 #endif // BOOST_MATH_TOOLS_PRECISION_INCLUDED
0418