From 2940c257a640e7c0f40a457c513a1bc199c204a4 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Marek=20Ol=C5=A1=C3=A1k?= Date: Sun, 23 Sep 2018 12:57:51 -0400 Subject: [PATCH] util: import public domain code for integer division by a constant MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit 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. Reviewed-by: Jason Ekstrand Reviewed-by: Marek Olšák --- src/util/Makefile.sources | 2 + src/util/fast_idiv_by_const.c | 224 ++++++++++++++++++++++++++++++++++ src/util/fast_idiv_by_const.h | 137 +++++++++++++++++++++ src/util/meson.build | 2 + 4 files changed, 365 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 5b1548c733c..e8558f561f8 100644 --- a/src/util/Makefile.sources +++ b/src/util/Makefile.sources @@ -11,6 +11,8 @@ MESA_UTIL_FILES := \ 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 \ diff --git a/src/util/fast_idiv_by_const.c b/src/util/fast_idiv_by_const.c new file mode 100644 index 00000000000..0bc9b60878b --- /dev/null +++ b/src/util/fast_idiv_by_const.c @@ -0,0 +1,224 @@ +/* + * 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 +#include + +/* 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; + + /* 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 00000000000..ac10cf79ba8 --- /dev/null +++ b/src/util/fast_idiv_by_const.h @@ -0,0 +1,137 @@ +/* + * 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 +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +/* 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); + +#ifdef __cplusplus +} /* extern C */ +#endif + +#endif /* FAST_IDIV_BY_CONST_H */ diff --git a/src/util/meson.build b/src/util/meson.build index 9a99d60c158..cdbad98e7cb 100644 --- a/src/util/meson.build +++ b/src/util/meson.build @@ -35,6 +35,8 @@ files_mesa_util = files( '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', -- 2.30.2