Implement of randint()

BIND

PRNG engine: arcfour, with tick in [0; 0xffff]

Code:

static isc_uint16_t
dispatch_arc4uniformrandom(dns_dispatchmgr_t *mgr, isc_uint16_t upper_bound) {
        isc_uint16_t min, r;
        /* The caller must hold the manager lock. */

        if (upper_bound < 2)
                return (0);

        /*
         * Ensure the range of random numbers [min, 0xffff] be a multiple of
         * upper_bound and contain at least a half of the 16 bit range.
         */

        if (upper_bound > 0x8000)
                min = 1 + ~upper_bound; /* 0x8000 - upper_bound */
        else
                min = (isc_uint16_t)(0x10000 % (isc_uint32_t)upper_bound);

        /*
         * This could theoretically loop forever but each retry has
         * p > 0.5 (worst case, usually far better) of selecting a
         * number inside the range we need, so it should rarely need
         * to re-roll.
         */
        for (;;) {
                r = dispatch_arc4random(mgr);
                if (r >= min)
                        break;
        }

        return (r % upper_bound);
}

GSL

Code:

extern inline unsigned long int
gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n)
{
  unsigned long int offset = r->type->min;
  unsigned long int range = r->type->max - offset;
  unsigned long int scale;
  unsigned long int k;

  if (n > range || n == 0)
    {
      GSL_ERROR_VAL ("invalid n, either 0 or exceeds maximum value of generator",
                     GSL_EINVAL, 0) ;
    }

  scale = range / n;

  do
    {
      k = (((r->type->get) (r->state)) - offset) / scale;
    }
  while (k >= n);

  return k;
}

glib

PRNG engine: Mersenne Twister, with tick in [0; 0xffffffff]

Code:

/**
 * g_rand_int_range:
 * @rand_: a #GRand.
 * @begin: lower closed bound of the interval.
 * @end: upper open bound of the interval.
 *
 * Returns the next random #gint32 from @rand_ equally distributed over
 * the range [@begin..@end-1].
 *
 * Return value: A random number.
 **/
gint32
g_rand_int_range (GRand* rand, gint32 begin, gint32 end)
{
  guint32 dist = end - begin;
  guint32 random;

  g_return_val_if_fail (rand != NULL, begin);
  g_return_val_if_fail (end > begin, begin);

  switch (get_random_version ())
    {
    case 20:
      if (dist <= 0x10000L) /* 2^16 */
        {
          /* This method, which only calls g_rand_int once is only good
           * for (end - begin) <= 2^16, because we only have 32 bits set
           * from the one call to g_rand_int (). */

          /* we are using (trans + trans * trans), because g_rand_int only
           * covers [0..2^32-1] and thus g_rand_int * trans only covers
           * [0..1-2^-32], but the biggest double < 1 is 1-2^-52.
           */

          gdouble double_rand = g_rand_int (rand) *
            (G_RAND_DOUBLE_TRANSFORM +
             G_RAND_DOUBLE_TRANSFORM * G_RAND_DOUBLE_TRANSFORM);

          random = (gint32) (double_rand * dist);
        }
      else
        {
          /* Now we use g_rand_double_range (), which will set 52 bits for
             us, so that it is safe to round and still get a decent
             distribution */
          random = (gint32) g_rand_double_range (rand, 0, dist);
        }
      break;
    case 22:
      if (dist == 0)
        random = 0;
      else
        {
          /* maxvalue is set to the predecessor of the greatest
           * multiple of dist less or equal 2^32. */
          guint32 maxvalue;
          if (dist <= 0x80000000u) /* 2^31 */
            {
              /* maxvalue = 2^32 - 1 - (2^32 % dist) */
              guint32 leftover = (0x80000000u % dist) * 2;
              if (leftover >= dist) leftover -= dist;
              maxvalue = 0xffffffffu - leftover;
            }
          else
            maxvalue = dist - 1;

          do
            random = g_rand_int (rand);
          while (random > maxvalue);

          random %= dist;
        }
      break;
    default:
      random = 0;                /* Quiet GCC */
      g_assert_not_reached ();
    }

  return begin + random;
}

with:

/* transform [0..2^32] -> [0..1] */
#define G_RAND_DOUBLE_TRANSFORM 2.3283064365386962890625e-10

/**
 * g_rand_double:
 * @rand_: a #GRand.
 *
 * Returns the next random #gdouble from @rand_ equally distributed over
 * the range [0..1).
 *
 * Return value: A random number.
 **/
gdouble
g_rand_double (GRand* rand)
{
  /* We set all 52 bits after the point for this, not only the first
     32. Thats why we need two calls to g_rand_int */
  gdouble retval = g_rand_int (rand) * G_RAND_DOUBLE_TRANSFORM;
  retval = (retval + g_rand_int (rand)) * G_RAND_DOUBLE_TRANSFORM;

  /* The following might happen due to very bad rounding luck, but
   * actually this should be more than rare, we just try again then */
  if (retval >= 1.0)
    return g_rand_double (rand);

  return retval;
}

/**
 * g_rand_double_range:
 * @rand_: a #GRand.
 * @begin: lower closed bound of the interval.
 * @end: upper open bound of the interval.
 *
 * Returns the next random #gdouble from @rand_ equally distributed over
 * the range [@begin..@end).
 *
 * Return value: A random number.
 **/
gdouble
g_rand_double_range (GRand* rand, gdouble begin, gdouble end)
{
  return g_rand_double (rand) * (end - begin) + begin;
}