TensorIntDiv.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2014 Benoit Steiner <benoit.steiner.goog@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_CXX11_TENSOR_TENSOR_INTDIV_H
11 #define EIGEN_CXX11_TENSOR_TENSOR_INTDIV_H
12 
13 
14 namespace Eigen {
15 
29 namespace internal {
30 
31 namespace {
32  // Note: result is undefined if val == 0
33  template <typename T>
34  EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE int count_leading_zeros(const T val)
35  {
36 #ifdef __CUDA_ARCH__
37  return (sizeof(T) == 8) ? __clzll(val) : __clz(val);
38 #elif EIGEN_COMP_MSVC
39  DWORD leading_zeros = 0;
40  if (sizeof(T) == 8) {
41  _BitScanReverse64(&leading_zero, val);
42  }
43  else {
44  _BitScanReverse(&leading_zero, val);
45  }
46  return leading_zeros;
47 #else
48  return (sizeof(T) == 8) ?
49  __builtin_clzl(static_cast<uint64_t>(val)) :
50  __builtin_clz(static_cast<uint32_t>(val));
51 #endif
52  }
53 
54  template <typename T>
55  struct UnsignedTraits {
56  typedef typename conditional<sizeof(T) == 8, uint64_t, uint32_t>::type type;
57  };
58 
59  template <typename T>
60  struct DividerTraits {
61  typedef typename UnsignedTraits<T>::type type;
62  static const int N = sizeof(T) * 8;
63  };
64 
65  template <typename T>
66  EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint32_t muluh(const uint32_t a, const T b) {
67 #if defined(__CUDA_ARCH__)
68  return __umulhi(a, b);
69 #else
70  return (static_cast<uint64_t>(a) * b) >> 32;
71 #endif
72  }
73 
74  template <typename T>
75  EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint64_t muluh(const uint64_t a, const T b) {
76 #if defined(__CUDA_ARCH__)
77  return __umul64hi(a, b);
78 #elif defined(__SIZEOF_INT128__)
79  __uint128_t v = static_cast<__uint128_t>(a) * static_cast<__uint128_t>(b);
80  return static_cast<uint64_t>(v >> 64);
81 #else
82  return (TensorUInt128<static_val<0>, uint64_t>(a) * TensorUInt128<static_val<0>, uint64_t>(b)).upper();
83 #endif
84  }
85 
86  template <int N, typename T>
87  struct DividerHelper {
88  static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint32_t computeMultiplier(const int log_div, const T divider) {
89  EIGEN_STATIC_ASSERT(N == 32, YOU_MADE_A_PROGRAMMING_MISTAKE);
90  return static_cast<uint32_t>((static_cast<uint64_t>(1) << (N+log_div)) / divider - (static_cast<uint64_t>(1) << N) + 1);
91  }
92  };
93 
94  template <typename T>
95  struct DividerHelper<64, T> {
96  static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint64_t computeMultiplier(const int log_div, const T divider) {
97 #if defined(__SIZEOF_INT128__) && !defined(__CUDA_ARCH__)
98  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);
99 #else
100  const uint64_t shift = 1ULL << log_div;
101  TensorUInt128<uint64_t, uint64_t> result = (TensorUInt128<uint64_t, static_val<0> >(shift, 0) / TensorUInt128<static_val<0>, uint64_t>(divider) - TensorUInt128<static_val<1>, static_val<0> >(1, 0) + TensorUInt128<static_val<0>, static_val<1> >(1));
102  return static_cast<uint64_t>(result);
103 #endif
104  }
105  };
106 }
107 
108 
109 template <typename T, bool div_gt_one = false>
110 struct TensorIntDivisor {
111  public:
112  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() {
113  multiplier = 0;
114  shift1 = 0;
115  shift2 = 0;
116  }
117 
118  // Must have 0 < divider < 2^31. This is relaxed to
119  // 0 < divider < 2^63 when using 64-bit indices on platforms that support
120  // the __uint128_t type.
121  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor(const T divider) {
122  const int N = DividerTraits<T>::N;
123  eigen_assert(static_cast<typename UnsignedTraits<T>::type>(divider) < NumTraits<UnsignedType>::highest()/2);
124  eigen_assert(divider > 0);
125 
126  // fast ln2
127  const int leading_zeros = count_leading_zeros(static_cast<UnsignedType>(divider));
128  int log_div = N - leading_zeros;
129  // if divider is a power of two then log_div is 1 more than it should be.
130  if ((static_cast<typename UnsignedTraits<T>::type>(1) << (log_div-1)) == static_cast<typename UnsignedTraits<T>::type>(divider))
131  log_div--;
132 
133  multiplier = DividerHelper<N, T>::computeMultiplier(log_div, divider);
134  shift1 = log_div > 1 ? 1 : log_div;
135  shift2 = log_div > 1 ? log_div-1 : 0;
136  }
137 
138  // Must have 0 <= numerator. On platforms that dont support the __uint128_t
139  // type numerator should also be less than 2^32-1.
140  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T divide(const T numerator) const {
141  eigen_assert(static_cast<typename UnsignedTraits<T>::type>(numerator) < NumTraits<UnsignedType>::highest()/2);
142  //eigen_assert(numerator >= 0); // this is implicitly asserted by the line above
143 
144  UnsignedType t1 = muluh(multiplier, numerator);
145  UnsignedType t = (static_cast<UnsignedType>(numerator) - t1) >> shift1;
146  return (t1 + t) >> shift2;
147  }
148 
149  private:
150  typedef typename DividerTraits<T>::type UnsignedType;
151  UnsignedType multiplier;
152  int32_t shift1;
153  int32_t shift2;
154 };
155 
156 
157 // Optimized version for signed 32 bit integers.
158 // Derived from Hacker's Delight.
159 // Only works for divisors strictly greater than one
160 template <>
161 class TensorIntDivisor<int32_t, true> {
162  public:
163  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() {
164  magic = 0;
165  shift = 0;
166  }
167  // Must have 2 <= divider
168  EIGEN_DEVICE_FUNC TensorIntDivisor(int32_t divider) {
169  eigen_assert(divider >= 2);
170  calcMagic(divider);
171  }
172 
173  EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE int divide(const int32_t n) const {
174 #ifdef __CUDA_ARCH__
175  return (__umulhi(magic, n) >> shift);
176 #else
177  uint64_t v = static_cast<uint64_t>(magic) * static_cast<uint64_t>(n);
178  return (static_cast<uint32_t>(v >> 32) >> shift);
179 #endif
180  }
181 
182 private:
183  // Compute the magic numbers. See Hacker's Delight section 10 for an in
184  // depth explanation.
185  EIGEN_DEVICE_FUNC void calcMagic(int32_t d) {
186  const unsigned two31 = 0x80000000; // 2**31.
187  unsigned ad = d;
188  unsigned t = two31 + (ad >> 31);
189  unsigned anc = t - 1 - t%ad; // Absolute value of nc.
190  int p = 31; // Init. p.
191  unsigned q1 = two31/anc; // Init. q1 = 2**p/|nc|.
192  unsigned r1 = two31 - q1*anc; // Init. r1 = rem(2**p, |nc|).
193  unsigned q2 = two31/ad; // Init. q2 = 2**p/|d|.
194  unsigned r2 = two31 - q2*ad; // Init. r2 = rem(2**p, |d|).
195  unsigned delta = 0;
196  do {
197  p = p + 1;
198  q1 = 2*q1; // Update q1 = 2**p/|nc|.
199  r1 = 2*r1; // Update r1 = rem(2**p, |nc|).
200  if (r1 >= anc) { // (Must be an unsigned
201  q1 = q1 + 1; // comparison here).
202  r1 = r1 - anc;}
203  q2 = 2*q2; // Update q2 = 2**p/|d|.
204  r2 = 2*r2; // Update r2 = rem(2**p, |d|).
205  if (r2 >= ad) { // (Must be an unsigned
206  q2 = q2 + 1; // comparison here).
207  r2 = r2 - ad;}
208  delta = ad - r2;
209  } while (q1 < delta || (q1 == delta && r1 == 0));
210 
211  magic = (unsigned)(q2 + 1);
212  shift = p - 32;
213  }
214 
215  uint32_t magic;
216  int32_t shift;
217 };
218 
219 
220 template <typename T, bool div_gt_one>
221 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator / (const T& numerator, const TensorIntDivisor<T, div_gt_one>& divisor) {
222  return divisor.divide(numerator);
223 }
224 
225 
226 } // end namespace internal
227 } // end namespace Eigen
228 
229 #endif // EIGEN_CXX11_TENSOR_TENSOR_INTDIV_H
Namespace containing all symbols from the Eigen library.
Definition: CXX11Meta.h:13