TensorIntDiv.h
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
 All Classes Functions Variables Typedefs Enumerator