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