| /* |
| * 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> |
| |
| struct util_fast_udiv_info |
| util_compute_fast_udiv_info(uint64_t D, unsigned num_bits, unsigned UINT_BITS) |
| { |
| /* The numerator must fit in a uint64_t */ |
| assert(num_bits > 0 && num_bits <= UINT_BITS); |
| assert(D != 0); |
| |
| /* The eventual result */ |
| struct util_fast_udiv_info result; |
| |
| if (util_is_power_of_two_or_zero64(D)) { |
| unsigned div_shift = util_logbase2_64(D); |
| |
| if (div_shift) { |
| /* Dividing by a power of two. */ |
| result.multiplier = 1ull << (UINT_BITS - div_shift); |
| result.pre_shift = 0; |
| result.post_shift = 0; |
| result.increment = 0; |
| return result; |
| } else { |
| /* Dividing by 1. */ |
| /* Assuming: floor((num + 1) * (2^32 - 1) / 2^32) = num */ |
| result.multiplier = UINT_BITS == 64 ? UINT64_MAX : |
| (1ull << UINT_BITS) - 1; |
| result.pre_shift = 0; |
| result.post_shift = 0; |
| result.increment = 1; |
| return result; |
| } |
| } |
| |
| /* 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 uint64_t initial_power_of_2 = (uint64_t)1 << (UINT_BITS-1); |
| |
| /* The remainder and quotient of our power of 2 divided by d */ |
| uint64_t quotient = initial_power_of_2 / D; |
| uint64_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 */ |
| uint64_t down_multiplier = 0; |
| unsigned down_exponent = 0; |
| int has_magic_down = 0; |
| |
| /* Compute ceil(log_2 D) */ |
| ceil_log_2_D = 0; |
| uint64_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) <= ((uint64_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 <= ((uint64_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; |
| uint64_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, |
| UINT_BITS); |
| /* expect no increment or pre_shift in this path */ |
| assert(result.increment == 0 && result.pre_shift == 0); |
| result.pre_shift = pre_shift; |
| } |
| return result; |
| } |
| |
| static inline int64_t |
| sign_extend(int64_t x, unsigned SINT_BITS) |
| { |
| return (int64_t)((uint64_t)x << (64 - SINT_BITS)) >> (64 - SINT_BITS); |
| } |
| |
| struct util_fast_sdiv_info |
| util_compute_fast_sdiv_info(int64_t D, unsigned SINT_BITS) |
| { |
| /* 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; |
| |
| /* Absolute value of D (we know D is not the most negative value since |
| * that's a power of 2) |
| */ |
| const uint64_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 uint64_t initial_power_of_2 = (uint64_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 uint64_t tmp = initial_power_of_2 + (D < 0); |
| const uint64_t abs_test_numer = tmp - 1 - tmp % abs_d; |
| |
| /* Initialize our quotients and remainders (q1, r1, q2, r2 in Warren) */ |
| uint64_t quotient1 = initial_power_of_2 / abs_test_numer; |
| uint64_t remainder1 = initial_power_of_2 % abs_test_numer; |
| uint64_t quotient2 = initial_power_of_2 / abs_d; |
| uint64_t remainder2 = initial_power_of_2 % abs_d; |
| uint64_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 = sign_extend(quotient2 + 1, SINT_BITS); |
| if (D < 0) result.multiplier = -result.multiplier; |
| result.shift = exponent - SINT_BITS; |
| return result; |
| } |