<p>Without loss of generality, the problem of generating random integers on [a, b] can be reduced to the problem of generating random integers on [0, s). The state of the art for generating random integers on a bounded range from a uniform PRNG is represented by the following recent publication:</p>

<p>Daniel Lemire,"Fast Random Integer Generation in an Interval." <em>ACM Trans. Model. Comput. Simul.</em> 29, 1, Article 3 (January 2019) (<a href="https://arxiv.org/abs/1805.10941" rel="nofollow noreferrer">ArXiv draft</a>)</p>

<p>Lemire shows that his algorithm provides unbiased results, and motivated by the growing popularity of very fast high-quality PRNGs such as Melissa O'Neill's <a href="http://www.pcg-random.org/" rel="nofollow noreferrer">PCG generators</a>, shows how to the results can be computed fast, avoiding slow division operations almost all of the time. </p>

<p>An exemplary ISO-C implementation of his algorithm is shown in <code>randint()</code> below. Here I demonstrate it in conjunction with George Marsaglia's older <a href="https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J" rel="nofollow noreferrer">KISS64</a> PRNG. For performance reasons, the required 64×64→128 bit unsigned multiplication is typically best implemented via machine-specific intrinsics or inline assembly that map directly to appropriate hardware instructions.</p>

<pre class="lang-c prettyprint-override"><code>#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;stdint.h&gt;

/* PRNG state */
typedef struct Prng_T *Prng_T;
/* Returns uniformly distributed integers in [0, 2**64-1] */
uint64_t random64 (Prng_T);
/* Multiplies two 64-bit factors into a 128-bit product */
void umul64wide (uint64_t, uint64_t, uint64_t *, uint64_t *);

/* Generate in bias-free manner a random integer in [0, s) with Lemire's fast
   algorithm that uses integer division only rarely. s must be in [0, 2**64-1].

   Daniel Lemire, "Fast Random Integer Generation in an Interval," ACM Trans.
   Model. Comput. Simul. 29, 1, Article 3 (January 2019)
*/
uint64_t randint (Prng_T prng, uint64_t s) 
{
    uint64_t x, h, l, t;
    x = random64 (prng);
    umul64wide (x, s, &amp;h, &amp;l);
    if (l &lt; s) {
        t = (0 - s) % s;
        while (l &lt; t) {
            x = random64 (prng);
            umul64wide (x, s, &amp;h, &amp;l);
        }
    }
    return h;
}

#define X86_INLINE_ASM (0)

/* Multiply two 64-bit unsigned integers into a 128 bit unsined product. Return
   the least significant 64 bist of the product to the location pointed to by
   lo, and the most signfiicant 64 bits of the product to the location pointed
   to by hi.
*/
void umul64wide (uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo)
{
#if X86_INLINE_ASM
    uint64_t l, h;
    __asm__ (
        "movq  %2, %%rax;\n\t"  // rax = a
        "mulq  %3;\n\t"         // rdx:rax = a * b
        "movq  %%rax, %0;\n\t"  // l = (a * b)&lt;31:0&gt;
        "movq  %%rdx, %1;\n\t"  // h = (a * b)&lt;63:32&gt;
        : "=r"(l), "=r"(h)
        : "r"(a), "r"(b)
        : "%rax", "%rdx");
    *lo = l;
    *hi = h;
#else // X86_INLINE_ASM
    uint64_t a_lo = (uint64_t)(uint32_t)a;
    uint64_t a_hi = a &gt;&gt; 32;
    uint64_t b_lo = (uint64_t)(uint32_t)b;
    uint64_t b_hi = b &gt;&gt; 32;

    uint64_t p0 = a_lo * b_lo;
    uint64_t p1 = a_lo * b_hi;
    uint64_t p2 = a_hi * b_lo;
    uint64_t p3 = a_hi * b_hi;

    uint32_t cy = (uint32_t)(((p0 &gt;&gt; 32) + (uint32_t)p1 + (uint32_t)p2) &gt;&gt; 32);

    *lo = p0 + (p1 &lt;&lt; 32) + (p2 &lt;&lt; 32);
    *hi = p3 + (p1 &gt;&gt; 32) + (p2 &gt;&gt; 32) + cy;
#endif // X86_INLINE_ASM
}

