1 // libdivide.h
2 // Copyright 2010 - 2018 ridiculous_fish
3 //
4 // libdivide is dual-licensed under the Boost or zlib licenses.
5 // You may use libdivide under the terms of either of these.
6 // See LICENSE.txt for more details.
7
8 #ifndef LIBDIVIDE_H
9 #define LIBDIVIDE_H
10
11 #if defined(_MSC_VER)
12 // disable warning C4146: unary minus operator applied to
13 // unsigned type, result still unsigned
14 #pragma warning(disable: 4146)
15 #define LIBDIVIDE_VC
16 #endif
17
18 #ifdef __cplusplus
19 #include <cstdlib>
20 #include <cstdio>
21 #else
22 #include <stdlib.h>
23 #include <stdio.h>
24 #endif
25
26 #include <stdint.h>
27
28 #if defined(LIBDIVIDE_USE_SSE2)
29 #include <emmintrin.h>
30 #endif
31
32 #if defined(LIBDIVIDE_VC)
33 #include <intrin.h>
34 #endif
35
36 #ifndef __has_builtin
37 #define __has_builtin(x) 0 // Compatibility with non-clang compilers.
38 #endif
39
40 #if defined(__SIZEOF_INT128__)
41 #define HAS_INT128_T
42 #endif
43
44 #if defined(__x86_64__) || defined(_WIN64) || defined(_M_X64)
45 #define LIBDIVIDE_IS_X86_64
46 #endif
47
48 #if defined(__i386__)
49 #define LIBDIVIDE_IS_i386
50 #endif
51
52 #if defined(__GNUC__) || defined(__clang__)
53 #define LIBDIVIDE_GCC_STYLE_ASM
54 #endif
55
56 #if defined(__cplusplus) || defined(LIBDIVIDE_VC)
57 #define LIBDIVIDE_FUNCTION __FUNCTION__
58 #else
59 #define LIBDIVIDE_FUNCTION __func__
60 #endif
61
62 #define LIBDIVIDE_ERROR(msg) \
63 do { \
64 fprintf(stderr, "libdivide.h:%d: %s(): Error: %s\n", \
65 __LINE__, LIBDIVIDE_FUNCTION, msg); \
66 exit(-1); \
67 } while (0)
68
69 #if defined(LIBDIVIDE_ASSERTIONS_ON)
70 #define LIBDIVIDE_ASSERT(x) \
71 do { \
72 if (!(x)) { \
73 fprintf(stderr, "libdivide.h:%d: %s(): Assertion failed: %s\n", \
74 __LINE__, LIBDIVIDE_FUNCTION, #x); \
75 exit(-1); \
76 } \
77 } while (0)
78 #else
79 #define LIBDIVIDE_ASSERT(x)
80 #endif
81
82 // libdivide may use the pmuldq (vector signed 32x32->64 mult instruction)
83 // which is in SSE 4.1. However, signed multiplication can be emulated
84 // efficiently with unsigned multiplication, and SSE 4.1 is currently rare, so
85 // it is OK to not turn this on.
86 #ifdef LIBDIVIDE_USE_SSE4_1
87 #include <smmintrin.h>
88 #endif
89
90 #ifdef __cplusplus
91 // We place libdivide within the libdivide namespace, and that goes in an
92 // anonymous namespace so that the functions are only visible to files that
93 // #include this header and don't get external linkage. At least that's the
94 // theory.
95 namespace {
96 namespace libdivide {
97 #endif
98
99 // Explanation of "more" field: bit 6 is whether to use shift path. If we are
100 // using the shift path, bit 7 is whether the divisor is negative in the signed
101 // case; in the unsigned case it is 0. Bits 0-4 is shift value (for shift
102 // path or mult path). In 32 bit case, bit 5 is always 0. We use bit 7 as the
103 // "negative divisor indicator" so that we can use sign extension to
104 // efficiently go to a full-width -1.
105 //
106 // u32: [0-4] shift value
107 // [5] ignored
108 // [6] add indicator
109 // [7] shift path
110 //
111 // s32: [0-4] shift value
112 // [5] shift path
113 // [6] add indicator
114 // [7] indicates negative divisor
115 //
116 // u64: [0-5] shift value
117 // [6] add indicator
118 // [7] shift path
119 //
120 // s64: [0-5] shift value
121 // [6] add indicator
122 // [7] indicates negative divisor
123 // magic number of 0 indicates shift path (we ran out of bits!)
124 //
125 // In s32 and s64 branchfree modes, the magic number is negated according to
126 // whether the divisor is negated. In branchfree strategy, it is not negated.
127
128 enum {
129 LIBDIVIDE_32_SHIFT_MASK = 0x1F,
130 LIBDIVIDE_64_SHIFT_MASK = 0x3F,
131 LIBDIVIDE_ADD_MARKER = 0x40,
132 LIBDIVIDE_U32_SHIFT_PATH = 0x80,
133 LIBDIVIDE_U64_SHIFT_PATH = 0x80,
134 LIBDIVIDE_S32_SHIFT_PATH = 0x20,
135 LIBDIVIDE_NEGATIVE_DIVISOR = 0x80
136 };
137
138 // pack divider structs to prevent compilers from padding.
139 // This reduces memory usage by up to 43% when using a large
140 // array of libdivide dividers and improves performance
141 // by up to 10% because of reduced memory bandwidth.
142 #pragma pack(push, 1)
143
144 struct libdivide_u32_t {
145 uint32_t magic;
146 uint8_t more;
147 };
148
149 struct libdivide_s32_t {
150 int32_t magic;
151 uint8_t more;
152 };
153
154 struct libdivide_u64_t {
155 uint64_t magic;
156 uint8_t more;
157 };
158
159 struct libdivide_s64_t {
160 int64_t magic;
161 uint8_t more;
162 };
163
164 struct libdivide_u32_branchfree_t {
165 uint32_t magic;
166 uint8_t more;
167 };
168
169 struct libdivide_s32_branchfree_t {
170 int32_t magic;
171 uint8_t more;
172 };
173
174 struct libdivide_u64_branchfree_t {
175 uint64_t magic;
176 uint8_t more;
177 };
178
179 struct libdivide_s64_branchfree_t {
180 int64_t magic;
181 uint8_t more;
182 };
183
184 #pragma pack(pop)
185
186 #ifndef LIBDIVIDE_API
187 #ifdef __cplusplus
188 // In C++, we don't want our public functions to be static, because
189 // they are arguments to templates and static functions can't do that.
190 // They get internal linkage through virtue of the anonymous namespace.
191 // In C, they should be static.
192 #define LIBDIVIDE_API
193 #else
194 #define LIBDIVIDE_API static inline
195 #endif
196 #endif
197
198 LIBDIVIDE_API struct libdivide_s32_t libdivide_s32_gen(int32_t y);
199 LIBDIVIDE_API struct libdivide_u32_t libdivide_u32_gen(uint32_t y);
200 LIBDIVIDE_API struct libdivide_s64_t libdivide_s64_gen(int64_t y);
201 LIBDIVIDE_API struct libdivide_u64_t libdivide_u64_gen(uint64_t y);
202
203 LIBDIVIDE_API struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t y);
204 LIBDIVIDE_API struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t y);
205 LIBDIVIDE_API struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t y);
206 LIBDIVIDE_API struct libdivide_u64_branchfree_t libdivide_u64_branchfree_gen(uint64_t y);
207
208 LIBDIVIDE_API int32_t libdivide_s32_do(int32_t numer, const struct libdivide_s32_t *denom);
209 LIBDIVIDE_API uint32_t libdivide_u32_do(uint32_t numer, const struct libdivide_u32_t *denom);
210 LIBDIVIDE_API int64_t libdivide_s64_do(int64_t numer, const struct libdivide_s64_t *denom);
211 LIBDIVIDE_API uint64_t libdivide_u64_do(uint64_t y, const struct libdivide_u64_t *denom);
212
213 LIBDIVIDE_API int32_t libdivide_s32_branchfree_do(int32_t numer, const struct libdivide_s32_branchfree_t *denom);
214 LIBDIVIDE_API uint32_t libdivide_u32_branchfree_do(uint32_t numer, const struct libdivide_u32_branchfree_t *denom);
215 LIBDIVIDE_API int64_t libdivide_s64_branchfree_do(int64_t numer, const struct libdivide_s64_branchfree_t *denom);
216 LIBDIVIDE_API uint64_t libdivide_u64_branchfree_do(uint64_t y, const struct libdivide_u64_branchfree_t *denom);
217
218 LIBDIVIDE_API int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom);
219 LIBDIVIDE_API uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom);
220 LIBDIVIDE_API int64_t libdivide_s64_recover(const struct libdivide_s64_t *denom);
221 LIBDIVIDE_API uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom);
222
223 LIBDIVIDE_API int32_t libdivide_s32_branchfree_recover(const struct libdivide_s32_branchfree_t *denom);
224 LIBDIVIDE_API uint32_t libdivide_u32_branchfree_recover(const struct libdivide_u32_branchfree_t *denom);
225 LIBDIVIDE_API int64_t libdivide_s64_branchfree_recover(const struct libdivide_s64_branchfree_t *denom);
226 LIBDIVIDE_API uint64_t libdivide_u64_branchfree_recover(const struct libdivide_u64_branchfree_t *denom);
227
228 LIBDIVIDE_API int libdivide_u32_get_algorithm(const struct libdivide_u32_t *denom);
229 LIBDIVIDE_API uint32_t libdivide_u32_do_alg0(uint32_t numer, const struct libdivide_u32_t *denom);
230 LIBDIVIDE_API uint32_t libdivide_u32_do_alg1(uint32_t numer, const struct libdivide_u32_t *denom);
231 LIBDIVIDE_API uint32_t libdivide_u32_do_alg2(uint32_t numer, const struct libdivide_u32_t *denom);
232
233 LIBDIVIDE_API int libdivide_u64_get_algorithm(const struct libdivide_u64_t *denom);
234 LIBDIVIDE_API uint64_t libdivide_u64_do_alg0(uint64_t numer, const struct libdivide_u64_t *denom);
235 LIBDIVIDE_API uint64_t libdivide_u64_do_alg1(uint64_t numer, const struct libdivide_u64_t *denom);
236 LIBDIVIDE_API uint64_t libdivide_u64_do_alg2(uint64_t numer, const struct libdivide_u64_t *denom);
237
238 LIBDIVIDE_API int libdivide_s32_get_algorithm(const struct libdivide_s32_t *denom);
239 LIBDIVIDE_API int32_t libdivide_s32_do_alg0(int32_t numer, const struct libdivide_s32_t *denom);
240 LIBDIVIDE_API int32_t libdivide_s32_do_alg1(int32_t numer, const struct libdivide_s32_t *denom);
241 LIBDIVIDE_API int32_t libdivide_s32_do_alg2(int32_t numer, const struct libdivide_s32_t *denom);
242 LIBDIVIDE_API int32_t libdivide_s32_do_alg3(int32_t numer, const struct libdivide_s32_t *denom);
243 LIBDIVIDE_API int32_t libdivide_s32_do_alg4(int32_t numer, const struct libdivide_s32_t *denom);
244
245 LIBDIVIDE_API int libdivide_s64_get_algorithm(const struct libdivide_s64_t *denom);
246 LIBDIVIDE_API int64_t libdivide_s64_do_alg0(int64_t numer, const struct libdivide_s64_t *denom);
247 LIBDIVIDE_API int64_t libdivide_s64_do_alg1(int64_t numer, const struct libdivide_s64_t *denom);
248 LIBDIVIDE_API int64_t libdivide_s64_do_alg2(int64_t numer, const struct libdivide_s64_t *denom);
249 LIBDIVIDE_API int64_t libdivide_s64_do_alg3(int64_t numer, const struct libdivide_s64_t *denom);
250 LIBDIVIDE_API int64_t libdivide_s64_do_alg4(int64_t numer, const struct libdivide_s64_t *denom);
251
252 #if defined(LIBDIVIDE_USE_SSE2)
253
254 LIBDIVIDE_API __m128i libdivide_u32_do_vector(__m128i numers, const struct libdivide_u32_t *denom);
255 LIBDIVIDE_API __m128i libdivide_s32_do_vector(__m128i numers, const struct libdivide_s32_t *denom);
256 LIBDIVIDE_API __m128i libdivide_u64_do_vector(__m128i numers, const struct libdivide_u64_t *denom);
257 LIBDIVIDE_API __m128i libdivide_s64_do_vector(__m128i numers, const struct libdivide_s64_t *denom);
258
259 LIBDIVIDE_API __m128i libdivide_u32_do_vector_alg0(__m128i numers, const struct libdivide_u32_t *denom);
260 LIBDIVIDE_API __m128i libdivide_u32_do_vector_alg1(__m128i numers, const struct libdivide_u32_t *denom);
261 LIBDIVIDE_API __m128i libdivide_u32_do_vector_alg2(__m128i numers, const struct libdivide_u32_t *denom);
262
263 LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg0(__m128i numers, const struct libdivide_s32_t *denom);
264 LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg1(__m128i numers, const struct libdivide_s32_t *denom);
265 LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg2(__m128i numers, const struct libdivide_s32_t *denom);
266 LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg3(__m128i numers, const struct libdivide_s32_t *denom);
267 LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg4(__m128i numers, const struct libdivide_s32_t *denom);
268
269 LIBDIVIDE_API __m128i libdivide_u64_do_vector_alg0(__m128i numers, const struct libdivide_u64_t *denom);
270 LIBDIVIDE_API __m128i libdivide_u64_do_vector_alg1(__m128i numers, const struct libdivide_u64_t *denom);
271 LIBDIVIDE_API __m128i libdivide_u64_do_vector_alg2(__m128i numers, const struct libdivide_u64_t *denom);
272
273 LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg0(__m128i numers, const struct libdivide_s64_t *denom);
274 LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg1(__m128i numers, const struct libdivide_s64_t *denom);
275 LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg2(__m128i numers, const struct libdivide_s64_t *denom);
276 LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg3(__m128i numers, const struct libdivide_s64_t *denom);
277 LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg4(__m128i numers, const struct libdivide_s64_t *denom);
278
279 LIBDIVIDE_API __m128i libdivide_u32_branchfree_do_vector(__m128i numers, const struct libdivide_u32_branchfree_t *denom);
280 LIBDIVIDE_API __m128i libdivide_s32_branchfree_do_vector(__m128i numers, const struct libdivide_s32_branchfree_t *denom);
281 LIBDIVIDE_API __m128i libdivide_u64_branchfree_do_vector(__m128i numers, const struct libdivide_u64_branchfree_t *denom);
282 LIBDIVIDE_API __m128i libdivide_s64_branchfree_do_vector(__m128i numers, const struct libdivide_s64_branchfree_t *denom);
283
284 #endif
285
286 //////// Internal Utility Functions
287
libdivide__mullhi_u32(uint32_t x,uint32_t y)288 static inline uint32_t libdivide__mullhi_u32(uint32_t x, uint32_t y) {
289 uint64_t xl = x, yl = y;
290 uint64_t rl = xl * yl;
291 return (uint32_t)(rl >> 32);
292 }
293
libdivide__mullhi_u64(uint64_t x,uint64_t y)294 static uint64_t libdivide__mullhi_u64(uint64_t x, uint64_t y) {
295 #if defined(LIBDIVIDE_VC) && defined(LIBDIVIDE_IS_X86_64)
296 return __umulh(x, y);
297 #elif defined(HAS_INT128_T)
298 __uint128_t xl = x, yl = y;
299 __uint128_t rl = xl * yl;
300 return (uint64_t)(rl >> 64);
301 #else
302 // full 128 bits are x0 * y0 + (x0 * y1 << 32) + (x1 * y0 << 32) + (x1 * y1 << 64)
303 uint32_t mask = 0xFFFFFFFF;
304 uint32_t x0 = (uint32_t)(x & mask);
305 uint32_t x1 = (uint32_t)(x >> 32);
306 uint32_t y0 = (uint32_t)(y & mask);
307 uint32_t y1 = (uint32_t)(y >> 32);
308 uint32_t x0y0_hi = libdivide__mullhi_u32(x0, y0);
309 uint64_t x0y1 = x0 * (uint64_t)y1;
310 uint64_t x1y0 = x1 * (uint64_t)y0;
311 uint64_t x1y1 = x1 * (uint64_t)y1;
312 uint64_t temp = x1y0 + x0y0_hi;
313 uint64_t temp_lo = temp & mask;
314 uint64_t temp_hi = temp >> 32;
315
316 return x1y1 + temp_hi + ((temp_lo + x0y1) >> 32);
317 #endif
318 }
319
libdivide__mullhi_s64(int64_t x,int64_t y)320 static inline int64_t libdivide__mullhi_s64(int64_t x, int64_t y) {
321 #if defined(LIBDIVIDE_VC) && defined(LIBDIVIDE_IS_X86_64)
322 return __mulh(x, y);
323 #elif defined(HAS_INT128_T)
324 __int128_t xl = x, yl = y;
325 __int128_t rl = xl * yl;
326 return (int64_t)(rl >> 64);
327 #else
328 // full 128 bits are x0 * y0 + (x0 * y1 << 32) + (x1 * y0 << 32) + (x1 * y1 << 64)
329 uint32_t mask = 0xFFFFFFFF;
330 uint32_t x0 = (uint32_t)(x & mask);
331 uint32_t y0 = (uint32_t)(y & mask);
332 int32_t x1 = (int32_t)(x >> 32);
333 int32_t y1 = (int32_t)(y >> 32);
334 uint32_t x0y0_hi = libdivide__mullhi_u32(x0, y0);
335 int64_t t = x1 * (int64_t)y0 + x0y0_hi;
336 int64_t w1 = x0 * (int64_t)y1 + (t & mask);
337
338 return x1 * (int64_t)y1 + (t >> 32) + (w1 >> 32);
339 #endif
340 }
341
342 #if defined(LIBDIVIDE_USE_SSE2)
343
libdivide__u64_to_m128(uint64_t x)344 static inline __m128i libdivide__u64_to_m128(uint64_t x) {
345 #if defined(LIBDIVIDE_VC) && !defined(_WIN64)
346 // 64 bit windows doesn't seem to have an implementation of any of these
347 // load intrinsics, and 32 bit Visual C++ crashes
348 _declspec(align(16)) uint64_t temp[2] = {x, x};
349 return _mm_load_si128((const __m128i*)temp);
350 #else
351 // everyone else gets it right
352 return _mm_set1_epi64x(x);
353 #endif
354 }
355
libdivide_get_FFFFFFFF00000000(void)356 static inline __m128i libdivide_get_FFFFFFFF00000000(void) {
357 // returns the same as _mm_set1_epi64(0xFFFFFFFF00000000ULL)
358 // without touching memory.
359 // optimizes to pcmpeqd on OS X
360 __m128i result = _mm_set1_epi8(-1);
361 return _mm_slli_epi64(result, 32);
362 }
363
libdivide_get_00000000FFFFFFFF(void)364 static inline __m128i libdivide_get_00000000FFFFFFFF(void) {
365 // returns the same as _mm_set1_epi64(0x00000000FFFFFFFFULL)
366 // without touching memory.
367 // optimizes to pcmpeqd on OS X
368 __m128i result = _mm_set1_epi8(-1);
369 result = _mm_srli_epi64(result, 32);
370 return result;
371 }
372
libdivide_s64_signbits(__m128i v)373 static inline __m128i libdivide_s64_signbits(__m128i v) {
374 // we want to compute v >> 63, that is, _mm_srai_epi64(v, 63). But there
375 // is no 64 bit shift right arithmetic instruction in SSE2. So we have to
376 // fake it by first duplicating the high 32 bit values, and then using a 32
377 // bit shift. Another option would be to use _mm_srli_epi64(v, 63) and
378 // then subtract that from 0, but that approach appears to be substantially
379 // slower for unknown reasons
380 __m128i hiBitsDuped = _mm_shuffle_epi32(v, _MM_SHUFFLE(3, 3, 1, 1));
381 __m128i signBits = _mm_srai_epi32(hiBitsDuped, 31);
382 return signBits;
383 }
384
385 // Returns an __m128i whose low 32 bits are equal to amt and has zero elsewhere.
libdivide_u32_to_m128i(uint32_t amt)386 static inline __m128i libdivide_u32_to_m128i(uint32_t amt) {
387 return _mm_set_epi32(0, 0, 0, amt);
388 }
389
libdivide_s64_shift_right_vector(__m128i v,int amt)390 static inline __m128i libdivide_s64_shift_right_vector(__m128i v, int amt) {
391 // implementation of _mm_sra_epi64. Here we have two 64 bit values which
392 // are shifted right to logically become (64 - amt) values, and are then
393 // sign extended from a (64 - amt) bit number.
394 const int b = 64 - amt;
395 __m128i m = libdivide__u64_to_m128(1ULL << (b - 1));
396 __m128i x = _mm_srl_epi64(v, libdivide_u32_to_m128i(amt));
397 __m128i result = _mm_sub_epi64(_mm_xor_si128(x, m), m); // result = x^m - m
398 return result;
399 }
400
401 // Here, b is assumed to contain one 32 bit value repeated four times.
402 // If it did not, the function would not work.
libdivide__mullhi_u32_flat_vector(__m128i a,__m128i b)403 static inline __m128i libdivide__mullhi_u32_flat_vector(__m128i a, __m128i b) {
404 __m128i hi_product_0Z2Z = _mm_srli_epi64(_mm_mul_epu32(a, b), 32);
405 __m128i a1X3X = _mm_srli_epi64(a, 32);
406 __m128i mask = libdivide_get_FFFFFFFF00000000();
407 __m128i hi_product_Z1Z3 = _mm_and_si128(_mm_mul_epu32(a1X3X, b), mask);
408 return _mm_or_si128(hi_product_0Z2Z, hi_product_Z1Z3); // = hi_product_0123
409 }
410
411 // Here, y is assumed to contain one 64 bit value repeated twice.
libdivide_mullhi_u64_flat_vector(__m128i x,__m128i y)412 static inline __m128i libdivide_mullhi_u64_flat_vector(__m128i x, __m128i y) {
413 // full 128 bits are x0 * y0 + (x0 * y1 << 32) + (x1 * y0 << 32) + (x1 * y1 << 64)
414 __m128i mask = libdivide_get_00000000FFFFFFFF();
415 // x0 is low half of 2 64 bit values, x1 is high half in low slots
416 __m128i x0 = _mm_and_si128(x, mask);
417 __m128i x1 = _mm_srli_epi64(x, 32);
418 __m128i y0 = _mm_and_si128(y, mask);
419 __m128i y1 = _mm_srli_epi64(y, 32);
420 // x0 happens to have the low half of the two 64 bit values in 32 bit slots
421 // 0 and 2, so _mm_mul_epu32 computes their full product, and then we shift
422 // right by 32 to get just the high values
423 __m128i x0y0_hi = _mm_srli_epi64(_mm_mul_epu32(x0, y0), 32);
424 __m128i x0y1 = _mm_mul_epu32(x0, y1);
425 __m128i x1y0 = _mm_mul_epu32(x1, y0);
426 __m128i x1y1 = _mm_mul_epu32(x1, y1);
427 __m128i temp = _mm_add_epi64(x1y0, x0y0_hi);
428 __m128i temp_lo = _mm_and_si128(temp, mask);
429 __m128i temp_hi = _mm_srli_epi64(temp, 32);
430 temp_lo = _mm_srli_epi64(_mm_add_epi64(temp_lo, x0y1), 32);
431 temp_hi = _mm_add_epi64(x1y1, temp_hi);
432
433 return _mm_add_epi64(temp_lo, temp_hi);
434 }
435
436 // y is one 64 bit value repeated twice
libdivide_mullhi_s64_flat_vector(__m128i x,__m128i y)437 static inline __m128i libdivide_mullhi_s64_flat_vector(__m128i x, __m128i y) {
438 __m128i p = libdivide_mullhi_u64_flat_vector(x, y);
439 __m128i t1 = _mm_and_si128(libdivide_s64_signbits(x), y);
440 p = _mm_sub_epi64(p, t1);
441 __m128i t2 = _mm_and_si128(libdivide_s64_signbits(y), x);
442 p = _mm_sub_epi64(p, t2);
443 return p;
444 }
445
446 #ifdef LIBDIVIDE_USE_SSE4_1
447
448 // b is one 32 bit value repeated four times.
libdivide_mullhi_s32_flat_vector(__m128i a,__m128i b)449 static inline __m128i libdivide_mullhi_s32_flat_vector(__m128i a, __m128i b) {
450 __m128i hi_product_0Z2Z = _mm_srli_epi64(_mm_mul_epi32(a, b), 32);
451 __m128i a1X3X = _mm_srli_epi64(a, 32);
452 __m128i mask = libdivide_get_FFFFFFFF00000000();
453 __m128i hi_product_Z1Z3 = _mm_and_si128(_mm_mul_epi32(a1X3X, b), mask);
454 return _mm_or_si128(hi_product_0Z2Z, hi_product_Z1Z3); // = hi_product_0123
455 }
456
457 #else
458
459 // SSE2 does not have a signed multiplication instruction, but we can convert
460 // unsigned to signed pretty efficiently. Again, b is just a 32 bit value
461 // repeated four times.
libdivide_mullhi_s32_flat_vector(__m128i a,__m128i b)462 static inline __m128i libdivide_mullhi_s32_flat_vector(__m128i a, __m128i b) {
463 __m128i p = libdivide__mullhi_u32_flat_vector(a, b);
464 __m128i t1 = _mm_and_si128(_mm_srai_epi32(a, 31), b); // t1 = (a >> 31) & y, arithmetic shift
465 __m128i t2 = _mm_and_si128(_mm_srai_epi32(b, 31), a);
466 p = _mm_sub_epi32(p, t1);
467 p = _mm_sub_epi32(p, t2);
468 return p;
469 }
470
471 #endif // LIBDIVIDE_USE_SSE4_1
472
473 #endif // LIBDIVIDE_USE_SSE2
474
libdivide__count_leading_zeros32(uint32_t val)475 static inline int32_t libdivide__count_leading_zeros32(uint32_t val) {
476 #if defined(__GNUC__) || __has_builtin(__builtin_clz)
477 // Fast way to count leading zeros
478 return __builtin_clz(val);
479 #elif defined(LIBDIVIDE_VC)
480 unsigned long result;
481 if (_BitScanReverse(&result, val)) {
482 return 31 - result;
483 }
484 return 0;
485 #else
486 int32_t result = 0;
487 uint32_t hi = 1U << 31;
488
489 while (~val & hi) {
490 hi >>= 1;
491 result++;
492 }
493 return result;
494 #endif
495 }
496
libdivide__count_leading_zeros64(uint64_t val)497 static inline int32_t libdivide__count_leading_zeros64(uint64_t val) {
498 #if defined(__GNUC__) || __has_builtin(__builtin_clzll)
499 // Fast way to count leading zeros
500 return __builtin_clzll(val);
501 #elif defined(LIBDIVIDE_VC) && defined(_WIN64)
502 unsigned long result;
503 if (_BitScanReverse64(&result, val)) {
504 return 63 - result;
505 }
506 return 0;
507 #else
508 uint32_t hi = val >> 32;
509 uint32_t lo = val & 0xFFFFFFFF;
510 if (hi != 0) return libdivide__count_leading_zeros32(hi);
511 return 32 + libdivide__count_leading_zeros32(lo);
512 #endif
513 }
514
515 #if (defined(LIBDIVIDE_IS_i386) || defined(LIBDIVIDE_IS_X86_64)) && \
516 defined(LIBDIVIDE_GCC_STYLE_ASM)
517
518 // libdivide_64_div_32_to_32: divides a 64 bit uint {u1, u0} by a 32 bit
519 // uint {v}. The result must fit in 32 bits.
520 // Returns the quotient directly and the remainder in *r
libdivide_64_div_32_to_32(uint32_t u1,uint32_t u0,uint32_t v,uint32_t * r)521 static uint32_t libdivide_64_div_32_to_32(uint32_t u1, uint32_t u0, uint32_t v, uint32_t *r) {
522 uint32_t result;
523 __asm__("divl %[v]"
524 : "=a"(result), "=d"(*r)
525 : [v] "r"(v), "a"(u0), "d"(u1)
526 );
527 return result;
528 }
529
530 #else
531
libdivide_64_div_32_to_32(uint32_t u1,uint32_t u0,uint32_t v,uint32_t * r)532 static uint32_t libdivide_64_div_32_to_32(uint32_t u1, uint32_t u0, uint32_t v, uint32_t *r) {
533 uint64_t n = (((uint64_t)u1) << 32) | u0;
534 uint32_t result = (uint32_t)(n / v);
535 *r = (uint32_t)(n - result * (uint64_t)v);
536 return result;
537 }
538
539 #endif
540
541 #if defined(LIBDIVIDE_IS_X86_64) && \
542 defined(LIBDIVIDE_GCC_STYLE_ASM)
543
libdivide_128_div_64_to_64(uint64_t u1,uint64_t u0,uint64_t v,uint64_t * r)544 static uint64_t libdivide_128_div_64_to_64(uint64_t u1, uint64_t u0, uint64_t v, uint64_t *r) {
545 // u0 -> rax
546 // u1 -> rdx
547 // divq
548 uint64_t result;
549 __asm__("divq %[v]"
550 : "=a"(result), "=d"(*r)
551 : [v] "r"(v), "a"(u0), "d"(u1)
552 );
553 return result;
554 }
555
556 #else
557
558 // Code taken from Hacker's Delight:
559 // http://www.hackersdelight.org/HDcode/divlu.c.
560 // License permits inclusion here per:
561 // http://www.hackersdelight.org/permissions.htm
562
libdivide_128_div_64_to_64(uint64_t u1,uint64_t u0,uint64_t v,uint64_t * r)563 static uint64_t libdivide_128_div_64_to_64(uint64_t u1, uint64_t u0, uint64_t v, uint64_t *r) {
564 const uint64_t b = (1ULL << 32); // Number base (16 bits)
565 uint64_t un1, un0; // Norm. dividend LSD's
566 uint64_t vn1, vn0; // Norm. divisor digits
567 uint64_t q1, q0; // Quotient digits
568 uint64_t un64, un21, un10; // Dividend digit pairs
569 uint64_t rhat; // A remainder
570 int s; // Shift amount for norm
571
572 // If overflow, set rem. to an impossible value,
573 // and return the largest possible quotient
574 if (u1 >= v) {
575 if (r != NULL)
576 *r = (uint64_t) -1;
577 return (uint64_t) -1;
578 }
579
580 // count leading zeros
581 s = libdivide__count_leading_zeros64(v);
582 if (s > 0) {
583 // Normalize divisor
584 v = v << s;
585 un64 = (u1 << s) | ((u0 >> (64 - s)) & (-s >> 31));
586 un10 = u0 << s; // Shift dividend left
587 } else {
588 // Avoid undefined behavior
589 un64 = u1 | u0;
590 un10 = u0;
591 }
592
593 // Break divisor up into two 32-bit digits
594 vn1 = v >> 32;
595 vn0 = v & 0xFFFFFFFF;
596
597 // Break right half of dividend into two digits
598 un1 = un10 >> 32;
599 un0 = un10 & 0xFFFFFFFF;
600
601 // Compute the first quotient digit, q1
602 q1 = un64 / vn1;
603 rhat = un64 - q1 * vn1;
604
605 while (q1 >= b || q1 * vn0 > b * rhat + un1) {
606 q1 = q1 - 1;
607 rhat = rhat + vn1;
608 if (rhat >= b)
609 break;
610 }
611
612 // Multiply and subtract
613 un21 = un64 * b + un1 - q1 * v;
614
615 // Compute the second quotient digit
616 q0 = un21 / vn1;
617 rhat = un21 - q0 * vn1;
618
619 while (q0 >= b || q0 * vn0 > b * rhat + un0) {
620 q0 = q0 - 1;
621 rhat = rhat + vn1;
622 if (rhat >= b)
623 break;
624 }
625
626 // If remainder is wanted, return it
627 if (r != NULL)
628 *r = (un21 * b + un0 - q0 * v) >> s;
629
630 return q1 * b + q0;
631 }
632
633 #endif
634
635 // Bitshift a u128 in place, left (signed_shift > 0) or right (signed_shift < 0)
libdivide_u128_shift(uint64_t * u1,uint64_t * u0,int32_t signed_shift)636 static inline void libdivide_u128_shift(uint64_t *u1, uint64_t *u0, int32_t signed_shift)
637 {
638 if (signed_shift > 0) {
639 uint32_t shift = signed_shift;
640 *u1 <<= shift;
641 *u1 |= *u0 >> (64 - shift);
642 *u0 <<= shift;
643 } else {
644 uint32_t shift = -signed_shift;
645 *u0 >>= shift;
646 *u0 |= *u1 << (64 - shift);
647 *u1 >>= shift;
648 }
649 }
650
651 // Computes a 128 / 128 -> 64 bit division, with a 128 bit remainder.
libdivide_128_div_128_to_64(uint64_t u_hi,uint64_t u_lo,uint64_t v_hi,uint64_t v_lo,uint64_t * r_hi,uint64_t * r_lo)652 static uint64_t libdivide_128_div_128_to_64(uint64_t u_hi, uint64_t u_lo, uint64_t v_hi, uint64_t v_lo, uint64_t *r_hi, uint64_t *r_lo) {
653 #if defined(HAS_INT128_T)
654 __uint128_t ufull = u_hi;
655 __uint128_t vfull = v_hi;
656 ufull = (ufull << 64) | u_lo;
657 vfull = (vfull << 64) | v_lo;
658 uint64_t res = (uint64_t)(ufull / vfull);
659 __uint128_t remainder = ufull - (vfull * res);
660 *r_lo = (uint64_t)remainder;
661 *r_hi = (uint64_t)(remainder >> 64);
662 return res;
663 #else
664 // Adapted from "Unsigned Doubleword Division" in Hacker's Delight
665 // We want to compute u / v
666 typedef struct { uint64_t hi; uint64_t lo; } u128_t;
667 u128_t u = {u_hi, u_lo};
668 u128_t v = {v_hi, v_lo};
669
670 if (v.hi == 0) {
671 // divisor v is a 64 bit value, so we just need one 128/64 division
672 // Note that we are simpler than Hacker's Delight here, because we know
673 // the quotient fits in 64 bits whereas Hacker's Delight demands a full
674 // 128 bit quotient
675 *r_hi = 0;
676 return libdivide_128_div_64_to_64(u.hi, u.lo, v.lo, r_lo);
677 }
678 // Here v >= 2**64
679 // We know that v.hi != 0, so count leading zeros is OK
680 // We have 0 <= n <= 63
681 uint32_t n = libdivide__count_leading_zeros64(v.hi);
682
683 // Normalize the divisor so its MSB is 1
684 u128_t v1t = v;
685 libdivide_u128_shift(&v1t.hi, &v1t.lo, n);
686 uint64_t v1 = v1t.hi; // i.e. v1 = v1t >> 64
687
688 // To ensure no overflow
689 u128_t u1 = u;
690 libdivide_u128_shift(&u1.hi, &u1.lo, -1);
691
692 // Get quotient from divide unsigned insn.
693 uint64_t rem_ignored;
694 uint64_t q1 = libdivide_128_div_64_to_64(u1.hi, u1.lo, v1, &rem_ignored);
695
696 // Undo normalization and division of u by 2.
697 u128_t q0 = {0, q1};
698 libdivide_u128_shift(&q0.hi, &q0.lo, n);
699 libdivide_u128_shift(&q0.hi, &q0.lo, -63);
700
701 // Make q0 correct or too small by 1
702 // Equivalent to `if (q0 != 0) q0 = q0 - 1;`
703 if (q0.hi != 0 || q0.lo != 0) {
704 q0.hi -= (q0.lo == 0); // borrow
705 q0.lo -= 1;
706 }
707
708 // Now q0 is correct.
709 // Compute q0 * v as q0v
710 // = (q0.hi << 64 + q0.lo) * (v.hi << 64 + v.lo)
711 // = (q0.hi * v.hi << 128) + (q0.hi * v.lo << 64) +
712 // (q0.lo * v.hi << 64) + q0.lo * v.lo)
713 // Each term is 128 bit
714 // High half of full product (upper 128 bits!) are dropped
715 u128_t q0v = {0, 0};
716 q0v.hi = q0.hi*v.lo + q0.lo*v.hi + libdivide__mullhi_u64(q0.lo, v.lo);
717 q0v.lo = q0.lo*v.lo;
718
719 // Compute u - q0v as u_q0v
720 // This is the remainder
721 u128_t u_q0v = u;
722 u_q0v.hi -= q0v.hi + (u.lo < q0v.lo); // second term is borrow
723 u_q0v.lo -= q0v.lo;
724
725 // Check if u_q0v >= v
726 // This checks if our remainder is larger than the divisor
727 if ((u_q0v.hi > v.hi) ||
728 (u_q0v.hi == v.hi && u_q0v.lo >= v.lo)) {
729 // Increment q0
730 q0.lo += 1;
731 q0.hi += (q0.lo == 0); // carry
732
733 // Subtract v from remainder
734 u_q0v.hi -= v.hi + (u_q0v.lo < v.lo);
735 u_q0v.lo -= v.lo;
736 }
737
738 *r_hi = u_q0v.hi;
739 *r_lo = u_q0v.lo;
740
741 LIBDIVIDE_ASSERT(q0.hi == 0);
742 return q0.lo;
743 #endif
744 }
745
746 ////////// UINT32
747
libdivide_internal_u32_gen(uint32_t d,int branchfree)748 static inline struct libdivide_u32_t libdivide_internal_u32_gen(uint32_t d, int branchfree) {
749 if (d == 0) {
750 LIBDIVIDE_ERROR("divider must be != 0");
751 }
752
753 struct libdivide_u32_t result;
754 uint32_t floor_log_2_d = 31 - libdivide__count_leading_zeros32(d);
755 if ((d & (d - 1)) == 0) {
756 // Power of 2
757 if (! branchfree) {
758 result.magic = 0;
759 result.more = floor_log_2_d | LIBDIVIDE_U32_SHIFT_PATH;
760 } else {
761 // We want a magic number of 2**32 and a shift of floor_log_2_d
762 // but one of the shifts is taken up by LIBDIVIDE_ADD_MARKER,
763 // so we subtract 1 from the shift
764 result.magic = 0;
765 result.more = (floor_log_2_d-1) | LIBDIVIDE_ADD_MARKER;
766 }
767 } else {
768 uint8_t more;
769 uint32_t rem, proposed_m;
770 proposed_m = libdivide_64_div_32_to_32(1U << floor_log_2_d, 0, d, &rem);
771
772 LIBDIVIDE_ASSERT(rem > 0 && rem < d);
773 const uint32_t e = d - rem;
774
775 // This power works if e < 2**floor_log_2_d.
776 if (!branchfree && (e < (1U << floor_log_2_d))) {
777 // This power works
778 more = floor_log_2_d;
779 } else {
780 // We have to use the general 33-bit algorithm. We need to compute
781 // (2**power) / d. However, we already have (2**(power-1))/d and
782 // its remainder. By doubling both, and then correcting the
783 // remainder, we can compute the larger division.
784 // don't care about overflow here - in fact, we expect it
785 proposed_m += proposed_m;
786 const uint32_t twice_rem = rem + rem;
787 if (twice_rem >= d || twice_rem < rem) proposed_m += 1;
788 more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;
789 }
790 result.magic = 1 + proposed_m;
791 result.more = more;
792 // result.more's shift should in general be ceil_log_2_d. But if we
793 // used the smaller power, we subtract one from the shift because we're
794 // using the smaller power. If we're using the larger power, we
795 // subtract one from the shift because it's taken care of by the add
796 // indicator. So floor_log_2_d happens to be correct in both cases.
797 }
798 return result;
799 }
800
libdivide_u32_gen(uint32_t d)801 struct libdivide_u32_t libdivide_u32_gen(uint32_t d) {
802 return libdivide_internal_u32_gen(d, 0);
803 }
804
libdivide_u32_branchfree_gen(uint32_t d)805 struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t d) {
806 if (d == 1) {
807 LIBDIVIDE_ERROR("branchfree divider must be != 1");
808 }
809 struct libdivide_u32_t tmp = libdivide_internal_u32_gen(d, 1);
810 struct libdivide_u32_branchfree_t ret = {tmp.magic, (uint8_t)(tmp.more & LIBDIVIDE_32_SHIFT_MASK)};
811 return ret;
812 }
813
libdivide_u32_do(uint32_t numer,const struct libdivide_u32_t * denom)814 uint32_t libdivide_u32_do(uint32_t numer, const struct libdivide_u32_t *denom) {
815 uint8_t more = denom->more;
816 if (more & LIBDIVIDE_U32_SHIFT_PATH) {
817 return numer >> (more & LIBDIVIDE_32_SHIFT_MASK);
818 }
819 else {
820 uint32_t q = libdivide__mullhi_u32(denom->magic, numer);
821 if (more & LIBDIVIDE_ADD_MARKER) {
822 uint32_t t = ((numer - q) >> 1) + q;
823 return t >> (more & LIBDIVIDE_32_SHIFT_MASK);
824 }
825 else {
826 // all upper bits are 0 - don't need to mask them off
827 return q >> more;
828 }
829 }
830 }
831
libdivide_u32_recover(const struct libdivide_u32_t * denom)832 uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom) {
833 uint8_t more = denom->more;
834 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
835 if (more & LIBDIVIDE_U32_SHIFT_PATH) {
836 return 1U << shift;
837 } else if (!(more & LIBDIVIDE_ADD_MARKER)) {
838 // We compute q = n/d = n*m / 2^(32 + shift)
839 // Therefore we have d = 2^(32 + shift) / m
840 // We need to ceil it.
841 // We know d is not a power of 2, so m is not a power of 2,
842 // so we can just add 1 to the floor
843 uint32_t hi_dividend = 1U << shift;
844 uint32_t rem_ignored;
845 return 1 + libdivide_64_div_32_to_32(hi_dividend, 0, denom->magic, &rem_ignored);
846 } else {
847 // Here we wish to compute d = 2^(32+shift+1)/(m+2^32).
848 // Notice (m + 2^32) is a 33 bit number. Use 64 bit division for now
849 // Also note that shift may be as high as 31, so shift + 1 will
850 // overflow. So we have to compute it as 2^(32+shift)/(m+2^32), and
851 // then double the quotient and remainder.
852 uint64_t half_n = 1ULL << (32 + shift);
853 uint64_t d = (1ULL << 32) | denom->magic;
854 // Note that the quotient is guaranteed <= 32 bits, but the remainder
855 // may need 33!
856 uint32_t half_q = (uint32_t)(half_n / d);
857 uint64_t rem = half_n % d;
858 // We computed 2^(32+shift)/(m+2^32)
859 // Need to double it, and then add 1 to the quotient if doubling th
860 // remainder would increase the quotient.
861 // Note that rem<<1 cannot overflow, since rem < d and d is 33 bits
862 uint32_t full_q = half_q + half_q + ((rem<<1) >= d);
863
864 // We rounded down in gen unless we're a power of 2 (i.e. in branchfree case)
865 // We can detect that by looking at m. If m zero, we're a power of 2
866 return full_q + (denom->magic != 0);
867 }
868 }
869
libdivide_u32_branchfree_recover(const struct libdivide_u32_branchfree_t * denom)870 uint32_t libdivide_u32_branchfree_recover(const struct libdivide_u32_branchfree_t *denom) {
871 struct libdivide_u32_t denom_u32 = {denom->magic, (uint8_t)(denom->more | LIBDIVIDE_ADD_MARKER)};
872 return libdivide_u32_recover(&denom_u32);
873 }
874
libdivide_u32_get_algorithm(const struct libdivide_u32_t * denom)875 int libdivide_u32_get_algorithm(const struct libdivide_u32_t *denom) {
876 uint8_t more = denom->more;
877 if (more & LIBDIVIDE_U32_SHIFT_PATH) return 0;
878 else if (!(more & LIBDIVIDE_ADD_MARKER)) return 1;
879 else return 2;
880 }
881
libdivide_u32_do_alg0(uint32_t numer,const struct libdivide_u32_t * denom)882 uint32_t libdivide_u32_do_alg0(uint32_t numer, const struct libdivide_u32_t *denom) {
883 return numer >> (denom->more & LIBDIVIDE_32_SHIFT_MASK);
884 }
885
libdivide_u32_do_alg1(uint32_t numer,const struct libdivide_u32_t * denom)886 uint32_t libdivide_u32_do_alg1(uint32_t numer, const struct libdivide_u32_t *denom) {
887 uint32_t q = libdivide__mullhi_u32(denom->magic, numer);
888 return q >> denom->more;
889 }
890
libdivide_u32_do_alg2(uint32_t numer,const struct libdivide_u32_t * denom)891 uint32_t libdivide_u32_do_alg2(uint32_t numer, const struct libdivide_u32_t *denom) {
892 // denom->add != 0
893 uint32_t q = libdivide__mullhi_u32(denom->magic, numer);
894 uint32_t t = ((numer - q) >> 1) + q;
895 // Note that this mask is typically free. Only the low bits are meaningful
896 // to a shift, so compilers can optimize out this AND.
897 return t >> (denom->more & LIBDIVIDE_32_SHIFT_MASK);
898 }
899
900 // same as algo 2
libdivide_u32_branchfree_do(uint32_t numer,const struct libdivide_u32_branchfree_t * denom)901 uint32_t libdivide_u32_branchfree_do(uint32_t numer, const struct libdivide_u32_branchfree_t *denom) {
902 uint32_t q = libdivide__mullhi_u32(denom->magic, numer);
903 uint32_t t = ((numer - q) >> 1) + q;
904 return t >> denom->more;
905 }
906
907 #if defined(LIBDIVIDE_USE_SSE2)
908
libdivide_u32_do_vector(__m128i numers,const struct libdivide_u32_t * denom)909 __m128i libdivide_u32_do_vector(__m128i numers, const struct libdivide_u32_t *denom) {
910 uint8_t more = denom->more;
911 if (more & LIBDIVIDE_U32_SHIFT_PATH) {
912 return _mm_srl_epi32(numers, libdivide_u32_to_m128i(more & LIBDIVIDE_32_SHIFT_MASK));
913 }
914 else {
915 __m128i q = libdivide__mullhi_u32_flat_vector(numers, _mm_set1_epi32(denom->magic));
916 if (more & LIBDIVIDE_ADD_MARKER) {
917 // uint32_t t = ((numer - q) >> 1) + q;
918 // return t >> denom->shift;
919 __m128i t = _mm_add_epi32(_mm_srli_epi32(_mm_sub_epi32(numers, q), 1), q);
920 return _mm_srl_epi32(t, libdivide_u32_to_m128i(more & LIBDIVIDE_32_SHIFT_MASK));
921
922 }
923 else {
924 // q >> denom->shift
925 return _mm_srl_epi32(q, libdivide_u32_to_m128i(more));
926 }
927 }
928 }
929
libdivide_u32_do_vector_alg0(__m128i numers,const struct libdivide_u32_t * denom)930 __m128i libdivide_u32_do_vector_alg0(__m128i numers, const struct libdivide_u32_t *denom) {
931 return _mm_srl_epi32(numers, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_32_SHIFT_MASK));
932 }
933
libdivide_u32_do_vector_alg1(__m128i numers,const struct libdivide_u32_t * denom)934 __m128i libdivide_u32_do_vector_alg1(__m128i numers, const struct libdivide_u32_t *denom) {
935 __m128i q = libdivide__mullhi_u32_flat_vector(numers, _mm_set1_epi32(denom->magic));
936 return _mm_srl_epi32(q, libdivide_u32_to_m128i(denom->more));
937 }
938
libdivide_u32_do_vector_alg2(__m128i numers,const struct libdivide_u32_t * denom)939 __m128i libdivide_u32_do_vector_alg2(__m128i numers, const struct libdivide_u32_t *denom) {
940 __m128i q = libdivide__mullhi_u32_flat_vector(numers, _mm_set1_epi32(denom->magic));
941 __m128i t = _mm_add_epi32(_mm_srli_epi32(_mm_sub_epi32(numers, q), 1), q);
942 return _mm_srl_epi32(t, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_32_SHIFT_MASK));
943 }
944
945 // same as algo 2
libdivide_u32_branchfree_do_vector(__m128i numers,const struct libdivide_u32_branchfree_t * denom)946 LIBDIVIDE_API __m128i libdivide_u32_branchfree_do_vector(__m128i numers, const struct libdivide_u32_branchfree_t *denom) {
947 __m128i q = libdivide__mullhi_u32_flat_vector(numers, _mm_set1_epi32(denom->magic));
948 __m128i t = _mm_add_epi32(_mm_srli_epi32(_mm_sub_epi32(numers, q), 1), q);
949 return _mm_srl_epi32(t, libdivide_u32_to_m128i(denom->more));
950 }
951
952 #endif
953
954 /////////// UINT64
955
libdivide_internal_u64_gen(uint64_t d,int branchfree)956 static inline struct libdivide_u64_t libdivide_internal_u64_gen(uint64_t d, int branchfree) {
957 if (d == 0) {
958 LIBDIVIDE_ERROR("divider must be != 0");
959 }
960
961 struct libdivide_u64_t result;
962 uint32_t floor_log_2_d = 63 - libdivide__count_leading_zeros64(d);
963 if ((d & (d - 1)) == 0) {
964 // Power of 2
965 if (! branchfree) {
966 result.magic = 0;
967 result.more = floor_log_2_d | LIBDIVIDE_U64_SHIFT_PATH;
968 } else {
969 // We want a magic number of 2**64 and a shift of floor_log_2_d
970 // but one of the shifts is taken up by LIBDIVIDE_ADD_MARKER,
971 // so we subtract 1 from the shift
972 result.magic = 0;
973 result.more = (floor_log_2_d-1) | LIBDIVIDE_ADD_MARKER;
974 }
975 } else {
976 uint64_t proposed_m, rem;
977 uint8_t more;
978 // (1 << (64 + floor_log_2_d)) / d
979 proposed_m = libdivide_128_div_64_to_64(1ULL << floor_log_2_d, 0, d, &rem);
980
981 LIBDIVIDE_ASSERT(rem > 0 && rem < d);
982 const uint64_t e = d - rem;
983
984 // This power works if e < 2**floor_log_2_d.
985 if (!branchfree && e < (1ULL << floor_log_2_d)) {
986 // This power works
987 more = floor_log_2_d;
988 } else {
989 // We have to use the general 65-bit algorithm. We need to compute
990 // (2**power) / d. However, we already have (2**(power-1))/d and
991 // its remainder. By doubling both, and then correcting the
992 // remainder, we can compute the larger division.
993 // don't care about overflow here - in fact, we expect it
994 proposed_m += proposed_m;
995 const uint64_t twice_rem = rem + rem;
996 if (twice_rem >= d || twice_rem < rem) proposed_m += 1;
997 more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;
998 }
999 result.magic = 1 + proposed_m;
1000 result.more = more;
1001 // result.more's shift should in general be ceil_log_2_d. But if we
1002 // used the smaller power, we subtract one from the shift because we're
1003 // using the smaller power. If we're using the larger power, we
1004 // subtract one from the shift because it's taken care of by the add
1005 // indicator. So floor_log_2_d happens to be correct in both cases,
1006 // which is why we do it outside of the if statement.
1007 }
1008 return result;
1009 }
1010
libdivide_u64_gen(uint64_t d)1011 struct libdivide_u64_t libdivide_u64_gen(uint64_t d) {
1012 return libdivide_internal_u64_gen(d, 0);
1013 }
1014
libdivide_u64_branchfree_gen(uint64_t d)1015 struct libdivide_u64_branchfree_t libdivide_u64_branchfree_gen(uint64_t d) {
1016 if (d == 1) {
1017 LIBDIVIDE_ERROR("branchfree divider must be != 1");
1018 }
1019 struct libdivide_u64_t tmp = libdivide_internal_u64_gen(d, 1);
1020 struct libdivide_u64_branchfree_t ret = {tmp.magic, (uint8_t)(tmp.more & LIBDIVIDE_64_SHIFT_MASK)};
1021 return ret;
1022 }
1023
libdivide_u64_do(uint64_t numer,const struct libdivide_u64_t * denom)1024 uint64_t libdivide_u64_do(uint64_t numer, const struct libdivide_u64_t *denom) {
1025 uint8_t more = denom->more;
1026 if (more & LIBDIVIDE_U64_SHIFT_PATH) {
1027 return numer >> (more & LIBDIVIDE_64_SHIFT_MASK);
1028 }
1029 else {
1030 uint64_t q = libdivide__mullhi_u64(denom->magic, numer);
1031 if (more & LIBDIVIDE_ADD_MARKER) {
1032 uint64_t t = ((numer - q) >> 1) + q;
1033 return t >> (more & LIBDIVIDE_64_SHIFT_MASK);
1034 }
1035 else {
1036 // all upper bits are 0 - don't need to mask them off
1037 return q >> more;
1038 }
1039 }
1040 }
1041
libdivide_u64_recover(const struct libdivide_u64_t * denom)1042 uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom) {
1043 uint8_t more = denom->more;
1044 uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
1045 if (more & LIBDIVIDE_U64_SHIFT_PATH) {
1046 return 1ULL << shift;
1047 } else if (!(more & LIBDIVIDE_ADD_MARKER)) {
1048 // We compute q = n/d = n*m / 2^(64 + shift)
1049 // Therefore we have d = 2^(64 + shift) / m
1050 // We need to ceil it.
1051 // We know d is not a power of 2, so m is not a power of 2,
1052 // so we can just add 1 to the floor
1053 uint64_t hi_dividend = 1ULL << shift;
1054 uint64_t rem_ignored;
1055 return 1 + libdivide_128_div_64_to_64(hi_dividend, 0, denom->magic, &rem_ignored);
1056 } else {
1057 // Here we wish to compute d = 2^(64+shift+1)/(m+2^64).
1058 // Notice (m + 2^64) is a 65 bit number. This gets hairy. See
1059 // libdivide_u32_recover for more on what we do here.
1060 // TODO: do something better than 128 bit math
1061
1062 // Hack: if d is not a power of 2, this is a 128/128->64 divide
1063 // If d is a power of 2, this may be a bigger divide
1064 // However we can optimize that easily
1065 if (denom->magic == 0) {
1066 // 2^(64 + shift + 1) / (2^64) == 2^(shift + 1)
1067 return 1ULL << (shift + 1);
1068 }
1069
1070 // Full n is a (potentially) 129 bit value
1071 // half_n is a 128 bit value
1072 // Compute the hi half of half_n. Low half is 0.
1073 uint64_t half_n_hi = 1ULL << shift, half_n_lo = 0;
1074 // d is a 65 bit value. The high bit is always set to 1.
1075 const uint64_t d_hi = 1, d_lo = denom->magic;
1076 // Note that the quotient is guaranteed <= 64 bits,
1077 // but the remainder may need 65!
1078 uint64_t r_hi, r_lo;
1079 uint64_t half_q = libdivide_128_div_128_to_64(half_n_hi, half_n_lo, d_hi, d_lo, &r_hi, &r_lo);
1080 // We computed 2^(64+shift)/(m+2^64)
1081 // Double the remainder ('dr') and check if that is larger than d
1082 // Note that d is a 65 bit value, so r1 is small and so r1 + r1 cannot
1083 // overflow
1084 uint64_t dr_lo = r_lo + r_lo;
1085 uint64_t dr_hi = r_hi + r_hi + (dr_lo < r_lo); // last term is carry
1086 int dr_exceeds_d = (dr_hi > d_hi) || (dr_hi == d_hi && dr_lo >= d_lo);
1087 uint64_t full_q = half_q + half_q + (dr_exceeds_d ? 1 : 0);
1088 return full_q + 1;
1089 }
1090 }
1091
libdivide_u64_branchfree_recover(const struct libdivide_u64_branchfree_t * denom)1092 uint64_t libdivide_u64_branchfree_recover(const struct libdivide_u64_branchfree_t *denom) {
1093 struct libdivide_u64_t denom_u64 = {denom->magic, (uint8_t)(denom->more | LIBDIVIDE_ADD_MARKER)};
1094 return libdivide_u64_recover(&denom_u64);
1095 }
1096
libdivide_u64_get_algorithm(const struct libdivide_u64_t * denom)1097 int libdivide_u64_get_algorithm(const struct libdivide_u64_t *denom) {
1098 uint8_t more = denom->more;
1099 if (more & LIBDIVIDE_U64_SHIFT_PATH) return 0;
1100 else if (!(more & LIBDIVIDE_ADD_MARKER)) return 1;
1101 else return 2;
1102 }
1103
libdivide_u64_do_alg0(uint64_t numer,const struct libdivide_u64_t * denom)1104 uint64_t libdivide_u64_do_alg0(uint64_t numer, const struct libdivide_u64_t *denom) {
1105 return numer >> (denom->more & LIBDIVIDE_64_SHIFT_MASK);
1106 }
1107
libdivide_u64_do_alg1(uint64_t numer,const struct libdivide_u64_t * denom)1108 uint64_t libdivide_u64_do_alg1(uint64_t numer, const struct libdivide_u64_t *denom) {
1109 uint64_t q = libdivide__mullhi_u64(denom->magic, numer);
1110 return q >> denom->more;
1111 }
1112
libdivide_u64_do_alg2(uint64_t numer,const struct libdivide_u64_t * denom)1113 uint64_t libdivide_u64_do_alg2(uint64_t numer, const struct libdivide_u64_t *denom) {
1114 uint64_t q = libdivide__mullhi_u64(denom->magic, numer);
1115 uint64_t t = ((numer - q) >> 1) + q;
1116 return t >> (denom->more & LIBDIVIDE_64_SHIFT_MASK);
1117 }
1118
1119 // same as alg 2
libdivide_u64_branchfree_do(uint64_t numer,const struct libdivide_u64_branchfree_t * denom)1120 uint64_t libdivide_u64_branchfree_do(uint64_t numer, const struct libdivide_u64_branchfree_t *denom) {
1121 uint64_t q = libdivide__mullhi_u64(denom->magic, numer);
1122 uint64_t t = ((numer - q) >> 1) + q;
1123 return t >> denom->more;
1124 }
1125
1126 #if defined(LIBDIVIDE_USE_SSE2)
1127
libdivide_u64_do_vector(__m128i numers,const struct libdivide_u64_t * denom)1128 __m128i libdivide_u64_do_vector(__m128i numers, const struct libdivide_u64_t *denom) {
1129 uint8_t more = denom->more;
1130 if (more & LIBDIVIDE_U64_SHIFT_PATH) {
1131 return _mm_srl_epi64(numers, libdivide_u32_to_m128i(more & LIBDIVIDE_64_SHIFT_MASK));
1132 }
1133 else {
1134 __m128i q = libdivide_mullhi_u64_flat_vector(numers, libdivide__u64_to_m128(denom->magic));
1135 if (more & LIBDIVIDE_ADD_MARKER) {
1136 // uint32_t t = ((numer - q) >> 1) + q;
1137 // return t >> denom->shift;
1138 __m128i t = _mm_add_epi64(_mm_srli_epi64(_mm_sub_epi64(numers, q), 1), q);
1139 return _mm_srl_epi64(t, libdivide_u32_to_m128i(more & LIBDIVIDE_64_SHIFT_MASK));
1140 }
1141 else {
1142 // q >> denom->shift
1143 return _mm_srl_epi64(q, libdivide_u32_to_m128i(more));
1144 }
1145 }
1146 }
1147
libdivide_u64_do_vector_alg0(__m128i numers,const struct libdivide_u64_t * denom)1148 __m128i libdivide_u64_do_vector_alg0(__m128i numers, const struct libdivide_u64_t *denom) {
1149 return _mm_srl_epi64(numers, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_64_SHIFT_MASK));
1150 }
1151
libdivide_u64_do_vector_alg1(__m128i numers,const struct libdivide_u64_t * denom)1152 __m128i libdivide_u64_do_vector_alg1(__m128i numers, const struct libdivide_u64_t *denom) {
1153 __m128i q = libdivide_mullhi_u64_flat_vector(numers, libdivide__u64_to_m128(denom->magic));
1154 return _mm_srl_epi64(q, libdivide_u32_to_m128i(denom->more));
1155 }
1156
libdivide_u64_do_vector_alg2(__m128i numers,const struct libdivide_u64_t * denom)1157 __m128i libdivide_u64_do_vector_alg2(__m128i numers, const struct libdivide_u64_t *denom) {
1158 __m128i q = libdivide_mullhi_u64_flat_vector(numers, libdivide__u64_to_m128(denom->magic));
1159 __m128i t = _mm_add_epi64(_mm_srli_epi64(_mm_sub_epi64(numers, q), 1), q);
1160 return _mm_srl_epi64(t, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_64_SHIFT_MASK));
1161 }
1162
libdivide_u64_branchfree_do_vector(__m128i numers,const struct libdivide_u64_branchfree_t * denom)1163 __m128i libdivide_u64_branchfree_do_vector(__m128i numers, const struct libdivide_u64_branchfree_t *denom) {
1164 __m128i q = libdivide_mullhi_u64_flat_vector(numers, libdivide__u64_to_m128(denom->magic));
1165 __m128i t = _mm_add_epi64(_mm_srli_epi64(_mm_sub_epi64(numers, q), 1), q);
1166 return _mm_srl_epi64(t, libdivide_u32_to_m128i(denom->more));
1167 }
1168
1169 #endif
1170
1171 /////////// SINT32
1172
libdivide__mullhi_s32(int32_t x,int32_t y)1173 static inline int32_t libdivide__mullhi_s32(int32_t x, int32_t y) {
1174 int64_t xl = x, yl = y;
1175 int64_t rl = xl * yl;
1176 // needs to be arithmetic shift
1177 return (int32_t)(rl >> 32);
1178 }
1179
libdivide_internal_s32_gen(int32_t d,int branchfree)1180 static inline struct libdivide_s32_t libdivide_internal_s32_gen(int32_t d, int branchfree) {
1181 if (d == 0) {
1182 LIBDIVIDE_ERROR("divider must be != 0");
1183 }
1184
1185 struct libdivide_s32_t result;
1186
1187 // If d is a power of 2, or negative a power of 2, we have to use a shift.
1188 // This is especially important because the magic algorithm fails for -1.
1189 // To check if d is a power of 2 or its inverse, it suffices to check
1190 // whether its absolute value has exactly one bit set. This works even for
1191 // INT_MIN, because abs(INT_MIN) == INT_MIN, and INT_MIN has one bit set
1192 // and is a power of 2.
1193 uint32_t ud = (uint32_t)d;
1194 uint32_t absD = (d < 0) ? -ud : ud;
1195 uint32_t floor_log_2_d = 31 - libdivide__count_leading_zeros32(absD);
1196 // check if exactly one bit is set,
1197 // don't care if absD is 0 since that's divide by zero
1198 if ((absD & (absD - 1)) == 0) {
1199 // Branchfree and normal paths are exactly the same
1200 result.magic = 0;
1201 result.more = floor_log_2_d | (d < 0 ? LIBDIVIDE_NEGATIVE_DIVISOR : 0) | LIBDIVIDE_S32_SHIFT_PATH;
1202 } else {
1203 LIBDIVIDE_ASSERT(floor_log_2_d >= 1);
1204
1205 uint8_t more;
1206 // the dividend here is 2**(floor_log_2_d + 31), so the low 32 bit word
1207 // is 0 and the high word is floor_log_2_d - 1
1208 uint32_t rem, proposed_m;
1209 proposed_m = libdivide_64_div_32_to_32(1U << (floor_log_2_d - 1), 0, absD, &rem);
1210 const uint32_t e = absD - rem;
1211
1212 // We are going to start with a power of floor_log_2_d - 1.
1213 // This works if works if e < 2**floor_log_2_d.
1214 if (!branchfree && e < (1U << floor_log_2_d)) {
1215 // This power works
1216 more = floor_log_2_d - 1;
1217 } else {
1218 // We need to go one higher. This should not make proposed_m
1219 // overflow, but it will make it negative when interpreted as an
1220 // int32_t.
1221 proposed_m += proposed_m;
1222 const uint32_t twice_rem = rem + rem;
1223 if (twice_rem >= absD || twice_rem < rem) proposed_m += 1;
1224 more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;
1225 }
1226
1227 proposed_m += 1;
1228 int32_t magic = (int32_t)proposed_m;
1229
1230 // Mark if we are negative. Note we only negate the magic number in the
1231 // branchfull case.
1232 if (d < 0) {
1233 more |= LIBDIVIDE_NEGATIVE_DIVISOR;
1234 if (!branchfree) {
1235 magic = -magic;
1236 }
1237 }
1238
1239 result.more = more;
1240 result.magic = magic;
1241 }
1242 return result;
1243 }
1244
libdivide_s32_gen(int32_t d)1245 LIBDIVIDE_API struct libdivide_s32_t libdivide_s32_gen(int32_t d) {
1246 return libdivide_internal_s32_gen(d, 0);
1247 }
1248
libdivide_s32_branchfree_gen(int32_t d)1249 LIBDIVIDE_API struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t d) {
1250 if (d == 1) {
1251 LIBDIVIDE_ERROR("branchfree divider must be != 1");
1252 }
1253 if (d == -1) {
1254 LIBDIVIDE_ERROR("branchfree divider must be != -1");
1255 }
1256 struct libdivide_s32_t tmp = libdivide_internal_s32_gen(d, 1);
1257 struct libdivide_s32_branchfree_t result = {tmp.magic, tmp.more};
1258 return result;
1259 }
1260
libdivide_s32_do(int32_t numer,const struct libdivide_s32_t * denom)1261 int32_t libdivide_s32_do(int32_t numer, const struct libdivide_s32_t *denom) {
1262 uint8_t more = denom->more;
1263 uint32_t sign = (int8_t)more >> 7;
1264 if (more & LIBDIVIDE_S32_SHIFT_PATH) {
1265 uint8_t shifter = more & LIBDIVIDE_32_SHIFT_MASK;
1266 uint32_t uq = (uint32_t)(numer + ((numer >> 31) & ((1U << shifter) - 1)));
1267 int32_t q = (int32_t)uq;
1268 q = q >> shifter;
1269 q = (q ^ sign) - sign;
1270 return q;
1271 } else {
1272 uint32_t uq = (uint32_t)libdivide__mullhi_s32(denom->magic, numer);
1273 if (more & LIBDIVIDE_ADD_MARKER) {
1274 // must be arithmetic shift and then sign extend
1275 int32_t sign = (int8_t)more >> 7;
1276 // q += (more < 0 ? -numer : numer), casts to avoid UB
1277 uq += ((uint32_t)numer ^ sign) - sign;
1278 }
1279 int32_t q = (int32_t)uq;
1280 q >>= more & LIBDIVIDE_32_SHIFT_MASK;
1281 q += (q < 0);
1282 return q;
1283 }
1284 }
1285
libdivide_s32_branchfree_do(int32_t numer,const struct libdivide_s32_branchfree_t * denom)1286 int32_t libdivide_s32_branchfree_do(int32_t numer, const struct libdivide_s32_branchfree_t *denom) {
1287 uint8_t more = denom->more;
1288 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
1289 // must be arithmetic shift and then sign extend
1290 int32_t sign = (int8_t)more >> 7;
1291 int32_t magic = denom->magic;
1292 int32_t q = libdivide__mullhi_s32(magic, numer);
1293 q += numer;
1294
1295 // If q is non-negative, we have nothing to do
1296 // If q is negative, we want to add either (2**shift)-1 if d is a power of
1297 // 2, or (2**shift) if it is not a power of 2
1298 uint32_t is_power_of_2 = !!(more & LIBDIVIDE_S32_SHIFT_PATH);
1299 uint32_t q_sign = (uint32_t)(q >> 31);
1300 q += q_sign & ((1 << shift) - is_power_of_2);
1301
1302 // Now arithmetic right shift
1303 q >>= shift;
1304
1305 // Negate if needed
1306 q = (q ^ sign) - sign;
1307
1308 return q;
1309 }
1310
libdivide_s32_recover(const struct libdivide_s32_t * denom)1311 int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom) {
1312 uint8_t more = denom->more;
1313 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
1314 if (more & LIBDIVIDE_S32_SHIFT_PATH) {
1315 uint32_t absD = 1U << shift;
1316 if (more & LIBDIVIDE_NEGATIVE_DIVISOR) {
1317 absD = -absD;
1318 }
1319 return (int32_t)absD;
1320 } else {
1321 // Unsigned math is much easier
1322 // We negate the magic number only in the branchfull case, and we don't
1323 // know which case we're in. However we have enough information to
1324 // determine the correct sign of the magic number. The divisor was
1325 // negative if LIBDIVIDE_NEGATIVE_DIVISOR is set. If ADD_MARKER is set,
1326 // the magic number's sign is opposite that of the divisor.
1327 // We want to compute the positive magic number.
1328 int negative_divisor = (more & LIBDIVIDE_NEGATIVE_DIVISOR);
1329 int magic_was_negated = (more & LIBDIVIDE_ADD_MARKER)
1330 ? denom->magic > 0 : denom->magic < 0;
1331
1332 // Handle the power of 2 case (including branchfree)
1333 if (denom->magic == 0) {
1334 int32_t result = 1 << shift;
1335 return negative_divisor ? -result : result;
1336 }
1337
1338 uint32_t d = (uint32_t)(magic_was_negated ? -denom->magic : denom->magic);
1339 uint64_t n = 1ULL << (32 + shift); // this shift cannot exceed 30
1340 uint32_t q = (uint32_t)(n / d);
1341 int32_t result = (int32_t)q;
1342 result += 1;
1343 return negative_divisor ? -result : result;
1344 }
1345 }
1346
libdivide_s32_branchfree_recover(const struct libdivide_s32_branchfree_t * denom)1347 int32_t libdivide_s32_branchfree_recover(const struct libdivide_s32_branchfree_t *denom) {
1348 return libdivide_s32_recover((const struct libdivide_s32_t *)denom);
1349 }
1350
libdivide_s32_get_algorithm(const struct libdivide_s32_t * denom)1351 int libdivide_s32_get_algorithm(const struct libdivide_s32_t *denom) {
1352 uint8_t more = denom->more;
1353 int positiveDivisor = !(more & LIBDIVIDE_NEGATIVE_DIVISOR);
1354 if (more & LIBDIVIDE_S32_SHIFT_PATH) return (positiveDivisor ? 0 : 1);
1355 else if (more & LIBDIVIDE_ADD_MARKER) return (positiveDivisor ? 2 : 3);
1356 else return 4;
1357 }
1358
libdivide_s32_do_alg0(int32_t numer,const struct libdivide_s32_t * denom)1359 int32_t libdivide_s32_do_alg0(int32_t numer, const struct libdivide_s32_t *denom) {
1360 uint8_t shifter = denom->more & LIBDIVIDE_32_SHIFT_MASK;
1361 int32_t q = numer + ((numer >> 31) & ((1U << shifter) - 1));
1362 return q >> shifter;
1363 }
1364
libdivide_s32_do_alg1(int32_t numer,const struct libdivide_s32_t * denom)1365 int32_t libdivide_s32_do_alg1(int32_t numer, const struct libdivide_s32_t *denom) {
1366 uint8_t shifter = denom->more & LIBDIVIDE_32_SHIFT_MASK;
1367 int32_t q = numer + ((numer >> 31) & ((1U << shifter) - 1));
1368 return - (q >> shifter);
1369 }
1370
libdivide_s32_do_alg2(int32_t numer,const struct libdivide_s32_t * denom)1371 int32_t libdivide_s32_do_alg2(int32_t numer, const struct libdivide_s32_t *denom) {
1372 int32_t q = libdivide__mullhi_s32(denom->magic, numer);
1373 q += numer;
1374 q >>= denom->more & LIBDIVIDE_32_SHIFT_MASK;
1375 q += (q < 0);
1376 return q;
1377 }
1378
libdivide_s32_do_alg3(int32_t numer,const struct libdivide_s32_t * denom)1379 int32_t libdivide_s32_do_alg3(int32_t numer, const struct libdivide_s32_t *denom) {
1380 int32_t q = libdivide__mullhi_s32(denom->magic, numer);
1381 q -= numer;
1382 q >>= denom->more & LIBDIVIDE_32_SHIFT_MASK;
1383 q += (q < 0);
1384 return q;
1385 }
1386
libdivide_s32_do_alg4(int32_t numer,const struct libdivide_s32_t * denom)1387 int32_t libdivide_s32_do_alg4(int32_t numer, const struct libdivide_s32_t *denom) {
1388 int32_t q = libdivide__mullhi_s32(denom->magic, numer);
1389 q >>= denom->more & LIBDIVIDE_32_SHIFT_MASK;
1390 q += (q < 0);
1391 return q;
1392 }
1393
1394 #if defined(LIBDIVIDE_USE_SSE2)
1395
libdivide_s32_do_vector(__m128i numers,const struct libdivide_s32_t * denom)1396 __m128i libdivide_s32_do_vector(__m128i numers, const struct libdivide_s32_t *denom) {
1397 uint8_t more = denom->more;
1398 if (more & LIBDIVIDE_S32_SHIFT_PATH) {
1399 uint32_t shifter = more & LIBDIVIDE_32_SHIFT_MASK;
1400 __m128i roundToZeroTweak = _mm_set1_epi32((1U << shifter) - 1); // could use _mm_srli_epi32 with an all -1 register
1401 __m128i q = _mm_add_epi32(numers, _mm_and_si128(_mm_srai_epi32(numers, 31), roundToZeroTweak)); //q = numer + ((numer >> 31) & roundToZeroTweak);
1402 q = _mm_sra_epi32(q, libdivide_u32_to_m128i(shifter)); // q = q >> shifter
1403 __m128i shiftMask = _mm_set1_epi32((int32_t)((int8_t)more >> 7)); // set all bits of shift mask = to the sign bit of more
1404 q = _mm_sub_epi32(_mm_xor_si128(q, shiftMask), shiftMask); // q = (q ^ shiftMask) - shiftMask;
1405 return q;
1406 }
1407 else {
1408 __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(denom->magic));
1409 if (more & LIBDIVIDE_ADD_MARKER) {
1410 __m128i sign = _mm_set1_epi32((int32_t)(int8_t)more >> 7); // must be arithmetic shift
1411 q = _mm_add_epi32(q, _mm_sub_epi32(_mm_xor_si128(numers, sign), sign)); // q += ((numer ^ sign) - sign);
1412 }
1413 q = _mm_sra_epi32(q, libdivide_u32_to_m128i(more & LIBDIVIDE_32_SHIFT_MASK)); // q >>= shift
1414 q = _mm_add_epi32(q, _mm_srli_epi32(q, 31)); // q += (q < 0)
1415 return q;
1416 }
1417 }
1418
libdivide_s32_do_vector_alg0(__m128i numers,const struct libdivide_s32_t * denom)1419 __m128i libdivide_s32_do_vector_alg0(__m128i numers, const struct libdivide_s32_t *denom) {
1420 uint8_t shifter = denom->more & LIBDIVIDE_32_SHIFT_MASK;
1421 __m128i roundToZeroTweak = _mm_set1_epi32((1U << shifter) - 1);
1422 __m128i q = _mm_add_epi32(numers, _mm_and_si128(_mm_srai_epi32(numers, 31), roundToZeroTweak));
1423 return _mm_sra_epi32(q, libdivide_u32_to_m128i(shifter));
1424 }
1425
libdivide_s32_do_vector_alg1(__m128i numers,const struct libdivide_s32_t * denom)1426 __m128i libdivide_s32_do_vector_alg1(__m128i numers, const struct libdivide_s32_t *denom) {
1427 uint8_t shifter = denom->more & LIBDIVIDE_32_SHIFT_MASK;
1428 __m128i roundToZeroTweak = _mm_set1_epi32((1U << shifter) - 1);
1429 __m128i q = _mm_add_epi32(numers, _mm_and_si128(_mm_srai_epi32(numers, 31), roundToZeroTweak));
1430 return _mm_sub_epi32(_mm_setzero_si128(), _mm_sra_epi32(q, libdivide_u32_to_m128i(shifter)));
1431 }
1432
libdivide_s32_do_vector_alg2(__m128i numers,const struct libdivide_s32_t * denom)1433 __m128i libdivide_s32_do_vector_alg2(__m128i numers, const struct libdivide_s32_t *denom) {
1434 __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(denom->magic));
1435 q = _mm_add_epi32(q, numers);
1436 q = _mm_sra_epi32(q, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_32_SHIFT_MASK));
1437 q = _mm_add_epi32(q, _mm_srli_epi32(q, 31));
1438 return q;
1439 }
1440
libdivide_s32_do_vector_alg3(__m128i numers,const struct libdivide_s32_t * denom)1441 __m128i libdivide_s32_do_vector_alg3(__m128i numers, const struct libdivide_s32_t *denom) {
1442 __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(denom->magic));
1443 q = _mm_sub_epi32(q, numers);
1444 q = _mm_sra_epi32(q, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_32_SHIFT_MASK));
1445 q = _mm_add_epi32(q, _mm_srli_epi32(q, 31));
1446 return q;
1447 }
1448
libdivide_s32_do_vector_alg4(__m128i numers,const struct libdivide_s32_t * denom)1449 __m128i libdivide_s32_do_vector_alg4(__m128i numers, const struct libdivide_s32_t *denom) {
1450 uint8_t more = denom->more;
1451 __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(denom->magic));
1452 q = _mm_sra_epi32(q, libdivide_u32_to_m128i(more & LIBDIVIDE_32_SHIFT_MASK)); //q >>= shift
1453 q = _mm_add_epi32(q, _mm_srli_epi32(q, 31)); // q += (q < 0)
1454 return q;
1455 }
1456
libdivide_s32_branchfree_do_vector(__m128i numers,const struct libdivide_s32_branchfree_t * denom)1457 __m128i libdivide_s32_branchfree_do_vector(__m128i numers, const struct libdivide_s32_branchfree_t *denom) {
1458 int32_t magic = denom->magic;
1459 uint8_t more = denom->more;
1460 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
1461 // must be arithmetic shift
1462 __m128i sign = _mm_set1_epi32((int32_t)(int8_t)more >> 7);
1463
1464 // libdivide__mullhi_s32(numers, magic);
1465 __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(magic));
1466 q = _mm_add_epi32(q, numers); // q += numers
1467
1468 // If q is non-negative, we have nothing to do
1469 // If q is negative, we want to add either (2**shift)-1 if d is a power of
1470 // 2, or (2**shift) if it is not a power of 2
1471 uint32_t is_power_of_2 = (magic == 0);
1472 __m128i q_sign = _mm_srai_epi32(q, 31); // q_sign = q >> 31
1473 __m128i mask = _mm_set1_epi32((1 << shift) - is_power_of_2);
1474 q = _mm_add_epi32(q, _mm_and_si128(q_sign, mask)); // q = q + (q_sign & mask)
1475 q = _mm_srai_epi32(q, shift); //q >>= shift
1476 q = _mm_sub_epi32(_mm_xor_si128(q, sign), sign); // q = (q ^ sign) - sign
1477 return q;
1478 }
1479
1480 #endif
1481
1482 ///////////// SINT64
1483
libdivide_internal_s64_gen(int64_t d,int branchfree)1484 static inline struct libdivide_s64_t libdivide_internal_s64_gen(int64_t d, int branchfree) {
1485 if (d == 0) {
1486 LIBDIVIDE_ERROR("divider must be != 0");
1487 }
1488
1489 struct libdivide_s64_t result;
1490
1491 // If d is a power of 2, or negative a power of 2, we have to use a shift.
1492 // This is especially important because the magic algorithm fails for -1.
1493 // To check if d is a power of 2 or its inverse, it suffices to check
1494 // whether its absolute value has exactly one bit set. This works even for
1495 // INT_MIN, because abs(INT_MIN) == INT_MIN, and INT_MIN has one bit set
1496 // and is a power of 2.
1497 uint64_t ud = (uint64_t)d;
1498 uint64_t absD = (d < 0) ? -ud : ud;
1499 uint32_t floor_log_2_d = 63 - libdivide__count_leading_zeros64(absD);
1500 // check if exactly one bit is set,
1501 // don't care if absD is 0 since that's divide by zero
1502 if ((absD & (absD - 1)) == 0) {
1503 // Branchfree and non-branchfree cases are the same
1504 result.magic = 0;
1505 result.more = floor_log_2_d | (d < 0 ? LIBDIVIDE_NEGATIVE_DIVISOR : 0);
1506 } else {
1507 // the dividend here is 2**(floor_log_2_d + 63), so the low 64 bit word
1508 // is 0 and the high word is floor_log_2_d - 1
1509 uint8_t more;
1510 uint64_t rem, proposed_m;
1511 proposed_m = libdivide_128_div_64_to_64(1ULL << (floor_log_2_d - 1), 0, absD, &rem);
1512 const uint64_t e = absD - rem;
1513
1514 // We are going to start with a power of floor_log_2_d - 1.
1515 // This works if works if e < 2**floor_log_2_d.
1516 if (!branchfree && e < (1ULL << floor_log_2_d)) {
1517 // This power works
1518 more = floor_log_2_d - 1;
1519 } else {
1520 // We need to go one higher. This should not make proposed_m
1521 // overflow, but it will make it negative when interpreted as an
1522 // int32_t.
1523 proposed_m += proposed_m;
1524 const uint64_t twice_rem = rem + rem;
1525 if (twice_rem >= absD || twice_rem < rem) proposed_m += 1;
1526 // note that we only set the LIBDIVIDE_NEGATIVE_DIVISOR bit if we
1527 // also set ADD_MARKER this is an annoying optimization that
1528 // enables algorithm #4 to avoid the mask. However we always set it
1529 // in the branchfree case
1530 more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;
1531 }
1532 proposed_m += 1;
1533 int64_t magic = (int64_t)proposed_m;
1534
1535 // Mark if we are negative
1536 if (d < 0) {
1537 more |= LIBDIVIDE_NEGATIVE_DIVISOR;
1538 if (!branchfree) {
1539 magic = -magic;
1540 }
1541 }
1542
1543 result.more = more;
1544 result.magic = magic;
1545 }
1546 return result;
1547 }
1548
libdivide_s64_gen(int64_t d)1549 struct libdivide_s64_t libdivide_s64_gen(int64_t d) {
1550 return libdivide_internal_s64_gen(d, 0);
1551 }
1552
libdivide_s64_branchfree_gen(int64_t d)1553 struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t d) {
1554 if (d == 1) {
1555 LIBDIVIDE_ERROR("branchfree divider must be != 1");
1556 }
1557 if (d == -1) {
1558 LIBDIVIDE_ERROR("branchfree divider must be != -1");
1559 }
1560 struct libdivide_s64_t tmp = libdivide_internal_s64_gen(d, 1);
1561 struct libdivide_s64_branchfree_t ret = {tmp.magic, tmp.more};
1562 return ret;
1563 }
1564
libdivide_s64_do(int64_t numer,const struct libdivide_s64_t * denom)1565 int64_t libdivide_s64_do(int64_t numer, const struct libdivide_s64_t *denom) {
1566 uint8_t more = denom->more;
1567 int64_t magic = denom->magic;
1568 if (magic == 0) { //shift path
1569 uint32_t shifter = more & LIBDIVIDE_64_SHIFT_MASK;
1570 uint64_t uq = (uint64_t)numer + ((numer >> 63) & ((1ULL << shifter) - 1));
1571 int64_t q = (int64_t)uq;
1572 q = q >> shifter;
1573 // must be arithmetic shift and then sign-extend
1574 int64_t shiftMask = (int8_t)more >> 7;
1575 q = (q ^ shiftMask) - shiftMask;
1576 return q;
1577 } else {
1578 uint64_t uq = (uint64_t)libdivide__mullhi_s64(magic, numer);
1579 if (more & LIBDIVIDE_ADD_MARKER) {
1580 // must be arithmetic shift and then sign extend
1581 int64_t sign = (int8_t)more >> 7;
1582 uq += ((uint64_t)numer ^ sign) - sign;
1583 }
1584 int64_t q = (int64_t)uq;
1585 q >>= more & LIBDIVIDE_64_SHIFT_MASK;
1586 q += (q < 0);
1587 return q;
1588 }
1589 }
1590
libdivide_s64_branchfree_do(int64_t numer,const struct libdivide_s64_branchfree_t * denom)1591 int64_t libdivide_s64_branchfree_do(int64_t numer, const struct libdivide_s64_branchfree_t *denom) {
1592 uint8_t more = denom->more;
1593 uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
1594 // must be arithmetic shift and then sign extend
1595 int64_t sign = (int8_t)more >> 7;
1596 int64_t magic = denom->magic;
1597 int64_t q = libdivide__mullhi_s64(magic, numer);
1598 q += numer;
1599
1600 // If q is non-negative, we have nothing to do.
1601 // If q is negative, we want to add either (2**shift)-1 if d is a power of
1602 // 2, or (2**shift) if it is not a power of 2.
1603 uint32_t is_power_of_2 = (magic == 0);
1604 uint64_t q_sign = (uint64_t)(q >> 63);
1605 q += q_sign & ((1ULL << shift) - is_power_of_2);
1606
1607 // Arithmetic right shift
1608 q >>= shift;
1609
1610 // Negate if needed
1611 q = (q ^ sign) - sign;
1612 return q;
1613 }
1614
libdivide_s64_recover(const struct libdivide_s64_t * denom)1615 int64_t libdivide_s64_recover(const struct libdivide_s64_t *denom) {
1616 uint8_t more = denom->more;
1617 uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
1618 if (denom->magic == 0) { // shift path
1619 uint64_t absD = 1ULL << shift;
1620 if (more & LIBDIVIDE_NEGATIVE_DIVISOR) {
1621 absD = -absD;
1622 }
1623 return (int64_t)absD;
1624 } else {
1625 // Unsigned math is much easier
1626 int negative_divisor = (more & LIBDIVIDE_NEGATIVE_DIVISOR);
1627 int magic_was_negated = (more & LIBDIVIDE_ADD_MARKER)
1628 ? denom->magic > 0 : denom->magic < 0;
1629
1630 uint64_t d = (uint64_t)(magic_was_negated ? -denom->magic : denom->magic);
1631 uint64_t n_hi = 1ULL << shift, n_lo = 0;
1632 uint64_t rem_ignored;
1633 uint64_t q = libdivide_128_div_64_to_64(n_hi, n_lo, d, &rem_ignored);
1634 int64_t result = (int64_t)(q + 1);
1635 if (negative_divisor) {
1636 result = -result;
1637 }
1638 return result;
1639 }
1640 }
1641
libdivide_s64_branchfree_recover(const struct libdivide_s64_branchfree_t * denom)1642 int64_t libdivide_s64_branchfree_recover(const struct libdivide_s64_branchfree_t *denom) {
1643 return libdivide_s64_recover((const struct libdivide_s64_t *)denom);
1644 }
1645
libdivide_s64_get_algorithm(const struct libdivide_s64_t * denom)1646 int libdivide_s64_get_algorithm(const struct libdivide_s64_t *denom) {
1647 uint8_t more = denom->more;
1648 int positiveDivisor = !(more & LIBDIVIDE_NEGATIVE_DIVISOR);
1649 if (denom->magic == 0) return (positiveDivisor ? 0 : 1); // shift path
1650 else if (more & LIBDIVIDE_ADD_MARKER) return (positiveDivisor ? 2 : 3);
1651 else return 4;
1652 }
1653
libdivide_s64_do_alg0(int64_t numer,const struct libdivide_s64_t * denom)1654 int64_t libdivide_s64_do_alg0(int64_t numer, const struct libdivide_s64_t *denom) {
1655 uint32_t shifter = denom->more & LIBDIVIDE_64_SHIFT_MASK;
1656 int64_t q = numer + ((numer >> 63) & ((1ULL << shifter) - 1));
1657 return q >> shifter;
1658 }
1659
libdivide_s64_do_alg1(int64_t numer,const struct libdivide_s64_t * denom)1660 int64_t libdivide_s64_do_alg1(int64_t numer, const struct libdivide_s64_t *denom) {
1661 // denom->shifter != -1 && demo->shiftMask != 0
1662 uint32_t shifter = denom->more & LIBDIVIDE_64_SHIFT_MASK;
1663 int64_t q = numer + ((numer >> 63) & ((1ULL << shifter) - 1));
1664 return - (q >> shifter);
1665 }
1666
libdivide_s64_do_alg2(int64_t numer,const struct libdivide_s64_t * denom)1667 int64_t libdivide_s64_do_alg2(int64_t numer, const struct libdivide_s64_t *denom) {
1668 int64_t q = libdivide__mullhi_s64(denom->magic, numer);
1669 q += numer;
1670 q >>= denom->more & LIBDIVIDE_64_SHIFT_MASK;
1671 q += (q < 0);
1672 return q;
1673 }
1674
libdivide_s64_do_alg3(int64_t numer,const struct libdivide_s64_t * denom)1675 int64_t libdivide_s64_do_alg3(int64_t numer, const struct libdivide_s64_t *denom) {
1676 int64_t q = libdivide__mullhi_s64(denom->magic, numer);
1677 q -= numer;
1678 q >>= denom->more & LIBDIVIDE_64_SHIFT_MASK;
1679 q += (q < 0);
1680 return q;
1681 }
1682
libdivide_s64_do_alg4(int64_t numer,const struct libdivide_s64_t * denom)1683 int64_t libdivide_s64_do_alg4(int64_t numer, const struct libdivide_s64_t *denom) {
1684 int64_t q = libdivide__mullhi_s64(denom->magic, numer);
1685 q >>= denom->more & LIBDIVIDE_64_SHIFT_MASK;
1686 q += (q < 0);
1687 return q;
1688 }
1689
1690 #if defined(LIBDIVIDE_USE_SSE2)
1691
libdivide_s64_do_vector(__m128i numers,const struct libdivide_s64_t * denom)1692 __m128i libdivide_s64_do_vector(__m128i numers, const struct libdivide_s64_t *denom) {
1693 uint8_t more = denom->more;
1694 int64_t magic = denom->magic;
1695 if (magic == 0) { // shift path
1696 uint32_t shifter = more & LIBDIVIDE_64_SHIFT_MASK;
1697 __m128i roundToZeroTweak = libdivide__u64_to_m128((1ULL << shifter) - 1);
1698 __m128i q = _mm_add_epi64(numers, _mm_and_si128(libdivide_s64_signbits(numers), roundToZeroTweak)); // q = numer + ((numer >> 63) & roundToZeroTweak);
1699 q = libdivide_s64_shift_right_vector(q, shifter); // q = q >> shifter
1700 __m128i shiftMask = _mm_set1_epi32((int32_t)((int8_t)more >> 7));
1701 q = _mm_sub_epi64(_mm_xor_si128(q, shiftMask), shiftMask); // q = (q ^ shiftMask) - shiftMask;
1702 return q;
1703 }
1704 else {
1705 __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide__u64_to_m128(magic));
1706 if (more & LIBDIVIDE_ADD_MARKER) {
1707 __m128i sign = _mm_set1_epi32((int32_t)((int8_t)more >> 7)); // must be arithmetic shift
1708 q = _mm_add_epi64(q, _mm_sub_epi64(_mm_xor_si128(numers, sign), sign)); // q += ((numer ^ sign) - sign);
1709 }
1710 // q >>= denom->mult_path.shift
1711 q = libdivide_s64_shift_right_vector(q, more & LIBDIVIDE_64_SHIFT_MASK);
1712 q = _mm_add_epi64(q, _mm_srli_epi64(q, 63)); // q += (q < 0)
1713 return q;
1714 }
1715 }
1716
libdivide_s64_do_vector_alg0(__m128i numers,const struct libdivide_s64_t * denom)1717 __m128i libdivide_s64_do_vector_alg0(__m128i numers, const struct libdivide_s64_t *denom) {
1718 uint32_t shifter = denom->more & LIBDIVIDE_64_SHIFT_MASK;
1719 __m128i roundToZeroTweak = libdivide__u64_to_m128((1ULL << shifter) - 1);
1720 __m128i q = _mm_add_epi64(numers, _mm_and_si128(libdivide_s64_signbits(numers), roundToZeroTweak));
1721 q = libdivide_s64_shift_right_vector(q, shifter);
1722 return q;
1723 }
1724
libdivide_s64_do_vector_alg1(__m128i numers,const struct libdivide_s64_t * denom)1725 __m128i libdivide_s64_do_vector_alg1(__m128i numers, const struct libdivide_s64_t *denom) {
1726 uint32_t shifter = denom->more & LIBDIVIDE_64_SHIFT_MASK;
1727 __m128i roundToZeroTweak = libdivide__u64_to_m128((1ULL << shifter) - 1);
1728 __m128i q = _mm_add_epi64(numers, _mm_and_si128(libdivide_s64_signbits(numers), roundToZeroTweak));
1729 q = libdivide_s64_shift_right_vector(q, shifter);
1730 return _mm_sub_epi64(_mm_setzero_si128(), q);
1731 }
1732
libdivide_s64_do_vector_alg2(__m128i numers,const struct libdivide_s64_t * denom)1733 __m128i libdivide_s64_do_vector_alg2(__m128i numers, const struct libdivide_s64_t *denom) {
1734 __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide__u64_to_m128(denom->magic));
1735 q = _mm_add_epi64(q, numers);
1736 q = libdivide_s64_shift_right_vector(q, denom->more & LIBDIVIDE_64_SHIFT_MASK);
1737 q = _mm_add_epi64(q, _mm_srli_epi64(q, 63)); // q += (q < 0)
1738 return q;
1739 }
1740
libdivide_s64_do_vector_alg3(__m128i numers,const struct libdivide_s64_t * denom)1741 __m128i libdivide_s64_do_vector_alg3(__m128i numers, const struct libdivide_s64_t *denom) {
1742 __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide__u64_to_m128(denom->magic));
1743 q = _mm_sub_epi64(q, numers);
1744 q = libdivide_s64_shift_right_vector(q, denom->more & LIBDIVIDE_64_SHIFT_MASK);
1745 q = _mm_add_epi64(q, _mm_srli_epi64(q, 63)); // q += (q < 0)
1746 return q;
1747 }
1748
libdivide_s64_do_vector_alg4(__m128i numers,const struct libdivide_s64_t * denom)1749 __m128i libdivide_s64_do_vector_alg4(__m128i numers, const struct libdivide_s64_t *denom) {
1750 __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide__u64_to_m128(denom->magic));
1751 q = libdivide_s64_shift_right_vector(q, denom->more & LIBDIVIDE_64_SHIFT_MASK);
1752 q = _mm_add_epi64(q, _mm_srli_epi64(q, 63));
1753 return q;
1754 }
1755
libdivide_s64_branchfree_do_vector(__m128i numers,const struct libdivide_s64_branchfree_t * denom)1756 __m128i libdivide_s64_branchfree_do_vector(__m128i numers, const struct libdivide_s64_branchfree_t *denom) {
1757 int64_t magic = denom->magic;
1758 uint8_t more = denom->more;
1759 uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
1760 // must be arithmetic shift
1761 __m128i sign = _mm_set1_epi32((int32_t)(int8_t)more >> 7);
1762
1763 // libdivide__mullhi_s64(numers, magic);
1764 __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide__u64_to_m128(magic));
1765 q = _mm_add_epi64(q, numers); // q += numers
1766
1767 // If q is non-negative, we have nothing to do.
1768 // If q is negative, we want to add either (2**shift)-1 if d is a power of
1769 // 2, or (2**shift) if it is not a power of 2.
1770 uint32_t is_power_of_2 = (magic == 0);
1771 __m128i q_sign = libdivide_s64_signbits(q); // q_sign = q >> 63
1772 __m128i mask = libdivide__u64_to_m128((1ULL << shift) - is_power_of_2);
1773 q = _mm_add_epi64(q, _mm_and_si128(q_sign, mask)); // q = q + (q_sign & mask)
1774 q = libdivide_s64_shift_right_vector(q, shift); // q >>= shift
1775 q = _mm_sub_epi64(_mm_xor_si128(q, sign), sign); // q = (q ^ sign) - sign
1776 return q;
1777 }
1778
1779 #endif
1780
1781 /////////// C++ stuff
1782
1783 #ifdef __cplusplus
1784
1785 // Our divider struct is templated on both a type (like uint64_t) and an
1786 // algorithm index. BRANCHFULL is the default algorithm, BRANCHFREE is the
1787 // branchfree variant, and the indexed variants are for unswitching.
1788 enum {
1789 BRANCHFULL = -1,
1790 BRANCHFREE = -2,
1791 ALGORITHM0 = 0,
1792 ALGORITHM1 = 1,
1793 ALGORITHM2 = 2,
1794 ALGORITHM3 = 3,
1795 ALGORITHM4 = 4
1796 };
1797
1798 namespace libdivide_internal {
1799
1800 #if defined(LIBDIVIDE_USE_SSE2)
1801 #define MAYBE_VECTOR(X) X
1802 #define MAYBE_VECTOR_PARAM(X) __m128i vector_func(__m128i, const X *)
1803 #else
1804 #define MAYBE_VECTOR(X) 0
1805 #define MAYBE_VECTOR_PARAM(X) int unused
1806 #endif
1807
1808 // The following convenience macros are used to build a type of the base
1809 // divider class and give it as template arguments the C functions
1810 // related to the macro name and the macro type paramaters.
1811
1812 #define BRANCHFULL_DIVIDER(INT, TYPE) \
1813 typedef base<INT, \
1814 libdivide_##TYPE##_t, \
1815 libdivide_##TYPE##_gen, \
1816 libdivide_##TYPE##_do, \
1817 MAYBE_VECTOR(libdivide_##TYPE##_do_vector)>
1818
1819 #define BRANCHFREE_DIVIDER(INT, TYPE) \
1820 typedef base<INT, \
1821 libdivide_##TYPE##_branchfree_t, \
1822 libdivide_##TYPE##_branchfree_gen, \
1823 libdivide_##TYPE##_branchfree_do, \
1824 MAYBE_VECTOR(libdivide_##TYPE##_branchfree_do_vector)>
1825
1826 #define ALGORITHM_DIVIDER(INT, TYPE, ALGO) \
1827 typedef base<INT, \
1828 libdivide_##TYPE##_t, \
1829 libdivide_##TYPE##_gen, \
1830 libdivide_##TYPE##_do_##ALGO, \
1831 MAYBE_VECTOR(libdivide_##TYPE##_do_vector_##ALGO)>
1832
1833 #define CRASH_DIVIDER(INT, TYPE) \
1834 typedef base<INT, \
1835 libdivide_##TYPE##_t, \
1836 libdivide_##TYPE##_gen, \
1837 libdivide_##TYPE##_crash, \
1838 MAYBE_VECTOR(libdivide_##TYPE##_crash_vector)>
1839
1840 // Base divider, provides storage for the actual divider.
1841 // @IntType: e.g. uint32_t
1842 // @DenomType: e.g. libdivide_u32_t
1843 // @gen_func(): e.g. libdivide_u32_gen
1844 // @do_func(): e.g. libdivide_u32_do
1845 // @MAYBE_VECTOR_PARAM: e.g. libdivide_u32_do_vector
1846 template<typename IntType,
1847 typename DenomType,
1848 DenomType gen_func(IntType),
1849 IntType do_func(IntType, const DenomType *),
1850 MAYBE_VECTOR_PARAM(DenomType)>
1851 struct base {
1852 // Storage for the actual divider
1853 DenomType denom;
1854
1855 // Constructor that takes a divisor value, and applies the gen function
basebase1856 base(IntType d) : denom(gen_func(d)) { }
1857
1858 // Default constructor to allow uninitialized uses in e.g. arrays
basebase1859 base() {}
1860
1861 // Needed for unswitch
basebase1862 base(const DenomType& d) : denom(d) { }
1863
perform_dividebase1864 IntType perform_divide(IntType val) const {
1865 return do_func(val, &denom);
1866 }
1867
1868 #if defined(LIBDIVIDE_USE_SSE2)
perform_divide_vectorbase1869 __m128i perform_divide_vector(__m128i val) const {
1870 return vector_func(val, &denom);
1871 }
1872 #endif
1873 };
1874
1875 // Functions that will never be called but are required to be able
1876 // to use unswitch in C++ template code. Unsigned has fewer algorithms
1877 // than signed i.e. alg3 and alg4 are not defined for unsigned. In
1878 // order to make templates compile we need to define unsigned alg3 and
1879 // alg4 as crash functions.
libdivide_u32_crash(uint32_t,const libdivide_u32_t *)1880 uint32_t libdivide_u32_crash(uint32_t, const libdivide_u32_t *) { exit(-1); }
libdivide_u64_crash(uint64_t,const libdivide_u64_t *)1881 uint64_t libdivide_u64_crash(uint64_t, const libdivide_u64_t *) { exit(-1); }
1882
1883 #if defined(LIBDIVIDE_USE_SSE2)
libdivide_u32_crash_vector(__m128i,const libdivide_u32_t *)1884 __m128i libdivide_u32_crash_vector(__m128i, const libdivide_u32_t *) { exit(-1); }
libdivide_u64_crash_vector(__m128i,const libdivide_u64_t *)1885 __m128i libdivide_u64_crash_vector(__m128i, const libdivide_u64_t *) { exit(-1); }
1886 #endif
1887
1888 template<typename T, int ALGO> struct dispatcher { };
1889
1890 // Templated dispatch using partial specialization
1891 template<> struct dispatcher<int32_t, BRANCHFULL> { BRANCHFULL_DIVIDER(int32_t, s32) divider; };
1892 template<> struct dispatcher<int32_t, BRANCHFREE> { BRANCHFREE_DIVIDER(int32_t, s32) divider; };
1893 template<> struct dispatcher<int32_t, ALGORITHM0> { ALGORITHM_DIVIDER(int32_t, s32, alg0) divider; };
1894 template<> struct dispatcher<int32_t, ALGORITHM1> { ALGORITHM_DIVIDER(int32_t, s32, alg1) divider; };
1895 template<> struct dispatcher<int32_t, ALGORITHM2> { ALGORITHM_DIVIDER(int32_t, s32, alg2) divider; };
1896 template<> struct dispatcher<int32_t, ALGORITHM3> { ALGORITHM_DIVIDER(int32_t, s32, alg3) divider; };
1897 template<> struct dispatcher<int32_t, ALGORITHM4> { ALGORITHM_DIVIDER(int32_t, s32, alg4) divider; };
1898
1899 template<> struct dispatcher<uint32_t, BRANCHFULL> { BRANCHFULL_DIVIDER(uint32_t, u32) divider; };
1900 template<> struct dispatcher<uint32_t, BRANCHFREE> { BRANCHFREE_DIVIDER(uint32_t, u32) divider; };
1901 template<> struct dispatcher<uint32_t, ALGORITHM0> { ALGORITHM_DIVIDER(uint32_t, u32, alg0) divider; };
1902 template<> struct dispatcher<uint32_t, ALGORITHM1> { ALGORITHM_DIVIDER(uint32_t, u32, alg1) divider; };
1903 template<> struct dispatcher<uint32_t, ALGORITHM2> { ALGORITHM_DIVIDER(uint32_t, u32, alg2) divider; };
1904 template<> struct dispatcher<uint32_t, ALGORITHM3> { CRASH_DIVIDER(uint32_t, u32) divider; };
1905 template<> struct dispatcher<uint32_t, ALGORITHM4> { CRASH_DIVIDER(uint32_t, u32) divider; };
1906
1907 template<> struct dispatcher<int64_t, BRANCHFULL> { BRANCHFULL_DIVIDER(int64_t, s64) divider; };
1908 template<> struct dispatcher<int64_t, BRANCHFREE> { BRANCHFREE_DIVIDER(int64_t, s64) divider; };
1909 template<> struct dispatcher<int64_t, ALGORITHM0> { ALGORITHM_DIVIDER (int64_t, s64, alg0) divider; };
1910 template<> struct dispatcher<int64_t, ALGORITHM1> { ALGORITHM_DIVIDER (int64_t, s64, alg1) divider; };
1911 template<> struct dispatcher<int64_t, ALGORITHM2> { ALGORITHM_DIVIDER (int64_t, s64, alg2) divider; };
1912 template<> struct dispatcher<int64_t, ALGORITHM3> { ALGORITHM_DIVIDER (int64_t, s64, alg3) divider; };
1913 template<> struct dispatcher<int64_t, ALGORITHM4> { ALGORITHM_DIVIDER (int64_t, s64, alg4) divider; };
1914
1915 template<> struct dispatcher<uint64_t, BRANCHFULL> { BRANCHFULL_DIVIDER(uint64_t, u64) divider; };
1916 template<> struct dispatcher<uint64_t, BRANCHFREE> { BRANCHFREE_DIVIDER(uint64_t, u64) divider; };
1917 template<> struct dispatcher<uint64_t, ALGORITHM0> { ALGORITHM_DIVIDER(uint64_t, u64, alg0) divider; };
1918 template<> struct dispatcher<uint64_t, ALGORITHM1> { ALGORITHM_DIVIDER(uint64_t, u64, alg1) divider; };
1919 template<> struct dispatcher<uint64_t, ALGORITHM2> { ALGORITHM_DIVIDER(uint64_t, u64, alg2) divider; };
1920 template<> struct dispatcher<uint64_t, ALGORITHM3> { CRASH_DIVIDER(uint64_t, u64) divider; };
1921 template<> struct dispatcher<uint64_t, ALGORITHM4> { CRASH_DIVIDER(uint64_t, u64) divider; };
1922
1923 // Overloads that don't depend on the algorithm
1924 inline int32_t recover(const libdivide_s32_t *s) { return libdivide_s32_recover(s); }
1925 inline uint32_t recover(const libdivide_u32_t *s) { return libdivide_u32_recover(s); }
1926 inline int64_t recover(const libdivide_s64_t *s) { return libdivide_s64_recover(s); }
1927 inline uint64_t recover(const libdivide_u64_t *s) { return libdivide_u64_recover(s); }
1928
1929 inline int32_t recover(const libdivide_s32_branchfree_t *s) { return libdivide_s32_branchfree_recover(s); }
1930 inline uint32_t recover(const libdivide_u32_branchfree_t *s) { return libdivide_u32_branchfree_recover(s); }
1931 inline int64_t recover(const libdivide_s64_branchfree_t *s) { return libdivide_s64_branchfree_recover(s); }
1932 inline uint64_t recover(const libdivide_u64_branchfree_t *s) { return libdivide_u64_branchfree_recover(s); }
1933
1934 inline int get_algorithm(const libdivide_s32_t *s) { return libdivide_s32_get_algorithm(s); }
1935 inline int get_algorithm(const libdivide_u32_t *s) { return libdivide_u32_get_algorithm(s); }
1936 inline int get_algorithm(const libdivide_s64_t *s) { return libdivide_s64_get_algorithm(s); }
1937 inline int get_algorithm(const libdivide_u64_t *s) { return libdivide_u64_get_algorithm(s); }
1938
1939 // Fallback for branchfree variants, which do not support unswitching
1940 template<typename T> int get_algorithm(const T *) { return -1; }
1941 }
1942
1943 // This is the main divider class for use by the user (C++ API).
1944 // The divider itself is stored in the div variable who's
1945 // type is chosen by the dispatcher based on the template paramaters.
1946 template<typename T, int ALGO = BRANCHFULL>
1947 class divider
1948 {
1949 private:
1950 // Here's the actual divider
1951 typedef typename libdivide_internal::dispatcher<T, ALGO>::divider div_t;
1952 div_t div;
1953
1954 // unswitch() friend declaration
1955 template<int NEW_ALGO, typename S>
1956 friend divider<S, NEW_ALGO> unswitch(const divider<S, BRANCHFULL> & d);
1957
1958 // Constructor used by the unswitch friend
1959 divider(const div_t& denom) : div(denom) { }
1960
1961 public:
1962 // Ordinary constructor that takes the divisor as a parameter
1963 divider(T n) : div(n) { }
1964
1965 // Default constructor. We leave this deliberately undefined so that
1966 // creating an array of divider and then initializing them
1967 // doesn't slow us down.
1968 divider() { }
1969
1970 // Divides the parameter by the divisor, returning the quotient
1971 T perform_divide(T val) const {
1972 return div.perform_divide(val);
1973 }
1974
1975 // Recovers the divisor that was used to initialize the divider
1976 T recover_divisor() const {
1977 return libdivide_internal::recover(&div.denom);
1978 }
1979
1980 #if defined(LIBDIVIDE_USE_SSE2)
1981 // Treats the vector as either two or four packed values (depending on the
1982 // size), and divides each of them by the divisor,
1983 // returning the packed quotients.
1984 __m128i perform_divide_vector(__m128i val) const {
1985 return div.perform_divide_vector(val);
1986 }
1987 #endif
1988
1989 // Returns the index of algorithm, for use in the unswitch function. Does
1990 // not apply to branchfree variant.
1991 // Returns the algorithm for unswitching.
1992 int get_algorithm() const {
1993 return libdivide_internal::get_algorithm(&div.denom);
1994 }
1995
1996 bool operator==(const divider<T, ALGO>& him) const {
1997 return div.denom.magic == him.div.denom.magic &&
1998 div.denom.more == him.div.denom.more;
1999 }
2000
2001 bool operator!=(const divider<T, ALGO>& him) const {
2002 return !(*this == him);
2003 }
2004 };
2005
2006 #if __cplusplus >= 201103L || \
2007 (defined(_MSC_VER) && _MSC_VER >= 1800)
2008
2009 // libdivdie::branchfree_divider<T>
2010 template <typename T>
2011 using branchfree_divider = divider<T, BRANCHFREE>;
2012
2013 #endif
2014
2015 // Returns a divider specialized for the given algorithm
2016 template<int NEW_ALGO, typename T>
2017 divider<T, NEW_ALGO> unswitch(const divider<T, BRANCHFULL>& d) {
2018 return divider<T, NEW_ALGO>(d.div.denom);
2019 }
2020
2021 // Overload of the / operator for scalar division
2022 template<typename int_type, int ALGO>
2023 int_type operator/(int_type numer, const divider<int_type, ALGO>& denom) {
2024 return denom.perform_divide(numer);
2025 }
2026
2027 // Overload of the /= operator for scalar division
2028 template<typename int_type, int ALGO>
2029 int_type operator/=(int_type& numer, const divider<int_type, ALGO>& denom) {
2030 numer = denom.perform_divide(numer);
2031 return numer;
2032 }
2033
2034 #if defined(LIBDIVIDE_USE_SSE2)
2035
2036 // Overload of the / operator for vector division
2037 template<typename int_type, int ALGO>
2038 __m128i operator/(__m128i numer, const divider<int_type, ALGO>& denom) {
2039 return denom.perform_divide_vector(numer);
2040 }
2041
2042 // Overload of the /= operator for vector division
2043 template<typename int_type, int ALGO>
2044 __m128i operator/=(__m128i& numer, const divider<int_type, ALGO>& denom) {
2045 numer = denom.perform_divide_vector(numer);
2046 return numer;
2047 }
2048
2049 #endif
2050
2051 } // namespace libdivide
2052 } // anonymous namespace
2053
2054 #endif // __cplusplus
2055
2056 #endif // LIBDIVIDE_H
2057