On Tue, Oct 2, 2018, 1:15 PM Jason Ekstrand <ja...@jlekstrand.net> wrote:
> Reading through things in a bit more detail, I do believe that importing > this version in some form would be better than using mine for a number of > reasons: > > * It is better optimized for signed integers > * The struct of division factors is much better than what I did. (I did > consider a struct and discarded the idea; I was wrong). > * Computation of the division factors doesn't involve N*2-bit > multiplication > * The round-up algorithm here results in significantly better code than > the N+1-bit round-down. > * I trust ridiculousfish to get this right more than I trust myself > > That said, I have a few caveats on merging this as-is: > > * I would like to see some unit tests. I already spent the time to write > some; they just have to be ported. > * It needs to be adjusted to handle 64-bit integers (right now, it > appears to only work for num_bits <= 32) > * We shouldn't define uint_t and sint_t in a header > > How do you want to proceed? > I don't have a plan. Anything that works for you would be OK with me, so if you wanna just rework it according to you, that's fine. Changing the types is tricky. Template code in a C header included several times would work. C++ templates would be ideal. What's your timeframe for this? Mine is certainly more than a month. Marek > --Jason > > > On Sun, Sep 23, 2018 at 11:58 AM Marek Olšák <mar...@gmail.com> wrote: > >> From: Marek Olšák <marek.ol...@amd.com> >> >> Compilers can use this to generate optimal code for integer division >> by a constant. >> >> Additionally, an unsigned division by a uniform that is constant but not >> known at compile time can still be optimized by passing 2-4 division >> factors to the shader as uniforms and executing one of the fast_udiv* >> variants. The signed division algorithm doesn't have this capability. >> --- >> src/util/Makefile.sources | 2 + >> src/util/fast_idiv_by_const.c | 245 >> ++++++++++++++++++++++++++++++++++++++++++ >> src/util/fast_idiv_by_const.h | 173 +++++++++++++++++++++++++++++ >> src/util/meson.build | 2 + >> 4 files changed, 422 insertions(+) >> create mode 100644 src/util/fast_idiv_by_const.c >> create mode 100644 src/util/fast_idiv_by_const.h >> >> diff --git a/src/util/Makefile.sources b/src/util/Makefile.sources >> index b562d6c..f741b2a 100644 >> --- a/src/util/Makefile.sources >> +++ b/src/util/Makefile.sources >> @@ -3,20 +3,22 @@ MESA_UTIL_FILES := \ >> bitscan.h \ >> bitset.h \ >> build_id.c \ >> build_id.h \ >> crc32.c \ >> crc32.h \ >> debug.c \ >> debug.h \ >> disk_cache.c \ >> disk_cache.h \ >> + fast_idiv_by_const.c \ >> + fast_idiv_by_const.h \ >> format_r11g11b10f.h \ >> format_rgb9e5.h \ >> format_srgb.h \ >> futex.h \ >> half_float.c \ >> half_float.h \ >> hash_table.c \ >> hash_table.h \ >> list.h \ >> macros.h \ >> diff --git a/src/util/fast_idiv_by_const.c b/src/util/fast_idiv_by_const.c >> new file mode 100644 >> index 0000000..f247b66 >> --- /dev/null >> +++ b/src/util/fast_idiv_by_const.c >> @@ -0,0 +1,245 @@ >> +/* >> + * Copyright © 2018 Advanced Micro Devices, Inc. >> + * >> + * Permission is hereby granted, free of charge, to any person obtaining >> a >> + * copy of this software and associated documentation files (the >> "Software"), >> + * to deal in the Software without restriction, including without >> limitation >> + * the rights to use, copy, modify, merge, publish, distribute, >> sublicense, >> + * and/or sell copies of the Software, and to permit persons to whom the >> + * Software is furnished to do so, subject to the following conditions: >> + * >> + * The above copyright notice and this permission notice (including the >> next >> + * paragraph) shall be included in all copies or substantial portions of >> the >> + * Software. >> + * >> + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, >> EXPRESS OR >> + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF >> MERCHANTABILITY, >> + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT >> SHALL >> + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR >> OTHER >> + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, >> ARISING >> + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER >> DEALINGS >> + * IN THE SOFTWARE. >> + */ >> + >> +/* Imported from: >> + * >> https://raw.githubusercontent.com/ridiculousfish/libdivide/master/divide_by_constants_codegen_reference.c >> + * Paper: >> + * >> http://ridiculousfish.com/files/faster_unsigned_division_by_constants.pdf >> + * >> + * The author, ridiculous_fish, wrote: >> + * >> + * ''Reference implementations of computing and using the "magic number" >> + * approach to dividing by constants, including codegen instructions. >> + * The unsigned division incorporates the "round down" optimization >> per >> + * ridiculous_fish. >> + * >> + * This is free and unencumbered software. Any copyright is dedicated >> + * to the Public Domain.'' >> + */ >> + >> +#include "fast_idiv_by_const.h" >> +#include "u_math.h" >> +#include <limits.h> >> +#include <assert.h> >> + >> +/* uint_t and sint_t can be replaced by different integer types and the >> code >> + * will work as-is. The only requirement is that sizeof(uintN) == >> sizeof(intN). >> + */ >> + >> +struct util_fast_udiv_info >> +util_compute_fast_udiv_info(uint_t D, unsigned num_bits) >> +{ >> + /* The numerator must fit in a uint_t */ >> + assert(num_bits > 0 && num_bits <= sizeof(uint_t) * CHAR_BIT); >> + assert(D != 0); >> + >> + /* The eventual result */ >> + struct util_fast_udiv_info result; >> + >> + if (util_is_power_of_two_nonzero(D)) { >> + unsigned div_shift = util_logbase2(D); >> + >> + if (div_shift) { >> + /* Dividing by a power of two. */ >> + result.multiplier = 1 << 31; >> + result.pre_shift = 0; >> + result.post_shift = div_shift - 1; >> + result.increment = 0; >> + return result; >> + } else { >> + /* Dividing by 1. */ >> + /* Assuming: floor((num + 1) * (2^32 - 1) / 2^32) = num */ >> + result.multiplier = UINT_MAX; >> + result.pre_shift = 0; >> + result.post_shift = 0; >> + result.increment = 1; >> + return result; >> + } >> + } >> + >> + /* Bits in a uint_t */ >> + const unsigned UINT_BITS = sizeof(uint_t) * CHAR_BIT; >> + >> + /* The extra shift implicit in the difference between UINT_BITS and >> num_bits >> + */ >> + const unsigned extra_shift = UINT_BITS - num_bits; >> + >> + /* The initial power of 2 is one less than the first one that can >> possibly >> + * work. >> + */ >> + const uint_t initial_power_of_2 = (uint_t)1 << (UINT_BITS-1); >> + >> + /* The remainder and quotient of our power of 2 divided by d */ >> + uint_t quotient = initial_power_of_2 / D; >> + uint_t remainder = initial_power_of_2 % D; >> + >> + /* ceil(log_2 D) */ >> + unsigned ceil_log_2_D; >> + >> + /* The magic info for the variant "round down" algorithm */ >> + uint_t down_multiplier = 0; >> + unsigned down_exponent = 0; >> + int has_magic_down = 0; >> + >> + /* Compute ceil(log_2 D) */ >> + ceil_log_2_D = 0; >> + uint_t tmp; >> + for (tmp = D; tmp > 0; tmp >>= 1) >> + ceil_log_2_D += 1; >> + >> + >> + /* Begin a loop that increments the exponent, until we find a power >> of 2 >> + * that works. >> + */ >> + unsigned exponent; >> + for (exponent = 0; ; exponent++) { >> + /* Quotient and remainder is from previous exponent; compute it >> for this >> + * exponent. >> + */ >> + if (remainder >= D - remainder) { >> + /* Doubling remainder will wrap around D */ >> + quotient = quotient * 2 + 1; >> + remainder = remainder * 2 - D; >> + } else { >> + /* Remainder will not wrap */ >> + quotient = quotient * 2; >> + remainder = remainder * 2; >> + } >> + >> + /* We're done if this exponent works for the round_up algorithm. >> + * Note that exponent may be larger than the maximum shift >> supported, >> + * so the check for >= ceil_log_2_D is critical. >> + */ >> + if ((exponent + extra_shift >= ceil_log_2_D) || >> + (D - remainder) <= ((uint_t)1 << (exponent + extra_shift))) >> + break; >> + >> + /* Set magic_down if we have not set it yet and this exponent >> works for >> + * the round_down algorithm >> + */ >> + if (!has_magic_down && >> + remainder <= ((uint_t)1 << (exponent + extra_shift))) { >> + has_magic_down = 1; >> + down_multiplier = quotient; >> + down_exponent = exponent; >> + } >> + } >> + >> + if (exponent < ceil_log_2_D) { >> + /* magic_up is efficient */ >> + result.multiplier = quotient + 1; >> + result.pre_shift = 0; >> + result.post_shift = exponent; >> + result.increment = 0; >> + } else if (D & 1) { >> + /* Odd divisor, so use magic_down, which must have been set */ >> + assert(has_magic_down); >> + result.multiplier = down_multiplier; >> + result.pre_shift = 0; >> + result.post_shift = down_exponent; >> + result.increment = 1; >> + } else { >> + /* Even divisor, so use a prefix-shifted dividend */ >> + unsigned pre_shift = 0; >> + uint_t shifted_D = D; >> + while ((shifted_D & 1) == 0) { >> + shifted_D >>= 1; >> + pre_shift += 1; >> + } >> + result = util_compute_fast_udiv_info(shifted_D, num_bits - >> pre_shift); >> + /* expect no increment or pre_shift in this path */ >> + assert(result.increment == 0 && result.pre_shift == 0); >> + result.pre_shift = pre_shift; >> + } >> + return result; >> +} >> + >> +struct util_fast_sdiv_info >> +util_compute_fast_sdiv_info(sint_t D) >> +{ >> + /* D must not be zero. */ >> + assert(D != 0); >> + /* The result is not correct for these divisors. */ >> + assert(D != 1 && D != -1); >> + >> + /* Our result */ >> + struct util_fast_sdiv_info result; >> + >> + /* Bits in an sint_t */ >> + const unsigned SINT_BITS = sizeof(sint_t) * CHAR_BIT; >> + >> + /* Absolute value of D (we know D is not the most negative value since >> + * that's a power of 2) >> + */ >> + const uint_t abs_d = (D < 0 ? -D : D); >> + >> + /* The initial power of 2 is one less than the first one that can >> possibly >> + * work */ >> + /* "two31" in Warren */ >> + unsigned exponent = SINT_BITS - 1; >> + const uint_t initial_power_of_2 = (uint_t)1 << exponent; >> + >> + /* Compute the absolute value of our "test numerator," >> + * which is the largest dividend whose remainder with d is d-1. >> + * This is called anc in Warren. >> + */ >> + const uint_t tmp = initial_power_of_2 + (D < 0); >> + const uint_t abs_test_numer = tmp - 1 - tmp % abs_d; >> + >> + /* Initialize our quotients and remainders (q1, r1, q2, r2 in Warren) >> */ >> + uint_t quotient1 = initial_power_of_2 / abs_test_numer; >> + uint_t remainder1 = initial_power_of_2 % abs_test_numer; >> + uint_t quotient2 = initial_power_of_2 / abs_d; >> + uint_t remainder2 = initial_power_of_2 % abs_d; >> + uint_t delta; >> + >> + /* Begin our loop */ >> + do { >> + /* Update the exponent */ >> + exponent++; >> + >> + /* Update quotient1 and remainder1 */ >> + quotient1 *= 2; >> + remainder1 *= 2; >> + if (remainder1 >= abs_test_numer) { >> + quotient1 += 1; >> + remainder1 -= abs_test_numer; >> + } >> + >> + /* Update quotient2 and remainder2 */ >> + quotient2 *= 2; >> + remainder2 *= 2; >> + if (remainder2 >= abs_d) { >> + quotient2 += 1; >> + remainder2 -= abs_d; >> + } >> + >> + /* Keep going as long as (2**exponent) / abs_d <= delta */ >> + delta = abs_d - remainder2; >> + } while (quotient1 < delta || (quotient1 == delta && remainder1 == >> 0)); >> + >> + result.multiplier = quotient2 + 1; >> + if (D < 0) result.multiplier = -result.multiplier; >> + result.shift = exponent - SINT_BITS; >> + return result; >> +} >> diff --git a/src/util/fast_idiv_by_const.h b/src/util/fast_idiv_by_const.h >> new file mode 100644 >> index 0000000..e8debbf >> --- /dev/null >> +++ b/src/util/fast_idiv_by_const.h >> @@ -0,0 +1,173 @@ >> +/* >> + * Copyright © 2018 Advanced Micro Devices, Inc. >> + * >> + * Permission is hereby granted, free of charge, to any person obtaining >> a >> + * copy of this software and associated documentation files (the >> "Software"), >> + * to deal in the Software without restriction, including without >> limitation >> + * the rights to use, copy, modify, merge, publish, distribute, >> sublicense, >> + * and/or sell copies of the Software, and to permit persons to whom the >> + * Software is furnished to do so, subject to the following conditions: >> + * >> + * The above copyright notice and this permission notice (including the >> next >> + * paragraph) shall be included in all copies or substantial portions of >> the >> + * Software. >> + * >> + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, >> EXPRESS OR >> + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF >> MERCHANTABILITY, >> + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT >> SHALL >> + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR >> OTHER >> + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, >> ARISING >> + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER >> DEALINGS >> + * IN THE SOFTWARE. >> + */ >> + >> +#ifndef FAST_IDIV_BY_CONST_H >> +#define FAST_IDIV_BY_CONST_H >> + >> +/* Imported from: >> + * >> https://raw.githubusercontent.com/ridiculousfish/libdivide/master/divide_by_constants_codegen_reference.c >> + */ >> + >> +#include <inttypes.h> >> +#include <limits.h> >> +#include <assert.h> >> + >> +/* You can set these to different types to get different precision. */ >> +typedef int32_t sint_t; >> +typedef uint32_t uint_t; >> + >> +/* Computes "magic info" for performing signed division by a fixed >> integer D. >> + * The type 'sint_t' is assumed to be defined as a signed integer type >> large >> + * enough to hold both the dividend and the divisor. >> + * Here >> is arithmetic (signed) shift, and >>> is logical shift. >> + * >> + * To emit code for n/d, rounding towards zero, use the following >> sequence: >> + * >> + * m = compute_signed_magic_info(D) >> + * emit("result = (m.multiplier * n) >> SINT_BITS"); >> + * if d > 0 and m.multiplier < 0: emit("result += n") >> + * if d < 0 and m.multiplier > 0: emit("result -= n") >> + * if m.post_shift > 0: emit("result >>= m.shift") >> + * emit("result += (result < 0)") >> + * >> + * The shifts by SINT_BITS may be "free" if the high half of the full >> multiply >> + * is put in a separate register. >> + * >> + * The final add can of course be implemented via the sign bit, e.g. >> + * result += (result >>> (SINT_BITS - 1)) >> + * or >> + * result -= (result >> (SINT_BITS - 1)) >> + * >> + * This code is heavily indebted to Hacker's Delight by Henry Warren. >> + * See http://www.hackersdelight.org/HDcode/magic.c.txt >> + * Used with permission from >> http://www.hackersdelight.org/permissions.htm >> + */ >> + >> +struct util_fast_sdiv_info { >> + sint_t multiplier; /* the "magic number" multiplier */ >> + unsigned shift; /* shift for the dividend after multiplying */ >> +}; >> + >> +struct util_fast_sdiv_info >> +util_compute_fast_sdiv_info(sint_t D); >> + >> +/* Computes "magic info" for performing unsigned division by a fixed >> positive >> + * integer D. The type 'uint_t' is assumed to be defined as an unsigned >> + * integer type large enough to hold both the dividend and the divisor. >> + * num_bits can be set appropriately if n is known to be smaller than >> + * the largest uint_t; if this is not known then pass >> + * "(sizeof(uint_t) * CHAR_BIT)" for num_bits. >> + * >> + * Assume we have a hardware register of width UINT_BITS, a known >> constant D >> + * which is not zero and not a power of 2, and a variable n of width >> num_bits >> + * (which may be up to UINT_BITS). To emit code for n/d, use one of the >> two >> + * following sequences (here >>> refers to a logical bitshift): >> + * >> + * m = compute_unsigned_magic_info(D, num_bits) >> + * if m.pre_shift > 0: emit("n >>>= m.pre_shift") >> + * if m.increment: emit("n = saturated_increment(n)") >> + * emit("result = (m.multiplier * n) >>> UINT_BITS") >> + * if m.post_shift > 0: emit("result >>>= m.post_shift") >> + * >> + * or >> + * >> + * m = compute_unsigned_magic_info(D, num_bits) >> + * if m.pre_shift > 0: emit("n >>>= m.pre_shift") >> + * emit("result = m.multiplier * n") >> + * if m.increment: emit("result = result + m.multiplier") >> + * emit("result >>>= UINT_BITS") >> + * if m.post_shift > 0: emit("result >>>= m.post_shift") >> + * >> + * The shifts by UINT_BITS may be "free" if the high half of the full >> multiply >> + * is put in a separate register. >> + * >> + * saturated_increment(n) means "increment n unless it would wrap to 0," >> i.e. >> + * if n == (1 << UINT_BITS)-1: result = n >> + * else: result = n+1 >> + * A common way to implement this is with the carry bit. For example, on >> x86: >> + * add 1 >> + * sbb 0 >> + * >> + * Some invariants: >> + * 1: At least one of pre_shift and increment is zero >> + * 2: multiplier is never zero >> + * >> + * This code incorporates the "round down" optimization per >> ridiculous_fish. >> + */ >> + >> +struct util_fast_udiv_info { >> + uint_t multiplier; /* the "magic number" multiplier */ >> + unsigned pre_shift; /* shift for the dividend before multiplying */ >> + unsigned post_shift; /* shift for the dividend after multiplying */ >> + int increment; /* 0 or 1; if set then increment the numerator, using >> one of >> + the two strategies */ >> +}; >> + >> +struct util_fast_udiv_info >> +util_compute_fast_udiv_info(uint_t D, unsigned num_bits); >> + >> +/* Below are possible options for dividing by a uniform in a shader where >> + * the divisor is constant but not known at compile time. >> + */ >> + >> +/* Full version. */ >> +static inline unsigned >> +fast_udiv(unsigned n, struct util_fast_udiv_info info) >> +{ >> + n = n >> info.pre_shift; >> + /* For non-power-of-two divisors, use a 32-bit ADD that clamps to >> UINT_MAX. */ >> + n = (((uint64_t)n + info.increment) * info.multiplier) >> 32; >> + n = n >> info.post_shift; >> + return n; >> +} >> + >> +/* A little more efficient version if n != UINT_MAX, i.e. no unsigned >> + * wraparound in the computation. >> + */ >> +static inline unsigned >> +fast_udiv_nuw(unsigned n, struct util_fast_udiv_info info) >> +{ >> + assert(n != UINT_MAX); >> + n = n >> info.pre_shift; >> + n = n + info.increment; >> + n = ((uint64_t)n * info.multiplier) >> 32; >> + n = n >> info.post_shift; >> + return n; >> +} >> + >> +/* Even faster version but both operands must be 31-bit unsigned integers >> + * and the divisor must be greater than 1. >> + * >> + * info must be computed with num_bits == 31. >> + */ >> +static inline unsigned >> +fast_udiv_u31_d_not_one(unsigned n, struct util_fast_udiv_info info) >> +{ >> + assert(info.pre_shift == 0); >> + assert(info.increment == 0); >> + n = ((uint64_t)n * info.multiplier) >> 32; >> + n = n >> info.post_shift; >> + return n; >> +} >> + >> +#endif >> diff --git a/src/util/meson.build b/src/util/meson.build >> index 027bc5b..ebaeb47 100644 >> --- a/src/util/meson.build >> +++ b/src/util/meson.build >> @@ -27,20 +27,22 @@ files_mesa_util = files( >> 'bitscan.h', >> 'bitset.h', >> 'build_id.c', >> 'build_id.h', >> 'crc32.c', >> 'crc32.h', >> 'debug.c', >> 'debug.h', >> 'disk_cache.c', >> 'disk_cache.h', >> + 'fast_idiv_by_const.c', >> + 'fast_idiv_by_const.h', >> 'format_r11g11b10f.h', >> 'format_rgb9e5.h', >> 'format_srgb.h', >> 'futex.h', >> 'half_float.c', >> 'half_float.h', >> 'hash_table.c', >> 'hash_table.h', >> 'list.h', >> 'macros.h', >> -- >> 2.7.4 >> >> _______________________________________________ >> mesa-dev mailing list >> mesa-dev@lists.freedesktop.org >> https://lists.freedesktop.org/mailman/listinfo/mesa-dev >> >
_______________________________________________ mesa-dev mailing list mesa-dev@lists.freedesktop.org https://lists.freedesktop.org/mailman/listinfo/mesa-dev
