precision.hpp 13 KB

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