precision.hpp 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395
  1. // Copyright John Maddock 2005-2006.
  2. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. #ifndef BOOST_MATH_TOOLS_PRECISION_INCLUDED
  6. #define BOOST_MATH_TOOLS_PRECISION_INCLUDED
  7. #ifdef _MSC_VER
  8. #pragma once
  9. #endif
  10. #include <boost/math/tools/assert.hpp>
  11. #include <boost/math/policies/policy.hpp>
  12. #include <type_traits>
  13. #include <limits>
  14. #include <climits>
  15. #include <cmath>
  16. #include <cstdint>
  17. #include <cfloat> // LDBL_MANT_DIG
  18. namespace boost{ namespace math
  19. {
  20. namespace tools
  21. {
  22. // If T is not specialized, the functions digits, max_value and min_value,
  23. // all get synthesised automatically from std::numeric_limits.
  24. // However, if numeric_limits is not specialised for type RealType,
  25. // for example with NTL::RR type, then you will get a compiler error
  26. // when code tries to use these functions, unless you explicitly specialise them.
  27. // For example if the precision of RealType varies at runtime,
  28. // then numeric_limits support may not be appropriate,
  29. // see boost/math/tools/ntl.hpp for examples like
  30. // template <> NTL::RR max_value<NTL::RR> ...
  31. // See Conceptual Requirements for Real Number Types.
  32. template <class T>
  33. inline constexpr int digits(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) noexcept
  34. {
  35. static_assert( ::std::numeric_limits<T>::is_specialized, "Type T must be specialized");
  36. static_assert( ::std::numeric_limits<T>::radix == 2 || ::std::numeric_limits<T>::radix == 10, "Type T must have a radix of 2 or 10");
  37. return std::numeric_limits<T>::radix == 2
  38. ? std::numeric_limits<T>::digits
  39. : ((std::numeric_limits<T>::digits + 1) * 1000L) / 301L;
  40. }
  41. template <class T>
  42. inline constexpr T max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(std::is_floating_point<T>::value)
  43. {
  44. static_assert( ::std::numeric_limits<T>::is_specialized, "Type T must be specialized");
  45. return (std::numeric_limits<T>::max)();
  46. } // Also used as a finite 'infinite' value for - and +infinity, for example:
  47. // -max_value<double> = -1.79769e+308, max_value<double> = 1.79769e+308.
  48. template <class T>
  49. inline constexpr T min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(std::is_floating_point<T>::value)
  50. {
  51. static_assert( ::std::numeric_limits<T>::is_specialized, "Type T must be specialized");
  52. return (std::numeric_limits<T>::min)();
  53. }
  54. namespace detail{
  55. //
  56. // Logarithmic limits come next, note that although
  57. // we can compute these from the log of the max value
  58. // that is not in general thread safe (if we cache the value)
  59. // so it's better to specialise these:
  60. //
  61. // For type float first:
  62. //
  63. template <class T>
  64. 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)
  65. {
  66. return 88.0f;
  67. }
  68. template <class T>
  69. 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)
  70. {
  71. return -87.0f;
  72. }
  73. //
  74. // Now double:
  75. //
  76. template <class T>
  77. 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)
  78. {
  79. return 709.0;
  80. }
  81. template <class T>
  82. 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)
  83. {
  84. return -708.0;
  85. }
  86. //
  87. // 80 and 128-bit long doubles:
  88. //
  89. template <class T>
  90. 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)
  91. {
  92. return 11356.0L;
  93. }
  94. template <class T>
  95. 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)
  96. {
  97. return -11355.0L;
  98. }
  99. template <class T>
  100. inline T log_max_value(const std::integral_constant<int, 0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
  101. {
  102. BOOST_MATH_STD_USING
  103. #ifdef __SUNPRO_CC
  104. static const T m = boost::math::tools::max_value<T>();
  105. static const T val = log(m);
  106. #else
  107. static const T val = log(boost::math::tools::max_value<T>());
  108. #endif
  109. return val;
  110. }
  111. template <class T>
  112. inline T log_min_value(const std::integral_constant<int, 0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
  113. {
  114. BOOST_MATH_STD_USING
  115. #ifdef __SUNPRO_CC
  116. static const T m = boost::math::tools::min_value<T>();
  117. static const T val = log(m);
  118. #else
  119. static const T val = log(boost::math::tools::min_value<T>());
  120. #endif
  121. return val;
  122. }
  123. template <class T>
  124. inline constexpr T epsilon(const std::true_type& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(std::is_floating_point<T>::value)
  125. {
  126. return std::numeric_limits<T>::epsilon();
  127. }
  128. #if defined(__GNUC__) && ((LDBL_MANT_DIG == 106) || (__LDBL_MANT_DIG__ == 106))
  129. template <>
  130. 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)
  131. {
  132. // numeric_limits on Darwin (and elsewhere) tells lies here:
  133. // the issue is that long double on a few platforms is
  134. // really a "double double" which has a non-contiguous
  135. // mantissa: 53 bits followed by an unspecified number of
  136. // zero bits, followed by 53 more bits. Thus the apparent
  137. // precision of the type varies depending where it's been.
  138. // Set epsilon to the value that a 106 bit fixed mantissa
  139. // type would have, as that will give us sensible behaviour everywhere.
  140. //
  141. // This static assert fails for some unknown reason, so
  142. // disabled for now...
  143. // static_assert(std::numeric_limits<long double>::digits == 106);
  144. return 2.4651903288156618919116517665087e-32L;
  145. }
  146. #endif
  147. template <class T>
  148. inline T epsilon(const std::false_type& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
  149. {
  150. // Note: don't cache result as precision may vary at runtime:
  151. BOOST_MATH_STD_USING // for ADL of std names
  152. return ldexp(static_cast<T>(1), 1-policies::digits<T, policies::policy<> >());
  153. }
  154. template <class T>
  155. struct log_limit_traits
  156. {
  157. typedef typename std::conditional<
  158. (std::numeric_limits<T>::radix == 2) &&
  159. (std::numeric_limits<T>::max_exponent == 128
  160. || std::numeric_limits<T>::max_exponent == 1024
  161. || std::numeric_limits<T>::max_exponent == 16384),
  162. std::integral_constant<int, (std::numeric_limits<T>::max_exponent > INT_MAX ? INT_MAX : static_cast<int>(std::numeric_limits<T>::max_exponent))>,
  163. std::integral_constant<int, 0>
  164. >::type tag_type;
  165. static constexpr bool value = (tag_type::value != 0);
  166. static_assert(::std::numeric_limits<T>::is_specialized || !value, "Type T must be specialized or equal to 0");
  167. };
  168. template <class T, bool b> struct log_limit_noexcept_traits_imp : public log_limit_traits<T> {};
  169. template <class T> struct log_limit_noexcept_traits_imp<T, false> : public std::integral_constant<bool, false> {};
  170. template <class T>
  171. struct log_limit_noexcept_traits : public log_limit_noexcept_traits_imp<T, std::is_floating_point<T>::value> {};
  172. } // namespace detail
  173. #ifdef _MSC_VER
  174. #pragma warning(push)
  175. #pragma warning(disable:4309)
  176. #endif
  177. template <class T>
  178. inline constexpr T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
  179. {
  180. #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
  181. return detail::log_max_value<T>(typename detail::log_limit_traits<T>::tag_type());
  182. #else
  183. BOOST_MATH_ASSERT(::std::numeric_limits<T>::is_specialized);
  184. BOOST_MATH_STD_USING
  185. static const T val = log((std::numeric_limits<T>::max)());
  186. return val;
  187. #endif
  188. }
  189. template <class T>
  190. inline constexpr T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
  191. {
  192. #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
  193. return detail::log_min_value<T>(typename detail::log_limit_traits<T>::tag_type());
  194. #else
  195. BOOST_MATH_ASSERT(::std::numeric_limits<T>::is_specialized);
  196. BOOST_MATH_STD_USING
  197. static const T val = log((std::numeric_limits<T>::min)());
  198. return val;
  199. #endif
  200. }
  201. #ifdef _MSC_VER
  202. #pragma warning(pop)
  203. #endif
  204. template <class T>
  205. inline constexpr T epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) noexcept(std::is_floating_point<T>::value)
  206. {
  207. #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
  208. return detail::epsilon<T>(std::integral_constant<bool, ::std::numeric_limits<T>::is_specialized>());
  209. #else
  210. return ::std::numeric_limits<T>::is_specialized ?
  211. detail::epsilon<T>(std::true_type()) :
  212. detail::epsilon<T>(std::false_type());
  213. #endif
  214. }
  215. namespace detail{
  216. template <class T>
  217. inline constexpr T root_epsilon_imp(const std::integral_constant<int, 24>&) noexcept(std::is_floating_point<T>::value)
  218. {
  219. return static_cast<T>(0.00034526698300124390839884978618400831996329879769945L);
  220. }
  221. template <class T>
  222. inline constexpr T root_epsilon_imp(const T*, const std::integral_constant<int, 53>&) noexcept(std::is_floating_point<T>::value)
  223. {
  224. return static_cast<T>(0.1490116119384765625e-7L);
  225. }
  226. template <class T>
  227. inline constexpr T root_epsilon_imp(const T*, const std::integral_constant<int, 64>&) noexcept(std::is_floating_point<T>::value)
  228. {
  229. return static_cast<T>(0.32927225399135962333569506281281311031656150598474e-9L);
  230. }
  231. template <class T>
  232. inline constexpr T root_epsilon_imp(const T*, const std::integral_constant<int, 113>&) noexcept(std::is_floating_point<T>::value)
  233. {
  234. return static_cast<T>(0.1387778780781445675529539585113525390625e-16L);
  235. }
  236. template <class T, class Tag>
  237. inline T root_epsilon_imp(const T*, const Tag&)
  238. {
  239. BOOST_MATH_STD_USING
  240. static const T r_eps = sqrt(tools::epsilon<T>());
  241. return r_eps;
  242. }
  243. template <class T>
  244. inline T root_epsilon_imp(const T*, const std::integral_constant<int, 0>&)
  245. {
  246. BOOST_MATH_STD_USING
  247. return sqrt(tools::epsilon<T>());
  248. }
  249. template <class T>
  250. inline constexpr T cbrt_epsilon_imp(const std::integral_constant<int, 24>&) noexcept(std::is_floating_point<T>::value)
  251. {
  252. return static_cast<T>(0.0049215666011518482998719164346805794944150447839903L);
  253. }
  254. template <class T>
  255. inline constexpr T cbrt_epsilon_imp(const T*, const std::integral_constant<int, 53>&) noexcept(std::is_floating_point<T>::value)
  256. {
  257. return static_cast<T>(6.05545445239333906078989272793696693569753008995e-6L);
  258. }
  259. template <class T>
  260. inline constexpr T cbrt_epsilon_imp(const T*, const std::integral_constant<int, 64>&) noexcept(std::is_floating_point<T>::value)
  261. {
  262. return static_cast<T>(4.76837158203125e-7L);
  263. }
  264. template <class T>
  265. inline constexpr T cbrt_epsilon_imp(const T*, const std::integral_constant<int, 113>&) noexcept(std::is_floating_point<T>::value)
  266. {
  267. return static_cast<T>(5.7749313854154005630396773604745549542403508090496e-12L);
  268. }
  269. template <class T, class Tag>
  270. inline T cbrt_epsilon_imp(const T*, const Tag&)
  271. {
  272. BOOST_MATH_STD_USING;
  273. static const T cbrt_eps = pow(tools::epsilon<T>(), T(1) / 3);
  274. return cbrt_eps;
  275. }
  276. template <class T>
  277. inline T cbrt_epsilon_imp(const T*, const std::integral_constant<int, 0>&)
  278. {
  279. BOOST_MATH_STD_USING;
  280. return pow(tools::epsilon<T>(), T(1) / 3);
  281. }
  282. template <class T>
  283. inline constexpr T forth_root_epsilon_imp(const T*, const std::integral_constant<int, 24>&) noexcept(std::is_floating_point<T>::value)
  284. {
  285. return static_cast<T>(0.018581361171917516667460937040007436176452688944747L);
  286. }
  287. template <class T>
  288. inline constexpr T forth_root_epsilon_imp(const T*, const std::integral_constant<int, 53>&) noexcept(std::is_floating_point<T>::value)
  289. {
  290. return static_cast<T>(0.0001220703125L);
  291. }
  292. template <class T>
  293. inline constexpr T forth_root_epsilon_imp(const T*, const std::integral_constant<int, 64>&) noexcept(std::is_floating_point<T>::value)
  294. {
  295. return static_cast<T>(0.18145860519450699870567321328132261891067079047605e-4L);
  296. }
  297. template <class T>
  298. inline constexpr T forth_root_epsilon_imp(const T*, const std::integral_constant<int, 113>&) noexcept(std::is_floating_point<T>::value)
  299. {
  300. return static_cast<T>(0.37252902984619140625e-8L);
  301. }
  302. template <class T, class Tag>
  303. inline T forth_root_epsilon_imp(const T*, const Tag&)
  304. {
  305. BOOST_MATH_STD_USING
  306. static const T r_eps = sqrt(sqrt(tools::epsilon<T>()));
  307. return r_eps;
  308. }
  309. template <class T>
  310. inline T forth_root_epsilon_imp(const T*, const std::integral_constant<int, 0>&)
  311. {
  312. BOOST_MATH_STD_USING
  313. return sqrt(sqrt(tools::epsilon<T>()));
  314. }
  315. template <class T>
  316. struct root_epsilon_traits
  317. {
  318. 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;
  319. static constexpr bool has_noexcept = (tag_type::value == 113) || (tag_type::value == 64) || (tag_type::value == 53) || (tag_type::value == 24);
  320. };
  321. }
  322. template <class T>
  323. inline constexpr T root_epsilon() noexcept(std::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)
  324. {
  325. return detail::root_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());
  326. }
  327. template <class T>
  328. inline constexpr T cbrt_epsilon() noexcept(std::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)
  329. {
  330. return detail::cbrt_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());
  331. }
  332. template <class T>
  333. inline constexpr T forth_root_epsilon() noexcept(std::is_floating_point<T>::value && detail::root_epsilon_traits<T>::has_noexcept)
  334. {
  335. return detail::forth_root_epsilon_imp(static_cast<T const*>(nullptr), typename detail::root_epsilon_traits<T>::tag_type());
  336. }
  337. } // namespace tools
  338. } // namespace math
  339. } // namespace boost
  340. #endif // BOOST_MATH_TOOLS_PRECISION_INCLUDED