src/util/fast_idiv_by_const.c

   1 /*
   2  * Copyright © 2018 Advanced Micro Devices, Inc.
   3  *
   4  * Permission is hereby granted, free of charge, to any person obtaining a
   5  * copy of this software and associated documentation files (the "Software"),
   6  * to deal in the Software without restriction, including without limitation
   7  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
   8  * and/or sell copies of the Software, and to permit persons to whom the
   9  * Software is furnished to do so, subject to the following conditions:
  10  *
  11  * The above copyright notice and this permission notice (including the next
  12  * paragraph) shall be included in all copies or substantial portions of the
  13  * Software.
  14  *
  15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  16  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  17  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  18  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  19  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  20  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
  21  * IN THE SOFTWARE.
  22  */
  23
  24 /* Imported from:
  25  *   https://raw.githubusercontent.com/ridiculousfish/libdivide/master/divide_by_constants_codegen_reference.c
  26  * Paper:
  27  *   http://ridiculousfish.com/files/faster_unsigned_division_by_constants.pdf
  28  *
  29  * The author, ridiculous_fish, wrote:
  30  *
  31  *  ''Reference implementations of computing and using the "magic number"
  32  *    approach to dividing by constants, including codegen instructions.
  33  *    The unsigned division incorporates the "round down" optimization per
  34  *    ridiculous_fish.
  35  *
  36  *    This is free and unencumbered software. Any copyright is dedicated
  37  *    to the Public Domain.''
  38  */
  39
  40 #include "fast_idiv_by_const.h"
  41 #include "u_math.h"
  42 #include <limits.h>
  43 #include <assert.h>
  44
  45 struct util_fast_udiv_info
  46 util_compute_fast_udiv_info(uint64_t D, unsigned num_bits, unsigned UINT_BITS)
  47 {
  48    /* The numerator must fit in a uint64_t */
  49    assert(num_bits > 0 && num_bits <= UINT_BITS);
  50    assert(D != 0);
  51
  52    /* The eventual result */
  53    struct util_fast_udiv_info result;
  54
  55    if (util_is_power_of_two_or_zero64(D)) {
  56       unsigned div_shift = util_logbase2_64(D);
  57
  58       if (div_shift) {
  59          /* Dividing by a power of two. */
  60          result.multiplier = 1ull << (UINT_BITS - div_shift);
  61          result.pre_shift = 0;
  62          result.post_shift = 0;
  63          result.increment = 0;
  64          return result;
  65       } else {
  66          /* Dividing by 1. */
  67          /* Assuming: floor((num + 1) * (2^32 - 1) / 2^32) = num */
  68          result.multiplier = UINT_BITS == 64 ? UINT64_MAX :
  69                                                (1ull << UINT_BITS) - 1;
  70          result.pre_shift = 0;
  71          result.post_shift = 0;
  72          result.increment = 1;
  73          return result;
  74       }
  75    }
  76
  77    /* The extra shift implicit in the difference between UINT_BITS and num_bits
  78     */
  79    const unsigned extra_shift = UINT_BITS - num_bits;
  80
  81    /* The initial power of 2 is one less than the first one that can possibly
  82     * work.
  83     */
  84    const uint64_t initial_power_of_2 = (uint64_t)1 << (UINT_BITS-1);
  85
  86    /* The remainder and quotient of our power of 2 divided by d */
  87    uint64_t quotient = initial_power_of_2 / D;
  88    uint64_t remainder = initial_power_of_2 % D;
  89
  90    /* ceil(log_2 D) */
  91    unsigned ceil_log_2_D;
  92
  93    /* The magic info for the variant "round down" algorithm */
  94    uint64_t down_multiplier = 0;
  95    unsigned down_exponent = 0;
  96    int has_magic_down = 0;
  97
  98    /* Compute ceil(log_2 D) */
  99    ceil_log_2_D = 0;
 100    uint64_t tmp;
 101    for (tmp = D; tmp > 0; tmp >>= 1)
 102       ceil_log_2_D += 1;
 103
 104
 105    /* Begin a loop that increments the exponent, until we find a power of 2
 106     * that works.
 107     */
 108    unsigned exponent;
 109    for (exponent = 0; ; exponent++) {
 110       /* Quotient and remainder is from previous exponent; compute it for this
 111        * exponent.
 112        */
 113       if (remainder >= D - remainder) {
 114          /* Doubling remainder will wrap around D */
 115          quotient = quotient * 2 + 1;
 116          remainder = remainder * 2 - D;
 117       } else {
 118          /* Remainder will not wrap */
 119          quotient = quotient * 2;
 120          remainder = remainder * 2;
 121       }
 122
 123       /* We're done if this exponent works for the round_up algorithm.
 124        * Note that exponent may be larger than the maximum shift supported,
 125        * so the check for >= ceil_log_2_D is critical.
 126        */
 127       if ((exponent + extra_shift >= ceil_log_2_D) ||
 128           (D - remainder) <= ((uint64_t)1 << (exponent + extra_shift)))
 129          break;
 130
 131       /* Set magic_down if we have not set it yet and this exponent works for
 132        * the round_down algorithm
 133        */
 134       if (!has_magic_down &&
 135           remainder <= ((uint64_t)1 << (exponent + extra_shift))) {
 136          has_magic_down = 1;
 137          down_multiplier = quotient;
 138          down_exponent = exponent;
 139       }
 140    }
 141
 142    if (exponent < ceil_log_2_D) {
 143       /* magic_up is efficient */
 144       result.multiplier = quotient + 1;
 145       result.pre_shift = 0;
 146       result.post_shift = exponent;
 147       result.increment = 0;
 148    } else if (D & 1) {
 149       /* Odd divisor, so use magic_down, which must have been set */
 150       assert(has_magic_down);
 151       result.multiplier = down_multiplier;
 152       result.pre_shift = 0;
 153       result.post_shift = down_exponent;
 154       result.increment = 1;
 155    } else {
 156       /* Even divisor, so use a prefix-shifted dividend */
 157       unsigned pre_shift = 0;
 158       uint64_t shifted_D = D;
 159       while ((shifted_D & 1) == 0) {
 160          shifted_D >>= 1;
 161          pre_shift += 1;
 162       }
 163       result = util_compute_fast_udiv_info(shifted_D, num_bits - pre_shift,
 164                                            UINT_BITS);
 165       /* expect no increment or pre_shift in this path */
 166       assert(result.increment == 0 && result.pre_shift == 0);
 167       result.pre_shift = pre_shift;
 168    }
 169    return result;
 170 }
 171
 172 static inline int64_t
 173 sign_extend(int64_t x, unsigned SINT_BITS)
 174 {
 175    return (int64_t)((uint64_t)x << (64 - SINT_BITS)) >> (64 - SINT_BITS);
 176 }
 177
 178 struct util_fast_sdiv_info
 179 util_compute_fast_sdiv_info(int64_t D, unsigned SINT_BITS)
 180 {
 181    /* D must not be zero. */
 182    assert(D != 0);
 183    /* The result is not correct for these divisors. */
 184    assert(D != 1 && D != -1);
 185
 186    /* Our result */
 187    struct util_fast_sdiv_info result;
 188
 189    /* Absolute value of D (we know D is not the most negative value since
 190     * that's a power of 2)
 191     */
 192    const uint64_t abs_d = (D < 0 ? -D : D);
 193
 194    /* The initial power of 2 is one less than the first one that can possibly
 195     * work */
 196    /* "two31" in Warren */
 197    unsigned exponent = SINT_BITS - 1;
 198    const uint64_t initial_power_of_2 = (uint64_t)1 << exponent;
 199
 200    /* Compute the absolute value of our "test numerator,"
 201     * which is the largest dividend whose remainder with d is d-1.
 202     * This is called anc in Warren.
 203     */
 204    const uint64_t tmp = initial_power_of_2 + (D < 0);
 205    const uint64_t abs_test_numer = tmp - 1 - tmp % abs_d;
 206
 207    /* Initialize our quotients and remainders (q1, r1, q2, r2 in Warren) */
 208    uint64_t quotient1 = initial_power_of_2 / abs_test_numer;
 209    uint64_t remainder1 = initial_power_of_2 % abs_test_numer;
 210    uint64_t quotient2 = initial_power_of_2 / abs_d;
 211    uint64_t remainder2 = initial_power_of_2 % abs_d;
 212    uint64_t delta;
 213
 214    /* Begin our loop */
 215    do {
 216       /* Update the exponent */
 217       exponent++;
 218
 219       /* Update quotient1 and remainder1 */
 220       quotient1 *= 2;
 221       remainder1 *= 2;
 222       if (remainder1 >= abs_test_numer) {
 223          quotient1 += 1;
 224          remainder1 -= abs_test_numer;
 225       }
 226
 227       /* Update quotient2 and remainder2 */
 228       quotient2 *= 2;
 229       remainder2 *= 2;
 230       if (remainder2 >= abs_d) {
 231          quotient2 += 1;
 232          remainder2 -= abs_d;
 233       }
 234
 235       /* Keep going as long as (2**exponent) / abs_d <= delta */
 236       delta = abs_d - remainder2;
 237    } while (quotient1 < delta || (quotient1 == delta && remainder1 == 0));
 238
 239    result.multiplier = sign_extend(quotient2 + 1, SINT_BITS);
 240    if (D < 0) result.multiplier = -result.multiplier;
 241    result.shift = exponent - SINT_BITS;
 242    return result;
 243 }