codekingpro/portable-devtools
114k
1// libdivide.h - Optimized integer division
2// https://libdivide.com
3//
4// Copyright (C) 2010 - 2019 ridiculous_fish, <libdivide@ridiculousfish.com>
5// Copyright (C) 2016 - 2019 Kim Walisch, <kim.walisch@gmail.com>
6//
7// libdivide is dual-licensed under the Boost or zlib licenses.
8// You may use libdivide under the terms of either of these.
9// See LICENSE.txt for more details.
10
11#ifndef NUMPY_CORE_INCLUDE_NUMPY_LIBDIVIDE_LIBDIVIDE_H_
12#define NUMPY_CORE_INCLUDE_NUMPY_LIBDIVIDE_LIBDIVIDE_H_
13
14#define LIBDIVIDE_VERSION "3.0"
15#define LIBDIVIDE_VERSION_MAJOR 3
16#define LIBDIVIDE_VERSION_MINOR 0
17
18#include <stdint.h>
19
20#if defined(__cplusplus)
21 #include <cstdlib>
22 #include <cstdio>
23 #include <type_traits>
24#else
25 #include <stdlib.h>
26 #include <stdio.h>
27#endif
28
29#if defined(LIBDIVIDE_AVX512)
30 #include <immintrin.h>
31#elif defined(LIBDIVIDE_AVX2)
32 #include <immintrin.h>
33#elif defined(LIBDIVIDE_SSE2)
34 #include <emmintrin.h>
35#endif
36
37#if defined(_MSC_VER)
38 #include <intrin.h>
39 // disable warning C4146: unary minus operator applied
40 // to unsigned type, result still unsigned
41 #pragma warning(disable: 4146)
42 #define LIBDIVIDE_VC
43#endif
44
45#if !defined(__has_builtin)
46 #define __has_builtin(x) 0
47#endif
48
49#if defined(__SIZEOF_INT128__)
50 #define HAS_INT128_T
51 // clang-cl on Windows does not yet support 128-bit division
52 #if !(defined(__clang__) && defined(LIBDIVIDE_VC))
53 #define HAS_INT128_DIV
54 #endif
55#endif
56
57#if defined(__x86_64__) || defined(_M_X64)
58 #define LIBDIVIDE_X86_64
59#endif
60
61#if defined(__i386__)
62 #define LIBDIVIDE_i386
63#endif
64
65#if defined(__GNUC__) || defined(__clang__)
66 #define LIBDIVIDE_GCC_STYLE_ASM
67#endif
68
69#if defined(__cplusplus) || defined(LIBDIVIDE_VC)
70 #define LIBDIVIDE_FUNCTION __FUNCTION__
71#else
72 #define LIBDIVIDE_FUNCTION __func__
73#endif
74
75#define LIBDIVIDE_ERROR(msg) \
76 do { \
77 fprintf(stderr, "libdivide.h:%d: %s(): Error: %s\n", \
78 __LINE__, LIBDIVIDE_FUNCTION, msg); \
79 abort(); \
80 } while (0)
81
82#if defined(LIBDIVIDE_ASSERTIONS_ON)
83 #define LIBDIVIDE_ASSERT(x) \
84 do { \
85 if (!(x)) { \
86 fprintf(stderr, "libdivide.h:%d: %s(): Assertion failed: %s\n", \
87 __LINE__, LIBDIVIDE_FUNCTION, #x); \
88 abort(); \
89 } \
90 } while (0)
91#else
92 #define LIBDIVIDE_ASSERT(x)
93#endif
94
95#ifdef __cplusplus
96namespace libdivide {
97#endif
98
99// pack divider structs to prevent compilers from padding.
100// This reduces memory usage by up to 43% when using a large
101// array of libdivide dividers and improves performance
102// by up to 10% because of reduced memory bandwidth.
103#pragma pack(push, 1)
104
105struct libdivide_u32_t {
106 uint32_t magic;
107 uint8_t more;
108};
109
110struct libdivide_s32_t {
111 int32_t magic;
112 uint8_t more;
113};
114
115struct libdivide_u64_t {
116 uint64_t magic;
117 uint8_t more;
118};
119
120struct libdivide_s64_t {
121 int64_t magic;
122 uint8_t more;
123};
124
125struct libdivide_u32_branchfree_t {
126 uint32_t magic;
127 uint8_t more;
128};
129
130struct libdivide_s32_branchfree_t {
131 int32_t magic;
132 uint8_t more;
133};
134
135struct libdivide_u64_branchfree_t {
136 uint64_t magic;
137 uint8_t more;
138};
139
140struct libdivide_s64_branchfree_t {
141 int64_t magic;
142 uint8_t more;
143};
144
145#pragma pack(pop)
146
147// Explanation of the "more" field:
148//
149// * Bits 0-5 is the shift value (for shift path or mult path).
150// * Bit 6 is the add indicator for mult path.
151// * Bit 7 is set if the divisor is negative. We use bit 7 as the negative
152// divisor indicator so that we can efficiently use sign extension to
153// create a bitmask with all bits set to 1 (if the divisor is negative)
154// or 0 (if the divisor is positive).
155//
156// u32: [0-4] shift value
157// [5] ignored
158// [6] add indicator
159// magic number of 0 indicates shift path
160//
161// s32: [0-4] shift value
162// [5] ignored
163// [6] add indicator
164// [7] indicates negative divisor
165// magic number of 0 indicates shift path
166//
167// u64: [0-5] shift value
168// [6] add indicator
169// magic number of 0 indicates shift path
170//
171// s64: [0-5] shift value
172// [6] add indicator
173// [7] indicates negative divisor
174// magic number of 0 indicates shift path
175//
176// In s32 and s64 branchfree modes, the magic number is negated according to
177// whether the divisor is negated. In branchfree strategy, it is not negated.
178
179enum {
180 LIBDIVIDE_32_SHIFT_MASK = 0x1F,
181 LIBDIVIDE_64_SHIFT_MASK = 0x3F,
182 LIBDIVIDE_ADD_MARKER = 0x40,
183 LIBDIVIDE_NEGATIVE_DIVISOR = 0x80
184};
185
186static inline struct libdivide_s32_t libdivide_s32_gen(int32_t d);
187static inline struct libdivide_u32_t libdivide_u32_gen(uint32_t d);
188static inline struct libdivide_s64_t libdivide_s64_gen(int64_t d);
189static inline struct libdivide_u64_t libdivide_u64_gen(uint64_t d);
190
191static inline struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t d);
192static inline struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t d);
193static inline struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t d);
194static inline struct libdivide_u64_branchfree_t libdivide_u64_branchfree_gen(uint64_t d);
195
196static inline int32_t libdivide_s32_do(int32_t numer, const struct libdivide_s32_t *denom);
197static inline uint32_t libdivide_u32_do(uint32_t numer, const struct libdivide_u32_t *denom);
198static inline int64_t libdivide_s64_do(int64_t numer, const struct libdivide_s64_t *denom);
199static inline uint64_t libdivide_u64_do(uint64_t numer, const struct libdivide_u64_t *denom);
200
201static inline int32_t libdivide_s32_branchfree_do(int32_t numer, const struct libdivide_s32_branchfree_t *denom);
202static inline uint32_t libdivide_u32_branchfree_do(uint32_t numer, const struct libdivide_u32_branchfree_t *denom);
203static inline int64_t libdivide_s64_branchfree_do(int64_t numer, const struct libdivide_s64_branchfree_t *denom);
204static inline uint64_t libdivide_u64_branchfree_do(uint64_t numer, const struct libdivide_u64_branchfree_t *denom);
205
206static inline int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom);
207static inline uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom);
208static inline int64_t libdivide_s64_recover(const struct libdivide_s64_t *denom);
209static inline uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom);
210
211static inline int32_t libdivide_s32_branchfree_recover(const struct libdivide_s32_branchfree_t *denom);
212static inline uint32_t libdivide_u32_branchfree_recover(const struct libdivide_u32_branchfree_t *denom);
213static inline int64_t libdivide_s64_branchfree_recover(const struct libdivide_s64_branchfree_t *denom);
214static inline uint64_t libdivide_u64_branchfree_recover(const struct libdivide_u64_branchfree_t *denom);
215
216//////// Internal Utility Functions
217
218static inline uint32_t libdivide_mullhi_u32(uint32_t x, uint32_t y) {
219 uint64_t xl = x, yl = y;
220 uint64_t rl = xl * yl;
221 return (uint32_t)(rl >> 32);
222}
223
224static inline int32_t libdivide_mullhi_s32(int32_t x, int32_t y) {
225 int64_t xl = x, yl = y;
226 int64_t rl = xl * yl;
227 // needs to be arithmetic shift
228 return (int32_t)(rl >> 32);
229}
230
231static inline uint64_t libdivide_mullhi_u64(uint64_t x, uint64_t y) {
232#if defined(LIBDIVIDE_VC) && \
233 defined(LIBDIVIDE_X86_64)
234 return __umulh(x, y);
235#elif defined(HAS_INT128_T)
236 __uint128_t xl = x, yl = y;
237 __uint128_t rl = xl * yl;
238 return (uint64_t)(rl >> 64);
239#else
240 // full 128 bits are x0 * y0 + (x0 * y1 << 32) + (x1 * y0 << 32) + (x1 * y1 << 64)
241 uint32_t mask = 0xFFFFFFFF;
242 uint32_t x0 = (uint32_t)(x & mask);
243 uint32_t x1 = (uint32_t)(x >> 32);
244 uint32_t y0 = (uint32_t)(y & mask);
245 uint32_t y1 = (uint32_t)(y >> 32);
246 uint32_t x0y0_hi = libdivide_mullhi_u32(x0, y0);
247 uint64_t x0y1 = x0 * (uint64_t)y1;
248 uint64_t x1y0 = x1 * (uint64_t)y0;
249 uint64_t x1y1 = x1 * (uint64_t)y1;
250 uint64_t temp = x1y0 + x0y0_hi;
251 uint64_t temp_lo = temp & mask;
252 uint64_t temp_hi = temp >> 32;
253
254 return x1y1 + temp_hi + ((temp_lo + x0y1) >> 32);
255#endif
256}
257
258static inline int64_t libdivide_mullhi_s64(int64_t x, int64_t y) {
259#if defined(LIBDIVIDE_VC) && \
260 defined(LIBDIVIDE_X86_64)
261 return __mulh(x, y);
262#elif defined(HAS_INT128_T)
263 __int128_t xl = x, yl = y;
264 __int128_t rl = xl * yl;
265 return (int64_t)(rl >> 64);
266#else
267 // full 128 bits are x0 * y0 + (x0 * y1 << 32) + (x1 * y0 << 32) + (x1 * y1 << 64)
268 uint32_t mask = 0xFFFFFFFF;
269 uint32_t x0 = (uint32_t)(x & mask);
270 uint32_t y0 = (uint32_t)(y & mask);
271 int32_t x1 = (int32_t)(x >> 32);
272 int32_t y1 = (int32_t)(y >> 32);
273 uint32_t x0y0_hi = libdivide_mullhi_u32(x0, y0);
274 int64_t t = x1 * (int64_t)y0 + x0y0_hi;
275 int64_t w1 = x0 * (int64_t)y1 + (t & mask);
276
277 return x1 * (int64_t)y1 + (t >> 32) + (w1 >> 32);
278#endif
279}
280
281static inline int32_t libdivide_count_leading_zeros32(uint32_t val) {
282#if defined(__GNUC__) || \
283 __has_builtin(__builtin_clz)
284 // Fast way to count leading zeros
285 return __builtin_clz(val);
286#elif defined(LIBDIVIDE_VC)
287 unsigned long result;
288 if (_BitScanReverse(&result, val)) {
289 return 31 - result;
290 }
291 return 0;
292#else
293 if (val == 0)
294 return 32;
295 int32_t result = 8;
296 uint32_t hi = 0xFFU << 24;
297 while ((val & hi) == 0) {
298 hi >>= 8;
299 result += 8;
300 }
301 while (val & hi) {
302 result -= 1;
303 hi <<= 1;
304 }
305 return result;
306#endif
307}
308
309static inline int32_t libdivide_count_leading_zeros64(uint64_t val) {
310#if defined(__GNUC__) || \
311 __has_builtin(__builtin_clzll)
312 // Fast way to count leading zeros
313 return __builtin_clzll(val);
314#elif defined(LIBDIVIDE_VC) && defined(_WIN64)
315 unsigned long result;
316 if (_BitScanReverse64(&result, val)) {
317 return 63 - result;
318 }
319 return 0;
320#else
321 uint32_t hi = val >> 32;
322 uint32_t lo = val & 0xFFFFFFFF;
323 if (hi != 0) return libdivide_count_leading_zeros32(hi);
324 return 32 + libdivide_count_leading_zeros32(lo);
325#endif
326}
327
328// libdivide_64_div_32_to_32: divides a 64-bit uint {u1, u0} by a 32-bit
329// uint {v}. The result must fit in 32 bits.
330// Returns the quotient directly and the remainder in *r
331static inline uint32_t libdivide_64_div_32_to_32(uint32_t u1, uint32_t u0, uint32_t v, uint32_t *r) {
332#if (defined(LIBDIVIDE_i386) || defined(LIBDIVIDE_X86_64)) && \
333 defined(LIBDIVIDE_GCC_STYLE_ASM)
334 uint32_t result;
335 __asm__("divl %[v]"
336 : "=a"(result), "=d"(*r)
337 : [v] "r"(v), "a"(u0), "d"(u1)
338 );
339 return result;
340#else
341 uint64_t n = ((uint64_t)u1 << 32) | u0;
342 uint32_t result = (uint32_t)(n / v);
343 *r = (uint32_t)(n - result * (uint64_t)v);
344 return result;
345#endif
346}
347
348// libdivide_128_div_64_to_64: divides a 128-bit uint {u1, u0} by a 64-bit
349// uint {v}. The result must fit in 64 bits.
350// Returns the quotient directly and the remainder in *r
351static uint64_t libdivide_128_div_64_to_64(uint64_t u1, uint64_t u0, uint64_t v, uint64_t *r) {
352#if defined(LIBDIVIDE_X86_64) && \
353 defined(LIBDIVIDE_GCC_STYLE_ASM)
354 uint64_t result;
355 __asm__("divq %[v]"
356 : "=a"(result), "=d"(*r)
357 : [v] "r"(v), "a"(u0), "d"(u1)
358 );
359 return result;
360#elif defined(HAS_INT128_T) && \
361 defined(HAS_INT128_DIV)
362 __uint128_t n = ((__uint128_t)u1 << 64) | u0;
363 uint64_t result = (uint64_t)(n / v);
364 *r = (uint64_t)(n - result * (__uint128_t)v);
365 return result;
366#else
367 // Code taken from Hacker's Delight:
368 // http://www.hackersdelight.org/HDcode/divlu.c.
369 // License permits inclusion here per:
370 // http://www.hackersdelight.org/permissions.htm
371
372 const uint64_t b = (1ULL << 32); // Number base (32 bits)
373 uint64_t un1, un0; // Norm. dividend LSD's
374 uint64_t vn1, vn0; // Norm. divisor digits
375 uint64_t q1, q0; // Quotient digits
376 uint64_t un64, un21, un10; // Dividend digit pairs
377 uint64_t rhat; // A remainder
378 int32_t s; // Shift amount for norm
379
380 // If overflow, set rem. to an impossible value,
381 // and return the largest possible quotient
382 if (u1 >= v) {
383 *r = (uint64_t) -1;
384 return (uint64_t) -1;
385 }
386
387 // count leading zeros
388 s = libdivide_count_leading_zeros64(v);
389 if (s > 0) {
390 // Normalize divisor
391 v = v << s;
392 un64 = (u1 << s) | (u0 >> (64 - s));
393 un10 = u0 << s; // Shift dividend left
394 } else {
395 // Avoid undefined behavior of (u0 >> 64).
396 // The behavior is undefined if the right operand is
397 // negative, or greater than or equal to the length
398 // in bits of the promoted left operand.
399 un64 = u1;
400 un10 = u0;
401 }
402
403 // Break divisor up into two 32-bit digits
404 vn1 = v >> 32;
405 vn0 = v & 0xFFFFFFFF;
406
407 // Break right half of dividend into two digits
408 un1 = un10 >> 32;
409 un0 = un10 & 0xFFFFFFFF;
410
411 // Compute the first quotient digit, q1
412 q1 = un64 / vn1;
413 rhat = un64 - q1 * vn1;
414
415 while (q1 >= b || q1 * vn0 > b * rhat + un1) {
416 q1 = q1 - 1;
417 rhat = rhat + vn1;
418 if (rhat >= b)
419 break;
420 }
421
422 // Multiply and subtract
423 un21 = un64 * b + un1 - q1 * v;
424
425 // Compute the second quotient digit
426 q0 = un21 / vn1;
427 rhat = un21 - q0 * vn1;
428
429 while (q0 >= b || q0 * vn0 > b * rhat + un0) {
430 q0 = q0 - 1;
431 rhat = rhat + vn1;
432 if (rhat >= b)
433 break;
434 }
435
436 *r = (un21 * b + un0 - q0 * v) >> s;
437 return q1 * b + q0;
438#endif
439}
440
441// Bitshift a u128 in place, left (signed_shift > 0) or right (signed_shift < 0)
442static inline void libdivide_u128_shift(uint64_t *u1, uint64_t *u0, int32_t signed_shift) {
443 if (signed_shift > 0) {
444 uint32_t shift = signed_shift;
445 *u1 <<= shift;
446 *u1 |= *u0 >> (64 - shift);
447 *u0 <<= shift;
448 }
449 else if (signed_shift < 0) {
450 uint32_t shift = -signed_shift;
451 *u0 >>= shift;
452 *u0 |= *u1 << (64 - shift);
453 *u1 >>= shift;
454 }
455}
456
457// Computes a 128 / 128 -> 64 bit division, with a 128 bit remainder.
458static 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) {
459#if defined(HAS_INT128_T) && \
460 defined(HAS_INT128_DIV)
461 __uint128_t ufull = u_hi;
462 __uint128_t vfull = v_hi;
463 ufull = (ufull << 64) | u_lo;
464 vfull = (vfull << 64) | v_lo;
465 uint64_t res = (uint64_t)(ufull / vfull);
466 __uint128_t remainder = ufull - (vfull * res);
467 *r_lo = (uint64_t)remainder;
468 *r_hi = (uint64_t)(remainder >> 64);
469 return res;
470#else
471 // Adapted from "Unsigned Doubleword Division" in Hacker's Delight
472 // We want to compute u / v
473 typedef struct { uint64_t hi; uint64_t lo; } u128_t;
474 u128_t u = {u_hi, u_lo};
475 u128_t v = {v_hi, v_lo};
476
477 if (v.hi == 0) {
478 // divisor v is a 64 bit value, so we just need one 128/64 division
479 // Note that we are simpler than Hacker's Delight here, because we know
480 // the quotient fits in 64 bits whereas Hacker's Delight demands a full
481 // 128 bit quotient
482 *r_hi = 0;
483 return libdivide_128_div_64_to_64(u.hi, u.lo, v.lo, r_lo);
484 }
485 // Here v >= 2**64
486 // We know that v.hi != 0, so count leading zeros is OK
487 // We have 0 <= n <= 63
488 uint32_t n = libdivide_count_leading_zeros64(v.hi);
489
490 // Normalize the divisor so its MSB is 1
491 u128_t v1t = v;
492 libdivide_u128_shift(&v1t.hi, &v1t.lo, n);
493 uint64_t v1 = v1t.hi; // i.e. v1 = v1t >> 64
494
495 // To ensure no overflow
496 u128_t u1 = u;
497 libdivide_u128_shift(&u1.hi, &u1.lo, -1);
498
499 // Get quotient from divide unsigned insn.
500 uint64_t rem_ignored;
501 uint64_t q1 = libdivide_128_div_64_to_64(u1.hi, u1.lo, v1, &rem_ignored);
502
503 // Undo normalization and division of u by 2.
504 u128_t q0 = {0, q1};
505 libdivide_u128_shift(&q0.hi, &q0.lo, n);
506 libdivide_u128_shift(&q0.hi, &q0.lo, -63);
507
508 // Make q0 correct or too small by 1
509 // Equivalent to `if (q0 != 0) q0 = q0 - 1;`
510 if (q0.hi != 0 || q0.lo != 0) {
511 q0.hi -= (q0.lo == 0); // borrow
512 q0.lo -= 1;
513 }
514
515 // Now q0 is correct.
516 // Compute q0 * v as q0v
517 // = (q0.hi << 64 + q0.lo) * (v.hi << 64 + v.lo)
518 // = (q0.hi * v.hi << 128) + (q0.hi * v.lo << 64) +
519 // (q0.lo * v.hi << 64) + q0.lo * v.lo)
520 // Each term is 128 bit
521 // High half of full product (upper 128 bits!) are dropped
522 u128_t q0v = {0, 0};
523 q0v.hi = q0.hi*v.lo + q0.lo*v.hi + libdivide_mullhi_u64(q0.lo, v.lo);
524 q0v.lo = q0.lo*v.lo;
525
526 // Compute u - q0v as u_q0v
527 // This is the remainder
528 u128_t u_q0v = u;
529 u_q0v.hi -= q0v.hi + (u.lo < q0v.lo); // second term is borrow
530 u_q0v.lo -= q0v.lo;
531
532 // Check if u_q0v >= v
533 // This checks if our remainder is larger than the divisor
534 if ((u_q0v.hi > v.hi) ||
535 (u_q0v.hi == v.hi && u_q0v.lo >= v.lo)) {
536 // Increment q0
537 q0.lo += 1;
538 q0.hi += (q0.lo == 0); // carry
539
540 // Subtract v from remainder
541 u_q0v.hi -= v.hi + (u_q0v.lo < v.lo);
542 u_q0v.lo -= v.lo;
543 }
544
545 *r_hi = u_q0v.hi;
546 *r_lo = u_q0v.lo;
547
548 LIBDIVIDE_ASSERT(q0.hi == 0);
549 return q0.lo;
550#endif
551}
552
553////////// UINT32
554
555static inline struct libdivide_u32_t libdivide_internal_u32_gen(uint32_t d, int branchfree) {
556 if (d == 0) {
557 LIBDIVIDE_ERROR("divider must be != 0");
558 }
559
560 struct libdivide_u32_t result;
561 uint32_t floor_log_2_d = 31 - libdivide_count_leading_zeros32(d);
562
563 // Power of 2
564 if ((d & (d - 1)) == 0) {
565 // We need to subtract 1 from the shift value in case of an unsigned
566 // branchfree divider because there is a hardcoded right shift by 1
567 // in its division algorithm. Because of this we also need to add back
568 // 1 in its recovery algorithm.
569 result.magic = 0;
570 result.more = (uint8_t)(floor_log_2_d - (branchfree != 0));
571 } else {
572 uint8_t more;
573 uint32_t rem, proposed_m;
574 proposed_m = libdivide_64_div_32_to_32(1U << floor_log_2_d, 0, d, &rem);
575
576 LIBDIVIDE_ASSERT(rem > 0 && rem < d);
577 const uint32_t e = d - rem;
578
579 // This power works if e < 2**floor_log_2_d.
580 if (!branchfree && (e < (1U << floor_log_2_d))) {
581 // This power works
582 more = floor_log_2_d;
583 } else {
584 // We have to use the general 33-bit algorithm. We need to compute
585 // (2**power) / d. However, we already have (2**(power-1))/d and
586 // its remainder. By doubling both, and then correcting the
587 // remainder, we can compute the larger division.
588 // don't care about overflow here - in fact, we expect it
589 proposed_m += proposed_m;
590 const uint32_t twice_rem = rem + rem;
591 if (twice_rem >= d || twice_rem < rem) proposed_m += 1;
592 more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;
593 }
594 result.magic = 1 + proposed_m;
595 result.more = more;
596 // result.more's shift should in general be ceil_log_2_d. But if we
597 // used the smaller power, we subtract one from the shift because we're
598 // using the smaller power. If we're using the larger power, we
599 // subtract one from the shift because it's taken care of by the add
600 // indicator. So floor_log_2_d happens to be correct in both cases.
601 }
602 return result;
603}
604
605struct libdivide_u32_t libdivide_u32_gen(uint32_t d) {
606 return libdivide_internal_u32_gen(d, 0);
607}
608
609struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t d) {
610 if (d == 1) {
611 LIBDIVIDE_ERROR("branchfree divider must be != 1");
612 }
613 struct libdivide_u32_t tmp = libdivide_internal_u32_gen(d, 1);
614 struct libdivide_u32_branchfree_t ret = {tmp.magic, (uint8_t)(tmp.more & LIBDIVIDE_32_SHIFT_MASK)};
615 return ret;
616}
617
618uint32_t libdivide_u32_do(uint32_t numer, const struct libdivide_u32_t *denom) {
619 uint8_t more = denom->more;
620 if (!denom->magic) {
621 return numer >> more;
622 }
623 else {
624 uint32_t q = libdivide_mullhi_u32(denom->magic, numer);
625 if (more & LIBDIVIDE_ADD_MARKER) {
626 uint32_t t = ((numer - q) >> 1) + q;
627 return t >> (more & LIBDIVIDE_32_SHIFT_MASK);
628 }
629 else {
630 // All upper bits are 0,
631 // don't need to mask them off.
632 return q >> more;
633 }
634 }
635}
636
637uint32_t libdivide_u32_branchfree_do(uint32_t numer, const struct libdivide_u32_branchfree_t *denom) {
638 uint32_t q = libdivide_mullhi_u32(denom->magic, numer);
639 uint32_t t = ((numer - q) >> 1) + q;
640 return t >> denom->more;
641}
642
643uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom) {
644 uint8_t more = denom->more;
645 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
646
647 if (!denom->magic) {
648 return 1U << shift;
649 } else if (!(more & LIBDIVIDE_ADD_MARKER)) {
650 // We compute q = n/d = n*m / 2^(32 + shift)
651 // Therefore we have d = 2^(32 + shift) / m
652 // We need to ceil it.
653 // We know d is not a power of 2, so m is not a power of 2,
654 // so we can just add 1 to the floor
655 uint32_t hi_dividend = 1U << shift;
656 uint32_t rem_ignored;
657 return 1 + libdivide_64_div_32_to_32(hi_dividend, 0, denom->magic, &rem_ignored);
658 } else {
659 // Here we wish to compute d = 2^(32+shift+1)/(m+2^32).
660 // Notice (m + 2^32) is a 33 bit number. Use 64 bit division for now
661 // Also note that shift may be as high as 31, so shift + 1 will
662 // overflow. So we have to compute it as 2^(32+shift)/(m+2^32), and
663 // then double the quotient and remainder.
664 uint64_t half_n = 1ULL << (32 + shift);
665 uint64_t d = (1ULL << 32) | denom->magic;
666 // Note that the quotient is guaranteed <= 32 bits, but the remainder
667 // may need 33!
668 uint32_t half_q = (uint32_t)(half_n / d);
669 uint64_t rem = half_n % d;
670 // We computed 2^(32+shift)/(m+2^32)
671 // Need to double it, and then add 1 to the quotient if doubling th
672 // remainder would increase the quotient.
673 // Note that rem<<1 cannot overflow, since rem < d and d is 33 bits
674 uint32_t full_q = half_q + half_q + ((rem<<1) >= d);
675
676 // We rounded down in gen (hence +1)
677 return full_q + 1;
678 }
679}
680
681uint32_t libdivide_u32_branchfree_recover(const struct libdivide_u32_branchfree_t *denom) {
682 uint8_t more = denom->more;
683 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
684
685 if (!denom->magic) {
686 return 1U << (shift + 1);
687 } else {
688 // Here we wish to compute d = 2^(32+shift+1)/(m+2^32).
689 // Notice (m + 2^32) is a 33 bit number. Use 64 bit division for now
690 // Also note that shift may be as high as 31, so shift + 1 will
691 // overflow. So we have to compute it as 2^(32+shift)/(m+2^32), and
692 // then double the quotient and remainder.
693 uint64_t half_n = 1ULL << (32 + shift);
694 uint64_t d = (1ULL << 32) | denom->magic;
695 // Note that the quotient is guaranteed <= 32 bits, but the remainder
696 // may need 33!
697 uint32_t half_q = (uint32_t)(half_n / d);
698 uint64_t rem = half_n % d;
699 // We computed 2^(32+shift)/(m+2^32)
700 // Need to double it, and then add 1 to the quotient if doubling th
701 // remainder would increase the quotient.
702 // Note that rem<<1 cannot overflow, since rem < d and d is 33 bits
703 uint32_t full_q = half_q + half_q + ((rem<<1) >= d);
704
705 // We rounded down in gen (hence +1)
706 return full_q + 1;
707 }
708}
709
710/////////// UINT64
711
712static inline struct libdivide_u64_t libdivide_internal_u64_gen(uint64_t d, int branchfree) {
713 if (d == 0) {
714 LIBDIVIDE_ERROR("divider must be != 0");
715 }
716
717 struct libdivide_u64_t result;
718 uint32_t floor_log_2_d = 63 - libdivide_count_leading_zeros64(d);
719
720 // Power of 2
721 if ((d & (d - 1)) == 0) {
722 // We need to subtract 1 from the shift value in case of an unsigned
723 // branchfree divider because there is a hardcoded right shift by 1
724 // in its division algorithm. Because of this we also need to add back
725 // 1 in its recovery algorithm.
726 result.magic = 0;
727 result.more = (uint8_t)(floor_log_2_d - (branchfree != 0));
728 } else {
729 uint64_t proposed_m, rem;
730 uint8_t more;
731 // (1 << (64 + floor_log_2_d)) / d
732 proposed_m = libdivide_128_div_64_to_64(1ULL << floor_log_2_d, 0, d, &rem);
733
734 LIBDIVIDE_ASSERT(rem > 0 && rem < d);
735 const uint64_t e = d - rem;
736
737 // This power works if e < 2**floor_log_2_d.
738 if (!branchfree && e < (1ULL << floor_log_2_d)) {
739 // This power works
740 more = floor_log_2_d;
741 } else {
742 // We have to use the general 65-bit algorithm. We need to compute
743 // (2**power) / d. However, we already have (2**(power-1))/d and
744 // its remainder. By doubling both, and then correcting the
745 // remainder, we can compute the larger division.
746 // don't care about overflow here - in fact, we expect it
747 proposed_m += proposed_m;
748 const uint64_t twice_rem = rem + rem;
749 if (twice_rem >= d || twice_rem < rem) proposed_m += 1;
750 more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;
751 }
752 result.magic = 1 + proposed_m;
753 result.more = more;
754 // result.more's shift should in general be ceil_log_2_d. But if we
755 // used the smaller power, we subtract one from the shift because we're
756 // using the smaller power. If we're using the larger power, we
757 // subtract one from the shift because it's taken care of by the add
758 // indicator. So floor_log_2_d happens to be correct in both cases,
759 // which is why we do it outside of the if statement.
760 }
761 return result;
762}
763
764struct libdivide_u64_t libdivide_u64_gen(uint64_t d) {
765 return libdivide_internal_u64_gen(d, 0);
766}
767
768struct libdivide_u64_branchfree_t libdivide_u64_branchfree_gen(uint64_t d) {
769 if (d == 1) {
770 LIBDIVIDE_ERROR("branchfree divider must be != 1");
771 }
772 struct libdivide_u64_t tmp = libdivide_internal_u64_gen(d, 1);
773 struct libdivide_u64_branchfree_t ret = {tmp.magic, (uint8_t)(tmp.more & LIBDIVIDE_64_SHIFT_MASK)};
774 return ret;
775}
776
777uint64_t libdivide_u64_do(uint64_t numer, const struct libdivide_u64_t *denom) {
778 uint8_t more = denom->more;
779 if (!denom->magic) {
780 return numer >> more;
781 }
782 else {
783 uint64_t q = libdivide_mullhi_u64(denom->magic, numer);
784 if (more & LIBDIVIDE_ADD_MARKER) {
785 uint64_t t = ((numer - q) >> 1) + q;
786 return t >> (more & LIBDIVIDE_64_SHIFT_MASK);
787 }
788 else {
789 // All upper bits are 0,
790 // don't need to mask them off.
791 return q >> more;
792 }
793 }
794}
795
796uint64_t libdivide_u64_branchfree_do(uint64_t numer, const struct libdivide_u64_branchfree_t *denom) {
797 uint64_t q = libdivide_mullhi_u64(denom->magic, numer);
798 uint64_t t = ((numer - q) >> 1) + q;
799 return t >> denom->more;
800}
801
802uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom) {
803 uint8_t more = denom->more;
804 uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
805
806 if (!denom->magic) {
807 return 1ULL << shift;
808 } else if (!(more & LIBDIVIDE_ADD_MARKER)) {
809 // We compute q = n/d = n*m / 2^(64 + shift)
810 // Therefore we have d = 2^(64 + shift) / m
811 // We need to ceil it.
812 // We know d is not a power of 2, so m is not a power of 2,
813 // so we can just add 1 to the floor
814 uint64_t hi_dividend = 1ULL << shift;
815 uint64_t rem_ignored;
816 return 1 + libdivide_128_div_64_to_64(hi_dividend, 0, denom->magic, &rem_ignored);
817 } else {
818 // Here we wish to compute d = 2^(64+shift+1)/(m+2^64).
819 // Notice (m + 2^64) is a 65 bit number. This gets hairy. See
820 // libdivide_u32_recover for more on what we do here.
821 // TODO: do something better than 128 bit math
822
823 // Full n is a (potentially) 129 bit value
824 // half_n is a 128 bit value
825 // Compute the hi half of half_n. Low half is 0.
826 uint64_t half_n_hi = 1ULL << shift, half_n_lo = 0;
827 // d is a 65 bit value. The high bit is always set to 1.
828 const uint64_t d_hi = 1, d_lo = denom->magic;
829 // Note that the quotient is guaranteed <= 64 bits,
830 // but the remainder may need 65!
831 uint64_t r_hi, r_lo;
832 uint64_t half_q = libdivide_128_div_128_to_64(half_n_hi, half_n_lo, d_hi, d_lo, &r_hi, &r_lo);
833 // We computed 2^(64+shift)/(m+2^64)
834 // Double the remainder ('dr') and check if that is larger than d
835 // Note that d is a 65 bit value, so r1 is small and so r1 + r1
836 // cannot overflow
837 uint64_t dr_lo = r_lo + r_lo;
838 uint64_t dr_hi = r_hi + r_hi + (dr_lo < r_lo); // last term is carry
839 int dr_exceeds_d = (dr_hi > d_hi) || (dr_hi == d_hi && dr_lo >= d_lo);
840 uint64_t full_q = half_q + half_q + (dr_exceeds_d ? 1 : 0);
841 return full_q + 1;
842 }
843}
844
845uint64_t libdivide_u64_branchfree_recover(const struct libdivide_u64_branchfree_t *denom) {
846 uint8_t more = denom->more;
847 uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
848
849 if (!denom->magic) {
850 return 1ULL << (shift + 1);
851 } else {
852 // Here we wish to compute d = 2^(64+shift+1)/(m+2^64).
853 // Notice (m + 2^64) is a 65 bit number. This gets hairy. See
854 // libdivide_u32_recover for more on what we do here.
855 // TODO: do something better than 128 bit math
856
857 // Full n is a (potentially) 129 bit value
858 // half_n is a 128 bit value
859 // Compute the hi half of half_n. Low half is 0.
860 uint64_t half_n_hi = 1ULL << shift, half_n_lo = 0;
861 // d is a 65 bit value. The high bit is always set to 1.
862 const uint64_t d_hi = 1, d_lo = denom->magic;
863 // Note that the quotient is guaranteed <= 64 bits,
864 // but the remainder may need 65!
865 uint64_t r_hi, r_lo;
866 uint64_t half_q = libdivide_128_div_128_to_64(half_n_hi, half_n_lo, d_hi, d_lo, &r_hi, &r_lo);
867 // We computed 2^(64+shift)/(m+2^64)
868 // Double the remainder ('dr') and check if that is larger than d
869 // Note that d is a 65 bit value, so r1 is small and so r1 + r1
870 // cannot overflow
871 uint64_t dr_lo = r_lo + r_lo;
872 uint64_t dr_hi = r_hi + r_hi + (dr_lo < r_lo); // last term is carry
873 int dr_exceeds_d = (dr_hi > d_hi) || (dr_hi == d_hi && dr_lo >= d_lo);
874 uint64_t full_q = half_q + half_q + (dr_exceeds_d ? 1 : 0);
875 return full_q + 1;
876 }
877}
878
879/////////// SINT32
880
881static inline struct libdivide_s32_t libdivide_internal_s32_gen(int32_t d, int branchfree) {
882 if (d == 0) {
883 LIBDIVIDE_ERROR("divider must be != 0");
884 }
885
886 struct libdivide_s32_t result;
887
888 // If d is a power of 2, or negative a power of 2, we have to use a shift.
889 // This is especially important because the magic algorithm fails for -1.
890 // To check if d is a power of 2 or its inverse, it suffices to check
891 // whether its absolute value has exactly one bit set. This works even for
892 // INT_MIN, because abs(INT_MIN) == INT_MIN, and INT_MIN has one bit set
893 // and is a power of 2.
894 uint32_t ud = (uint32_t)d;
895 uint32_t absD = (d < 0) ? -ud : ud;
896 uint32_t floor_log_2_d = 31 - libdivide_count_leading_zeros32(absD);
897 // check if exactly one bit is set,
898 // don't care if absD is 0 since that's divide by zero
899 if ((absD & (absD - 1)) == 0) {
900 // Branchfree and normal paths are exactly the same
901 result.magic = 0;
902 result.more = floor_log_2_d | (d < 0 ? LIBDIVIDE_NEGATIVE_DIVISOR : 0);
903 } else {
904 LIBDIVIDE_ASSERT(floor_log_2_d >= 1);
905
906 uint8_t more;
907 // the dividend here is 2**(floor_log_2_d + 31), so the low 32 bit word
908 // is 0 and the high word is floor_log_2_d - 1
909 uint32_t rem, proposed_m;
910 proposed_m = libdivide_64_div_32_to_32(1U << (floor_log_2_d - 1), 0, absD, &rem);
911 const uint32_t e = absD - rem;
912
913 // We are going to start with a power of floor_log_2_d - 1.
914 // This works if works if e < 2**floor_log_2_d.
915 if (!branchfree && e < (1U << floor_log_2_d)) {
916 // This power works
917 more = floor_log_2_d - 1;
918 } else {
919 // We need to go one higher. This should not make proposed_m
920 // overflow, but it will make it negative when interpreted as an
921 // int32_t.
922 proposed_m += proposed_m;
923 const uint32_t twice_rem = rem + rem;
924 if (twice_rem >= absD || twice_rem < rem) proposed_m += 1;
925 more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;
926 }
927
928 proposed_m += 1;
929 int32_t magic = (int32_t)proposed_m;
930
931 // Mark if we are negative. Note we only negate the magic number in the
932 // branchfull case.
933 if (d < 0) {
934 more |= LIBDIVIDE_NEGATIVE_DIVISOR;
935 if (!branchfree) {
936 magic = -magic;
937 }
938 }
939
940 result.more = more;
941 result.magic = magic;
942 }
943 return result;
944}
945
946struct libdivide_s32_t libdivide_s32_gen(int32_t d) {
947 return libdivide_internal_s32_gen(d, 0);
948}
949
950struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t d) {
951 struct libdivide_s32_t tmp = libdivide_internal_s32_gen(d, 1);
952 struct libdivide_s32_branchfree_t result = {tmp.magic, tmp.more};
953 return result;
954}
955
956int32_t libdivide_s32_do(int32_t numer, const struct libdivide_s32_t *denom) {
957 uint8_t more = denom->more;
958 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
959
960 if (!denom->magic) {
961 uint32_t sign = (int8_t)more >> 7;
962 uint32_t mask = (1U << shift) - 1;
963 uint32_t uq = numer + ((numer >> 31) & mask);
964 int32_t q = (int32_t)uq;
965 q >>= shift;
966 q = (q ^ sign) - sign;
967 return q;
968 } else {
969 uint32_t uq = (uint32_t)libdivide_mullhi_s32(denom->magic, numer);
970 if (more & LIBDIVIDE_ADD_MARKER) {
971 // must be arithmetic shift and then sign extend
972 int32_t sign = (int8_t)more >> 7;
973 // q += (more < 0 ? -numer : numer)
974 // cast required to avoid UB
975 uq += ((uint32_t)numer ^ sign) - sign;
976 }
977 int32_t q = (int32_t)uq;
978 q >>= shift;
979 q += (q < 0);
980 return q;
981 }
982}
983
984int32_t libdivide_s32_branchfree_do(int32_t numer, const struct libdivide_s32_branchfree_t *denom) {
985 uint8_t more = denom->more;
986 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
987 // must be arithmetic shift and then sign extend
988 int32_t sign = (int8_t)more >> 7;
989 int32_t magic = denom->magic;
990 int32_t q = libdivide_mullhi_s32(magic, numer);
991 q += numer;
992
993 // If q is non-negative, we have nothing to do
994 // If q is negative, we want to add either (2**shift)-1 if d is a power of
995 // 2, or (2**shift) if it is not a power of 2
996 uint32_t is_power_of_2 = (magic == 0);
997 uint32_t q_sign = (uint32_t)(q >> 31);
998 q += q_sign & ((1U << shift) - is_power_of_2);
999
1000 // Now arithmetic right shift
1001 q >>= shift;
1002 // Negate if needed
1003 q = (q ^ sign) - sign;
1004
1005 return q;
1006}
1007
1008int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom) {
1009 uint8_t more = denom->more;
1010 uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK;
1011 if (!denom->magic) {
1012 uint32_t absD = 1U << shift;
1013 if (more & LIBDIVIDE_NEGATIVE_DIVISOR) {
1014 absD = -absD;
1015 }
1016 return (int32_t)absD;
1017 } else {
1018 // Unsigned math is much easier
1019 // We negate the magic number only in the branchfull case, and we don't
1020 // know which case we're in. However we have enough information to
1021 // determine the correct sign of the magic number. The divisor was
1022 // negative if LIBDIVIDE_NEGATIVE_DIVISOR is set. If ADD_MARKER is set,
1023 // the magic number's sign is opposite that of the divisor.
1024 // We want to compute the positive magic number.
1025 int negative_divisor = (more & LIBDIVIDE_NEGATIVE_DIVISOR);
1026 int magic_was_negated = (more & LIBDIVIDE_ADD_MARKER)
1027 ? denom->magic > 0 : denom->magic < 0;
1028
1029 // Handle the power of 2 case (including branchfree)
1030 if (denom->magic == 0) {
1031 int32_t result = 1U << shift;
1032 return negative_divisor ? -result : result;
1033 }
1034
1035 uint32_t d = (uint32_t)(magic_was_negated ? -denom->magic : denom->magic);
1036 uint64_t n = 1ULL << (32 + shift); // this shift cannot exceed 30
1037 uint32_t q = (uint32_t)(n / d);
1038 int32_t result = (int32_t)q;
1039 result += 1;
1040 return negative_divisor ? -result : result;
1041 }
1042}
1043
1044int32_t libdivide_s32_branchfree_recover(const struct libdivide_s32_branchfree_t *denom) {
1045 return libdivide_s32_recover((const struct libdivide_s32_t *)denom);
1046}
1047
1048///////////// SINT64
1049
1050static inline struct libdivide_s64_t libdivide_internal_s64_gen(int64_t d, int branchfree) {
1051 if (d == 0) {
1052 LIBDIVIDE_ERROR("divider must be != 0");
1053 }
1054
1055 struct libdivide_s64_t result;
1056
1057 // If d is a power of 2, or negative a power of 2, we have to use a shift.
1058 // This is especially important because the magic algorithm fails for -1.
1059 // To check if d is a power of 2 or its inverse, it suffices to check
1060 // whether its absolute value has exactly one bit set. This works even for
1061 // INT_MIN, because abs(INT_MIN) == INT_MIN, and INT_MIN has one bit set
1062 // and is a power of 2.
1063 uint64_t ud = (uint64_t)d;
1064 uint64_t absD = (d < 0) ? -ud : ud;
1065 uint32_t floor_log_2_d = 63 - libdivide_count_leading_zeros64(absD);
1066 // check if exactly one bit is set,
1067 // don't care if absD is 0 since that's divide by zero
1068 if ((absD & (absD - 1)) == 0) {
1069 // Branchfree and non-branchfree cases are the same
1070 result.magic = 0;
1071 result.more = floor_log_2_d | (d < 0 ? LIBDIVIDE_NEGATIVE_DIVISOR : 0);
1072 } else {
1073 // the dividend here is 2**(floor_log_2_d + 63), so the low 64 bit word
1074 // is 0 and the high word is floor_log_2_d - 1
1075 uint8_t more;
1076 uint64_t rem, proposed_m;
1077 proposed_m = libdivide_128_div_64_to_64(1ULL << (floor_log_2_d - 1), 0, absD, &rem);
1078 const uint64_t e = absD - rem;
1079
1080 // We are going to start with a power of floor_log_2_d - 1.
1081 // This works if works if e < 2**floor_log_2_d.
1082 if (!branchfree && e < (1ULL << floor_log_2_d)) {
1083 // This power works
1084 more = floor_log_2_d - 1;
1085 } else {
1086 // We need to go one higher. This should not make proposed_m
1087 // overflow, but it will make it negative when interpreted as an
1088 // int32_t.
1089 proposed_m += proposed_m;
1090 const uint64_t twice_rem = rem + rem;
1091 if (twice_rem >= absD || twice_rem < rem) proposed_m += 1;
1092 // note that we only set the LIBDIVIDE_NEGATIVE_DIVISOR bit if we
1093 // also set ADD_MARKER this is an annoying optimization that
1094 // enables algorithm #4 to avoid the mask. However we always set it
1095 // in the branchfree case
1096 more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;
1097 }
1098 proposed_m += 1;
1099 int64_t magic = (int64_t)proposed_m;
1100
1101 // Mark if we are negative
1102 if (d < 0) {
1103 more |= LIBDIVIDE_NEGATIVE_DIVISOR;
1104 if (!branchfree) {
1105 magic = -magic;
1106 }
1107 }
1108
1109 result.more = more;
1110 result.magic = magic;
1111 }
1112 return result;
1113}
1114
1115struct libdivide_s64_t libdivide_s64_gen(int64_t d) {
1116 return libdivide_internal_s64_gen(d, 0);
1117}
1118
1119struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t d) {
1120 struct libdivide_s64_t tmp = libdivide_internal_s64_gen(d, 1);
1121 struct libdivide_s64_branchfree_t ret = {tmp.magic, tmp.more};
1122 return ret;
1123}
1124
1125int64_t libdivide_s64_do(int64_t numer, const struct libdivide_s64_t *denom) {
1126 uint8_t more = denom->more;
1127 uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
1128
1129 if (!denom->magic) { // shift path
1130 uint64_t mask = (1ULL << shift) - 1;
1131 uint64_t uq = numer + ((numer >> 63) & mask);
1132 int64_t q = (int64_t)uq;
1133 q >>= shift;
1134 // must be arithmetic shift and then sign-extend
1135 int64_t sign = (int8_t)more >> 7;
1136 q = (q ^ sign) - sign;
1137 return q;
1138 } else {
1139 uint64_t uq = (uint64_t)libdivide_mullhi_s64(denom->magic, numer);
1140 if (more & LIBDIVIDE_ADD_MARKER) {
1141 // must be arithmetic shift and then sign extend
1142 int64_t sign = (int8_t)more >> 7;
1143 // q += (more < 0 ? -numer : numer)
1144 // cast required to avoid UB
1145 uq += ((uint64_t)numer ^ sign) - sign;
1146 }
1147 int64_t q = (int64_t)uq;
1148 q >>= shift;
1149 q += (q < 0);
1150 return q;
1151 }
1152}
1153
1154int64_t libdivide_s64_branchfree_do(int64_t numer, const struct libdivide_s64_branchfree_t *denom) {
1155 uint8_t more = denom->more;
1156 uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
1157 // must be arithmetic shift and then sign extend
1158 int64_t sign = (int8_t)more >> 7;
1159 int64_t magic = denom->magic;
1160 int64_t q = libdivide_mullhi_s64(magic, numer);
1161 q += numer;
1162
1163 // If q is non-negative, we have nothing to do.
1164 // If q is negative, we want to add either (2**shift)-1 if d is a power of
1165 // 2, or (2**shift) if it is not a power of 2.
1166 uint64_t is_power_of_2 = (magic == 0);
1167 uint64_t q_sign = (uint64_t)(q >> 63);
1168 q += q_sign & ((1ULL << shift) - is_power_of_2);
1169
1170 // Arithmetic right shift
1171 q >>= shift;
1172 // Negate if needed
1173 q = (q ^ sign) - sign;
1174
1175 return q;
1176}
1177
1178int64_t libdivide_s64_recover(const struct libdivide_s64_t *denom) {
1179 uint8_t more = denom->more;
1180 uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK;
1181 if (denom->magic == 0) { // shift path
1182 uint64_t absD = 1ULL << shift;
1183 if (more & LIBDIVIDE_NEGATIVE_DIVISOR) {
1184 absD = -absD;
1185 }
1186 return (int64_t)absD;
1187 } else {
1188 // Unsigned math is much easier
1189 int negative_divisor = (more & LIBDIVIDE_NEGATIVE_DIVISOR);
1190 int magic_was_negated = (more & LIBDIVIDE_ADD_MARKER)
1191 ? denom->magic > 0 : denom->magic < 0;
1192
1193 uint64_t d = (uint64_t)(magic_was_negated ? -denom->magic : denom->magic);
1194 uint64_t n_hi = 1ULL << shift, n_lo = 0;
1195 uint64_t rem_ignored;
1196 uint64_t q = libdivide_128_div_64_to_64(n_hi, n_lo, d, &rem_ignored);
1197 int64_t result = (int64_t)(q + 1);
1198 if (negative_divisor) {
1199 result = -result;
1200 }