path: root/src/util/fast_idiv_by_const.c

/*
 * 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;
}