![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2014 Benoit Steiner <benoit.steiner.goog@gmail.com> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_CXX11_TENSOR_TENSOR_INTDIV_H 00011 #define EIGEN_CXX11_TENSOR_TENSOR_INTDIV_H 00012 00013 00014 namespace Eigen { 00015 00029 namespace internal { 00030 00031 namespace { 00032 00033 // Note: result is undefined if val == 0 00034 template <typename T> 00035 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE 00036 typename internal::enable_if<sizeof(T)==4,int>::type count_leading_zeros(const T val) 00037 { 00038 #ifdef __CUDA_ARCH__ 00039 return __clz(val); 00040 #elif EIGEN_COMP_MSVC 00041 unsigned long index; 00042 _BitScanReverse(&index, val); 00043 return 31 - index; 00044 #else 00045 EIGEN_STATIC_ASSERT(sizeof(unsigned long long) == 8, YOU_MADE_A_PROGRAMMING_MISTAKE); 00046 return __builtin_clz(static_cast<uint32_t>(val)); 00047 #endif 00048 } 00049 00050 template <typename T> 00051 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE 00052 typename internal::enable_if<sizeof(T)==8,int>::type count_leading_zeros(const T val) 00053 { 00054 #ifdef __CUDA_ARCH__ 00055 return __clzll(val); 00056 #elif EIGEN_COMP_MSVC && EIGEN_ARCH_x86_64 00057 unsigned long index; 00058 _BitScanReverse64(&index, val); 00059 return 63 - index; 00060 #elif EIGEN_COMP_MSVC 00061 // MSVC's _BitScanReverse64 is not available for 32bits builds. 00062 unsigned int lo = (unsigned int)(val&0xffffffff); 00063 unsigned int hi = (unsigned int)((val>>32)&0xffffffff); 00064 int n; 00065 if(hi==0) 00066 n = 32 + count_leading_zeros<unsigned int>(lo); 00067 else 00068 n = count_leading_zeros<unsigned int>(hi); 00069 return n; 00070 #else 00071 EIGEN_STATIC_ASSERT(sizeof(unsigned long long) == 8, YOU_MADE_A_PROGRAMMING_MISTAKE); 00072 return __builtin_clzll(static_cast<uint64_t>(val)); 00073 #endif 00074 } 00075 00076 template <typename T> 00077 struct UnsignedTraits { 00078 typedef typename conditional<sizeof(T) == 8, uint64_t, uint32_t>::type type; 00079 }; 00080 00081 template <typename T> 00082 struct DividerTraits { 00083 typedef typename UnsignedTraits<T>::type type; 00084 static const int N = sizeof(T) * 8; 00085 }; 00086 00087 template <typename T> 00088 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint32_t muluh(const uint32_t a, const T b) { 00089 #if defined(__CUDA_ARCH__) 00090 return __umulhi(a, b); 00091 #else 00092 return (static_cast<uint64_t>(a) * b) >> 32; 00093 #endif 00094 } 00095 00096 template <typename T> 00097 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint64_t muluh(const uint64_t a, const T b) { 00098 #if defined(__CUDA_ARCH__) 00099 return __umul64hi(a, b); 00100 #elif defined(__SIZEOF_INT128__) 00101 __uint128_t v = static_cast<__uint128_t>(a) * static_cast<__uint128_t>(b); 00102 return static_cast<uint64_t>(v >> 64); 00103 #else 00104 return (TensorUInt128<static_val<0>, uint64_t>(a) * TensorUInt128<static_val<0>, uint64_t>(b)).upper(); 00105 #endif 00106 } 00107 00108 template <int N, typename T> 00109 struct DividerHelper { 00110 static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint32_t computeMultiplier(const int log_div, const T divider) { 00111 EIGEN_STATIC_ASSERT(N == 32, YOU_MADE_A_PROGRAMMING_MISTAKE); 00112 return static_cast<uint32_t>((static_cast<uint64_t>(1) << (N+log_div)) / divider - (static_cast<uint64_t>(1) << N) + 1); 00113 } 00114 }; 00115 00116 template <typename T> 00117 struct DividerHelper<64, T> { 00118 static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint64_t computeMultiplier(const int log_div, const T divider) { 00119 #if defined(__SIZEOF_INT128__) && !defined(__CUDA_ARCH__) 00120 return static_cast<uint64_t>((static_cast<__uint128_t>(1) << (64+log_div)) / static_cast<__uint128_t>(divider) - (static_cast<__uint128_t>(1) << 64) + 1); 00121 #else 00122 const uint64_t shift = 1ULL << log_div; 00123 TensorUInt128<uint64_t, uint64_t> result = TensorUInt128<uint64_t, static_val<0> >(shift, 0) / TensorUInt128<static_val<0>, uint64_t>(divider) 00124 - TensorUInt128<static_val<1>, static_val<0> >(1, 0) 00125 + TensorUInt128<static_val<0>, static_val<1> >(1); 00126 return static_cast<uint64_t>(result); 00127 #endif 00128 } 00129 }; 00130 } 00131 00132 00133 template <typename T, bool div_gt_one = false> 00134 struct TensorIntDivisor { 00135 public: 00136 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() { 00137 multiplier = 0; 00138 shift1 = 0; 00139 shift2 = 0; 00140 } 00141 00142 // Must have 0 < divider < 2^31. This is relaxed to 00143 // 0 < divider < 2^63 when using 64-bit indices on platforms that support 00144 // the __uint128_t type. 00145 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor(const T divider) { 00146 const int N = DividerTraits<T>::N; 00147 eigen_assert(static_cast<typename UnsignedTraits<T>::type>(divider) < NumTraits<UnsignedType>::highest()/2); 00148 eigen_assert(divider > 0); 00149 00150 // fast ln2 00151 const int leading_zeros = count_leading_zeros(static_cast<UnsignedType>(divider)); 00152 int log_div = N - leading_zeros; 00153 // if divider is a power of two then log_div is 1 more than it should be. 00154 if ((static_cast<typename UnsignedTraits<T>::type>(1) << (log_div-1)) == static_cast<typename UnsignedTraits<T>::type>(divider)) 00155 log_div--; 00156 00157 multiplier = DividerHelper<N, T>::computeMultiplier(log_div, divider); 00158 shift1 = log_div > 1 ? 1 : log_div; 00159 shift2 = log_div > 1 ? log_div-1 : 0; 00160 } 00161 00162 // Must have 0 <= numerator. On platforms that dont support the __uint128_t 00163 // type numerator should also be less than 2^32-1. 00164 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T divide(const T numerator) const { 00165 eigen_assert(static_cast<typename UnsignedTraits<T>::type>(numerator) < NumTraits<UnsignedType>::highest()/2); 00166 //eigen_assert(numerator >= 0); // this is implicitly asserted by the line above 00167 00168 UnsignedType t1 = muluh(multiplier, numerator); 00169 UnsignedType t = (static_cast<UnsignedType>(numerator) - t1) >> shift1; 00170 return (t1 + t) >> shift2; 00171 } 00172 00173 private: 00174 typedef typename DividerTraits<T>::type UnsignedType; 00175 UnsignedType multiplier; 00176 int32_t shift1; 00177 int32_t shift2; 00178 }; 00179 00180 00181 // Optimized version for signed 32 bit integers. 00182 // Derived from Hacker's Delight. 00183 // Only works for divisors strictly greater than one 00184 template <> 00185 class TensorIntDivisor<int32_t, true> { 00186 public: 00187 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() { 00188 magic = 0; 00189 shift = 0; 00190 } 00191 // Must have 2 <= divider 00192 EIGEN_DEVICE_FUNC TensorIntDivisor(int32_t divider) { 00193 eigen_assert(divider >= 2); 00194 calcMagic(divider); 00195 } 00196 00197 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE int divide(const int32_t n) const { 00198 #ifdef __CUDA_ARCH__ 00199 return (__umulhi(magic, n) >> shift); 00200 #else 00201 uint64_t v = static_cast<uint64_t>(magic) * static_cast<uint64_t>(n); 00202 return (static_cast<uint32_t>(v >> 32) >> shift); 00203 #endif 00204 } 00205 00206 private: 00207 // Compute the magic numbers. See Hacker's Delight section 10 for an in 00208 // depth explanation. 00209 EIGEN_DEVICE_FUNC void calcMagic(int32_t d) { 00210 const unsigned two31 = 0x80000000; // 2**31. 00211 unsigned ad = d; 00212 unsigned t = two31 + (ad >> 31); 00213 unsigned anc = t - 1 - t%ad; // Absolute value of nc. 00214 int p = 31; // Init. p. 00215 unsigned q1 = two31/anc; // Init. q1 = 2**p/|nc|. 00216 unsigned r1 = two31 - q1*anc; // Init. r1 = rem(2**p, |nc|). 00217 unsigned q2 = two31/ad; // Init. q2 = 2**p/|d|. 00218 unsigned r2 = two31 - q2*ad; // Init. r2 = rem(2**p, |d|). 00219 unsigned delta = 0; 00220 do { 00221 p = p + 1; 00222 q1 = 2*q1; // Update q1 = 2**p/|nc|. 00223 r1 = 2*r1; // Update r1 = rem(2**p, |nc|). 00224 if (r1 >= anc) { // (Must be an unsigned 00225 q1 = q1 + 1; // comparison here). 00226 r1 = r1 - anc;} 00227 q2 = 2*q2; // Update q2 = 2**p/|d|. 00228 r2 = 2*r2; // Update r2 = rem(2**p, |d|). 00229 if (r2 >= ad) { // (Must be an unsigned 00230 q2 = q2 + 1; // comparison here). 00231 r2 = r2 - ad;} 00232 delta = ad - r2; 00233 } while (q1 < delta || (q1 == delta && r1 == 0)); 00234 00235 magic = (unsigned)(q2 + 1); 00236 shift = p - 32; 00237 } 00238 00239 uint32_t magic; 00240 int32_t shift; 00241 }; 00242 00243 00244 template <typename T, bool div_gt_one> 00245 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator / (const T& numerator, const TensorIntDivisor<T, div_gt_one>& divisor) { 00246 return divisor.divide(numerator); 00247 } 00248 00249 00250 } // end namespace internal 00251 } // end namespace Eigen 00252 00253 #endif // EIGEN_CXX11_TENSOR_TENSOR_INTDIV_H