<p>I think I can show the optimality of a scheme that I skteched in the question, the one that does rejection sampling but reuses the remaining entropy.</p>

<p><strong><em>Description of the scheme.</em></strong> Thinking as an algorithm, the scheme is the following. It has two parameters, an integer $n$, and a bound $m$, $0 \leq n &lt; m$, with the invariant that $n$ was drawn uniformly at random in $[0, m[$. (I'll rephrase everything to ranges of the form [0, n[ rather than $\{1, \ldots, n\}$ because it's simpler.)</p>

<p>At each step we do the following:</p>

<ul>
<li><p>if $m &lt; q$, then we need more entropy: let $n'$ be the result of calling the oracle that samples in $[0, p[$ uniformly, and recurse with $n := n p + n'$ (i.e., a combination of the current $n$ and the new draw) and $m := m p$ (a larger interval). It is clear that if $n$ is uniform in $[0, m[$ and $n'$ is uniform in $[0, p[$ then the combination $n p + n'$ is uniform in $[0, mp[$</p></li>
<li><p>if $m \geq q$, then let $M$ be the largest multiple of $q$ which is smaller than $m$. If $n$ is in $[0, M[$, then we return $n \text{ mod } q$ as the result: this is uniform. Otherwise, the remainder $n - M$ is uniform in the range $[0, m - M[$, so we recurse with $n := n - M$ and $q := q - M$.</p></li>
</ul>

<p>In pseudocode, reformulating the second point to substract $q$ multiple times instead of finding $M$, we get the following very simple algorithm:</p>

<pre><code>int algo(int n, int m) {
  // invariant: n is uniform in the range [0, m[
  if (m &lt; q)
    return algo(p*n + randint(p), m * p);  
  if (n &lt; q)
    return n;
  return algo(n-q, m-q);
}

algo(0, 1);
</code></pre>

<p><strong><em>Optimality.</em></strong> This algorithm is, I think, optimal. To see why, let me reformulate it in terms of decision trees. What the algorithm does is, at every level, if we have at least $q$ leaves then turn into leaves the largest multiple of $q$ that you can, and grow the other nodes; otherwise just grow all nodes. (This makes it clear that the scheme is uniform, because the labels in $[0, q[$ on the tree are used in exactly the same way.)</p>

<p>Let us consider the (generally infinite) tree $T$ obtained by such a scheme, and consider a hypothetical tree $T'$ which denotes a uniform strategy with a better expected score.</p>

<p>It is clear that the performance of the tree depends only on the number of internal nodes at every depth. In other words, if I can show that, at every level, the number of internal nodes of $T'$ is at least that of $T$, then $T'$ cannot be better than $T$.</p>

<p>But consider any depth $k$. It is clear that, in $T$, the number of internal nodes is precisely $p^k - M$ for $M$ the largest multiple of $q$ which is below $p^k$: let's write $M = mq$. Now if $T'$ does better than this and has less than this number of internal nodes, let us divide the $p^k$ possible outcomes at depth $k$ of $T'$ between the decided ones and the yet undecided ones: more than $mq$ must be decided (to be better than $T$) and they must be split between $q$ outcomes, so the most probable outcome $v$ at that level must have probability at least $(m+1)/p^k$. Now by maximality of $m$ we have $(m+1)q &gt; p^k$, so that $(m+1)/p^k &gt; 1/q$. But by uniformity the probability of $v$ overall (so, in particular, at depth $k$) should not be more than $1/q$. So $T'$ cannot be uniform: outcome $v$ has probability $&gt; 1/q$.</p>

<p>Hence, $T'$ never does better than $T$.</p>

<p><strong><em>Expected performance.</em></strong> To compute the expected performance of this optimal strategy, one can just write out an equation system based on which recursive calls are performed by the pseudocode above (with the indicated probability distribution). The value of the parameters is bounded, so presumably, for given values of $p$ and $q$, this should be solvable as a geometric series or something like that. (This can be interpreted as an equation system describing the expected number of draws required for a decision forest with $n$ root nodes, $n$ being the number of nodes at the current level.)</p>

<p><strong><em>Extending to multiple draws.</em></strong> One can generalize this scheme to draw multiple values, in the following way: when you have made a draw in $[0, q[$, instead of starting at a single-root decision tree for the next draw, start at a forest with $m$ nodes, the number of the outcomes with the same label at the level at which you decided, to use this remaining entropy. My conjecture would be that the resulting scheme would have the same performance as the $p$- and $q$-ary solution sketched in my question (so be optimal for multiple draws while concluding faster for a single draw). I have no proof of this fact, though.</p>

CC BY-SA 3.0
2015-10-08T10:45:43.687Z
2015-10-08T10:45:43.687Z
Antoine Amarilli 'a3nm'
<p>Let $p$ and $q$ be two positive integers. I have an oracle that can generate a uniform integer in $\{1, \ldots, p\}$, the integers thus produced being independent across oracle calls. My goal is to generate a random number in $\{1, \ldots, q\}$ with uniform probability, and to do so in a way that minimizes the expected number of oracle calls required.</p>

<p>My question is: <strong>What is an optimal algorithm to solve this problem?</strong> I expect this to be a standard problem, so I'd be interested in pointers to the name of this problem in the literature. A special case of this problem ($p = 5$ and $q = 7$) is a standard puzzle or interview question (see <a href="https://stackoverflow.com/questions/137783/expand-a-random-range-from-1-5-to-1-7/842343">here</a>, but I was unable to find pointers to a generalization, and I don't see any optimality proof in those answers. </p>

<p>Of course, in general, the problem will not be solvable with a bounded number of oracle calls. Take for instance $p = 2$ and $q = 3$, or more generally, cases where some prime divisor of $q$ is not a divisor of $p$, and it is clear that for any bound $N$, the number of possible outcomes when doing $N$ calls, i.e,. $p^N$, cannot be divided evenly in $q$ buckets for each possible value.</p>

<p>One natural approach outlined <a href="https://stackoverflow.com/a/891304/414272">here</a> is to see the oracle calls as generating a uniform real number in [0, 1] from its $p$-ary writing, and terminate whenever the first decimal in the $q$-ary writing of that number is certain: think of it as using the oracle to select a subdivision of size $1/p^k$ in the interval $[0, 1]$, and concluding whenever the interval that you fall in is included in one of the $q$ intervals $[i/q, (i+1)/q]$ for $0 \leq i &lt; q$. However, this approach does not seem optimal: for instance, for $p = q+1$, the approach only concludes in one call if the call returns 1, otherwise a second call is necessary. By contrast, a simpler rejection-based approach that calls the oracle, concludes if it falls in $\{1, \ldots, q\}$, and rejects and tries again if it falls in $q+1$, seems to have better performance.</p>

<p>One can try to generalize the rejection-based approach to the following: letting $k$ be minimal such that $p^k \geq q$, use $k$ oracle calls to get a uniform integer $n$ in $\{1, \ldots, p^k\}$. If $n$ falls in $\{1, \ldots, M\}$, for $M &gt; 0$ the smallest multiple of $q$ which is $\leq p^k$, return $n \text{ mod } q$, otherwise try again. However it would seem that the additional information of <em>which</em> excess value in $\{M+1, \ldots, p^k\}$ was obtained could be also used in the next attempt. I can imagine how to do this (see <a href="http://paste.fulltxt.net/1Kze5y" rel="nofollow noreferrer">here</a> for a proposed algorithm, which I think is correct), but I have no idea of how to prove the optimality of such a scheme.</p>

<p>I am also curious about the natural variation of this problem where you want to generate not a single random number in $\{1, \ldots, q\}$, but a stream of such numbers. What's the expected number of calls per number, and an optimal method in this case?</p>

<p>[Note that an alternative vision of the problem that does not involve randomness at all is to think of the outcome trees: each node is either a leaf with label in $\{1, \ldots, q\}$ corresponding to a decision, or an internal node with $p$ children corresponding to an oracle call. An algorithm is a (generally infinite) such tree where, for each $1 \leq l \leq q$, the <em>probability mass</em> of all the leaves labeled $l$ (i.e., the sum of the mass over leaves, where the mass of each leaf is $1/p^h$, where $h$ is its height) is equal to $1/q$; and its expected performance is the sum, over all leaves, of $h/p^h$ (probability mass of this leaf, times the number of oracle calls performed in that case).]</p>

CC BY-SA 3.0
2015-10-07T13:15:31.257Z
did you see this? http://www.prismmodelchecker.org/bibitem.php?key=KY76 here is why i think it is related http://math.stackexchange.com/a/1357/512
Radu GRIGore
CC BY-SA 3.0
2015-10-07T13:22:44.717Z
The math.SE thread is indeed related (it's the case $p=2$, $q=3$). The Knuth and Yao paper seems related, but it would then be specific to the case $p=2$? (and much more general in this case). It doesn't help that I can't locate a copy of that paper online. Thanks a lot for the suggestion!
Antoine Amarilli 'a3nm'
CC BY-SA 3.0
2015-10-07T14:15:39.150Z
One related reference (for a more general problem, where the input and output distribution are arbitrary): "Interval algorithm for random number generation", Hao and Hoshi. Their algorithm looks sort of like arithmetic coding.
usul
CC BY-SA 3.0
2015-10-07T19:50:07.400Z
That link you gave is for generating a stream of numbers, not a single value.  It's basically arithmetic encoding base q, so its optimality is hard to refute.
KWillets
CC BY-SA 3.0
2015-10-08T17:48:42.407Z
The product of k uniform random variables is not uniformly distributed.
KWillets
CC BY-SA 3.0
2015-10-08T21:44:23.827Z
@KWillets: sure, the product of uniform variables isn't uniform, but I'm not sure why do you make this point?
Antoine Amarilli 'a3nm'
CC BY-SA 3.0
2015-10-08T21:52:50.807Z
Ah, I misunderstood how you're using the k oracle calls.  You're concatenating the digits base p I think.
KWillets
CC BY-SA 3.0
2015-10-09T08:52:37.857Z
@KWillets: Yes, that is precisely what I do. (Or, in the scheme in my answer, I'm also using values that are not really in base $p$ as their bound is not $p$ but $p^k - M$.)
Antoine Amarilli 'a3nm'
CC BY-SA 3.0
2015-10-13T01:41:37.710Z
It seems this or similar problems has been approached on a number of SE sites. Here's another from [cs.SE](http://cs.stackexchange.com/q/2576/528).
Mark Hurd
CC BY-SA 3.0
2015-10-14T11:00:54.993Z
@a3nm: You can find a Haskell implementation of rejection sampling (and reusing remaining entropy) here: 
https://hackage.haskell.org/package/dice-entropy-conduit-1.0.0.1/docs/System-Random-Dice.html
Peter
CC BY-SA 3.0
2015-10-14T20:34:10.893Z
@Peter: Brilliant! It looks like the dRoll function is indeed doing exactly the algorithm I was sketching in my answer, plus the trick of reusing entropy when generating multiple values. Did you come up with it independently, or is it from some existing source? Do you know anything about optimality? Anyway, thanks for pointing me to this!
Antoine Amarilli 'a3nm'
<reference-request><randomness><pseudorandom-generators>
Generating uniform integers in a range from a random generator with another range
Antoine Amarilli 'a3nm'
# Generating uniform integers in a range from a random generator with another range

(Asked by: username_0 on 2015-10-07)

Let $p$ and $q$ be two positive integers. I have an oracle that can generate a uniform integer in $\{1, \ldots, p\}$, the integers thus produced being independent across oracle calls. My goal is to generate a random number in $\{1, \ldots, q\}$ with uniform probability, and to do so in a way that minimizes the expected number of oracle calls required.

My question is: **What is an optimal algorithm to solve this problem?** I expect this to be a standard problem, so I'd be interested in pointers to the name of this problem in the literature. A special case of this problem ($p = 5$ and $q = 7$) is a standard puzzle or interview question (see [here](https://stackoverflow.com/questions/137783/expand-a-random-range-from-1-5-to-1-7/842343), but I was unable to find pointers to a generalization, and I don't see any optimality proof in those answers.

Of course, in general, the problem will not be solvable with a bounded number of oracle calls. Take for instance $p = 2$ and $q = 3$, or more generally, cases where some prime divisor of $q$ is not a divisor of $p$, and it is clear that for any bound $N$, the number of possible outcomes when doing $N$ calls, i.e,. $p^N$, cannot be divided evenly in $q$ buckets for each possible value.

One natural approach outlined [here](https://stackoverflow.com/a/891304/414272) is to see the oracle calls as generating a uniform real number in [0, 1] from its $p$-ary writing, and terminate whenever the first decimal in the $q$-ary writing of that number is certain: think of it as using the oracle to select a subdivision of size $1/p^k$ in the interval $[0, 1]$, and concluding whenever the interval that you fall in is included in one of the $q$ intervals $[i/q, (i+1)/q]$ for $0 \leq i < q$. However, this approach does not seem optimal: for instance, for $p = q+1$, the approach only concludes in one call if the call returns 1, otherwise a second call is necessary. By contrast, a simpler rejection-based approach that calls the oracle, concludes if it falls in $\{1, \ldots, q\}$, and rejects and tries again if it falls in $q+1$, seems to have better performance.

One can try to generalize the rejection-based approach to the following: letting $k$ be minimal such that $p^k \geq q$, use $k$ oracle calls to get a uniform integer $n$ in $\{1, \ldots, p^k\}$. If $n$ falls in $\{1, \ldots, M\}$, for $M > 0$ the smallest multiple of $q$ which is $\leq p^k$, return $n \text{ mod } q$, otherwise try again. However it would seem that the additional information of *which* excess value in $\{M+1, \ldots, p^k\}$ was obtained could be also used in the next attempt. I can imagine how to do this (see [here](http://paste.fulltxt.net/1Kze5y) for a proposed algorithm, which I think is correct), but I have no idea of how to prove the optimality of such a scheme.

I am also curious about the natural variation of this problem where you want to generate not a single random number in $\{1, \ldots, q\}$, but a stream of such numbers. What's the expected number of calls per number, and an optimal method in this case?

[Note that an alternative vision of the problem that does not involve randomness at all is to think of the outcome trees: each node is either a leaf with label in $\{1, \ldots, q\}$ corresponding to a decision, or an internal node with $p$ children corresponding to an oracle call. An algorithm is a (generally infinite) such tree where, for each $1 \leq l \leq q$, the *probability mass* of all the leaves labeled $l$ (i.e., the sum of the mass over leaves, where the mass of each leaf is $1/p^h$, where $h$ is its height) is equal to $1/q$; and its expected performance is the sum, over all leaves, of $h/p^h$ (probability mass of this leaf, times the number of oracle calls performed in that case).]

**Tags:** reference-request, randomness, pseudorandom-generators
> **Comments:**
> **username_1**: did you see this? http://www.prismmodelchecker.org/bibitem.php?key=KY76 here is why i think it is related http://math.stackexchange.com/a/1357/512  
> **username_0**: The math.SE thread is indeed related (it's the case $p=2$, $q=3$). The Knuth and Yao paper seems related, but it would then be specific to the case $p=2$? (and much more general in this case). It doesn't help that I can't locate a copy of that paper online. Thanks a lot for the suggestion!  
> **username_2**: One related reference (for a more general problem, where the input and output distribution are arbitrary): "Interval algorithm for random number generation", Hao and Hoshi. Their algorithm looks sort of like arithmetic coding.  
> **username_3**: That link you gave is for generating a stream of numbers, not a single value.  It's basically arithmetic encoding base q, so its optimality is hard to refute.  
> **username_3**: The product of k uniform random variables is not uniformly distributed.  
> **username_0**: @KWillets: sure, the product of uniform variables isn't uniform, but I'm not sure why do you make this point?  
> **username_3**: Ah, I misunderstood how you're using the k oracle calls.  You're concatenating the digits base p I think.  
> **username_0**: @KWillets: Yes, that is precisely what I do. (Or, in the scheme in my answer, I'm also using values that are not really in base $p$ as their bound is not $p$ but $p^k - M$.)  
> **username_4**: It seems this or similar problems has been approached on a number of SE sites. Here's another from [cs.SE](http://cs.stackexchange.com/q/2576/528).  
> **username_5**: @a3nm: You can find a Haskell implementation of rejection sampling (and reusing remaining entropy) here:  https://hackage.haskell.org/package/dice-entropy-conduit-1.0.0.1/docs/System-Random-Dice.html  
> **username_0**: @Peter: Brilliant! It looks like the dRoll function is indeed doing exactly the algorithm I was sketching in my answer, plus the trick of reusing entropy when generating multiple values. Did you come up with it independently, or is it from some existing source? Do you know anything about optimality? Anyway, thanks for pointing me to this!  

---

## 1 answers

### Answer by username_0 on 2015-10-08. Score: 2 | Accepted answer

I think I can show the optimality of a scheme that I skteched in the question, the one that does rejection sampling but reuses the remaining entropy.
Description of the scheme. Thinking as an algorithm, the scheme is the following. It has two parameters, an integer $n$, and a bound $m$, $0 \leq n < m$, with the invariant that $n$ was drawn uniformly at random in $[0, m[$. (I'll rephrase everything to ranges of the form [0, n[ rather than $\{1, \ldots, n\}$ because it's simpler.)
At each step we do the following:

if $m < q$, then we need more entropy: let $n'$ be the result of calling the oracle that samples in $[0, p[$ uniformly, and recurse with $n := n p + n'$ (i.e., a combination of the current $n$ and the new draw) and $m := m p$ (a larger interval). It is clear that if $n$ is uniform in $[0, m[$ and $n'$ is uniform in $[0, p[$ then the combination $n p + n'$ is uniform in $[0, mp[$
if $m \geq q$, then let $M$ be the largest multiple of $q$ which is smaller than $m$. If $n$ is in $[0, M[$, then we return $n \text{ mod } q$ as the result: this is uniform. Otherwise, the remainder $n - M$ is uniform in the range $[0, m - M[$, so we recurse with $n := n - M$ and $q := q - M$.

In pseudocode, reformulating the second point to substract $q$ multiple times instead of finding $M$, we get the following very simple algorithm:
int algo(int n, int m) {
  // invariant: n is uniform in the range [0, m[
  if (m < q)
    return algo(p*n + randint(p), m * p);  
  if (n < q)
    return n;
  return algo(n-q, m-q);
}

algo(0, 1);

Optimality. This algorithm is, I think, optimal. To see why, let me reformulate it in terms of decision trees. What the algorithm does is, at every level, if we have at least $q$ leaves then turn into leaves the largest multiple of $q$ that you can, and grow the other nodes; otherwise just grow all nodes. (This makes it clear that the scheme is uniform, because the labels in $[0, q[$ on the tree are used in exactly the same way.)
Let us consider the (generally infinite) tree $T$ obtained by such a scheme, and consider a hypothetical tree $T'$ which denotes a uniform strategy with a better expected score.
It is clear that the performance of the tree depends only on the number of internal nodes at every depth. In other words, if I can show that, at every level, the number of internal nodes of $T'$ is at least that of $T$, then $T'$ cannot be better than $T$.
But consider any depth $k$. It is clear that, in $T$, the number of internal nodes is precisely $p^k - M$ for $M$ the largest multiple of $q$ which is below $p^k$: let's write $M = mq$. Now if $T'$ does better than this and has less than this number of internal nodes, let us divide the $p^k$ possible outcomes at depth $k$ of $T'$ between the decided ones and the yet undecided ones: more than $mq$ must be decided (to be better than $T$) and they must be split between $q$ outcomes, so the most probable outcome $v$ at that level must have probability at least $(m+1)/p^k$. Now by maximality of $m$ we have $(m+1)q > p^k$, so that $(m+1)/p^k > 1/q$. But by uniformity the probability of $v$ overall (so, in particular, at depth $k$) should not be more than $1/q$. So $T'$ cannot be uniform: outcome $v$ has probability $> 1/q$.
Hence, $T'$ never does better than $T$.
Expected performance. To compute the expected performance of this optimal strategy, one can just write out an equation system based on which recursive calls are performed by the pseudocode above (with the indicated probability distribution). The value of the parameters is bounded, so presumably, for given values of $p$ and $q$, this should be solvable as a geometric series or something like that. (This can be interpreted as an equation system describing the expected number of draws required for a decision forest with $n$ root nodes, $n$ being the number of nodes at the current level.)
Extending to multiple draws. One can generalize this scheme to draw multiple values, in the following way: when you have made a draw in $[0, q[$, instead of starting at a single-root decision tree for the next draw, start at a forest with $m$ nodes, the number of the outcomes with the same label at the level at which you decided, to use this remaining entropy. My conjecture would be that the resulting scheme would have the same performance as the $p$- and $q$-ary solution sketched in my question (so be optimal for multiple draws while concluding faster for a single draw). I have no proof of this fact, though.


---

