precision.hpp 16 KB

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