/* George Marsaglia's KISS64 generator, posted to comp.lang.c on 28 Feb 2009
   https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
*/
struct Prng_T {
    uint64_t x, c, y, z, t;
};

struct Prng_T kiss64 = {1234567890987654321ULL, 123456123456123456ULL,
                        362436362436362436ULL, 1066149217761810ULL, 0ULL};

/* KISS64 state equations */
#define MWC64 (kiss64-&gt;t = (kiss64-&gt;x &lt;&lt; 58) + kiss64-&gt;c,            \
               kiss64-&gt;c = (kiss64-&gt;x &gt;&gt; 6), kiss64-&gt;x += kiss64-&gt;t, \
               kiss64-&gt;c += (kiss64-&gt;x &lt; kiss64-&gt;t), kiss64-&gt;x)
#define XSH64 (kiss64-&gt;y ^= (kiss64-&gt;y &lt;&lt; 13), kiss64-&gt;y ^= (kiss64-&gt;y &gt;&gt; 17), \
               kiss64-&gt;y ^= (kiss64-&gt;y &lt;&lt; 43))
#define CNG64 (kiss64-&gt;z = 6906969069ULL * kiss64-&gt;z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
uint64_t random64 (Prng_T kiss64)
{
    return KISS64;
}

int main (void)
{
    int i;
    Prng_T state = &amp;kiss64;

    for (i = 0; i &lt; 1000; i++) {
        printf ("%llu\n", randint (state, 10));
    }
    return EXIT_SUCCESS;
}
</code></pre>
What is the optimal algorithm for generating an unbiased random integer within a range?
<p>In this StackOverflow question:</p>

<p><a href="https://stackoverflow.com/questions/5008804/generating-random-integer-from-a-range">Generating random integer from a range</a></p>

<p>the accepted answer suggests the following formula for generating a random integer in between given <code>min</code> and <code>max</code>, with <code>min</code> and <code>max</code> being included into the range:</p>

<pre><code>output = min + (rand() % (int)(max - min + 1))
</code></pre>

<p>But it also says that</p>

<blockquote>
  <p>This is still <em>slightly</em> biased towards lower numbers ... It's also
  possible to extend it so that it removes the bias.</p>
</blockquote>

<p>But it doesn't explain why it's biased towards lower numbers or how to remove the bias. So, the question is: is this the most optimal approach to generation of a random integer within a (signed) range while not relying on anything fancy, just <code>rand()</code> function, and in case if it is optimal, how to remove the bias?</p>

<p><strong>EDIT:</strong></p>

<p>I've just tested the <code>while</code>-loop algorithm suggested by @Joey against floating-point extrapolation:</p>

<pre><code>static const double s_invRandMax = 1.0/((double)RAND_MAX + 1.0);
return min + (int)(((double)(max + 1 - min))*rand()*s_invRandMax);
</code></pre>

<p>to see how much uniformly "balls" are "falling" into and are being distributed among a number of "buckets", one test for the floating-point extrapolation and another for the <code>while</code>-loop algorithm. But results turned out to be varying depending on the number of "balls" (and "buckets") so I couldn't easily pick a winner. The working code can be found at <a href="http://ideone.com/oPhqK" rel="noreferrer">this Ideone page</a>. For example, with 10 buckets and 100 balls the maximum deviation from the ideal probability among buckets is less for the floating-point extrapolation than for the <code>while</code>-loop algorithm (0.04 and 0.05 respectively) but with 1000 balls, the maximum deviation of the <code>while</code>-loop algorithm is lesser (0.024 and 0.011), and with 10000 balls, the floating-point extrapolation is again doing better (0.0034 and 0.0053), and so on without much of consistency. Thinking of the possibility that none of the algorithms consistently produces uniform distribution better than that of the other algorithm, makes me lean towards the floating-point extrapolation since it appears to perform faster than the <code>while</code>-loop algorithm. So is it fine to choose the floating-point extrapolation algorithm or my testings/conclusions are not completely correct?</p>