math.hpp 7.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290
  1. /////////////////////////////////////////////////////////////////////////////
  2. //
  3. // (C) Copyright Ion Gaztanaga 2014-2014
  4. //
  5. // Distributed under the Boost Software License, Version 1.0.
  6. // (See accompanying file LICENSE_1_0.txt or copy at
  7. // http://www.boost.org/LICENSE_1_0.txt)
  8. //
  9. // See http://www.boost.org/libs/intrusive for documentation.
  10. //
  11. /////////////////////////////////////////////////////////////////////////////
  12. #ifndef BOOST_INTRUSIVE_DETAIL_MATH_HPP
  13. #define BOOST_INTRUSIVE_DETAIL_MATH_HPP
  14. #ifndef BOOST_CONFIG_HPP
  15. # include <boost/config.hpp>
  16. #endif
  17. #if defined(BOOST_HAS_PRAGMA_ONCE)
  18. # pragma once
  19. #endif
  20. #include <cstddef>
  21. #include <climits>
  22. #include <boost/intrusive/detail/mpl.hpp>
  23. #include <cstring>
  24. namespace boost {
  25. namespace intrusive {
  26. namespace detail {
  27. ///////////////////////////
  28. // floor_log2 Dispatcher
  29. ////////////////////////////
  30. #if defined(_MSC_VER) && (_MSC_VER >= 1300)
  31. }}} //namespace boost::intrusive::detail
  32. //Use _BitScanReverseXX intrinsics
  33. #if defined(_M_X64) || defined(_M_AMD64) || defined(_M_IA64) //64 bit target
  34. #define BOOST_INTRUSIVE_BSR_INTRINSIC_64_BIT
  35. #endif
  36. #ifndef __INTRIN_H_ // Avoid including any windows system header
  37. #ifdef __cplusplus
  38. extern "C" {
  39. #endif // __cplusplus
  40. #if defined(BOOST_INTRUSIVE_BSR_INTRINSIC_64_BIT) //64 bit target
  41. unsigned char _BitScanReverse64(unsigned long *index, unsigned __int64 mask);
  42. #pragma intrinsic(_BitScanReverse64)
  43. #else //32 bit target
  44. unsigned char _BitScanReverse(unsigned long *index, unsigned long mask);
  45. #pragma intrinsic(_BitScanReverse)
  46. #endif
  47. #ifdef __cplusplus
  48. }
  49. #endif // __cplusplus
  50. #endif // __INTRIN_H_
  51. #ifdef BOOST_INTRUSIVE_BSR_INTRINSIC_64_BIT
  52. #define BOOST_INTRUSIVE_BSR_INTRINSIC _BitScanReverse64
  53. #undef BOOST_INTRUSIVE_BSR_INTRINSIC_64_BIT
  54. #else
  55. #define BOOST_INTRUSIVE_BSR_INTRINSIC _BitScanReverse
  56. #endif
  57. namespace boost {
  58. namespace intrusive {
  59. namespace detail {
  60. inline std::size_t floor_log2 (std::size_t x)
  61. {
  62. unsigned long log2;
  63. BOOST_INTRUSIVE_BSR_INTRINSIC( &log2, x );
  64. return static_cast<std::size_t>(log2);
  65. }
  66. #undef BOOST_INTRUSIVE_BSR_INTRINSIC
  67. #elif defined(__GNUC__) && ((__GNUC__ >= 4) || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4)) //GCC >=3.4
  68. //Compile-time error in case of missing specialization
  69. template<class Uint>
  70. struct builtin_clz_dispatch;
  71. #if defined(BOOST_HAS_LONG_LONG)
  72. template<>
  73. struct builtin_clz_dispatch< ::boost::ulong_long_type >
  74. {
  75. static ::boost::ulong_long_type call(::boost::ulong_long_type n)
  76. { return __builtin_clzll(n); }
  77. };
  78. #endif
  79. template<>
  80. struct builtin_clz_dispatch<unsigned long>
  81. {
  82. static unsigned long call(unsigned long n)
  83. { return __builtin_clzl(n); }
  84. };
  85. template<>
  86. struct builtin_clz_dispatch<unsigned int>
  87. {
  88. static unsigned int call(unsigned int n)
  89. { return __builtin_clz(n); }
  90. };
  91. inline std::size_t floor_log2(std::size_t n)
  92. {
  93. return sizeof(std::size_t)*CHAR_BIT - std::size_t(1) - builtin_clz_dispatch<std::size_t>::call(n);
  94. }
  95. #else //Portable methods
  96. ////////////////////////////
  97. // Generic method
  98. ////////////////////////////
  99. inline std::size_t floor_log2_get_shift(std::size_t n, true_ )//power of two size_t
  100. { return n >> 1; }
  101. inline std::size_t floor_log2_get_shift(std::size_t n, false_ )//non-power of two size_t
  102. { return (n >> 1) + ((n & 1u) & (n != 1)); }
  103. template<std::size_t N>
  104. inline std::size_t floor_log2 (std::size_t x, integral_constant<std::size_t, N>)
  105. {
  106. const std::size_t Bits = N;
  107. const bool Size_t_Bits_Power_2= !(Bits & (Bits-1));
  108. std::size_t n = x;
  109. std::size_t log2 = 0;
  110. std::size_t remaining_bits = Bits;
  111. std::size_t shift = floor_log2_get_shift(remaining_bits, bool_<Size_t_Bits_Power_2>());
  112. while(shift){
  113. std::size_t tmp = n >> shift;
  114. if (tmp){
  115. log2 += shift, n = tmp;
  116. }
  117. shift = floor_log2_get_shift(shift, bool_<Size_t_Bits_Power_2>());
  118. }
  119. return log2;
  120. }
  121. ////////////////////////////
  122. // DeBruijn method
  123. ////////////////////////////
  124. //Taken from:
  125. //http://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
  126. //Thanks to Desmond Hume
  127. inline std::size_t floor_log2 (std::size_t v, integral_constant<std::size_t, 32>)
  128. {
  129. static const int MultiplyDeBruijnBitPosition[32] =
  130. {
  131. 0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
  132. 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
  133. };
  134. v |= v >> 1;
  135. v |= v >> 2;
  136. v |= v >> 4;
  137. v |= v >> 8;
  138. v |= v >> 16;
  139. return MultiplyDeBruijnBitPosition[(std::size_t)(v * 0x07C4ACDDU) >> 27];
  140. }
  141. inline std::size_t floor_log2 (std::size_t v, integral_constant<std::size_t, 64>)
  142. {
  143. static const std::size_t MultiplyDeBruijnBitPosition[64] = {
  144. 63, 0, 58, 1, 59, 47, 53, 2,
  145. 60, 39, 48, 27, 54, 33, 42, 3,
  146. 61, 51, 37, 40, 49, 18, 28, 20,
  147. 55, 30, 34, 11, 43, 14, 22, 4,
  148. 62, 57, 46, 52, 38, 26, 32, 41,
  149. 50, 36, 17, 19, 29, 10, 13, 21,
  150. 56, 45, 25, 31, 35, 16, 9, 12,
  151. 44, 24, 15, 8, 23, 7, 6, 5};
  152. v |= v >> 1;
  153. v |= v >> 2;
  154. v |= v >> 4;
  155. v |= v >> 8;
  156. v |= v >> 16;
  157. v |= v >> 32;
  158. return MultiplyDeBruijnBitPosition[((std::size_t)((v - (v >> 1))*0x07EDD5E59A4E28C2ULL)) >> 58];
  159. }
  160. inline std::size_t floor_log2 (std::size_t x)
  161. {
  162. const std::size_t Bits = sizeof(std::size_t)*CHAR_BIT;
  163. return floor_log2(x, integral_constant<std::size_t, Bits>());
  164. }
  165. #endif
  166. //Thanks to Laurent de Soras in
  167. //http://www.flipcode.com/archives/Fast_log_Function.shtml
  168. inline float fast_log2 (float val)
  169. {
  170. float f = val;
  171. unsigned x;
  172. std::memcpy(&x, &val, sizeof(f));
  173. const int log_2 = int((x >> 23) & 255) - 128;
  174. x &= ~(unsigned(255u) << 23u);
  175. x += unsigned(127) << 23u;
  176. std::memcpy(&val, &x, sizeof(f));
  177. //1+log2(m), m ranging from 1 to 2
  178. //3rd degree polynomial keeping first derivate continuity.
  179. //For less precision the line can be commented out
  180. val = ((-1.f/3.f) * val + 2.f) * val - (2.f/3.f);
  181. return val + static_cast<float>(log_2);
  182. }
  183. inline bool is_pow2(std::size_t x)
  184. { return (x & (x-1)) == 0; }
  185. template<std::size_t N>
  186. struct static_is_pow2
  187. {
  188. static const bool value = (N & (N-1)) == 0;
  189. };
  190. inline std::size_t ceil_log2 (std::size_t x)
  191. {
  192. return static_cast<std::size_t>(!(is_pow2)(x)) + floor_log2(x);
  193. }
  194. inline std::size_t ceil_pow2 (std::size_t x)
  195. {
  196. return std::size_t(1u) << (ceil_log2)(x);
  197. }
  198. inline std::size_t previous_or_equal_pow2(std::size_t x)
  199. {
  200. return std::size_t(1u) << floor_log2(x);
  201. }
  202. template<class SizeType, std::size_t N>
  203. struct numbits_eq
  204. {
  205. static const bool value = sizeof(SizeType)*CHAR_BIT == N;
  206. };
  207. template<class SizeType, class Enabler = void >
  208. struct sqrt2_pow_max;
  209. template <class SizeType>
  210. struct sqrt2_pow_max<SizeType, typename voider<typename enable_if< numbits_eq<SizeType, 32> >::type>::type>
  211. {
  212. static const SizeType value = 0xb504f334;
  213. static const std::size_t pow = 31;
  214. };
  215. #ifndef BOOST_NO_INT64_T
  216. template <class SizeType>
  217. struct sqrt2_pow_max<SizeType, typename voider<typename enable_if< numbits_eq<SizeType, 64> >::type>::type>
  218. {
  219. static const SizeType value = 0xb504f333f9de6484ull;
  220. static const std::size_t pow = 63;
  221. };
  222. #endif //BOOST_NO_INT64_T
  223. // Returns floor(pow(sqrt(2), x * 2 + 1)).
  224. // Defined for X from 0 up to the number of bits in size_t minus 1.
  225. inline std::size_t sqrt2_pow_2xplus1 (std::size_t x)
  226. {
  227. const std::size_t value = (std::size_t)sqrt2_pow_max<std::size_t>::value;
  228. const std::size_t pow = (std::size_t)sqrt2_pow_max<std::size_t>::pow;
  229. return (value >> (pow - x)) + 1;
  230. }
  231. } //namespace detail
  232. } //namespace intrusive
  233. } //namespace boost
  234. #endif //BOOST_INTRUSIVE_DETAIL_MATH_HPP