Random floating point double in Inclusive Range

We can easily get random floating point numbers within a desired range [X,Y) (note that X is inclusive and Y is exclusive) with the function listed below since Math.random() (and most pseudorandom number generators, AFAIK) produce numbers in [0,1):

function randomInRange(min, max) {
  return Math.random() * (max-min) + min;
}
// Notice that we can get "min" exactly but never "max".

How can we get a random number in a desired range inclusive to both bounds, i.e. [X,Y]?

I suppose we could "increment" our value from Math.random() (or equivalent) by "rolling" the bits of an IEE-754 floating point double precision to put the maximum possible value at 1.0 exactly but that seems like a pain to get right, especially in languages poorly suited for bit manipulation. Is there an easier way?

(As an aside, why do random number generators produce numbers in [0,1) instead of [0,1]?)

[Edit] Please note that I have no need for this and I am fully aware that the distinction is pedantic. Just being curious and hoping for some interesting answers. Feel free to vote to close if this question is inappropriate.

I believe there is much better decision but this one should work :)

function randomInRange(min, max) {
  return Math.random() < 0.5 ? ((1-Math.random()) * (max-min) + min) : (Math.random() * (max-min) + min);
}

First off, there's a problem in your code: Try randomInRange(0,5e-324) or just enter Math.random()*5e-324 in your browser's JavaScript console.

Even without overflow/underflow/denorms, it's difficult to reason reliably about floating point ops. After a bit of digging, I can find a counterexample:

>>> a=1.0
>>> b=2**-54
>>> rand=a-2*b
>>> a
1.0
>>> b
5.551115123125783e-17
>>> rand
0.9999999999999999
>>> (a-b)*rand+b
1.0

It's easier to explain why this happens with a=2⁵³ and b=0.5: 2⁵³-1 is the next representable number down. The default rounding mode ("round to nearest even") rounds 2⁵³-0.5 up (because 2⁵³ is "even" [LSB = 0] and 2⁵³-1 is "odd" [LSB = 1]), so you subtract b and get 2⁵³, multiply to get 2⁵³-1, and add b to get 2⁵³ again.

To answer your second question: Because the underlying PRNG almost always generates a random number in the interval [0,2ⁿ-1], i.e. it generates random bits. It's very easy to pick a suitable n (the bits of precision in your floating point representation) and divide by 2ⁿ and get a predictable distribution. Note that there are some numbers in [0,1) that you will will never generate using this method (anything in (0,2^-53) with IEEE doubles).

It also means that you can do a[Math.floor(Math.random()*a.length)] and not worry about overflow (homework: In IEEE binary floating point, prove that b < 1 implies a*b < a for positive integer a).

The other nice thing is that you can think of each random output x as representing an interval [x,x+2^-53) (the not-so-nice thing is that the average value returned is slightly less than 0.5). If you return in [0,1], do you return the endpoints with the same probability as everything else, or should they only have half the probability because they only represent half the interval as everything else?

To answer the simpler question of returning a number in [0,1], the method below effectively generates an integer [0,2ⁿ] (by generating an integer in [0,2ⁿ⁺¹-1] and throwing it away if it's too big) and dividing by 2ⁿ:

function randominclusive() {
  // Generate a random "top bit". Is it set?
  while (Math.random() >= 0.5) {
    // Generate the rest of the random bits. Are they zero?
    // If so, then we've generated 2^n, and dividing by 2^n gives us 1.
    if (Math.random() == 0) { return 1.0; }
    // If not, generate a new random number.
  }
  // If the top bits are not set, just divide by 2^n.
  return Math.random();
}

The comments imply base 2, but I think the assumptions are thus:

• 0 and 1 should be returned equiprobably (i.e. the Math.random() doesn't make use of the closer spacing of floating point numbers near 0).
• Math.random() >= 0.5 with probability 1/2 (should be true for even bases)
• The underlying PRNG is good enough that we can do this.

Note that random numbers are always generated in pairs: the one in the while (a) is always followed by either the one in the if or the one at the end (b). It's fairly easy to verify that it's sensible by considering a PRNG that returns either 0 or 0.5:

• a=0   b=0  : return 0
• a=0   b=0.5: return 0.5
• a=0.5 b=0  : return 1
• a=0.5 b=0.5: loop

Problems:

• The assumptions might not be true. In particular, a common PRNG is to take the top 32 bits of a 48-bit LCG (Firefox and Java do this). To generate a double, you take 53 bits from two consecutive outputs and divide by 2⁵³, but some outputs are impossible (you can't generate 2⁵³ outputs with 48 bits of state!). I suspect some of them never return 0 (assuming single-threaded access), but I don't feel like checking Java's implementation right now.
• Math.random() is twice for every potential output as a consequence of needing to get the extra bit, but this places more constraints on the PRNG (requiring us to reason about four consecutive outputs of the above LCG).
• Math.random() is called on average about four times per output. A bit slow.
• It throws away results deterministically (assuming single-threaded access), so is pretty much guaranteed to reduce the output space.

My solution to this problem has always been to use the following in place of your upper bound.

Math.nextAfter(upperBound,upperBound+1)

or

upperBound + Double.MIN_VALUE

So your code would look like this:

double myRandomNum = Math.random() * Math.nextAfter(upperBound,upperBound+1) + lowerBound;

or

double myRandomNum = Math.random() * (upperBound + Double.MIN_VALUE) + lowerBound;

This simply increments your upper bound by the smallest double (Double.MIN_VALUE) so that your upper bound will be included as a possibility in the random calculation.

This is a good way to go about it because it does not skew the probabilities in favor of any one number.

The only case this wouldn't work is where your upper bound is equal to Double.MAX_VALUE