Eigen-Contrib  5.0.1
 
Loading...
Searching...
No Matches
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// SPDX-License-Identifier: MPL-2.0
10
11#ifndef EIGEN_TENSOR_TENSOR_INTDIV_H
12#define EIGEN_TENSOR_TENSOR_INTDIV_H
13
14// IWYU pragma: private
15#include "./InternalHeaderCheck.h"
16
17namespace Eigen {
18
19namespace internal {
20
21// Note: result is undefined if val == 0
22template <typename T>
23EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::enable_if_t<sizeof(T) == 4, int> count_leading_zeros(const T val) {
24#ifdef EIGEN_GPU_COMPILE_PHASE
25 return __clz(val);
26#elif defined(SYCL_DEVICE_ONLY)
27 return cl::sycl::clz(val);
28#elif EIGEN_COMP_MSVC
29 unsigned long index;
30 _BitScanReverse(&index, val);
31 return 31 - index;
32#else
33 EIGEN_STATIC_ASSERT(sizeof(unsigned long long) == 8, YOU_MADE_A_PROGRAMMING_MISTAKE);
34 return __builtin_clz(static_cast<uint32_t>(val));
35#endif
36}
37
38template <typename T>
39EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::enable_if_t<sizeof(T) == 8, int> count_leading_zeros(const T val) {
40#ifdef EIGEN_GPU_COMPILE_PHASE
41 return __clzll(val);
42#elif defined(SYCL_DEVICE_ONLY)
43 return static_cast<int>(cl::sycl::clz(val));
44#elif EIGEN_COMP_MSVC && EIGEN_ARCH_x86_64
45 unsigned long index;
46 _BitScanReverse64(&index, val);
47 return 63 - index;
48#elif EIGEN_COMP_MSVC
49 // MSVC's _BitScanReverse64 is not available for 32bits builds.
50 unsigned int lo = (unsigned int)(val & 0xffffffff);
51 unsigned int hi = (unsigned int)((val >> 32) & 0xffffffff);
52 int n;
53 if (hi == 0)
54 n = 32 + count_leading_zeros<unsigned int>(lo);
55 else
56 n = count_leading_zeros<unsigned int>(hi);
57 return n;
58#else
59 EIGEN_STATIC_ASSERT(sizeof(unsigned long long) == 8, YOU_MADE_A_PROGRAMMING_MISTAKE);
60 return __builtin_clzll(static_cast<uint64_t>(val));
61#endif
62}
63
64template <typename T>
65struct UnsignedTraits {
66 typedef std::conditional_t<sizeof(T) == 8, uint64_t, uint32_t> type;
67};
68
69template <typename T>
70struct DividerTraits {
71 typedef typename UnsignedTraits<T>::type type;
72 static constexpr int N = sizeof(T) * 8;
73};
74
75template <typename T>
76EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint32_t muluh(const uint32_t a, const T b) {
77#if defined(EIGEN_GPU_COMPILE_PHASE)
78 return __umulhi(a, b);
79#elif defined(SYCL_DEVICE_ONLY)
80 return cl::sycl::mul_hi(a, static_cast<uint32_t>(b));
81#else
82 return (static_cast<uint64_t>(a) * b) >> 32;
83#endif
84}
85
86template <typename T>
87EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint64_t muluh(const uint64_t a, const T b) {
88#if defined(EIGEN_GPU_COMPILE_PHASE)
89 return __umul64hi(a, b);
90#elif defined(SYCL_DEVICE_ONLY)
91 return cl::sycl::mul_hi(a, static_cast<uint64_t>(b));
92#elif EIGEN_COMP_MSVC && (EIGEN_ARCH_x86_64 || EIGEN_ARCH_ARM64)
93 return __umulh(a, static_cast<uint64_t>(b));
94#elif EIGEN_HAS_BUILTIN_INT128
95 __uint128_t v = static_cast<__uint128_t>(a) * static_cast<__uint128_t>(b);
96 return static_cast<uint64_t>(v >> 64);
97#else
98 return (TensorUInt128<static_val<0>, uint64_t>(a) * TensorUInt128<static_val<0>, uint64_t>(b)).upper();
99#endif
100}
101
102template <int N, typename T>
103struct DividerHelper {
104 static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint32_t computeMultiplier(const int log_div, const T divider) {
105 EIGEN_STATIC_ASSERT(N == 32, YOU_MADE_A_PROGRAMMING_MISTAKE);
106 return static_cast<uint32_t>((static_cast<uint64_t>(1) << (N + log_div)) / divider -
107 (static_cast<uint64_t>(1) << N) + 1);
108 }
109};
110
111template <typename T>
112struct DividerHelper<64, T> {
113 static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint64_t computeMultiplier(const int log_div, const T divider) {
114#if EIGEN_HAS_BUILTIN_INT128 && !defined(EIGEN_GPU_COMPILE_PHASE) && !defined(SYCL_DEVICE_ONLY)
115 return static_cast<uint64_t>((static_cast<__uint128_t>(1) << (64 + log_div)) / static_cast<__uint128_t>(divider) -
116 (static_cast<__uint128_t>(1) << 64) + 1);
117#else
118 const uint64_t shift = 1ULL << log_div;
119 TensorUInt128<uint64_t, uint64_t> result =
120 TensorUInt128<uint64_t, static_val<0> >(shift, 0) / TensorUInt128<static_val<0>, uint64_t>(divider) -
121 TensorUInt128<static_val<1>, static_val<0> >(1, 0) + TensorUInt128<static_val<0>, static_val<1> >(1);
122 return static_cast<uint64_t>(result);
123#endif
124 }
125};
126
138template <typename T, bool div_gt_one = false>
139struct TensorIntDivisor {
140 public:
141 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() = default;
142
143 // Must have 0 < divider < 2^31. This is relaxed to
144 // 0 < divider < 2^63 when using 64-bit indices on platforms that support
145 // the __uint128_t type.
146 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor(const T divider) {
147 const int N = DividerTraits<T>::N;
148 eigen_assert(static_cast<typename UnsignedTraits<T>::type>(divider) < NumTraits<UnsignedType>::highest() / 2);
149 eigen_assert(divider > 0);
150
151 // fast ln2
152 const int leading_zeros = count_leading_zeros(static_cast<UnsignedType>(divider));
153 int log_div = N - leading_zeros;
154 // if divider is a power of two then log_div is 1 more than it should be.
155 if ((static_cast<typename UnsignedTraits<T>::type>(1) << (log_div - 1)) ==
156 static_cast<typename UnsignedTraits<T>::type>(divider))
157 log_div--;
158
159 multiplier = DividerHelper<N, T>::computeMultiplier(log_div, divider);
160 shift1 = log_div > 1 ? 1 : log_div;
161 shift2 = log_div > 1 ? log_div - 1 : 0;
162 }
163
164 // Must have 0 <= numerator. On platforms that don't support the __uint128_t
165 // type numerator should also be less than 2^32-1.
166 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T divide(const T numerator) const {
167 eigen_assert(static_cast<typename UnsignedTraits<T>::type>(numerator) < NumTraits<UnsignedType>::highest() / 2);
168
169 UnsignedType t1 = muluh(multiplier, numerator);
170 UnsignedType t = (static_cast<UnsignedType>(numerator) - t1) >> shift1;
171 return (t1 + t) >> shift2;
172 }
173
174 private:
175 typedef typename DividerTraits<T>::type UnsignedType;
176 UnsignedType multiplier = 0;
177 int32_t shift1 = 0;
178 int32_t shift2 = 0;
179};
180
181// Optimized version for signed 32 bit integers.
182// Derived from Hacker's Delight.
183// Only works for divisors strictly greater than one
184template <>
185class TensorIntDivisor<int32_t, true> {
186 public:
187 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() = default;
188 // Must have 2 <= divider
189 EIGEN_DEVICE_FUNC TensorIntDivisor(int32_t divider) {
190 eigen_assert(divider >= 2);
191 calcMagic(divider);
192 }
193
194 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE int divide(const int32_t n) const {
195#ifdef EIGEN_GPU_COMPILE_PHASE
196 return __umulhi(magic, n) >> shift;
197#elif defined(SYCL_DEVICE_ONLY)
198 return cl::sycl::mul_hi(magic, static_cast<uint32_t>(n)) >> shift;
199#else
200 uint64_t v = static_cast<uint64_t>(magic) * static_cast<uint64_t>(n);
201 return static_cast<uint32_t>(v >> 32) >> shift;
202#endif
203 }
204
205 private:
206 // Compute the magic numbers. See Hacker's Delight section 10 for an in-depth
207 // explanation.
208 EIGEN_DEVICE_FUNC void calcMagic(int32_t d) {
209 const unsigned two31 = 0x80000000; // 2**31.
210 unsigned ad = d;
211 unsigned t = two31 + (ad >> 31);
212 unsigned anc = t - 1 - t % ad; // Absolute value of nc.
213 int p = 31; // Init. p.
214 unsigned q1 = two31 / anc; // Init. q1 = 2**p/|nc|.
215 unsigned r1 = two31 - q1 * anc; // Init. r1 = rem(2**p, |nc|).
216 unsigned q2 = two31 / ad; // Init. q2 = 2**p/|d|.
217 unsigned r2 = two31 - q2 * ad; // Init. r2 = rem(2**p, |d|).
218 unsigned delta = 0;
219 do {
220 p = p + 1;
221 q1 = 2 * q1; // Update q1 = 2**p/|nc|.
222 r1 = 2 * r1; // Update r1 = rem(2**p, |nc|).
223 if (r1 >= anc) { // (Must be an unsigned
224 q1 = q1 + 1; // comparison here).
225 r1 = r1 - anc;
226 }
227 q2 = 2 * q2; // Update q2 = 2**p/|d|.
228 r2 = 2 * r2; // Update r2 = rem(2**p, |d|).
229 if (r2 >= ad) { // (Must be an unsigned
230 q2 = q2 + 1; // comparison here).
231 r2 = r2 - ad;
232 }
233 delta = ad - r2;
234 } while (q1 < delta || (q1 == delta && r1 == 0));
235
236 magic = (unsigned)(q2 + 1);
237 shift = p - 32;
238 }
239
240 uint32_t magic = 0;
241 int32_t shift = 0;
242};
243
244template <typename T, bool div_gt_one>
245EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator/(const T& numerator, const TensorIntDivisor<T, div_gt_one>& divisor) {
246 return divisor.divide(numerator);
247}
248
249} // end namespace internal
250} // end namespace Eigen
251
252#endif // EIGEN_TENSOR_TENSOR_INTDIV_H
Eigen
Namespace containing all symbols from the Eigen library.