Generating unsigned, bounded random value using signed bounded random values

Question

Use case: I have a random number generator that can generate random numbers in a range that is limited by a number $x$, say $x$ ranges from $2$ to $2^{(n - 1)} - 1$ exclusive. In this case $n$ is $16$, i.e. we can generate a number from $0$ to $32767$ exclusive. Programmers will directly identify the upper bound to be a possitive, signed value of 16 bits.

Now I want to generate a bounded number from $0$ to $y$, where $y$ ranges from $2$ to $2^n - 1$, i.e. we can generate a number from $0$ to $65535$ exclusive. Yes, this is an unsigned value of 16 bits. Is there an algorithm that can do this efficiently? Loops are allowed; the random in the range $0$ to $x$ may be generated repeatedly.

If this is not possible, I have no problem allowing $0$ to $65533$ exclusive because $32766 * 2 = 65532$ and generating $65533$ and $65534$ could significantly increase the complexity.

All generated random values must have a uniform distribution.

Answer 1

I assume $R(x)$ is the original generator, returning $r$ uniformly distributed with $0\le r<x$ for $x<2^{n-1}$, as does Java's int nextInt(int) for $n=32$; and we want to extend that to $R'(y)$, returning $r$ uniformly distributed with $0\le r<y$ for $y<2^n$. $R$ and $R'$ should treat an argument less than $2$ in the same way (perhaps accept it and return $0$, or reject if somehow). Out of laziness I'll only handle $n$ even and $n>2$, which covers practical use cases.

If $y<2^{n-1}$, $R'(y)=R(x)$ is fine.

If $y=a\cdot b$ with $a>1$ and $b>1$, $R'(y)=b\cdot R(a)+R(b)$ is fine. That is useful when $y$ is known at compilation time and composite.

Else (when $y\ge2^{n-1}$ and prime), there is no method guaranteed to work with a finite number of invocations of $R$, and we must use some iterative method.

---

Out of my head (not tested), a possible line of action when $y$ is known only at runtime:

• if $y<2^{n-1}$, return $R(y)$.
• if $y$ is even, return $2\cdot R(y/2)+R(2)$.
• if $y=2^n-1$, return $(2^{n/2}-1)\cdot R(2^{n/2}+1)+R(2^{n/2}-1)$.
• repeat
  • let $r:=R((y+1)/2)$
  • if $R(2)\ne0$, return $r$.
  • if $r\ne0$, return $y-r$.

Notes:

1. With the assumption that $n$ is even, we can handle $y=2^n-1$ as $(2^{n/2}+1)\cdot(2^{n/2}-1)$.
2. When $R((y+1)/2)$ is computed, this is valid since $1<(y+1)/2<2^{n-2}$.
3. At each iteration of the loop (which occurs only for odd $y$), the value returned (if any) is correctly distributed because:
   • $r$ is returned with odds $1/2$, and then is uniform on $[0\dots(y-1)/2]$
   • $y-r$ is returned with odds $1/2-1/(y+1)$, and then is uniform on $[(y+1)/2\dots y-1]$
   • the "repeat" loop continues with odds $1/(y+1)$, which is at most $1/2^{n-2}$, and thus low for practical $n$.

---

We can lower the expected number of calls to $R$ (by $46\%$ for $n=16$ and $2^{n-1}\le y<2^n-1$) by accumulating $n-2$ uniform bits (and a sentinel) in a variable $b$ surviving across invocations, set to $0$ initially (in the C programming language, that could be static short b = 0; in the body of a function). It is however at the expense of increased complexity, and perhaps of thread-safety, and thus a bad idea in most practical situations! The algorithm (not tested) can become:

• if $y<2^{n-1}$, return $R(y)$.
• if $y=2^n-1$, return $(2^{n/2}-1)\cdot R(2^{n/2}+1)+R(2^{n/2}-1)$.
• repeat
  • let $b:=\begin{cases}R(2^{n-2})+2^{n-2}& \text{if }b<4\\\lfloor b/2\rfloor&\text{otherwise}\end{cases}$
  • let $r:=R(\lfloor(y+1)/2\rfloor)$
  • if $b$ is even, return $r$.
  • if $y$ is even, return $y-1-r$.
  • if $r\ne0$, return $y-r$.