419 lines · plain
1// Copyright John Maddock 2005-2006.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_PRECISION_INCLUDED7#define BOOST_MATH_TOOLS_PRECISION_INCLUDED8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/tools/config.hpp>14#include <boost/math/tools/assert.hpp>15#include <boost/math/tools/type_traits.hpp>16#include <boost/math/tools/numeric_limits.hpp>17#include <boost/math/policies/policy.hpp>18 19#ifndef BOOST_MATH_HAS_NVRTC20#include <type_traits>21#include <limits>22#include <climits>23#include <cmath>24#include <cstdint>25#include <cfloat> // LDBL_MANT_DIG26#endif27 28namespace boost{ namespace math29{30namespace tools31{32// If T is not specialized, the functions digits, max_value and min_value,33// all get synthesised automatically from std::numeric_limits.34// However, if numeric_limits is not specialised for type RealType,35// for example with NTL::RR type, then you will get a compiler error36// when code tries to use these functions, unless you explicitly specialise them.37 38// For example if the precision of RealType varies at runtime,39// then numeric_limits support may not be appropriate,40// see boost/math/tools/ntl.hpp for examples like41// template <> NTL::RR max_value<NTL::RR> ...42// See Conceptual Requirements for Real Number Types.43 44template <class T>45BOOST_MATH_GPU_ENABLED inline constexpr int digits(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) noexcept46{47 static_assert( ::boost::math::numeric_limits<T>::is_specialized, "Type T must be specialized");48 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");49 50 return boost::math::numeric_limits<T>::radix == 251 ? boost::math::numeric_limits<T>::digits52 : ((boost::math::numeric_limits<T>::digits + 1) * 1000L) / 301L;53}54 55template <class T>56BOOST_MATH_GPU_ENABLED inline constexpr T max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)57{58 static_assert( ::boost::math::numeric_limits<T>::is_specialized, "Type T must be specialized");59 return (boost::math::numeric_limits<T>::max)();60} // Also used as a finite 'infinite' value for - and +infinity, for example:61// -max_value<double> = -1.79769e+308, max_value<double> = 1.79769e+308.62 63template <class T>64BOOST_MATH_GPU_ENABLED inline constexpr T min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(boost::math::is_floating_point<T>::value)65{66 static_assert( ::boost::math::numeric_limits<T>::is_specialized, "Type T must be specialized");67 68 return (boost::math::numeric_limits<T>::min)();69}70 71namespace detail{72//73// Logarithmic limits come next, note that although74// we can compute these from the log of the max value75// that is not in general thread safe (if we cache the value)76// so it's better to specialise these:77//78// For type float first:79//80template <class T>81BOOST_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)82{83 return 88.0f;84}85 86template <class T>87BOOST_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)88{89 return -87.0f;90}91//92// Now double:93//94template <class T>95BOOST_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)96{97 return 709.0;98}99 100template <class T>101BOOST_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)102{103 return -708.0;104}105//106// 80 and 128-bit long doubles:107//108template <class T>109BOOST_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)110{111 return 11356.0L;112}113 114template <class T>115BOOST_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)116{117 return -11355.0L;118}119 120template <class T>121BOOST_MATH_GPU_ENABLED inline T log_max_value(const boost::math::integral_constant<int, 0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))122{123 BOOST_MATH_STD_USING124#ifdef __SUNPRO_CC125 static const T m = boost::math::tools::max_value<T>();126 static const T val = log(m);127#else128 static const T val = log(boost::math::tools::max_value<T>());129#endif130 return val;131}132 133template <class T>134BOOST_MATH_GPU_ENABLED inline T log_min_value(const boost::math::integral_constant<int, 0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))135{136 BOOST_MATH_STD_USING137#ifdef __SUNPRO_CC138 static const T m = boost::math::tools::min_value<T>();139 static const T val = log(m);140#else141 static const T val = log(boost::math::tools::min_value<T>());142#endif143 return val;144}145 146template <class T>147BOOST_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)148{149 return boost::math::numeric_limits<T>::epsilon();150}151 152#if defined(__GNUC__) && ((LDBL_MANT_DIG == 106) || (__LDBL_MANT_DIG__ == 106))153template <>154BOOST_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)155{156 // numeric_limits on Darwin (and elsewhere) tells lies here:157 // the issue is that long double on a few platforms is158 // really a "double double" which has a non-contiguous159 // mantissa: 53 bits followed by an unspecified number of160 // zero bits, followed by 53 more bits. Thus the apparent161 // precision of the type varies depending where it's been.162 // Set epsilon to the value that a 106 bit fixed mantissa163 // type would have, as that will give us sensible behaviour everywhere.164 //165 // This static assert fails for some unknown reason, so166 // disabled for now...167 // static_assert(std::numeric_limits<long double>::digits == 106);168 return 2.4651903288156618919116517665087e-32L;169}170#endif171 172template <class T>173BOOST_MATH_GPU_ENABLED inline T epsilon(const boost::math::false_type& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))174{175 // Note: don't cache result as precision may vary at runtime:176 BOOST_MATH_STD_USING // for ADL of std names177 return ldexp(static_cast<T>(1), 1-policies::digits<T, policies::policy<> >());178}179 180template <class T>181struct log_limit_traits182{183 typedef typename boost::math::conditional<184 (boost::math::numeric_limits<T>::radix == 2) &&185 (186 ( boost::math::numeric_limits<T>::max_exponent == 128187 || boost::math::numeric_limits<T>::max_exponent == 1024188 || boost::math::numeric_limits<T>::max_exponent == 16384189 )190 && (-boost::math::numeric_limits<T>::min_exponent10 + 1 == boost::math::numeric_limits<T>::max_exponent10)191 ),192 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))>,193 boost::math::integral_constant<int, 0>194 >::type tag_type;195 static constexpr bool value = (tag_type::value != 0);196 static_assert(::boost::math::numeric_limits<T>::is_specialized || !value, "Type T must be specialized or equal to 0");197};198 199template <class T, bool b> struct log_limit_noexcept_traits_imp : public log_limit_traits<T> {};200template <class T> struct log_limit_noexcept_traits_imp<T, false> : public boost::math::integral_constant<bool, false> {};201 202template <class T>203struct log_limit_noexcept_traits : public log_limit_noexcept_traits_imp<T, boost::math::is_floating_point<T>::value> {};204 205} // namespace detail206 207#ifdef _MSC_VER208#pragma warning(push)209#pragma warning(disable:4309)210#endif211 212template <class T>213BOOST_MATH_GPU_ENABLED inline T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)214{215#ifndef BOOST_MATH_HAS_NVRTC216 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS217 return detail::log_max_value<T>(typename detail::log_limit_traits<T>::tag_type());218 #else219 BOOST_MATH_ASSERT(::boost::math::numeric_limits<T>::is_specialized);220 BOOST_MATH_STD_USING221 static const T val = log((boost::math::numeric_limits<T>::max)());222 return val;223 #endif224#else225 return log((boost::math::numeric_limits<T>::max)());226#endif227}228 229template <class T>230BOOST_MATH_GPU_ENABLED inline T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)231{232#ifndef BOOST_MATH_HAS_NVRTC233 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS234 return detail::log_min_value<T>(typename detail::log_limit_traits<T>::tag_type());235 #else236 BOOST_MATH_ASSERT(::boost::math::numeric_limits<T>::is_specialized);237 BOOST_MATH_STD_USING238 static const T val = log((boost::math::numeric_limits<T>::min)());239 return val;240 #endif241#else242 return log((boost::math::numeric_limits<T>::min)());243#endif244}245 246#ifdef _MSC_VER247#pragma warning(pop)248#endif249 250template <class T>251BOOST_MATH_GPU_ENABLED constexpr T epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) noexcept(boost::math::is_floating_point<T>::value)252{253 // NVRTC does not like this dispatching method so we just skip to where we want to go254#ifndef BOOST_MATH_HAS_NVRTC255 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS256 return detail::epsilon<T>(boost::math::integral_constant<bool, ::boost::math::numeric_limits<T>::is_specialized>());257 #else258 return ::boost::math::numeric_limits<T>::is_specialized ?259 detail::epsilon<T>(boost::math::true_type()) :260 detail::epsilon<T>(boost::math::false_type());261 #endif262#else263 return boost::math::numeric_limits<T>::epsilon();264#endif265}266 267namespace detail{268 269template <class T>270BOOST_MATH_GPU_ENABLED inline constexpr T root_epsilon_imp(const boost::math::integral_constant<int, 24>&) noexcept(boost::math::is_floating_point<T>::value)271{272 return static_cast<T>(0.00034526698300124390839884978618400831996329879769945L);273}274 275template <class T>276BOOST_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)277{278 return static_cast<T>(0.1490116119384765625e-7L);279}280 281template <class T>282BOOST_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)283{284 return static_cast<T>(0.32927225399135962333569506281281311031656150598474e-9L);285}286 287template <class T>288BOOST_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)289{290 return static_cast<T>(0.1387778780781445675529539585113525390625e-16L);291}292 293template <class T, class Tag>294BOOST_MATH_GPU_ENABLED inline T root_epsilon_imp(const T*, const Tag&)295{296 BOOST_MATH_STD_USING297 BOOST_MATH_STATIC_LOCAL_VARIABLE const T r_eps = sqrt(tools::epsilon<T>());298 return r_eps;299}300 301template <class T>302BOOST_MATH_GPU_ENABLED inline T root_epsilon_imp(const T*, const boost::math::integral_constant<int, 0>&)303{304 BOOST_MATH_STD_USING305 return sqrt(tools::epsilon<T>());306}307 308template <class T>309BOOST_MATH_GPU_ENABLED inline constexpr T cbrt_epsilon_imp(const boost::math::integral_constant<int, 24>&) noexcept(boost::math::is_floating_point<T>::value)310{311 return static_cast<T>(0.0049215666011518482998719164346805794944150447839903L);312}313 314template <class T>315BOOST_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)316{317 return static_cast<T>(6.05545445239333906078989272793696693569753008995e-6L);318}319 320template <class T>321BOOST_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)322{323 return static_cast<T>(4.76837158203125e-7L);324}325 326template <class T>327BOOST_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)328{329 return static_cast<T>(5.7749313854154005630396773604745549542403508090496e-12L);330}331 332template <class T, class Tag>333BOOST_MATH_GPU_ENABLED inline T cbrt_epsilon_imp(const T*, const Tag&)334{335 BOOST_MATH_STD_USING;336 static const T cbrt_eps = pow(tools::epsilon<T>(), T(1) / 3);337 return cbrt_eps;338}339 340template <class T>341BOOST_MATH_GPU_ENABLED inline T cbrt_epsilon_imp(const T*, const boost::math::integral_constant<int, 0>&)342{343 BOOST_MATH_STD_USING;344 return pow(tools::epsilon<T>(), T(1) / 3);345}346 347template <class T>348BOOST_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)349{350 return static_cast<T>(0.018581361171917516667460937040007436176452688944747L);351}352 353template <class T>354BOOST_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)355{356 return static_cast<T>(0.0001220703125L);357}358 359template <class T>360BOOST_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)361{362 return static_cast<T>(0.18145860519450699870567321328132261891067079047605e-4L);363}364 365template <class T>366BOOST_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)367{368 return static_cast<T>(0.37252902984619140625e-8L);369}370 371template <class T, class Tag>372BOOST_MATH_GPU_ENABLED inline T forth_root_epsilon_imp(const T*, const Tag&)373{374 BOOST_MATH_STD_USING375 static const T r_eps = sqrt(sqrt(tools::epsilon<T>()));376 return r_eps;377}378 379template <class T>380BOOST_MATH_GPU_ENABLED inline T forth_root_epsilon_imp(const T*, const boost::math::integral_constant<int, 0>&)381{382 BOOST_MATH_STD_USING383 return sqrt(sqrt(tools::epsilon<T>()));384}385 386template <class T>387struct root_epsilon_traits388{389 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;390 static constexpr bool has_noexcept = (tag_type::value == 113) || (tag_type::value == 64) || (tag_type::value == 53) || (tag_type::value == 24);391};392 393}394 395template <class T>396BOOST_MATH_GPU_ENABLED inline constexpr T root_epsilon() noexcept(boost::math::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)397{398 return detail::root_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());399}400 401template <class T>402BOOST_MATH_GPU_ENABLED inline constexpr T cbrt_epsilon() noexcept(boost::math::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)403{404 return detail::cbrt_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());405}406 407template <class T>408BOOST_MATH_GPU_ENABLED inline constexpr T forth_root_epsilon() noexcept(boost::math::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)409{410 return detail::forth_root_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());411}412 413} // namespace tools414} // namespace math415} // namespace boost416 417#endif // BOOST_MATH_TOOLS_PRECISION_INCLUDED418 419