From 7c455c39560f753964c12863e1413ef603f846f4 Mon Sep 17 00:00:00 2001 From: Daniel Micay Date: Sun, 18 Aug 2019 06:56:52 -0400 Subject: [PATCH] update libdivide to 2.0 --- third_party/libdivide.h | 1775 +++++++++++++++++++-------------------- 1 file changed, 870 insertions(+), 905 deletions(-) diff --git a/third_party/libdivide.h b/third_party/libdivide.h index d1c3f8d..7969d86 100644 --- a/third_party/libdivide.h +++ b/third_party/libdivide.h @@ -1,5 +1,8 @@ -// libdivide.h -// Copyright 2010 - 2019 ridiculous_fish +// libdivide.h - Optimized integer division +// https://libdivide.com +// +// Copyright (C) 2010 - 2019 ridiculous_fish, +// Copyright (C) 2016 - 2019 Kim Walisch, // // libdivide is dual-licensed under the Boost or zlib licenses. // You may use libdivide under the terms of either of these. @@ -8,9 +11,9 @@ #ifndef LIBDIVIDE_H #define LIBDIVIDE_H -#define LIBDIVIDE_VERSION "1.1" -#define LIBDIVIDE_VERSION_MAJOR 1 -#define LIBDIVIDE_VERSION_MINOR 1 +#define LIBDIVIDE_VERSION "2.0" +#define LIBDIVIDE_VERSION_MAJOR 2 +#define LIBDIVIDE_VERSION_MINOR 0 #include @@ -22,16 +25,20 @@ #include #endif -#if defined(LIBDIVIDE_USE_SSE2) +#if defined(LIBDIVIDE_AVX512) + #include +#elif defined(LIBDIVIDE_AVX2) + #include +#elif defined(LIBDIVIDE_SSE2) #include #endif -// libdivide may use the pmuldq (vector signed 32x32->64 mult instruction) -// which is in SSE 4.1. However, signed multiplication can be emulated -// efficiently with unsigned multiplication, and SSE 4.1 is currently rare, -// so it is OK to not turn this on. -#if defined(LIBDIVIDE_USE_SSE4_1) - #include +#if defined(_MSC_VER) + #include + // disable warning C4146: unary minus operator applied + // to unsigned type, result still unsigned + #pragma warning(disable: 4146) + #define LIBDIVIDE_VC #endif #if !defined(__has_builtin) @@ -40,6 +47,10 @@ #if defined(__SIZEOF_INT128__) #define HAS_INT128_T + // clang-cl on Windows does not yet support 128-bit division + #if !(defined(__clang__) && defined(LIBDIVIDE_VC)) + #define HAS_INT128_DIV + #endif #endif #if defined(__x86_64__) || defined(_M_X64) @@ -60,26 +71,6 @@ #define LIBDIVIDE_FUNCTION __func__ #endif -#if defined(_MSC_VER) - #include - // disable warning C4146: unary minus operator applied - // to unsigned type, result still unsigned - #pragma warning(disable: 4146) - #define LIBDIVIDE_VC - - // _udiv128() is available in Visual Studio 2019 - // or later on the x64 CPU architecture - #if defined(LIBDIVIDE_X86_64) && _MSC_VER >= 1920 - #if !defined(__has_include) - #include - #define LIBDIVIDE_VC_UDIV128 - #elif __has_include() - #include - #define LIBDIVIDE_VC_UDIV128 - #endif - #endif -#endif - #define LIBDIVIDE_ERROR(msg) \ do { \ fprintf(stderr, "libdivide.h:%d: %s(): Error: %s\n", \ @@ -101,53 +92,9 @@ #endif #ifdef __cplusplus -// We place libdivide within the libdivide namespace, and that goes in an -// anonymous namespace so that the functions are only visible to files that -// #include this header and don't get external linkage. At least that's the -// theory. -namespace { namespace libdivide { #endif -// Explanation of "more" field: bit 6 is whether to use shift path. If we are -// using the shift path, bit 7 is whether the divisor is negative in the signed -// case; in the unsigned case it is 0. Bits 0-4 is shift value (for shift -// path or mult path). In 32 bit case, bit 5 is always 0. We use bit 7 as the -// "negative divisor indicator" so that we can use sign extension to -// efficiently go to a full-width -1. -// -// u32: [0-4] shift value -// [5] ignored -// [6] add indicator -// [7] shift path -// -// s32: [0-4] shift value -// [5] shift path -// [6] add indicator -// [7] indicates negative divisor -// -// u64: [0-5] shift value -// [6] add indicator -// [7] shift path -// -// s64: [0-5] shift value -// [6] add indicator -// [7] indicates negative divisor -// magic number of 0 indicates shift path (we ran out of bits!) -// -// In s32 and s64 branchfree modes, the magic number is negated according to -// whether the divisor is negated. In branchfree strategy, it is not negated. - -enum { - LIBDIVIDE_32_SHIFT_MASK = 0x1F, - LIBDIVIDE_64_SHIFT_MASK = 0x3F, - LIBDIVIDE_ADD_MARKER = 0x40, - LIBDIVIDE_U32_SHIFT_PATH = 0x80, - LIBDIVIDE_U64_SHIFT_PATH = 0x80, - LIBDIVIDE_S32_SHIFT_PATH = 0x20, - LIBDIVIDE_NEGATIVE_DIVISOR = 0x80 -}; - // pack divider structs to prevent compilers from padding. // This reduces memory usage by up to 43% when using a large // array of libdivide dividers and improves performance @@ -196,71 +143,74 @@ struct libdivide_s64_branchfree_t { #pragma pack(pop) -#ifndef LIBDIVIDE_API - #ifdef __cplusplus - // In C++, we don't want our public functions to be static, because - // they are arguments to templates and static functions can't do that. - // They get internal linkage through virtue of the anonymous namespace. - // In C, they should be static. - #define LIBDIVIDE_API - #else - #define LIBDIVIDE_API static inline - #endif -#endif +// Explanation of the "more" field: +// +// * Bits 0-5 is the shift value (for shift path or mult path). +// * Bit 6 is the add indicator for mult path. +// * Bit 7 is set if the divisor is negative. We use bit 7 as the negative +// divisor indicator so that we can efficiently use sign extension to +// create a bitmask with all bits set to 1 (if the divisor is negative) +// or 0 (if the divisor is positive). +// +// u32: [0-4] shift value +// [5] ignored +// [6] add indicator +// magic number of 0 indicates shift path +// +// s32: [0-4] shift value +// [5] ignored +// [6] add indicator +// [7] indicates negative divisor +// magic number of 0 indicates shift path +// +// u64: [0-5] shift value +// [6] add indicator +// magic number of 0 indicates shift path +// +// s64: [0-5] shift value +// [6] add indicator +// [7] indicates negative divisor +// magic number of 0 indicates shift path +// +// In s32 and s64 branchfree modes, the magic number is negated according to +// whether the divisor is negated. In branchfree strategy, it is not negated. -LIBDIVIDE_API struct libdivide_s32_t libdivide_s32_gen(int32_t y); -LIBDIVIDE_API struct libdivide_u32_t libdivide_u32_gen(uint32_t y); -LIBDIVIDE_API struct libdivide_s64_t libdivide_s64_gen(int64_t y); -LIBDIVIDE_API struct libdivide_u64_t libdivide_u64_gen(uint64_t y); +enum { + LIBDIVIDE_32_SHIFT_MASK = 0x1F, + LIBDIVIDE_64_SHIFT_MASK = 0x3F, + LIBDIVIDE_ADD_MARKER = 0x40, + LIBDIVIDE_NEGATIVE_DIVISOR = 0x80 +}; -LIBDIVIDE_API struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t y); -LIBDIVIDE_API struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t y); -LIBDIVIDE_API struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t y); -LIBDIVIDE_API struct libdivide_u64_branchfree_t libdivide_u64_branchfree_gen(uint64_t y); +static inline struct libdivide_s32_t libdivide_s32_gen(int32_t d); +static inline struct libdivide_u32_t libdivide_u32_gen(uint32_t d); +static inline struct libdivide_s64_t libdivide_s64_gen(int64_t d); +static inline struct libdivide_u64_t libdivide_u64_gen(uint64_t d); -LIBDIVIDE_API int32_t libdivide_s32_do(int32_t numer, const struct libdivide_s32_t *denom); -LIBDIVIDE_API uint32_t libdivide_u32_do(uint32_t numer, const struct libdivide_u32_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_do(int64_t numer, const struct libdivide_s64_t *denom); -LIBDIVIDE_API uint64_t libdivide_u64_do(uint64_t y, const struct libdivide_u64_t *denom); +static inline struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t d); +static inline struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t d); +static inline struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t d); +static inline struct libdivide_u64_branchfree_t libdivide_u64_branchfree_gen(uint64_t d); -LIBDIVIDE_API int32_t libdivide_s32_branchfree_do(int32_t numer, const struct libdivide_s32_branchfree_t *denom); -LIBDIVIDE_API uint32_t libdivide_u32_branchfree_do(uint32_t numer, const struct libdivide_u32_branchfree_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_branchfree_do(int64_t numer, const struct libdivide_s64_branchfree_t *denom); -LIBDIVIDE_API uint64_t libdivide_u64_branchfree_do(uint64_t y, const struct libdivide_u64_branchfree_t *denom); +static inline int32_t libdivide_s32_do(int32_t numer, const struct libdivide_s32_t *denom); +static inline uint32_t libdivide_u32_do(uint32_t numer, const struct libdivide_u32_t *denom); +static inline int64_t libdivide_s64_do(int64_t numer, const struct libdivide_s64_t *denom); +static inline uint64_t libdivide_u64_do(uint64_t numer, const struct libdivide_u64_t *denom); -LIBDIVIDE_API int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom); -LIBDIVIDE_API uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_recover(const struct libdivide_s64_t *denom); -LIBDIVIDE_API uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom); +static inline int32_t libdivide_s32_branchfree_do(int32_t numer, const struct libdivide_s32_branchfree_t *denom); +static inline uint32_t libdivide_u32_branchfree_do(uint32_t numer, const struct libdivide_u32_branchfree_t *denom); +static inline int64_t libdivide_s64_branchfree_do(int64_t numer, const struct libdivide_s64_branchfree_t *denom); +static inline uint64_t libdivide_u64_branchfree_do(uint64_t numer, const struct libdivide_u64_branchfree_t *denom); -LIBDIVIDE_API int32_t libdivide_s32_branchfree_recover(const struct libdivide_s32_branchfree_t *denom); -LIBDIVIDE_API uint32_t libdivide_u32_branchfree_recover(const struct libdivide_u32_branchfree_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_branchfree_recover(const struct libdivide_s64_branchfree_t *denom); -LIBDIVIDE_API uint64_t libdivide_u64_branchfree_recover(const struct libdivide_u64_branchfree_t *denom); +static inline int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom); +static inline uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom); +static inline int64_t libdivide_s64_recover(const struct libdivide_s64_t *denom); +static inline uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom); -LIBDIVIDE_API int libdivide_u32_get_algorithm(const struct libdivide_u32_t *denom); -LIBDIVIDE_API uint32_t libdivide_u32_do_alg0(uint32_t numer, const struct libdivide_u32_t *denom); -LIBDIVIDE_API uint32_t libdivide_u32_do_alg1(uint32_t numer, const struct libdivide_u32_t *denom); -LIBDIVIDE_API uint32_t libdivide_u32_do_alg2(uint32_t numer, const struct libdivide_u32_t *denom); - -LIBDIVIDE_API int libdivide_u64_get_algorithm(const struct libdivide_u64_t *denom); -LIBDIVIDE_API uint64_t libdivide_u64_do_alg0(uint64_t numer, const struct libdivide_u64_t *denom); -LIBDIVIDE_API uint64_t libdivide_u64_do_alg1(uint64_t numer, const struct libdivide_u64_t *denom); -LIBDIVIDE_API uint64_t libdivide_u64_do_alg2(uint64_t numer, const struct libdivide_u64_t *denom); - -LIBDIVIDE_API int libdivide_s32_get_algorithm(const struct libdivide_s32_t *denom); -LIBDIVIDE_API int32_t libdivide_s32_do_alg0(int32_t numer, const struct libdivide_s32_t *denom); -LIBDIVIDE_API int32_t libdivide_s32_do_alg1(int32_t numer, const struct libdivide_s32_t *denom); -LIBDIVIDE_API int32_t libdivide_s32_do_alg2(int32_t numer, const struct libdivide_s32_t *denom); -LIBDIVIDE_API int32_t libdivide_s32_do_alg3(int32_t numer, const struct libdivide_s32_t *denom); -LIBDIVIDE_API int32_t libdivide_s32_do_alg4(int32_t numer, const struct libdivide_s32_t *denom); - -LIBDIVIDE_API int libdivide_s64_get_algorithm(const struct libdivide_s64_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_do_alg0(int64_t numer, const struct libdivide_s64_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_do_alg1(int64_t numer, const struct libdivide_s64_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_do_alg2(int64_t numer, const struct libdivide_s64_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_do_alg3(int64_t numer, const struct libdivide_s64_t *denom); -LIBDIVIDE_API int64_t libdivide_s64_do_alg4(int64_t numer, const struct libdivide_s64_t *denom); +static inline int32_t libdivide_s32_branchfree_recover(const struct libdivide_s32_branchfree_t *denom); +static inline uint32_t libdivide_u32_branchfree_recover(const struct libdivide_u32_branchfree_t *denom); +static inline int64_t libdivide_s64_branchfree_recover(const struct libdivide_s64_branchfree_t *denom); +static inline uint64_t libdivide_u64_branchfree_recover(const struct libdivide_u64_branchfree_t *denom); //////// Internal Utility Functions @@ -277,7 +227,7 @@ static inline int32_t libdivide_mullhi_s32(int32_t x, int32_t y) { return (int32_t)(rl >> 32); } -static uint64_t libdivide_mullhi_u64(uint64_t x, uint64_t y) { +static inline uint64_t libdivide_mullhi_u64(uint64_t x, uint64_t y) { #if defined(LIBDIVIDE_VC) && \ defined(LIBDIVIDE_X86_64) return __umulh(x, y); @@ -370,7 +320,7 @@ static inline int32_t libdivide_count_leading_zeros64(uint64_t val) { // libdivide_64_div_32_to_32: divides a 64-bit uint {u1, u0} by a 32-bit // uint {v}. The result must fit in 32 bits. // Returns the quotient directly and the remainder in *r -static uint32_t libdivide_64_div_32_to_32(uint32_t u1, uint32_t u0, uint32_t v, uint32_t *r) { +static inline uint32_t libdivide_64_div_32_to_32(uint32_t u1, uint32_t u0, uint32_t v, uint32_t *r) { #if (defined(LIBDIVIDE_i386) || defined(LIBDIVIDE_X86_64)) && \ defined(LIBDIVIDE_GCC_STYLE_ASM) uint32_t result; @@ -391,20 +341,16 @@ static uint32_t libdivide_64_div_32_to_32(uint32_t u1, uint32_t u0, uint32_t v, // uint {v}. The result must fit in 64 bits. // Returns the quotient directly and the remainder in *r static uint64_t libdivide_128_div_64_to_64(uint64_t u1, uint64_t u0, uint64_t v, uint64_t *r) { -#if defined(LIBDIVIDE_VC_UDIV128) - return _udiv128(u1, u0, v, r); - -#elif defined(LIBDIVIDE_X86_64) && \ - defined(LIBDIVIDE_GCC_STYLE_ASM) - +#if defined(LIBDIVIDE_X86_64) && \ + defined(LIBDIVIDE_GCC_STYLE_ASM) uint64_t result; __asm__("divq %[v]" : "=a"(result), "=d"(*r) : [v] "r"(v), "a"(u0), "d"(u1) ); return result; - -#elif defined(HAS_INT128_T) +#elif defined(HAS_INT128_T) && \ + defined(HAS_INT128_DIV) __uint128_t n = ((__uint128_t)u1 << 64) | u0; uint64_t result = (uint64_t)(n / v); *r = (uint64_t)(n - result * (__uint128_t)v); @@ -491,7 +437,8 @@ static inline void libdivide_u128_shift(uint64_t *u1, uint64_t *u0, int32_t sign *u1 <<= shift; *u1 |= *u0 >> (64 - shift); *u0 <<= shift; - } else { + } + else if (signed_shift < 0) { uint32_t shift = -signed_shift; *u0 >>= shift; *u0 |= *u1 << (64 - shift); @@ -501,7 +448,8 @@ static inline void libdivide_u128_shift(uint64_t *u1, uint64_t *u0, int32_t sign // Computes a 128 / 128 -> 64 bit division, with a 128 bit remainder. 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) { -#if defined(HAS_INT128_T) +#if defined(HAS_INT128_T) && \ + defined(HAS_INT128_DIV) __uint128_t ufull = u_hi; __uint128_t vfull = v_hi; ufull = (ufull << 64) | u_lo; @@ -603,18 +551,15 @@ static inline struct libdivide_u32_t libdivide_internal_u32_gen(uint32_t d, int struct libdivide_u32_t result; uint32_t floor_log_2_d = 31 - libdivide_count_leading_zeros32(d); + + // Power of 2 if ((d & (d - 1)) == 0) { - // Power of 2 - if (! branchfree) { - result.magic = 0; - result.more = (uint8_t)(floor_log_2_d | LIBDIVIDE_U32_SHIFT_PATH); - } else { - // We want a magic number of 2**32 and a shift of floor_log_2_d - // but one of the shifts is taken up by LIBDIVIDE_ADD_MARKER, - // so we subtract 1 from the shift - result.magic = 0; - result.more = (uint8_t)((floor_log_2_d-1) | LIBDIVIDE_ADD_MARKER); - } + // We need to subtract 1 from the shift value in case of an unsigned + // branchfree divider because there is a hardcoded right shift by 1 + // in its division algorithm. Because of this we also need to add back + // 1 in its recovery algorithm. + result.magic = 0; + result.more = (uint8_t)(floor_log_2_d - (branchfree != 0)); } else { uint8_t more; uint32_t rem, proposed_m; @@ -664,8 +609,8 @@ struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t d) { uint32_t libdivide_u32_do(uint32_t numer, const struct libdivide_u32_t *denom) { uint8_t more = denom->more; - if (more & LIBDIVIDE_U32_SHIFT_PATH) { - return numer >> (more & LIBDIVIDE_32_SHIFT_MASK); + if (!denom->magic) { + return numer >> more; } else { uint32_t q = libdivide_mullhi_u32(denom->magic, numer); @@ -681,10 +626,17 @@ uint32_t libdivide_u32_do(uint32_t numer, const struct libdivide_u32_t *denom) { } } +uint32_t libdivide_u32_branchfree_do(uint32_t numer, const struct libdivide_u32_branchfree_t *denom) { + uint32_t q = libdivide_mullhi_u32(denom->magic, numer); + uint32_t t = ((numer - q) >> 1) + q; + return t >> denom->more; +} + uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom) { uint8_t more = denom->more; uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK; - if (more & LIBDIVIDE_U32_SHIFT_PATH) { + + if (!denom->magic) { return 1U << shift; } else if (!(more & LIBDIVIDE_ADD_MARKER)) { // We compute q = n/d = n*m / 2^(32 + shift) @@ -713,50 +665,38 @@ uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom) { // Note that rem<<1 cannot overflow, since rem < d and d is 33 bits uint32_t full_q = half_q + half_q + ((rem<<1) >= d); - // We rounded down in gen unless we're a power of 2 (i.e. in branchfree case) - // We can detect that by looking at m. If m zero, we're a power of 2 - return full_q + (denom->magic != 0); + // We rounded down in gen (hence +1) + return full_q + 1; } } uint32_t libdivide_u32_branchfree_recover(const struct libdivide_u32_branchfree_t *denom) { - struct libdivide_u32_t denom_u32 = {denom->magic, (uint8_t)(denom->more | LIBDIVIDE_ADD_MARKER)}; - return libdivide_u32_recover(&denom_u32); -} - -int libdivide_u32_get_algorithm(const struct libdivide_u32_t *denom) { uint8_t more = denom->more; - if (more & LIBDIVIDE_U32_SHIFT_PATH) - return 0; - else if (!(more & LIBDIVIDE_ADD_MARKER)) - return 1; - else - return 2; -} + uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK; -uint32_t libdivide_u32_do_alg0(uint32_t numer, const struct libdivide_u32_t *denom) { - return numer >> (denom->more & LIBDIVIDE_32_SHIFT_MASK); -} + if (!denom->magic) { + return 1U << (shift + 1); + } else { + // Here we wish to compute d = 2^(32+shift+1)/(m+2^32). + // Notice (m + 2^32) is a 33 bit number. Use 64 bit division for now + // Also note that shift may be as high as 31, so shift + 1 will + // overflow. So we have to compute it as 2^(32+shift)/(m+2^32), and + // then double the quotient and remainder. + uint64_t half_n = 1ULL << (32 + shift); + uint64_t d = (1ULL << 32) | denom->magic; + // Note that the quotient is guaranteed <= 32 bits, but the remainder + // may need 33! + uint32_t half_q = (uint32_t)(half_n / d); + uint64_t rem = half_n % d; + // We computed 2^(32+shift)/(m+2^32) + // Need to double it, and then add 1 to the quotient if doubling th + // remainder would increase the quotient. + // Note that rem<<1 cannot overflow, since rem < d and d is 33 bits + uint32_t full_q = half_q + half_q + ((rem<<1) >= d); -uint32_t libdivide_u32_do_alg1(uint32_t numer, const struct libdivide_u32_t *denom) { - uint32_t q = libdivide_mullhi_u32(denom->magic, numer); - return q >> denom->more; -} - -uint32_t libdivide_u32_do_alg2(uint32_t numer, const struct libdivide_u32_t *denom) { - // denom->add != 0 - uint32_t q = libdivide_mullhi_u32(denom->magic, numer); - uint32_t t = ((numer - q) >> 1) + q; - // Note that this mask is typically free. Only the low bits are - // meaningful to a shift, so compilers can optimize this out. - return t >> (denom->more & LIBDIVIDE_32_SHIFT_MASK); -} - -// same as algo 2 -uint32_t libdivide_u32_branchfree_do(uint32_t numer, const struct libdivide_u32_branchfree_t *denom) { - uint32_t q = libdivide_mullhi_u32(denom->magic, numer); - uint32_t t = ((numer - q) >> 1) + q; - return t >> denom->more; + // We rounded down in gen (hence +1) + return full_q + 1; + } } /////////// UINT64 @@ -768,18 +708,15 @@ static inline struct libdivide_u64_t libdivide_internal_u64_gen(uint64_t d, int struct libdivide_u64_t result; uint32_t floor_log_2_d = 63 - libdivide_count_leading_zeros64(d); + + // Power of 2 if ((d & (d - 1)) == 0) { - // Power of 2 - if (! branchfree) { - result.magic = 0; - result.more = floor_log_2_d | LIBDIVIDE_U64_SHIFT_PATH; - } else { - // We want a magic number of 2**64 and a shift of floor_log_2_d - // but one of the shifts is taken up by LIBDIVIDE_ADD_MARKER, - // so we subtract 1 from the shift. - result.magic = 0; - result.more = (floor_log_2_d-1) | LIBDIVIDE_ADD_MARKER; - } + // We need to subtract 1 from the shift value in case of an unsigned + // branchfree divider because there is a hardcoded right shift by 1 + // in its division algorithm. Because of this we also need to add back + // 1 in its recovery algorithm. + result.magic = 0; + result.more = (uint8_t)(floor_log_2_d - (branchfree != 0)); } else { uint64_t proposed_m, rem; uint8_t more; @@ -831,8 +768,8 @@ struct libdivide_u64_branchfree_t libdivide_u64_branchfree_gen(uint64_t d) { uint64_t libdivide_u64_do(uint64_t numer, const struct libdivide_u64_t *denom) { uint8_t more = denom->more; - if (more & LIBDIVIDE_U64_SHIFT_PATH) { - return numer >> (more & LIBDIVIDE_64_SHIFT_MASK); + if (!denom->magic) { + return numer >> more; } else { uint64_t q = libdivide_mullhi_u64(denom->magic, numer); @@ -848,10 +785,17 @@ uint64_t libdivide_u64_do(uint64_t numer, const struct libdivide_u64_t *denom) { } } +uint64_t libdivide_u64_branchfree_do(uint64_t numer, const struct libdivide_u64_branchfree_t *denom) { + uint64_t q = libdivide_mullhi_u64(denom->magic, numer); + uint64_t t = ((numer - q) >> 1) + q; + return t >> denom->more; +} + uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom) { uint8_t more = denom->more; uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK; - if (more & LIBDIVIDE_U64_SHIFT_PATH) { + + if (!denom->magic) { return 1ULL << shift; } else if (!(more & LIBDIVIDE_ADD_MARKER)) { // We compute q = n/d = n*m / 2^(64 + shift) @@ -868,14 +812,6 @@ uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom) { // libdivide_u32_recover for more on what we do here. // TODO: do something better than 128 bit math - // Hack: if d is not a power of 2, this is a 128/128->64 divide - // If d is a power of 2, this may be a bigger divide - // However we can optimize that easily - if (denom->magic == 0) { - // 2^(64 + shift + 1) / (2^64) == 2^(shift + 1) - return 1ULL << (shift + 1); - } - // Full n is a (potentially) 129 bit value // half_n is a 128 bit value // Compute the hi half of half_n. Low half is 0. @@ -899,40 +835,37 @@ uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom) { } uint64_t libdivide_u64_branchfree_recover(const struct libdivide_u64_branchfree_t *denom) { - struct libdivide_u64_t denom_u64 = {denom->magic, (uint8_t)(denom->more | LIBDIVIDE_ADD_MARKER)}; - return libdivide_u64_recover(&denom_u64); -} - -int libdivide_u64_get_algorithm(const struct libdivide_u64_t *denom) { uint8_t more = denom->more; - if (more & LIBDIVIDE_U64_SHIFT_PATH) - return 0; - else if (!(more & LIBDIVIDE_ADD_MARKER)) - return 1; - else - return 2; -} + uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK; -uint64_t libdivide_u64_do_alg0(uint64_t numer, const struct libdivide_u64_t *denom) { - return numer >> (denom->more & LIBDIVIDE_64_SHIFT_MASK); -} + if (!denom->magic) { + return 1ULL << (shift + 1); + } else { + // Here we wish to compute d = 2^(64+shift+1)/(m+2^64). + // Notice (m + 2^64) is a 65 bit number. This gets hairy. See + // libdivide_u32_recover for more on what we do here. + // TODO: do something better than 128 bit math -uint64_t libdivide_u64_do_alg1(uint64_t numer, const struct libdivide_u64_t *denom) { - uint64_t q = libdivide_mullhi_u64(denom->magic, numer); - return q >> denom->more; -} - -uint64_t libdivide_u64_do_alg2(uint64_t numer, const struct libdivide_u64_t *denom) { - uint64_t q = libdivide_mullhi_u64(denom->magic, numer); - uint64_t t = ((numer - q) >> 1) + q; - return t >> (denom->more & LIBDIVIDE_64_SHIFT_MASK); -} - -// same as alg 2 -uint64_t libdivide_u64_branchfree_do(uint64_t numer, const struct libdivide_u64_branchfree_t *denom) { - uint64_t q = libdivide_mullhi_u64(denom->magic, numer); - uint64_t t = ((numer - q) >> 1) + q; - return t >> denom->more; + // Full n is a (potentially) 129 bit value + // half_n is a 128 bit value + // Compute the hi half of half_n. Low half is 0. + uint64_t half_n_hi = 1ULL << shift, half_n_lo = 0; + // d is a 65 bit value. The high bit is always set to 1. + const uint64_t d_hi = 1, d_lo = denom->magic; + // Note that the quotient is guaranteed <= 64 bits, + // but the remainder may need 65! + uint64_t r_hi, r_lo; + uint64_t half_q = libdivide_128_div_128_to_64(half_n_hi, half_n_lo, d_hi, d_lo, &r_hi, &r_lo); + // We computed 2^(64+shift)/(m+2^64) + // Double the remainder ('dr') and check if that is larger than d + // Note that d is a 65 bit value, so r1 is small and so r1 + r1 + // cannot overflow + uint64_t dr_lo = r_lo + r_lo; + uint64_t dr_hi = r_hi + r_hi + (dr_lo < r_lo); // last term is carry + int dr_exceeds_d = (dr_hi > d_hi) || (dr_hi == d_hi && dr_lo >= d_lo); + uint64_t full_q = half_q + half_q + (dr_exceeds_d ? 1 : 0); + return full_q + 1; + } } /////////// SINT32 @@ -958,7 +891,7 @@ static inline struct libdivide_s32_t libdivide_internal_s32_gen(int32_t d, int b if ((absD & (absD - 1)) == 0) { // Branchfree and normal paths are exactly the same result.magic = 0; - result.more = floor_log_2_d | (d < 0 ? LIBDIVIDE_NEGATIVE_DIVISOR : 0) | LIBDIVIDE_S32_SHIFT_PATH; + result.more = floor_log_2_d | (d < 0 ? LIBDIVIDE_NEGATIVE_DIVISOR : 0); } else { LIBDIVIDE_ASSERT(floor_log_2_d >= 1); @@ -1002,17 +935,11 @@ static inline struct libdivide_s32_t libdivide_internal_s32_gen(int32_t d, int b return result; } -LIBDIVIDE_API struct libdivide_s32_t libdivide_s32_gen(int32_t d) { +struct libdivide_s32_t libdivide_s32_gen(int32_t d) { return libdivide_internal_s32_gen(d, 0); } -LIBDIVIDE_API struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t d) { - if (d == 1) { - LIBDIVIDE_ERROR("branchfree divider must be != 1"); - } - if (d == -1) { - LIBDIVIDE_ERROR("branchfree divider must be != -1"); - } +struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t d) { struct libdivide_s32_t tmp = libdivide_internal_s32_gen(d, 1); struct libdivide_s32_branchfree_t result = {tmp.magic, tmp.more}; return result; @@ -1020,13 +947,14 @@ LIBDIVIDE_API struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int int32_t libdivide_s32_do(int32_t numer, const struct libdivide_s32_t *denom) { uint8_t more = denom->more; - if (more & LIBDIVIDE_S32_SHIFT_PATH) { + uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK; + + if (!denom->magic) { uint32_t sign = (int8_t)more >> 7; - uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK; uint32_t mask = (1U << shift) - 1; uint32_t uq = numer + ((numer >> 31) & mask); int32_t q = (int32_t)uq; - q = q >> shift; + q >>= shift; q = (q ^ sign) - sign; return q; } else { @@ -1034,11 +962,12 @@ int32_t libdivide_s32_do(int32_t numer, const struct libdivide_s32_t *denom) { if (more & LIBDIVIDE_ADD_MARKER) { // must be arithmetic shift and then sign extend int32_t sign = (int8_t)more >> 7; - // q += (more < 0 ? -numer : numer), casts to avoid UB + // q += (more < 0 ? -numer : numer) + // cast required to avoid UB uq += ((uint32_t)numer ^ sign) - sign; } int32_t q = (int32_t)uq; - q >>= more & LIBDIVIDE_32_SHIFT_MASK; + q >>= shift; q += (q < 0); return q; } @@ -1056,13 +985,12 @@ int32_t libdivide_s32_branchfree_do(int32_t numer, const struct libdivide_s32_br // If q is non-negative, we have nothing to do // If q is negative, we want to add either (2**shift)-1 if d is a power of // 2, or (2**shift) if it is not a power of 2 - uint32_t is_power_of_2 = !!(more & LIBDIVIDE_S32_SHIFT_PATH); + uint32_t is_power_of_2 = (magic == 0); uint32_t q_sign = (uint32_t)(q >> 31); - q += q_sign & ((1 << shift) - is_power_of_2); + q += q_sign & ((1U << shift) - is_power_of_2); // Now arithmetic right shift q >>= shift; - // Negate if needed q = (q ^ sign) - sign; @@ -1072,7 +1000,7 @@ int32_t libdivide_s32_branchfree_do(int32_t numer, const struct libdivide_s32_br int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom) { uint8_t more = denom->more; uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK; - if (more & LIBDIVIDE_S32_SHIFT_PATH) { + if (!denom->magic) { uint32_t absD = 1U << shift; if (more & LIBDIVIDE_NEGATIVE_DIVISOR) { absD = -absD; @@ -1092,7 +1020,7 @@ int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom) { // Handle the power of 2 case (including branchfree) if (denom->magic == 0) { - int32_t result = 1 << shift; + int32_t result = 1U << shift; return negative_divisor ? -result : result; } @@ -1109,54 +1037,6 @@ int32_t libdivide_s32_branchfree_recover(const struct libdivide_s32_branchfree_t return libdivide_s32_recover((const struct libdivide_s32_t *)denom); } -int libdivide_s32_get_algorithm(const struct libdivide_s32_t *denom) { - uint8_t more = denom->more; - int positiveDivisor = !(more & LIBDIVIDE_NEGATIVE_DIVISOR); - if (more & LIBDIVIDE_S32_SHIFT_PATH) - return (positiveDivisor ? 0 : 1); - else if (more & LIBDIVIDE_ADD_MARKER) - return (positiveDivisor ? 2 : 3); - else - return 4; -} - -int32_t libdivide_s32_do_alg0(int32_t numer, const struct libdivide_s32_t *denom) { - uint8_t shift = denom->more & LIBDIVIDE_32_SHIFT_MASK; - uint32_t mask = (1U << shift) - 1; - int32_t q = numer + ((numer >> 31) & mask); - return q >> shift; -} - -int32_t libdivide_s32_do_alg1(int32_t numer, const struct libdivide_s32_t *denom) { - uint8_t shift = denom->more & LIBDIVIDE_32_SHIFT_MASK; - uint32_t mask = (1U << shift) - 1; - int32_t q = numer + ((numer >> 31) & mask); - return - (q >> shift); -} - -int32_t libdivide_s32_do_alg2(int32_t numer, const struct libdivide_s32_t *denom) { - int32_t q = libdivide_mullhi_s32(denom->magic, numer); - q += numer; - q >>= denom->more & LIBDIVIDE_32_SHIFT_MASK; - q += (q < 0); - return q; -} - -int32_t libdivide_s32_do_alg3(int32_t numer, const struct libdivide_s32_t *denom) { - int32_t q = libdivide_mullhi_s32(denom->magic, numer); - q -= numer; - q >>= denom->more & LIBDIVIDE_32_SHIFT_MASK; - q += (q < 0); - return q; -} - -int32_t libdivide_s32_do_alg4(int32_t numer, const struct libdivide_s32_t *denom) { - int32_t q = libdivide_mullhi_s32(denom->magic, numer); - q >>= denom->more & LIBDIVIDE_32_SHIFT_MASK; - q += (q < 0); - return q; -} - ///////////// SINT64 static inline struct libdivide_s64_t libdivide_internal_s64_gen(int64_t d, int branchfree) { @@ -1229,12 +1109,6 @@ struct libdivide_s64_t libdivide_s64_gen(int64_t d) { } struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t d) { - if (d == 1) { - LIBDIVIDE_ERROR("branchfree divider must be != 1"); - } - if (d == -1) { - LIBDIVIDE_ERROR("branchfree divider must be != -1"); - } struct libdivide_s64_t tmp = libdivide_internal_s64_gen(d, 1); struct libdivide_s64_branchfree_t ret = {tmp.magic, tmp.more}; return ret; @@ -1242,26 +1116,28 @@ struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t d) { int64_t libdivide_s64_do(int64_t numer, const struct libdivide_s64_t *denom) { uint8_t more = denom->more; - int64_t magic = denom->magic; - if (magic == 0) { //shift path - uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + + if (!denom->magic) { // shift path uint64_t mask = (1ULL << shift) - 1; uint64_t uq = numer + ((numer >> 63) & mask); int64_t q = (int64_t)uq; - q = q >> shift; + q >>= shift; // must be arithmetic shift and then sign-extend - int64_t shiftMask = (int8_t)more >> 7; - q = (q ^ shiftMask) - shiftMask; + int64_t sign = (int8_t)more >> 7; + q = (q ^ sign) - sign; return q; } else { - uint64_t uq = (uint64_t)libdivide_mullhi_s64(magic, numer); + uint64_t uq = (uint64_t)libdivide_mullhi_s64(denom->magic, numer); if (more & LIBDIVIDE_ADD_MARKER) { // must be arithmetic shift and then sign extend int64_t sign = (int8_t)more >> 7; + // q += (more < 0 ? -numer : numer) + // cast required to avoid UB uq += ((uint64_t)numer ^ sign) - sign; } int64_t q = (int64_t)uq; - q >>= more & LIBDIVIDE_64_SHIFT_MASK; + q >>= shift; q += (q < 0); return q; } @@ -1269,7 +1145,7 @@ int64_t libdivide_s64_do(int64_t numer, const struct libdivide_s64_t *denom) { int64_t libdivide_s64_branchfree_do(int64_t numer, const struct libdivide_s64_branchfree_t *denom) { uint8_t more = denom->more; - uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK; // must be arithmetic shift and then sign extend int64_t sign = (int8_t)more >> 7; int64_t magic = denom->magic; @@ -1279,15 +1155,15 @@ int64_t libdivide_s64_branchfree_do(int64_t numer, const struct libdivide_s64_br // If q is non-negative, we have nothing to do. // If q is negative, we want to add either (2**shift)-1 if d is a power of // 2, or (2**shift) if it is not a power of 2. - uint32_t is_power_of_2 = (magic == 0); + uint64_t is_power_of_2 = (magic == 0); uint64_t q_sign = (uint64_t)(q >> 63); q += q_sign & ((1ULL << shift) - is_power_of_2); // Arithmetic right shift q >>= shift; - // Negate if needed q = (q ^ sign) - sign; + return q; } @@ -1322,203 +1198,528 @@ int64_t libdivide_s64_branchfree_recover(const struct libdivide_s64_branchfree_t return libdivide_s64_recover((const struct libdivide_s64_t *)denom); } -int libdivide_s64_get_algorithm(const struct libdivide_s64_t *denom) { - uint8_t more = denom->more; - int positiveDivisor = !(more & LIBDIVIDE_NEGATIVE_DIVISOR); - if (denom->magic == 0) - return (positiveDivisor ? 0 : 1); // shift path - else if (more & LIBDIVIDE_ADD_MARKER) - return (positiveDivisor ? 2 : 3); - else - return 4; -} +#if defined(LIBDIVIDE_AVX512) -int64_t libdivide_s64_do_alg0(int64_t numer, const struct libdivide_s64_t *denom) { - uint32_t shift = denom->more & LIBDIVIDE_64_SHIFT_MASK; - uint64_t mask = (1ULL << shift) - 1; - int64_t q = numer + ((numer >> 63) & mask); - return q >> shift; -} +static inline __m512i libdivide_u32_do_vector(__m512i numers, const struct libdivide_u32_t *denom); +static inline __m512i libdivide_s32_do_vector(__m512i numers, const struct libdivide_s32_t *denom); +static inline __m512i libdivide_u64_do_vector(__m512i numers, const struct libdivide_u64_t *denom); +static inline __m512i libdivide_s64_do_vector(__m512i numers, const struct libdivide_s64_t *denom); -int64_t libdivide_s64_do_alg1(int64_t numer, const struct libdivide_s64_t *denom) { - // denom->shift != -1 && demo->shiftMask != 0 - uint32_t shift = denom->more & LIBDIVIDE_64_SHIFT_MASK; - uint64_t mask = (1ULL << shift) - 1; - int64_t q = numer + ((numer >> 63) & mask); - return - (q >> shift); -} - -int64_t libdivide_s64_do_alg2(int64_t numer, const struct libdivide_s64_t *denom) { - int64_t q = libdivide_mullhi_s64(denom->magic, numer); - q += numer; - q >>= denom->more & LIBDIVIDE_64_SHIFT_MASK; - q += (q < 0); - return q; -} - -int64_t libdivide_s64_do_alg3(int64_t numer, const struct libdivide_s64_t *denom) { - int64_t q = libdivide_mullhi_s64(denom->magic, numer); - q -= numer; - q >>= denom->more & LIBDIVIDE_64_SHIFT_MASK; - q += (q < 0); - return q; -} - -int64_t libdivide_s64_do_alg4(int64_t numer, const struct libdivide_s64_t *denom) { - int64_t q = libdivide_mullhi_s64(denom->magic, numer); - q >>= denom->more & LIBDIVIDE_64_SHIFT_MASK; - q += (q < 0); - return q; -} - -#if defined(LIBDIVIDE_USE_SSE2) - -LIBDIVIDE_API __m128i libdivide_u32_do_vector(__m128i numers, const struct libdivide_u32_t *denom); -LIBDIVIDE_API __m128i libdivide_s32_do_vector(__m128i numers, const struct libdivide_s32_t *denom); -LIBDIVIDE_API __m128i libdivide_u64_do_vector(__m128i numers, const struct libdivide_u64_t *denom); -LIBDIVIDE_API __m128i libdivide_s64_do_vector(__m128i numers, const struct libdivide_s64_t *denom); - -LIBDIVIDE_API __m128i libdivide_u32_do_vector_alg0(__m128i numers, const struct libdivide_u32_t *denom); -LIBDIVIDE_API __m128i libdivide_u32_do_vector_alg1(__m128i numers, const struct libdivide_u32_t *denom); -LIBDIVIDE_API __m128i libdivide_u32_do_vector_alg2(__m128i numers, const struct libdivide_u32_t *denom); - -LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg0(__m128i numers, const struct libdivide_s32_t *denom); -LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg1(__m128i numers, const struct libdivide_s32_t *denom); -LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg2(__m128i numers, const struct libdivide_s32_t *denom); -LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg3(__m128i numers, const struct libdivide_s32_t *denom); -LIBDIVIDE_API __m128i libdivide_s32_do_vector_alg4(__m128i numers, const struct libdivide_s32_t *denom); - -LIBDIVIDE_API __m128i libdivide_u64_do_vector_alg0(__m128i numers, const struct libdivide_u64_t *denom); -LIBDIVIDE_API __m128i libdivide_u64_do_vector_alg1(__m128i numers, const struct libdivide_u64_t *denom); -LIBDIVIDE_API __m128i libdivide_u64_do_vector_alg2(__m128i numers, const struct libdivide_u64_t *denom); - -LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg0(__m128i numers, const struct libdivide_s64_t *denom); -LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg1(__m128i numers, const struct libdivide_s64_t *denom); -LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg2(__m128i numers, const struct libdivide_s64_t *denom); -LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg3(__m128i numers, const struct libdivide_s64_t *denom); -LIBDIVIDE_API __m128i libdivide_s64_do_vector_alg4(__m128i numers, const struct libdivide_s64_t *denom); - -LIBDIVIDE_API __m128i libdivide_u32_branchfree_do_vector(__m128i numers, const struct libdivide_u32_branchfree_t *denom); -LIBDIVIDE_API __m128i libdivide_s32_branchfree_do_vector(__m128i numers, const struct libdivide_s32_branchfree_t *denom); -LIBDIVIDE_API __m128i libdivide_u64_branchfree_do_vector(__m128i numers, const struct libdivide_u64_branchfree_t *denom); -LIBDIVIDE_API __m128i libdivide_s64_branchfree_do_vector(__m128i numers, const struct libdivide_s64_branchfree_t *denom); +static inline __m512i libdivide_u32_branchfree_do_vector(__m512i numers, const struct libdivide_u32_branchfree_t *denom); +static inline __m512i libdivide_s32_branchfree_do_vector(__m512i numers, const struct libdivide_s32_branchfree_t *denom); +static inline __m512i libdivide_u64_branchfree_do_vector(__m512i numers, const struct libdivide_u64_branchfree_t *denom); +static inline __m512i libdivide_s64_branchfree_do_vector(__m512i numers, const struct libdivide_s64_branchfree_t *denom); //////// Internal Utility Functions -static inline __m128i libdivide_u64_to_m128(uint64_t x) { - return _mm_set1_epi64x(x); +static inline __m512i libdivide_s64_signbits(__m512i v) {; + return _mm512_srai_epi64(v, 63); } -static inline __m128i libdivide_get_FFFFFFFF00000000(void) { - // returns the same as _mm_set1_epi64(0xFFFFFFFF00000000ULL) - // without touching memory. - // optimizes to pcmpeqd on OS X - __m128i result = _mm_set1_epi8(-1); - return _mm_slli_epi64(result, 32); +static inline __m512i libdivide_s64_shift_right_vector(__m512i v, int amt) { + return _mm512_srai_epi64(v, amt); } -static inline __m128i libdivide_get_00000000FFFFFFFF(void) { - // returns the same as _mm_set1_epi64(0x00000000FFFFFFFFULL) - // without touching memory. - // optimizes to pcmpeqd on OS X - __m128i result = _mm_set1_epi8(-1); - result = _mm_srli_epi64(result, 32); +// Here, b is assumed to contain one 32-bit value repeated. +static inline __m512i libdivide_mullhi_u32_vector(__m512i a, __m512i b) { + __m512i hi_product_0Z2Z = _mm512_srli_epi64(_mm512_mul_epu32(a, b), 32); + __m512i a1X3X = _mm512_srli_epi64(a, 32); + __m512i mask = _mm512_set_epi32(-1, 0, -1, 0, -1, 0, -1, 0, -1, 0, -1, 0, -1, 0, -1, 0); + __m512i hi_product_Z1Z3 = _mm512_and_si512(_mm512_mul_epu32(a1X3X, b), mask); + return _mm512_or_si512(hi_product_0Z2Z, hi_product_Z1Z3); +} + +// b is one 32-bit value repeated. +static inline __m512i libdivide_mullhi_s32_vector(__m512i a, __m512i b) { + __m512i hi_product_0Z2Z = _mm512_srli_epi64(_mm512_mul_epi32(a, b), 32); + __m512i a1X3X = _mm512_srli_epi64(a, 32); + __m512i mask = _mm512_set_epi32(-1, 0, -1, 0, -1, 0, -1, 0, -1, 0, -1, 0, -1, 0, -1, 0); + __m512i hi_product_Z1Z3 = _mm512_and_si512(_mm512_mul_epi32(a1X3X, b), mask); + return _mm512_or_si512(hi_product_0Z2Z, hi_product_Z1Z3); +} + +// Here, y is assumed to contain one 64-bit value repeated. +// https://stackoverflow.com/a/28827013 +static inline __m512i libdivide_mullhi_u64_vector(__m512i x, __m512i y) { + __m512i lomask = _mm512_set1_epi64(0xffffffff); + __m512i xh = _mm512_shuffle_epi32(x, (_MM_PERM_ENUM) 0xB1); + __m512i yh = _mm512_shuffle_epi32(y, (_MM_PERM_ENUM) 0xB1); + __m512i w0 = _mm512_mul_epu32(x, y); + __m512i w1 = _mm512_mul_epu32(x, yh); + __m512i w2 = _mm512_mul_epu32(xh, y); + __m512i w3 = _mm512_mul_epu32(xh, yh); + __m512i w0h = _mm512_srli_epi64(w0, 32); + __m512i s1 = _mm512_add_epi64(w1, w0h); + __m512i s1l = _mm512_and_si512(s1, lomask); + __m512i s1h = _mm512_srli_epi64(s1, 32); + __m512i s2 = _mm512_add_epi64(w2, s1l); + __m512i s2h = _mm512_srli_epi64(s2, 32); + __m512i hi = _mm512_add_epi64(w3, s1h); + hi = _mm512_add_epi64(hi, s2h); + + return hi; +} + +// y is one 64-bit value repeated. +static inline __m512i libdivide_mullhi_s64_vector(__m512i x, __m512i y) { + __m512i p = libdivide_mullhi_u64_vector(x, y); + __m512i t1 = _mm512_and_si512(libdivide_s64_signbits(x), y); + __m512i t2 = _mm512_and_si512(libdivide_s64_signbits(y), x); + p = _mm512_sub_epi64(p, t1); + p = _mm512_sub_epi64(p, t2); + return p; +} + +////////// UINT32 + +__m512i libdivide_u32_do_vector(__m512i numers, const struct libdivide_u32_t *denom) { + uint8_t more = denom->more; + if (!denom->magic) { + return _mm512_srli_epi32(numers, more); + } + else { + __m512i q = libdivide_mullhi_u32_vector(numers, _mm512_set1_epi32(denom->magic)); + if (more & LIBDIVIDE_ADD_MARKER) { + // uint32_t t = ((numer - q) >> 1) + q; + // return t >> denom->shift; + uint32_t shift = more & LIBDIVIDE_32_SHIFT_MASK; + __m512i t = _mm512_add_epi32(_mm512_srli_epi32(_mm512_sub_epi32(numers, q), 1), q); + return _mm512_srli_epi32(t, shift); + } + else { + return _mm512_srli_epi32(q, more); + } + } +} + +__m512i libdivide_u32_branchfree_do_vector(__m512i numers, const struct libdivide_u32_branchfree_t *denom) { + __m512i q = libdivide_mullhi_u32_vector(numers, _mm512_set1_epi32(denom->magic)); + __m512i t = _mm512_add_epi32(_mm512_srli_epi32(_mm512_sub_epi32(numers, q), 1), q); + return _mm512_srli_epi32(t, denom->more); +} + +////////// UINT64 + +__m512i libdivide_u64_do_vector(__m512i numers, const struct libdivide_u64_t *denom) { + uint8_t more = denom->more; + if (!denom->magic) { + return _mm512_srli_epi64(numers, more); + } + else { + __m512i q = libdivide_mullhi_u64_vector(numers, _mm512_set1_epi64(denom->magic)); + if (more & LIBDIVIDE_ADD_MARKER) { + // uint32_t t = ((numer - q) >> 1) + q; + // return t >> denom->shift; + uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + __m512i t = _mm512_add_epi64(_mm512_srli_epi64(_mm512_sub_epi64(numers, q), 1), q); + return _mm512_srli_epi64(t, shift); + } + else { + return _mm512_srli_epi64(q, more); + } + } +} + +__m512i libdivide_u64_branchfree_do_vector(__m512i numers, const struct libdivide_u64_branchfree_t *denom) { + __m512i q = libdivide_mullhi_u64_vector(numers, _mm512_set1_epi64(denom->magic)); + __m512i t = _mm512_add_epi64(_mm512_srli_epi64(_mm512_sub_epi64(numers, q), 1), q); + return _mm512_srli_epi64(t, denom->more); +} + +////////// SINT32 + +__m512i libdivide_s32_do_vector(__m512i numers, const struct libdivide_s32_t *denom) { + uint8_t more = denom->more; + if (!denom->magic) { + uint32_t shift = more & LIBDIVIDE_32_SHIFT_MASK; + uint32_t mask = (1U << shift) - 1; + __m512i roundToZeroTweak = _mm512_set1_epi32(mask); + // q = numer + ((numer >> 31) & roundToZeroTweak); + __m512i q = _mm512_add_epi32(numers, _mm512_and_si512(_mm512_srai_epi32(numers, 31), roundToZeroTweak)); + q = _mm512_srai_epi32(q, shift); + __m512i sign = _mm512_set1_epi32((int8_t)more >> 7); + // q = (q ^ sign) - sign; + q = _mm512_sub_epi32(_mm512_xor_si512(q, sign), sign); + return q; + } + else { + __m512i q = libdivide_mullhi_s32_vector(numers, _mm512_set1_epi32(denom->magic)); + if (more & LIBDIVIDE_ADD_MARKER) { + // must be arithmetic shift + __m512i sign = _mm512_set1_epi32((int8_t)more >> 7); + // q += ((numer ^ sign) - sign); + q = _mm512_add_epi32(q, _mm512_sub_epi32(_mm512_xor_si512(numers, sign), sign)); + } + // q >>= shift + q = _mm512_srai_epi32(q, more & LIBDIVIDE_32_SHIFT_MASK); + q = _mm512_add_epi32(q, _mm512_srli_epi32(q, 31)); // q += (q < 0) + return q; + } +} + +__m512i libdivide_s32_branchfree_do_vector(__m512i numers, const struct libdivide_s32_branchfree_t *denom) { + int32_t magic = denom->magic; + uint8_t more = denom->more; + uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK; + // must be arithmetic shift + __m512i sign = _mm512_set1_epi32((int8_t)more >> 7); + __m512i q = libdivide_mullhi_s32_vector(numers, _mm512_set1_epi32(magic)); + q = _mm512_add_epi32(q, numers); // q += numers + + // If q is non-negative, we have nothing to do + // If q is negative, we want to add either (2**shift)-1 if d is + // a power of 2, or (2**shift) if it is not a power of 2 + uint32_t is_power_of_2 = (magic == 0); + __m512i q_sign = _mm512_srai_epi32(q, 31); // q_sign = q >> 31 + __m512i mask = _mm512_set1_epi32((1U << shift) - is_power_of_2); + q = _mm512_add_epi32(q, _mm512_and_si512(q_sign, mask)); // q = q + (q_sign & mask) + q = _mm512_srai_epi32(q, shift); // q >>= shift + q = _mm512_sub_epi32(_mm512_xor_si512(q, sign), sign); // q = (q ^ sign) - sign + return q; +} + +////////// SINT64 + +__m512i libdivide_s64_do_vector(__m512i numers, const struct libdivide_s64_t *denom) { + uint8_t more = denom->more; + int64_t magic = denom->magic; + if (magic == 0) { // shift path + uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + uint64_t mask = (1ULL << shift) - 1; + __m512i roundToZeroTweak = _mm512_set1_epi64(mask); + // q = numer + ((numer >> 63) & roundToZeroTweak); + __m512i q = _mm512_add_epi64(numers, _mm512_and_si512(libdivide_s64_signbits(numers), roundToZeroTweak)); + q = libdivide_s64_shift_right_vector(q, shift); + __m512i sign = _mm512_set1_epi32((int8_t)more >> 7); + // q = (q ^ sign) - sign; + q = _mm512_sub_epi64(_mm512_xor_si512(q, sign), sign); + return q; + } + else { + __m512i q = libdivide_mullhi_s64_vector(numers, _mm512_set1_epi64(magic)); + if (more & LIBDIVIDE_ADD_MARKER) { + // must be arithmetic shift + __m512i sign = _mm512_set1_epi32((int8_t)more >> 7); + // q += ((numer ^ sign) - sign); + q = _mm512_add_epi64(q, _mm512_sub_epi64(_mm512_xor_si512(numers, sign), sign)); + } + // q >>= denom->mult_path.shift + q = libdivide_s64_shift_right_vector(q, more & LIBDIVIDE_64_SHIFT_MASK); + q = _mm512_add_epi64(q, _mm512_srli_epi64(q, 63)); // q += (q < 0) + return q; + } +} + +__m512i libdivide_s64_branchfree_do_vector(__m512i numers, const struct libdivide_s64_branchfree_t *denom) { + int64_t magic = denom->magic; + uint8_t more = denom->more; + uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + // must be arithmetic shift + __m512i sign = _mm512_set1_epi32((int8_t)more >> 7); + + // libdivide_mullhi_s64(numers, magic); + __m512i q = libdivide_mullhi_s64_vector(numers, _mm512_set1_epi64(magic)); + q = _mm512_add_epi64(q, numers); // q += numers + + // If q is non-negative, we have nothing to do. + // If q is negative, we want to add either (2**shift)-1 if d is + // a power of 2, or (2**shift) if it is not a power of 2. + uint32_t is_power_of_2 = (magic == 0); + __m512i q_sign = libdivide_s64_signbits(q); // q_sign = q >> 63 + __m512i mask = _mm512_set1_epi64((1ULL << shift) - is_power_of_2); + q = _mm512_add_epi64(q, _mm512_and_si512(q_sign, mask)); // q = q + (q_sign & mask) + q = libdivide_s64_shift_right_vector(q, shift); // q >>= shift + q = _mm512_sub_epi64(_mm512_xor_si512(q, sign), sign); // q = (q ^ sign) - sign + return q; +} + +#elif defined(LIBDIVIDE_AVX2) + +static inline __m256i libdivide_u32_do_vector(__m256i numers, const struct libdivide_u32_t *denom); +static inline __m256i libdivide_s32_do_vector(__m256i numers, const struct libdivide_s32_t *denom); +static inline __m256i libdivide_u64_do_vector(__m256i numers, const struct libdivide_u64_t *denom); +static inline __m256i libdivide_s64_do_vector(__m256i numers, const struct libdivide_s64_t *denom); + +static inline __m256i libdivide_u32_branchfree_do_vector(__m256i numers, const struct libdivide_u32_branchfree_t *denom); +static inline __m256i libdivide_s32_branchfree_do_vector(__m256i numers, const struct libdivide_s32_branchfree_t *denom); +static inline __m256i libdivide_u64_branchfree_do_vector(__m256i numers, const struct libdivide_u64_branchfree_t *denom); +static inline __m256i libdivide_s64_branchfree_do_vector(__m256i numers, const struct libdivide_s64_branchfree_t *denom); + +//////// Internal Utility Functions + +// Implementation of _mm256_srai_epi64(v, 63) (from AVX512). +static inline __m256i libdivide_s64_signbits(__m256i v) { + __m256i hiBitsDuped = _mm256_shuffle_epi32(v, _MM_SHUFFLE(3, 3, 1, 1)); + __m256i signBits = _mm256_srai_epi32(hiBitsDuped, 31); + return signBits; +} + +// Implementation of _mm256_srai_epi64 (from AVX512). +static inline __m256i libdivide_s64_shift_right_vector(__m256i v, int amt) { + const int b = 64 - amt; + __m256i m = _mm256_set1_epi64x(1ULL << (b - 1)); + __m256i x = _mm256_srli_epi64(v, amt); + __m256i result = _mm256_sub_epi64(_mm256_xor_si256(x, m), m); return result; } +// Here, b is assumed to contain one 32-bit value repeated. +static inline __m256i libdivide_mullhi_u32_vector(__m256i a, __m256i b) { + __m256i hi_product_0Z2Z = _mm256_srli_epi64(_mm256_mul_epu32(a, b), 32); + __m256i a1X3X = _mm256_srli_epi64(a, 32); + __m256i mask = _mm256_set_epi32(-1, 0, -1, 0, -1, 0, -1, 0); + __m256i hi_product_Z1Z3 = _mm256_and_si256(_mm256_mul_epu32(a1X3X, b), mask); + return _mm256_or_si256(hi_product_0Z2Z, hi_product_Z1Z3); +} + +// b is one 32-bit value repeated. +static inline __m256i libdivide_mullhi_s32_vector(__m256i a, __m256i b) { + __m256i hi_product_0Z2Z = _mm256_srli_epi64(_mm256_mul_epi32(a, b), 32); + __m256i a1X3X = _mm256_srli_epi64(a, 32); + __m256i mask = _mm256_set_epi32(-1, 0, -1, 0, -1, 0, -1, 0); + __m256i hi_product_Z1Z3 = _mm256_and_si256(_mm256_mul_epi32(a1X3X, b), mask); + return _mm256_or_si256(hi_product_0Z2Z, hi_product_Z1Z3); +} + +// Here, y is assumed to contain one 64-bit value repeated. +// https://stackoverflow.com/a/28827013 +static inline __m256i libdivide_mullhi_u64_vector(__m256i x, __m256i y) { + __m256i lomask = _mm256_set1_epi64x(0xffffffff); + __m256i xh = _mm256_shuffle_epi32(x, 0xB1); // x0l, x0h, x1l, x1h + __m256i yh = _mm256_shuffle_epi32(y, 0xB1); // y0l, y0h, y1l, y1h + __m256i w0 = _mm256_mul_epu32(x, y); // x0l*y0l, x1l*y1l + __m256i w1 = _mm256_mul_epu32(x, yh); // x0l*y0h, x1l*y1h + __m256i w2 = _mm256_mul_epu32(xh, y); // x0h*y0l, x1h*y0l + __m256i w3 = _mm256_mul_epu32(xh, yh); // x0h*y0h, x1h*y1h + __m256i w0h = _mm256_srli_epi64(w0, 32); + __m256i s1 = _mm256_add_epi64(w1, w0h); + __m256i s1l = _mm256_and_si256(s1, lomask); + __m256i s1h = _mm256_srli_epi64(s1, 32); + __m256i s2 = _mm256_add_epi64(w2, s1l); + __m256i s2h = _mm256_srli_epi64(s2, 32); + __m256i hi = _mm256_add_epi64(w3, s1h); + hi = _mm256_add_epi64(hi, s2h); + + return hi; +} + +// y is one 64-bit value repeated. +static inline __m256i libdivide_mullhi_s64_vector(__m256i x, __m256i y) { + __m256i p = libdivide_mullhi_u64_vector(x, y); + __m256i t1 = _mm256_and_si256(libdivide_s64_signbits(x), y); + __m256i t2 = _mm256_and_si256(libdivide_s64_signbits(y), x); + p = _mm256_sub_epi64(p, t1); + p = _mm256_sub_epi64(p, t2); + return p; +} + +////////// UINT32 + +__m256i libdivide_u32_do_vector(__m256i numers, const struct libdivide_u32_t *denom) { + uint8_t more = denom->more; + if (!denom->magic) { + return _mm256_srli_epi32(numers, more); + } + else { + __m256i q = libdivide_mullhi_u32_vector(numers, _mm256_set1_epi32(denom->magic)); + if (more & LIBDIVIDE_ADD_MARKER) { + // uint32_t t = ((numer - q) >> 1) + q; + // return t >> denom->shift; + uint32_t shift = more & LIBDIVIDE_32_SHIFT_MASK; + __m256i t = _mm256_add_epi32(_mm256_srli_epi32(_mm256_sub_epi32(numers, q), 1), q); + return _mm256_srli_epi32(t, shift); + } + else { + return _mm256_srli_epi32(q, more); + } + } +} + +__m256i libdivide_u32_branchfree_do_vector(__m256i numers, const struct libdivide_u32_branchfree_t *denom) { + __m256i q = libdivide_mullhi_u32_vector(numers, _mm256_set1_epi32(denom->magic)); + __m256i t = _mm256_add_epi32(_mm256_srli_epi32(_mm256_sub_epi32(numers, q), 1), q); + return _mm256_srli_epi32(t, denom->more); +} + +////////// UINT64 + +__m256i libdivide_u64_do_vector(__m256i numers, const struct libdivide_u64_t *denom) { + uint8_t more = denom->more; + if (!denom->magic) { + return _mm256_srli_epi64(numers, more); + } + else { + __m256i q = libdivide_mullhi_u64_vector(numers, _mm256_set1_epi64x(denom->magic)); + if (more & LIBDIVIDE_ADD_MARKER) { + // uint32_t t = ((numer - q) >> 1) + q; + // return t >> denom->shift; + uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + __m256i t = _mm256_add_epi64(_mm256_srli_epi64(_mm256_sub_epi64(numers, q), 1), q); + return _mm256_srli_epi64(t, shift); + } + else { + return _mm256_srli_epi64(q, more); + } + } +} + +__m256i libdivide_u64_branchfree_do_vector(__m256i numers, const struct libdivide_u64_branchfree_t *denom) { + __m256i q = libdivide_mullhi_u64_vector(numers, _mm256_set1_epi64x(denom->magic)); + __m256i t = _mm256_add_epi64(_mm256_srli_epi64(_mm256_sub_epi64(numers, q), 1), q); + return _mm256_srli_epi64(t, denom->more); +} + +////////// SINT32 + +__m256i libdivide_s32_do_vector(__m256i numers, const struct libdivide_s32_t *denom) { + uint8_t more = denom->more; + if (!denom->magic) { + uint32_t shift = more & LIBDIVIDE_32_SHIFT_MASK; + uint32_t mask = (1U << shift) - 1; + __m256i roundToZeroTweak = _mm256_set1_epi32(mask); + // q = numer + ((numer >> 31) & roundToZeroTweak); + __m256i q = _mm256_add_epi32(numers, _mm256_and_si256(_mm256_srai_epi32(numers, 31), roundToZeroTweak)); + q = _mm256_srai_epi32(q, shift); + __m256i sign = _mm256_set1_epi32((int8_t)more >> 7); + // q = (q ^ sign) - sign; + q = _mm256_sub_epi32(_mm256_xor_si256(q, sign), sign); + return q; + } + else { + __m256i q = libdivide_mullhi_s32_vector(numers, _mm256_set1_epi32(denom->magic)); + if (more & LIBDIVIDE_ADD_MARKER) { + // must be arithmetic shift + __m256i sign = _mm256_set1_epi32((int8_t)more >> 7); + // q += ((numer ^ sign) - sign); + q = _mm256_add_epi32(q, _mm256_sub_epi32(_mm256_xor_si256(numers, sign), sign)); + } + // q >>= shift + q = _mm256_srai_epi32(q, more & LIBDIVIDE_32_SHIFT_MASK); + q = _mm256_add_epi32(q, _mm256_srli_epi32(q, 31)); // q += (q < 0) + return q; + } +} + +__m256i libdivide_s32_branchfree_do_vector(__m256i numers, const struct libdivide_s32_branchfree_t *denom) { + int32_t magic = denom->magic; + uint8_t more = denom->more; + uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK; + // must be arithmetic shift + __m256i sign = _mm256_set1_epi32((int8_t)more >> 7); + __m256i q = libdivide_mullhi_s32_vector(numers, _mm256_set1_epi32(magic)); + q = _mm256_add_epi32(q, numers); // q += numers + + // If q is non-negative, we have nothing to do + // If q is negative, we want to add either (2**shift)-1 if d is + // a power of 2, or (2**shift) if it is not a power of 2 + uint32_t is_power_of_2 = (magic == 0); + __m256i q_sign = _mm256_srai_epi32(q, 31); // q_sign = q >> 31 + __m256i mask = _mm256_set1_epi32((1U << shift) - is_power_of_2); + q = _mm256_add_epi32(q, _mm256_and_si256(q_sign, mask)); // q = q + (q_sign & mask) + q = _mm256_srai_epi32(q, shift); // q >>= shift + q = _mm256_sub_epi32(_mm256_xor_si256(q, sign), sign); // q = (q ^ sign) - sign + return q; +} + +////////// SINT64 + +__m256i libdivide_s64_do_vector(__m256i numers, const struct libdivide_s64_t *denom) { + uint8_t more = denom->more; + int64_t magic = denom->magic; + if (magic == 0) { // shift path + uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + uint64_t mask = (1ULL << shift) - 1; + __m256i roundToZeroTweak = _mm256_set1_epi64x(mask); + // q = numer + ((numer >> 63) & roundToZeroTweak); + __m256i q = _mm256_add_epi64(numers, _mm256_and_si256(libdivide_s64_signbits(numers), roundToZeroTweak)); + q = libdivide_s64_shift_right_vector(q, shift); + __m256i sign = _mm256_set1_epi32((int8_t)more >> 7); + // q = (q ^ sign) - sign; + q = _mm256_sub_epi64(_mm256_xor_si256(q, sign), sign); + return q; + } + else { + __m256i q = libdivide_mullhi_s64_vector(numers, _mm256_set1_epi64x(magic)); + if (more & LIBDIVIDE_ADD_MARKER) { + // must be arithmetic shift + __m256i sign = _mm256_set1_epi32((int8_t)more >> 7); + // q += ((numer ^ sign) - sign); + q = _mm256_add_epi64(q, _mm256_sub_epi64(_mm256_xor_si256(numers, sign), sign)); + } + // q >>= denom->mult_path.shift + q = libdivide_s64_shift_right_vector(q, more & LIBDIVIDE_64_SHIFT_MASK); + q = _mm256_add_epi64(q, _mm256_srli_epi64(q, 63)); // q += (q < 0) + return q; + } +} + +__m256i libdivide_s64_branchfree_do_vector(__m256i numers, const struct libdivide_s64_branchfree_t *denom) { + int64_t magic = denom->magic; + uint8_t more = denom->more; + uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK; + // must be arithmetic shift + __m256i sign = _mm256_set1_epi32((int8_t)more >> 7); + + // libdivide_mullhi_s64(numers, magic); + __m256i q = libdivide_mullhi_s64_vector(numers, _mm256_set1_epi64x(magic)); + q = _mm256_add_epi64(q, numers); // q += numers + + // If q is non-negative, we have nothing to do. + // If q is negative, we want to add either (2**shift)-1 if d is + // a power of 2, or (2**shift) if it is not a power of 2. + uint32_t is_power_of_2 = (magic == 0); + __m256i q_sign = libdivide_s64_signbits(q); // q_sign = q >> 63 + __m256i mask = _mm256_set1_epi64x((1ULL << shift) - is_power_of_2); + q = _mm256_add_epi64(q, _mm256_and_si256(q_sign, mask)); // q = q + (q_sign & mask) + q = libdivide_s64_shift_right_vector(q, shift); // q >>= shift + q = _mm256_sub_epi64(_mm256_xor_si256(q, sign), sign); // q = (q ^ sign) - sign + return q; +} + +#elif defined(LIBDIVIDE_SSE2) + +static inline __m128i libdivide_u32_do_vector(__m128i numers, const struct libdivide_u32_t *denom); +static inline __m128i libdivide_s32_do_vector(__m128i numers, const struct libdivide_s32_t *denom); +static inline __m128i libdivide_u64_do_vector(__m128i numers, const struct libdivide_u64_t *denom); +static inline __m128i libdivide_s64_do_vector(__m128i numers, const struct libdivide_s64_t *denom); + +static inline __m128i libdivide_u32_branchfree_do_vector(__m128i numers, const struct libdivide_u32_branchfree_t *denom); +static inline __m128i libdivide_s32_branchfree_do_vector(__m128i numers, const struct libdivide_s32_branchfree_t *denom); +static inline __m128i libdivide_u64_branchfree_do_vector(__m128i numers, const struct libdivide_u64_branchfree_t *denom); +static inline __m128i libdivide_s64_branchfree_do_vector(__m128i numers, const struct libdivide_s64_branchfree_t *denom); + +//////// Internal Utility Functions + +// Implementation of _mm_srai_epi64(v, 63) (from AVX512). static inline __m128i libdivide_s64_signbits(__m128i v) { - // we want to compute v >> 63, that is, _mm_srai_epi64(v, 63). But there - // is no 64 bit shift right arithmetic instruction in SSE2. So we have to - // fake it by first duplicating the high 32 bit values, and then using a 32 - // bit shift. Another option would be to use _mm_srli_epi64(v, 63) and - // then subtract that from 0, but that approach appears to be substantially - // slower for unknown reasons __m128i hiBitsDuped = _mm_shuffle_epi32(v, _MM_SHUFFLE(3, 3, 1, 1)); __m128i signBits = _mm_srai_epi32(hiBitsDuped, 31); return signBits; } -// Returns an __m128i whose low 32 bits are equal to amt and has zero elsewhere. -static inline __m128i libdivide_u32_to_m128i(uint32_t amt) { - return _mm_set_epi32(0, 0, 0, amt); -} - +// Implementation of _mm_srai_epi64 (from AVX512). static inline __m128i libdivide_s64_shift_right_vector(__m128i v, int amt) { - // implementation of _mm_sra_epi64. Here we have two 64 bit values which - // are shifted right to logically become (64 - amt) values, and are then - // sign extended from a (64 - amt) bit number. const int b = 64 - amt; - __m128i m = libdivide_u64_to_m128(1ULL << (b - 1)); - __m128i x = _mm_srl_epi64(v, libdivide_u32_to_m128i(amt)); + __m128i m = _mm_set1_epi64x(1ULL << (b - 1)); + __m128i x = _mm_srli_epi64(v, amt); __m128i result = _mm_sub_epi64(_mm_xor_si128(x, m), m); return result; } -// Here, b is assumed to contain one 32 bit value repeated four times. -// If it did not, the function would not work. -static inline __m128i libdivide_mullhi_u32_flat_vector(__m128i a, __m128i b) { +// Here, b is assumed to contain one 32-bit value repeated. +static inline __m128i libdivide_mullhi_u32_vector(__m128i a, __m128i b) { __m128i hi_product_0Z2Z = _mm_srli_epi64(_mm_mul_epu32(a, b), 32); __m128i a1X3X = _mm_srli_epi64(a, 32); - __m128i mask = libdivide_get_FFFFFFFF00000000(); + __m128i mask = _mm_set_epi32(-1, 0, -1, 0); __m128i hi_product_Z1Z3 = _mm_and_si128(_mm_mul_epu32(a1X3X, b), mask); - // return hi_product_0123 return _mm_or_si128(hi_product_0Z2Z, hi_product_Z1Z3); } -// Here, y is assumed to contain one 64 bit value repeated twice. -static inline __m128i libdivide_mullhi_u64_flat_vector(__m128i x, __m128i y) { - // full 128 bits are x0 * y0 + (x0 * y1 << 32) + (x1 * y0 << 32) + (x1 * y1 << 64) - __m128i mask = libdivide_get_00000000FFFFFFFF(); - // x0 is low half of 2 64 bit values, x1 is high half in low slots - __m128i x0 = _mm_and_si128(x, mask); - __m128i x1 = _mm_srli_epi64(x, 32); - __m128i y0 = _mm_and_si128(y, mask); - __m128i y1 = _mm_srli_epi64(y, 32); - // x0 happens to have the low half of the two 64 bit values in 32 bit slots - // 0 and 2, so _mm_mul_epu32 computes their full product, and then we shift - // right by 32 to get just the high values - __m128i x0y0_hi = _mm_srli_epi64(_mm_mul_epu32(x0, y0), 32); - __m128i x0y1 = _mm_mul_epu32(x0, y1); - __m128i x1y0 = _mm_mul_epu32(x1, y0); - __m128i x1y1 = _mm_mul_epu32(x1, y1); - __m128i temp = _mm_add_epi64(x1y0, x0y0_hi); - __m128i temp_lo = _mm_and_si128(temp, mask); - __m128i temp_hi = _mm_srli_epi64(temp, 32); - temp_lo = _mm_srli_epi64(_mm_add_epi64(temp_lo, x0y1), 32); - temp_hi = _mm_add_epi64(x1y1, temp_hi); - - return _mm_add_epi64(temp_lo, temp_hi); -} - -// y is one 64 bit value repeated twice -static inline __m128i libdivide_mullhi_s64_flat_vector(__m128i x, __m128i y) { - __m128i p = libdivide_mullhi_u64_flat_vector(x, y); - __m128i t1 = _mm_and_si128(libdivide_s64_signbits(x), y); - p = _mm_sub_epi64(p, t1); - __m128i t2 = _mm_and_si128(libdivide_s64_signbits(y), x); - p = _mm_sub_epi64(p, t2); - return p; -} - -#ifdef LIBDIVIDE_USE_SSE4_1 - -// b is one 32 bit value repeated four times. -static inline __m128i libdivide_mullhi_s32_flat_vector(__m128i a, __m128i b) { - __m128i hi_product_0Z2Z = _mm_srli_epi64(_mm_mul_epi32(a, b), 32); - __m128i a1X3X = _mm_srli_epi64(a, 32); - __m128i mask = libdivide_get_FFFFFFFF00000000(); - __m128i hi_product_Z1Z3 = _mm_and_si128(_mm_mul_epi32(a1X3X, b), mask); - // return hi_product_0123 - return _mm_or_si128(hi_product_0Z2Z, hi_product_Z1Z3); -} - -#else - // SSE2 does not have a signed multiplication instruction, but we can convert // unsigned to signed pretty efficiently. Again, b is just a 32 bit value // repeated four times. -static inline __m128i libdivide_mullhi_s32_flat_vector(__m128i a, __m128i b) { - __m128i p = libdivide_mullhi_u32_flat_vector(a, b); +static inline __m128i libdivide_mullhi_s32_vector(__m128i a, __m128i b) { + __m128i p = libdivide_mullhi_u32_vector(a, b); // t1 = (a >> 31) & y, arithmetic shift __m128i t1 = _mm_and_si128(_mm_srai_epi32(a, 31), b); __m128i t2 = _mm_and_si128(_mm_srai_epi32(b, 31), a); @@ -1527,180 +1728,132 @@ static inline __m128i libdivide_mullhi_s32_flat_vector(__m128i a, __m128i b) { return p; } -#endif +// Here, y is assumed to contain one 64-bit value repeated. +// https://stackoverflow.com/a/28827013 +static inline __m128i libdivide_mullhi_u64_vector(__m128i x, __m128i y) { + __m128i lomask = _mm_set1_epi64x(0xffffffff); + __m128i xh = _mm_shuffle_epi32(x, 0xB1); // x0l, x0h, x1l, x1h + __m128i yh = _mm_shuffle_epi32(y, 0xB1); // y0l, y0h, y1l, y1h + __m128i w0 = _mm_mul_epu32(x, y); // x0l*y0l, x1l*y1l + __m128i w1 = _mm_mul_epu32(x, yh); // x0l*y0h, x1l*y1h + __m128i w2 = _mm_mul_epu32(xh, y); // x0h*y0l, x1h*y0l + __m128i w3 = _mm_mul_epu32(xh, yh); // x0h*y0h, x1h*y1h + __m128i w0h = _mm_srli_epi64(w0, 32); + __m128i s1 = _mm_add_epi64(w1, w0h); + __m128i s1l = _mm_and_si128(s1, lomask); + __m128i s1h = _mm_srli_epi64(s1, 32); + __m128i s2 = _mm_add_epi64(w2, s1l); + __m128i s2h = _mm_srli_epi64(s2, 32); + __m128i hi = _mm_add_epi64(w3, s1h); + hi = _mm_add_epi64(hi, s2h); + + return hi; +} + +// y is one 64-bit value repeated. +static inline __m128i libdivide_mullhi_s64_vector(__m128i x, __m128i y) { + __m128i p = libdivide_mullhi_u64_vector(x, y); + __m128i t1 = _mm_and_si128(libdivide_s64_signbits(x), y); + __m128i t2 = _mm_and_si128(libdivide_s64_signbits(y), x); + p = _mm_sub_epi64(p, t1); + p = _mm_sub_epi64(p, t2); + return p; +} ////////// UINT32 __m128i libdivide_u32_do_vector(__m128i numers, const struct libdivide_u32_t *denom) { uint8_t more = denom->more; - if (more & LIBDIVIDE_U32_SHIFT_PATH) { - uint32_t shift = more & LIBDIVIDE_32_SHIFT_MASK; - return _mm_srl_epi32(numers, libdivide_u32_to_m128i(shift)); + if (!denom->magic) { + return _mm_srli_epi32(numers, more); } else { - __m128i q = libdivide_mullhi_u32_flat_vector(numers, _mm_set1_epi32(denom->magic)); + __m128i q = libdivide_mullhi_u32_vector(numers, _mm_set1_epi32(denom->magic)); if (more & LIBDIVIDE_ADD_MARKER) { // uint32_t t = ((numer - q) >> 1) + q; // return t >> denom->shift; uint32_t shift = more & LIBDIVIDE_32_SHIFT_MASK; __m128i t = _mm_add_epi32(_mm_srli_epi32(_mm_sub_epi32(numers, q), 1), q); - return _mm_srl_epi32(t, libdivide_u32_to_m128i(shift)); + return _mm_srli_epi32(t, shift); } else { - // q >> denom->shift - return _mm_srl_epi32(q, libdivide_u32_to_m128i(more)); + return _mm_srli_epi32(q, more); } } } -__m128i libdivide_u32_do_vector_alg0(__m128i numers, const struct libdivide_u32_t *denom) { - return _mm_srl_epi32(numers, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_32_SHIFT_MASK)); -} - -__m128i libdivide_u32_do_vector_alg1(__m128i numers, const struct libdivide_u32_t *denom) { - __m128i q = libdivide_mullhi_u32_flat_vector(numers, _mm_set1_epi32(denom->magic)); - return _mm_srl_epi32(q, libdivide_u32_to_m128i(denom->more)); -} - -__m128i libdivide_u32_do_vector_alg2(__m128i numers, const struct libdivide_u32_t *denom) { - __m128i q = libdivide_mullhi_u32_flat_vector(numers, _mm_set1_epi32(denom->magic)); +__m128i libdivide_u32_branchfree_do_vector(__m128i numers, const struct libdivide_u32_branchfree_t *denom) { + __m128i q = libdivide_mullhi_u32_vector(numers, _mm_set1_epi32(denom->magic)); __m128i t = _mm_add_epi32(_mm_srli_epi32(_mm_sub_epi32(numers, q), 1), q); - return _mm_srl_epi32(t, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_32_SHIFT_MASK)); -} - -LIBDIVIDE_API __m128i libdivide_u32_branchfree_do_vector(__m128i numers, const struct libdivide_u32_branchfree_t *denom) { - __m128i q = libdivide_mullhi_u32_flat_vector(numers, _mm_set1_epi32(denom->magic)); - __m128i t = _mm_add_epi32(_mm_srli_epi32(_mm_sub_epi32(numers, q), 1), q); - return _mm_srl_epi32(t, libdivide_u32_to_m128i(denom->more)); + return _mm_srli_epi32(t, denom->more); } ////////// UINT64 __m128i libdivide_u64_do_vector(__m128i numers, const struct libdivide_u64_t *denom) { uint8_t more = denom->more; - if (more & LIBDIVIDE_U64_SHIFT_PATH) { - uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; - return _mm_srl_epi64(numers, libdivide_u32_to_m128i(shift)); + if (!denom->magic) { + return _mm_srli_epi64(numers, more); } else { - __m128i q = libdivide_mullhi_u64_flat_vector(numers, libdivide_u64_to_m128(denom->magic)); + __m128i q = libdivide_mullhi_u64_vector(numers, _mm_set1_epi64x(denom->magic)); if (more & LIBDIVIDE_ADD_MARKER) { // uint32_t t = ((numer - q) >> 1) + q; // return t >> denom->shift; uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; __m128i t = _mm_add_epi64(_mm_srli_epi64(_mm_sub_epi64(numers, q), 1), q); - return _mm_srl_epi64(t, libdivide_u32_to_m128i(shift)); + return _mm_srli_epi64(t, shift); } else { - // q >> denom->shift - return _mm_srl_epi64(q, libdivide_u32_to_m128i(more)); + return _mm_srli_epi64(q, more); } } } -__m128i libdivide_u64_do_vector_alg0(__m128i numers, const struct libdivide_u64_t *denom) { - return _mm_srl_epi64(numers, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_64_SHIFT_MASK)); -} - -__m128i libdivide_u64_do_vector_alg1(__m128i numers, const struct libdivide_u64_t *denom) { - __m128i q = libdivide_mullhi_u64_flat_vector(numers, libdivide_u64_to_m128(denom->magic)); - return _mm_srl_epi64(q, libdivide_u32_to_m128i(denom->more)); -} - -__m128i libdivide_u64_do_vector_alg2(__m128i numers, const struct libdivide_u64_t *denom) { - __m128i q = libdivide_mullhi_u64_flat_vector(numers, libdivide_u64_to_m128(denom->magic)); - __m128i t = _mm_add_epi64(_mm_srli_epi64(_mm_sub_epi64(numers, q), 1), q); - return _mm_srl_epi64(t, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_64_SHIFT_MASK)); -} - __m128i libdivide_u64_branchfree_do_vector(__m128i numers, const struct libdivide_u64_branchfree_t *denom) { - __m128i q = libdivide_mullhi_u64_flat_vector(numers, libdivide_u64_to_m128(denom->magic)); + __m128i q = libdivide_mullhi_u64_vector(numers, _mm_set1_epi64x(denom->magic)); __m128i t = _mm_add_epi64(_mm_srli_epi64(_mm_sub_epi64(numers, q), 1), q); - return _mm_srl_epi64(t, libdivide_u32_to_m128i(denom->more)); + return _mm_srli_epi64(t, denom->more); } ////////// SINT32 __m128i libdivide_s32_do_vector(__m128i numers, const struct libdivide_s32_t *denom) { uint8_t more = denom->more; - if (more & LIBDIVIDE_S32_SHIFT_PATH) { + if (!denom->magic) { uint32_t shift = more & LIBDIVIDE_32_SHIFT_MASK; uint32_t mask = (1U << shift) - 1; - // could use _mm_srli_epi32 with an all -1 register __m128i roundToZeroTweak = _mm_set1_epi32(mask); // q = numer + ((numer >> 31) & roundToZeroTweak); __m128i q = _mm_add_epi32(numers, _mm_and_si128(_mm_srai_epi32(numers, 31), roundToZeroTweak)); - q = _mm_sra_epi32(q, libdivide_u32_to_m128i(shift)); - // set all bits of shift mask = to the sign bit of more - __m128i shiftMask = _mm_set1_epi32((int32_t)((int8_t)more >> 7)); - // q = (q ^ shiftMask) - shiftMask; - q = _mm_sub_epi32(_mm_xor_si128(q, shiftMask), shiftMask); + q = _mm_srai_epi32(q, shift); + __m128i sign = _mm_set1_epi32((int8_t)more >> 7); + // q = (q ^ sign) - sign; + q = _mm_sub_epi32(_mm_xor_si128(q, sign), sign); return q; } else { - __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(denom->magic)); + __m128i q = libdivide_mullhi_s32_vector(numers, _mm_set1_epi32(denom->magic)); if (more & LIBDIVIDE_ADD_MARKER) { // must be arithmetic shift - __m128i sign = _mm_set1_epi32((int32_t)(int8_t)more >> 7); + __m128i sign = _mm_set1_epi32((int8_t)more >> 7); // q += ((numer ^ sign) - sign); q = _mm_add_epi32(q, _mm_sub_epi32(_mm_xor_si128(numers, sign), sign)); } // q >>= shift - q = _mm_sra_epi32(q, libdivide_u32_to_m128i(more & LIBDIVIDE_32_SHIFT_MASK)); + q = _mm_srai_epi32(q, more & LIBDIVIDE_32_SHIFT_MASK); q = _mm_add_epi32(q, _mm_srli_epi32(q, 31)); // q += (q < 0) return q; } } -__m128i libdivide_s32_do_vector_alg0(__m128i numers, const struct libdivide_s32_t *denom) { - uint8_t shift = denom->more & LIBDIVIDE_32_SHIFT_MASK; - uint32_t mask = (1U << shift) - 1; - __m128i roundToZeroTweak = _mm_set1_epi32(mask); - __m128i q = _mm_add_epi32(numers, _mm_and_si128(_mm_srai_epi32(numers, 31), roundToZeroTweak)); - return _mm_sra_epi32(q, libdivide_u32_to_m128i(shift)); -} - -__m128i libdivide_s32_do_vector_alg1(__m128i numers, const struct libdivide_s32_t *denom) { - uint8_t shift = denom->more & LIBDIVIDE_32_SHIFT_MASK; - uint32_t mask = (1U << shift) - 1; - __m128i roundToZeroTweak = _mm_set1_epi32(mask); - __m128i q = _mm_add_epi32(numers, _mm_and_si128(_mm_srai_epi32(numers, 31), roundToZeroTweak)); - return _mm_sub_epi32(_mm_setzero_si128(), _mm_sra_epi32(q, libdivide_u32_to_m128i(shift))); -} - -__m128i libdivide_s32_do_vector_alg2(__m128i numers, const struct libdivide_s32_t *denom) { - __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(denom->magic)); - q = _mm_add_epi32(q, numers); - q = _mm_sra_epi32(q, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_32_SHIFT_MASK)); - q = _mm_add_epi32(q, _mm_srli_epi32(q, 31)); - return q; -} - -__m128i libdivide_s32_do_vector_alg3(__m128i numers, const struct libdivide_s32_t *denom) { - __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(denom->magic)); - q = _mm_sub_epi32(q, numers); - q = _mm_sra_epi32(q, libdivide_u32_to_m128i(denom->more & LIBDIVIDE_32_SHIFT_MASK)); - q = _mm_add_epi32(q, _mm_srli_epi32(q, 31)); - return q; -} - -__m128i libdivide_s32_do_vector_alg4(__m128i numers, const struct libdivide_s32_t *denom) { - uint8_t more = denom->more; - __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(denom->magic)); - q = _mm_sra_epi32(q, libdivide_u32_to_m128i(more & LIBDIVIDE_32_SHIFT_MASK)); - q = _mm_add_epi32(q, _mm_srli_epi32(q, 31)); // q += (q < 0) - return q; -} - __m128i libdivide_s32_branchfree_do_vector(__m128i numers, const struct libdivide_s32_branchfree_t *denom) { int32_t magic = denom->magic; uint8_t more = denom->more; uint8_t shift = more & LIBDIVIDE_32_SHIFT_MASK; // must be arithmetic shift - __m128i sign = _mm_set1_epi32((int32_t)(int8_t)more >> 7); - - // libdivide_mullhi_s32(numers, magic); - __m128i q = libdivide_mullhi_s32_flat_vector(numers, _mm_set1_epi32(magic)); + __m128i sign = _mm_set1_epi32((int8_t)more >> 7); + __m128i q = libdivide_mullhi_s32_vector(numers, _mm_set1_epi32(magic)); q = _mm_add_epi32(q, numers); // q += numers // If q is non-negative, we have nothing to do @@ -1708,7 +1861,7 @@ __m128i libdivide_s32_branchfree_do_vector(__m128i numers, const struct libdivid // a power of 2, or (2**shift) if it is not a power of 2 uint32_t is_power_of_2 = (magic == 0); __m128i q_sign = _mm_srai_epi32(q, 31); // q_sign = q >> 31 - __m128i mask = _mm_set1_epi32((1 << shift) - is_power_of_2); + __m128i mask = _mm_set1_epi32((1U << shift) - is_power_of_2); q = _mm_add_epi32(q, _mm_and_si128(q_sign, mask)); // q = q + (q_sign & mask) q = _mm_srai_epi32(q, shift); // q >>= shift q = _mm_sub_epi32(_mm_xor_si128(q, sign), sign); // q = (q ^ sign) - sign @@ -1723,20 +1876,20 @@ __m128i libdivide_s64_do_vector(__m128i numers, const struct libdivide_s64_t *de if (magic == 0) { // shift path uint32_t shift = more & LIBDIVIDE_64_SHIFT_MASK; uint64_t mask = (1ULL << shift) - 1; - __m128i roundToZeroTweak = libdivide_u64_to_m128(mask); + __m128i roundToZeroTweak = _mm_set1_epi64x(mask); // q = numer + ((numer >> 63) & roundToZeroTweak); __m128i q = _mm_add_epi64(numers, _mm_and_si128(libdivide_s64_signbits(numers), roundToZeroTweak)); q = libdivide_s64_shift_right_vector(q, shift); - __m128i shiftMask = _mm_set1_epi32((int32_t)((int8_t)more >> 7)); - // q = (q ^ shiftMask) - shiftMask; - q = _mm_sub_epi64(_mm_xor_si128(q, shiftMask), shiftMask); + __m128i sign = _mm_set1_epi32((int8_t)more >> 7); + // q = (q ^ sign) - sign; + q = _mm_sub_epi64(_mm_xor_si128(q, sign), sign); return q; } else { - __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide_u64_to_m128(magic)); + __m128i q = libdivide_mullhi_s64_vector(numers, _mm_set1_epi64x(magic)); if (more & LIBDIVIDE_ADD_MARKER) { // must be arithmetic shift - __m128i sign = _mm_set1_epi32((int32_t)((int8_t)more >> 7)); + __m128i sign = _mm_set1_epi32((int8_t)more >> 7); // q += ((numer ^ sign) - sign); q = _mm_add_epi64(q, _mm_sub_epi64(_mm_xor_si128(numers, sign), sign)); } @@ -1747,56 +1900,15 @@ __m128i libdivide_s64_do_vector(__m128i numers, const struct libdivide_s64_t *de } } -__m128i libdivide_s64_do_vector_alg0(__m128i numers, const struct libdivide_s64_t *denom) { - uint32_t shift = denom->more & LIBDIVIDE_64_SHIFT_MASK; - uint64_t mask = (1ULL << shift) - 1; - __m128i roundToZeroTweak = libdivide_u64_to_m128(mask); - __m128i q = _mm_add_epi64(numers, _mm_and_si128(libdivide_s64_signbits(numers), roundToZeroTweak)); - q = libdivide_s64_shift_right_vector(q, shift); - return q; -} - -__m128i libdivide_s64_do_vector_alg1(__m128i numers, const struct libdivide_s64_t *denom) { - uint32_t shift = denom->more & LIBDIVIDE_64_SHIFT_MASK; - uint64_t mask = (1ULL << shift) - 1; - __m128i roundToZeroTweak = libdivide_u64_to_m128(mask); - __m128i q = _mm_add_epi64(numers, _mm_and_si128(libdivide_s64_signbits(numers), roundToZeroTweak)); - q = libdivide_s64_shift_right_vector(q, shift); - return _mm_sub_epi64(_mm_setzero_si128(), q); -} - -__m128i libdivide_s64_do_vector_alg2(__m128i numers, const struct libdivide_s64_t *denom) { - __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide_u64_to_m128(denom->magic)); - q = _mm_add_epi64(q, numers); - q = libdivide_s64_shift_right_vector(q, denom->more & LIBDIVIDE_64_SHIFT_MASK); - q = _mm_add_epi64(q, _mm_srli_epi64(q, 63)); // q += (q < 0) - return q; -} - -__m128i libdivide_s64_do_vector_alg3(__m128i numers, const struct libdivide_s64_t *denom) { - __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide_u64_to_m128(denom->magic)); - q = _mm_sub_epi64(q, numers); - q = libdivide_s64_shift_right_vector(q, denom->more & LIBDIVIDE_64_SHIFT_MASK); - q = _mm_add_epi64(q, _mm_srli_epi64(q, 63)); // q += (q < 0) - return q; -} - -__m128i libdivide_s64_do_vector_alg4(__m128i numers, const struct libdivide_s64_t *denom) { - __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide_u64_to_m128(denom->magic)); - q = libdivide_s64_shift_right_vector(q, denom->more & LIBDIVIDE_64_SHIFT_MASK); - q = _mm_add_epi64(q, _mm_srli_epi64(q, 63)); - return q; -} - __m128i libdivide_s64_branchfree_do_vector(__m128i numers, const struct libdivide_s64_branchfree_t *denom) { int64_t magic = denom->magic; uint8_t more = denom->more; uint8_t shift = more & LIBDIVIDE_64_SHIFT_MASK; // must be arithmetic shift - __m128i sign = _mm_set1_epi32((int32_t)(int8_t)more >> 7); + __m128i sign = _mm_set1_epi32((int8_t)more >> 7); // libdivide_mullhi_s64(numers, magic); - __m128i q = libdivide_mullhi_s64_flat_vector(numers, libdivide_u64_to_m128(magic)); + __m128i q = libdivide_mullhi_s64_vector(numers, _mm_set1_epi64x(magic)); q = _mm_add_epi64(q, numers); // q += numers // If q is non-negative, we have nothing to do. @@ -1804,7 +1916,7 @@ __m128i libdivide_s64_branchfree_do_vector(__m128i numers, const struct libdivid // a power of 2, or (2**shift) if it is not a power of 2. uint32_t is_power_of_2 = (magic == 0); __m128i q_sign = libdivide_s64_signbits(q); // q_sign = q >> 63 - __m128i mask = libdivide_u64_to_m128((1ULL << shift) - is_power_of_2); + __m128i mask = _mm_set1_epi64x((1ULL << shift) - is_power_of_2); q = _mm_add_epi64(q, _mm_and_si128(q_sign, mask)); // q = q + (q_sign & mask) q = libdivide_s64_shift_right_vector(q, shift); // q >>= shift q = _mm_sub_epi64(_mm_xor_si128(q, sign), sign); // q = (q ^ sign) - sign @@ -1817,290 +1929,143 @@ __m128i libdivide_s64_branchfree_do_vector(__m128i numers, const struct libdivid #ifdef __cplusplus -// Our divider struct is templated on both a type (like uint64_t) and an -// algorithm index. BRANCHFULL is the default algorithm, BRANCHFREE is the -// branchfree variant, and the indexed variants are for unswitching. +// The C++ divider class is templated on both an integer type +// (like uint64_t) and an algorithm type. +// * BRANCHFULL is the default algorithm type. +// * BRANCHFREE is the branchfree algorithm type. enum { - BRANCHFULL = -1, - BRANCHFREE = -2, - ALGORITHM0 = 0, - ALGORITHM1 = 1, - ALGORITHM2 = 2, - ALGORITHM3 = 3, - ALGORITHM4 = 4 + BRANCHFULL, + BRANCHFREE }; -namespace libdivide_internal { - -#if defined(LIBDIVIDE_USE_SSE2) && \ - defined(__GNUC__) && \ - __GNUC__ >= 6 - // Using vector functions as template arguments causes many - // -Wignored-attributes compiler warnings with GCC 9. - // These warnings can safely be turned off. - #pragma GCC diagnostic push - #pragma GCC diagnostic ignored "-Wignored-attributes" +#if defined(LIBDIVIDE_AVX512) + #define LIBDIVIDE_VECTOR_TYPE __m512i +#elif defined(LIBDIVIDE_AVX2) + #define LIBDIVIDE_VECTOR_TYPE __m256i +#elif defined(LIBDIVIDE_SSE2) + #define LIBDIVIDE_VECTOR_TYPE __m128i #endif -#if defined(LIBDIVIDE_USE_SSE2) - #define MAYBE_VECTOR(X) X - #define MAYBE_VECTOR_PARAM(X) __m128i vector_func(__m128i, const X *) +#if !defined(LIBDIVIDE_VECTOR_TYPE) + #define LIBDIVIDE_DIVIDE_VECTOR(ALGO) #else - #define MAYBE_VECTOR(X) 0 - #define MAYBE_VECTOR_PARAM(X) int unused + #define LIBDIVIDE_DIVIDE_VECTOR(ALGO) \ + LIBDIVIDE_VECTOR_TYPE divide(LIBDIVIDE_VECTOR_TYPE n) const { \ + return libdivide_##ALGO##_do_vector(n, &denom); \ + } #endif -// The following convenience macros are used to build a type of the base -// divider class and give it as template arguments the C functions -// related to the macro name and the macro type paramaters. - -#define BRANCHFULL_DIVIDER(INT, TYPE) \ - typedef base - -#define BRANCHFREE_DIVIDER(INT, TYPE) \ - typedef base - -#define ALGORITHM_DIVIDER(INT, TYPE, ALGO) \ - typedef base - -#define CRASH_DIVIDER(INT, TYPE) \ - typedef base - -// Base divider, provides storage for the actual divider. -// @T: e.g. uint32_t -// @DenomType: e.g. libdivide_u32_t -// @gen_func(): e.g. libdivide_u32_gen -// @do_func(): e.g. libdivide_u32_do -// @MAYBE_VECTOR_PARAM: e.g. libdivide_u32_do_vector -template -struct base { - // Storage for the actual divider - DenomType denom; - - // Constructor that takes a divisor value, and applies the gen function - base(T d) : denom(gen_func(d)) { } - - // Default constructor to allow uninitialized uses in e.g. arrays - base() {} - - // Needed for unswitch - base(const DenomType& d) : denom(d) { } - - T perform_divide(T val) const { - return do_func(val, &denom); +// The DISPATCHER_GEN() macro generates C++ methods (for the given integer +// and algorithm types) that redirect to libdivide's C API. +#define DISPATCHER_GEN(T, ALGO) \ + libdivide_##ALGO##_t denom; \ + dispatcher() { } \ + dispatcher(T d) \ + : denom(libdivide_##ALGO##_gen(d)) \ + { } \ + T divide(T n) const { \ + return libdivide_##ALGO##_do(n, &denom); \ + } \ + LIBDIVIDE_DIVIDE_VECTOR(ALGO) \ + T recover() const { \ + return libdivide_##ALGO##_recover(&denom); \ } -#if defined(LIBDIVIDE_USE_SSE2) - __m128i perform_divide_vector(__m128i val) const { - return vector_func(val, &denom); - } -#endif -}; - -// Functions that will never be called but are required to be able -// to use unswitch in C++ template code. Unsigned has fewer algorithms -// than signed i.e. alg3 and alg4 are not defined for unsigned. In -// order to make templates compile we need to define unsigned alg3 and -// alg4 as crash functions. -uint32_t libdivide_u32_crash(uint32_t, const libdivide_u32_t *) { exit(-1); } -uint64_t libdivide_u64_crash(uint64_t, const libdivide_u64_t *) { exit(-1); } - -#if defined(LIBDIVIDE_USE_SSE2) - __m128i libdivide_u32_crash_vector(__m128i, const libdivide_u32_t *) { exit(-1); } - __m128i libdivide_u64_crash_vector(__m128i, const libdivide_u64_t *) { exit(-1); } -#endif - +// The dispatcher selects a specific division algorithm for a given +// type and ALGO using partial template specialization. template struct dispatcher { }; -// Templated dispatch using partial specialization -template<> struct dispatcher { BRANCHFULL_DIVIDER(int32_t, s32) divider; }; -template<> struct dispatcher { BRANCHFREE_DIVIDER(int32_t, s32) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(int32_t, s32, alg0) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(int32_t, s32, alg1) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(int32_t, s32, alg2) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(int32_t, s32, alg3) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(int32_t, s32, alg4) divider; }; - -template<> struct dispatcher { BRANCHFULL_DIVIDER(uint32_t, u32) divider; }; -template<> struct dispatcher { BRANCHFREE_DIVIDER(uint32_t, u32) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(uint32_t, u32, alg0) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(uint32_t, u32, alg1) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(uint32_t, u32, alg2) divider; }; -template<> struct dispatcher { CRASH_DIVIDER(uint32_t, u32) divider; }; -template<> struct dispatcher { CRASH_DIVIDER(uint32_t, u32) divider; }; - -template<> struct dispatcher { BRANCHFULL_DIVIDER(int64_t, s64) divider; }; -template<> struct dispatcher { BRANCHFREE_DIVIDER(int64_t, s64) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER (int64_t, s64, alg0) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER (int64_t, s64, alg1) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER (int64_t, s64, alg2) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER (int64_t, s64, alg3) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER (int64_t, s64, alg4) divider; }; - -template<> struct dispatcher { BRANCHFULL_DIVIDER(uint64_t, u64) divider; }; -template<> struct dispatcher { BRANCHFREE_DIVIDER(uint64_t, u64) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(uint64_t, u64, alg0) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(uint64_t, u64, alg1) divider; }; -template<> struct dispatcher { ALGORITHM_DIVIDER(uint64_t, u64, alg2) divider; }; -template<> struct dispatcher { CRASH_DIVIDER(uint64_t, u64) divider; }; -template<> struct dispatcher { CRASH_DIVIDER(uint64_t, u64) divider; }; - -#if defined(LIBDIVIDE_USE_SSE2) && \ - defined(__GNUC__) && \ - __GNUC__ >= 6 - #pragma GCC diagnostic pop -#endif - -// Overloads that don't depend on the algorithm -inline int32_t recover(const libdivide_s32_t *s) { return libdivide_s32_recover(s); } -inline uint32_t recover(const libdivide_u32_t *s) { return libdivide_u32_recover(s); } -inline int64_t recover(const libdivide_s64_t *s) { return libdivide_s64_recover(s); } -inline uint64_t recover(const libdivide_u64_t *s) { return libdivide_u64_recover(s); } - -inline int32_t recover(const libdivide_s32_branchfree_t *s) { return libdivide_s32_branchfree_recover(s); } -inline uint32_t recover(const libdivide_u32_branchfree_t *s) { return libdivide_u32_branchfree_recover(s); } -inline int64_t recover(const libdivide_s64_branchfree_t *s) { return libdivide_s64_branchfree_recover(s); } -inline uint64_t recover(const libdivide_u64_branchfree_t *s) { return libdivide_u64_branchfree_recover(s); } - -inline int get_algorithm(const libdivide_s32_t *s) { return libdivide_s32_get_algorithm(s); } -inline int get_algorithm(const libdivide_u32_t *s) { return libdivide_u32_get_algorithm(s); } -inline int get_algorithm(const libdivide_s64_t *s) { return libdivide_s64_get_algorithm(s); } -inline int get_algorithm(const libdivide_u64_t *s) { return libdivide_u64_get_algorithm(s); } - -// Fallback for branchfree variants, which do not support unswitching -template int get_algorithm(const T *) { return -1; } - -} // namespace libdivide_internal +template<> struct dispatcher { DISPATCHER_GEN(int32_t, s32) }; +template<> struct dispatcher { DISPATCHER_GEN(int32_t, s32_branchfree) }; +template<> struct dispatcher { DISPATCHER_GEN(uint32_t, u32) }; +template<> struct dispatcher { DISPATCHER_GEN(uint32_t, u32_branchfree) }; +template<> struct dispatcher { DISPATCHER_GEN(int64_t, s64) }; +template<> struct dispatcher { DISPATCHER_GEN(int64_t, s64_branchfree) }; +template<> struct dispatcher { DISPATCHER_GEN(uint64_t, u64) }; +template<> struct dispatcher { DISPATCHER_GEN(uint64_t, u64_branchfree) }; // This is the main divider class for use by the user (C++ API). -// The divider itself is stored in the div variable who's -// type is chosen by the dispatcher based on the template paramaters. +// The actual division algorithm is selected using the dispatcher struct +// based on the integer and algorithm template parameters. template class divider { -private: - // Here's the actual divider - typedef typename libdivide_internal::dispatcher::divider div_t; - div_t div; - - // unswitch() friend declaration - template - friend divider unswitch(const divider & d); - - // Constructor used by the unswitch friend - divider(const div_t& denom) : div(denom) { } - public: - // Ordinary constructor that takes the divisor as a parameter - divider(T n) : div(n) { } - - // Default constructor. We leave this deliberately undefined so that - // creating an array of divider and then initializing them - // doesn't slow us down. + // We leave the default constructor empty so that creating + // an array of dividers and then initializing them + // later doesn't slow us down. divider() { } - // Divides the parameter by the divisor, returning the quotient - T perform_divide(T val) const { - return div.perform_divide(val); + // Constructor that takes the divisor as a parameter + divider(T d) : div(d) { } + + // Divides n by the divisor + T divide(T n) const { + return div.divide(n); } - // Recovers the divisor that was used to initialize the divider - T recover_divisor() const { - return libdivide_internal::recover(&div.denom); + // Recovers the divisor, returns the value that was + // used to initialize this divider object. + T recover() const { + return div.recover(); } -#if defined(LIBDIVIDE_USE_SSE2) - // Treats the vector as either two or four packed values (depending on the - // size), and divides each of them by the divisor, - // returning the packed quotients. - __m128i perform_divide_vector(__m128i val) const { - return div.perform_divide_vector(val); + bool operator==(const divider& other) const { + return div.denom.magic == other.denom.magic && + div.denom.more == other.denom.more; + } + + bool operator!=(const divider& other) const { + return !(*this == other); + } + +#if defined(LIBDIVIDE_VECTOR_TYPE) + // Treats the vector as packed integer values with the same type as + // the divider (e.g. s32, u32, s64, u64) and divides each of + // them by the divider, returning the packed quotients. + LIBDIVIDE_VECTOR_TYPE divide(LIBDIVIDE_VECTOR_TYPE n) const { + return div.divide(n); } #endif - - // Returns the index of algorithm, for use in the unswitch function. Does - // not apply to branchfree variant. - // Returns the algorithm for unswitching. - int get_algorithm() const { - return libdivide_internal::get_algorithm(&div.denom); - } - - bool operator==(const divider& him) const { - return div.denom.magic == him.div.denom.magic && - div.denom.more == him.div.denom.more; - } - - bool operator!=(const divider& him) const { - return !(*this == him); - } +private: + // Storage for the actual divisor + dispatcher div; }; +// Overload of operator / for scalar division +template +T operator/(T n, const divider& div) { + return div.divide(n); +} + +// Overload of operator /= for scalar division +template +T& operator/=(T& n, const divider& div) { + n = div.divide(n); + return n; +} + +#if defined(LIBDIVIDE_VECTOR_TYPE) + // Overload of operator / for vector division + template + LIBDIVIDE_VECTOR_TYPE operator/(LIBDIVIDE_VECTOR_TYPE n, const divider& div) { + return div.divide(n); + } + // Overload of operator /= for vector division + template + LIBDIVIDE_VECTOR_TYPE& operator/=(LIBDIVIDE_VECTOR_TYPE& n, const divider& div) { + n = div.divide(n); + return n; + } +#endif + #if __cplusplus >= 201103L || \ (defined(_MSC_VER) && _MSC_VER >= 1800) - -// libdivdie::branchfree_divider -template -using branchfree_divider = divider; - -#endif - -// Returns a divider specialized for the given algorithm -template -divider unswitch(const divider& d) { - return divider(d.div.denom); -} - -// Overload of the / operator for scalar division -template -T operator/(T numer, const divider& denom) { - return denom.perform_divide(numer); -} - -// Overload of the /= operator for scalar division -template -T operator/=(T& numer, const divider& denom) { - numer = denom.perform_divide(numer); - return numer; -} - -#if defined(LIBDIVIDE_USE_SSE2) - -// Overload of the / operator for vector division -template -__m128i operator/(__m128i numer, const divider& denom) { - return denom.perform_divide_vector(numer); -} - -// Overload of the /= operator for vector division -template -__m128i operator/=(__m128i& numer, const divider& denom) { - numer = denom.perform_divide_vector(numer); - return numer; -} - + // libdivdie::branchfree_divider + template + using branchfree_divider = divider; #endif } // namespace libdivide -} // anonymous namespace #endif // __cplusplus