src/util/fast_urem_by_const.h

   1 /*
   2  * Copyright © 2010 Valve Software
   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 #include <stdint.h>
  25
  26 /*
  27  * Code for fast 32-bit unsigned remainder, based off of "Faster Remainder by
  28  * Direct Computation: Applications to Compilers and Software Libraries,"
  29  * available at https://arxiv.org/pdf/1902.01961.pdf.
  30  *
  31  * util_fast_urem32(n, d, REMAINDER_MAGIC(d)) returns the same thing as
  32  * n % d for any unsigned n and d, however it compiles down to only a few
  33  * multiplications, so it should be faster than plain uint32_t modulo if the
  34  * same divisor is used many times.
  35  */
  36
  37 #define REMAINDER_MAGIC(divisor) \
  38    ((uint64_t) ~0ull / (divisor) + 1)
  39
  40 /*
  41  * Get bits 64-96 of a 32x64-bit multiply. If __int128_t is available, we use
  42  * it, which usually compiles down to one instruction on 64-bit architectures.
  43  * Otherwise on 32-bit architectures we usually get four instructions (one
  44  * 32x32->64 multiply, one 32x32->32 multiply, and one 64-bit add).
  45  */
  46
  47 static inline uint32_t
  48 _mul32by64_hi(uint32_t a, uint64_t b)
  49 {
  50 #ifdef HAVE_UINT128
  51    return ((__uint128_t) b * a) >> 64;
  52 #else
  53    /*
  54     * Let b = b0 + 2^32 * b1. Then a * b = a * b0 + 2^32 * a * b1. We would
  55     * have to do a 96-bit addition to get the full result, except that only
  56     * one term has non-zero lower 32 bits, which means that to get the high 32
  57     * bits, we only have to add the high 64 bits of each term. Unfortunately,
  58     * we have to do the 64-bit addition in case the low 32 bits overflow.
  59     */
  60    uint32_t b0 = (uint32_t) b;
  61    uint32_t b1 = b >> 32;
  62    return ((((uint64_t) a * b0) >> 32) + (uint64_t) a * b1) >> 32;
  63 #endif
  64 }
  65
  66 static inline uint32_t
  67 util_fast_urem32(uint32_t n, uint32_t d, uint64_t magic)
  68 {
  69    uint64_t lowbits = magic * n;
  70    uint32_t result = _mul32by64_hi(d, lowbits);
  71    assert(result == n % d);
  72    return result;
  73 }
  74