Function bounded_lemire_uint64

numpy/random/src/distributions/distributions.c:1190–1271

Static `Lemire rejection` function called by random_bounded_uint64(...) */

/* Static `Lemire rejection` function called by random_bounded_uint64(...) */
static inline uint64_t bounded_lemire_uint64(bitgen_t *bitgen_state,
                                             uint64_t rng) {
  /*
   * Uses Lemire's algorithm - https://arxiv.org/abs/1805.10941
   *
   * Note: `rng` should not be 0xFFFFFFFFFFFFFFFF. When this happens `rng_excl`
   * becomes zero.
   */
  const uint64_t rng_excl = rng + 1;

  assert(rng != 0xFFFFFFFFFFFFFFFFULL);

#if __SIZEOF_INT128__
  /* 128-bit uint available (e.g. GCC/clang). `m` is the __uint128_t scaled
   * integer. */
  __uint128_t m;
  uint64_t leftover;

  /* Generate a scaled random number. */
  m = ((__uint128_t)next_uint64(bitgen_state)) * rng_excl;

  /* Rejection sampling to remove any bias. */
  leftover = m & 0xFFFFFFFFFFFFFFFFULL;

  if (leftover < rng_excl) {
    /* `rng_excl` is a simple upper bound for `threshold`. */
    const uint64_t threshold = (UINT64_MAX - rng) % rng_excl;

    while (leftover < threshold) {
      m = ((__uint128_t)next_uint64(bitgen_state)) * rng_excl;
      leftover = m & 0xFFFFFFFFFFFFFFFFULL;
    }
  }

  return (m >> 64);
#else
  /* 128-bit uint NOT available (e.g. MSVS). `m1` is the upper 64-bits of the
   * scaled integer. */
  uint64_t m1;
  uint64_t x;
  uint64_t leftover;

  x = next_uint64(bitgen_state);

  /* Rejection sampling to remove any bias. */
  leftover = x * rng_excl; /* The lower 64-bits of the mult. */

  if (leftover < rng_excl) {
    /* `rng_excl` is a simple upper bound for `threshold`. */
    const uint64_t threshold = (UINT64_MAX - rng) % rng_excl;

    while (leftover < threshold) {
      x = next_uint64(bitgen_state);
      leftover = x * rng_excl;
    }
  }

#if defined(_MSC_VER) && defined(_WIN64)