random real number between 0 (inclusive) and 1 (inclusive)

[Piotr Kobzda]

Patricia Shanahan wrote:

[...]

There is a long used in implementation for computations, so equivalent
nextLong(long) will need 126 bits (or two longs) to do the similar --
not so straightforward reimplementation...

But that's the problem of the optimization only. Fortunately, do not
applies to the heart of the algorithm implemented there.

Tried that, and it seems that even special treatment of the power of 2
cases is not very complicated in implementation (64 bits is enough):

public long nextlong(long n) {
if (n <= 0)
throw new IllegalArgumentException("n must be positive");

if ((n & -n) == n) // i.e., n is a power of 2

long bits, val;
do {
bits = ((((long) next(32)) << 32) >>> 1) | next(31);
val = bits % n;
} while (bits - val + (n - 1) < 0);
return val;
}


piotr

 

[Hunter Gratzner]

Your ticket and the detector are in the mail, along
with the latest batch of purple pills and the, er, latex
novelty items you ordered. Glad, as always, to provide
service to my long-time customers!

Thanks for throwing in some free vitamin pills and a rubber duck. I'l
make sure the rubber duck doesn't get lost on the beaches of Costa
Rica.

 

[Eric Sosman]

Piotr Kobzda wrote On 11/05/07 14:28,:

[...]
So, what do you think about the following:

private double lastr = rand.nextDouble();

public double nextInclusiveDouble() {
double r = lastr + rand.nextDouble();
if (r > 1d) r -= 1d;
return lastr = r;
}

The probability of each double in range [0, 1] seams to equal. Is it
correct?

I don't know. [...]

... but I've been thinking about it, and I've found the
answer. Three answers, in fact: Yes, No, and Maybe.

YES: The first number returned is the sum of two uniformly
distributed values, possibly minus one. Call the two uniform
values X and Y, and think of them as coordinates in the unit
square. The returned value is x+y for any point in the lower
left triangle or x+y-1 in the upper right triangle. This value
will be less than some chosen r, 0<=r<=1, if (x,y) is in a small
triangle at the origin (where x+y<r) or in a band just above the
diagonal (where 1<x=y<1+r). The area of the small triangle is
0.5*r^2, and the area of the band is 0.5-0.5*(1-r)^2. Add those,
and the total area in which the returned value is <r comes to r,
which is precisely the definition of a uniformly distributed
variable. So lastr is uniform after the first use, and will
remain uniform for all subsequent uses: the lastr from one use
becomes the X of the next one, with nextDouble contributing a
new Y, and the same argument applies all over again.

NO: The above assumes sums and differences are computed with
perfect accuracy, but in fact a double has finite precision: 53
bits for Java (one implicit, for normalized values). To keep
the illustration easy, imagine that double had only five bits of
precision. Express X and Y as binary fractions .xxxxx and .yyyyy,
where the x's and y's are jumbles of 1's and 0's. If the sum
z=x+y is less than 1, it looks like .zzzzz and all is well. But
if it's 1 or greater, it's 1.zzzz: one bit has been rounded off.
Now subtract 1 from this, and you get .zzzz0 -- that is, after a
subtraction step, the lowest-order fraction bit is necessarily
zero. If x+y<1 the low-order bit is about equally likely to be
zero or one, but the other half of the time it is always zero;
all in all, the low-order bit is zero 75% of the time and one
only 25% of the time. (A comment in the source of nextDouble
says that early versions of that method had a related bug.)

MAYBE: A double has 53 bits of precision when stored, but
some JVM's use higher precision in intermediate computations if
strictfp is not specified. On such JVM's, the sum of x+y might
be 1.zzzz (rounded) or 1.zzzzz (exact, in a high-precision CPU
register), and in the latter case subtracting 1 will not introduce
a low-order 0. So you're at the mercy of the particular JVM and
of exactly how the JIT decides to compile the byte code.

"Go not to c.l.j.p. for advice, for they will say both Yes
and No -- and sometimes Maybe."

 

[Patricia Shanahan]

Eric Sosman wrote:
....

NO: The above assumes sums and differences are computed with
perfect accuracy, but in fact a double has finite precision: 53
bits for Java (one implicit, for normalized values). To keep
the illustration easy, imagine that double had only five bits of
precision. Express X and Y as binary fractions .xxxxx and .yyyyy,
where the x's and y's are jumbles of 1's and 0's. If the sum
z=x+y is less than 1, it looks like .zzzzz and all is well. But
if it's 1 or greater, it's 1.zzzz: one bit has been rounded off.
Now subtract 1 from this, and you get .zzzz0 -- that is, after a
subtraction step, the lowest-order fraction bit is necessarily
zero. If x+y<1 the low-order bit is about equally likely to be
zero or one, but the other half of the time it is always zero;
all in all, the low-order bit is zero 75% of the time and one
only 25% of the time. (A comment in the source of nextDouble
says that early versions of that method had a related bug.)

....

This problem, and the general difficulties of double rounding, could be
avoided by the rewriting of the problem I have previously suggested.
Think of it as the problem of getting a long in the closed range 0
through 2**53, converting it to double, and dividing it by 2**53.

The nextDouble result multiplied by 2**53 is a integer in the range 0
through 2**53 - 1, exactly convertible to long.

private static final long MULTIPLIER = 1L<<53;

private long lastr = (long)(rand.nextDouble()*MULTIPLIER);

public double nextInclusiveDouble() {
long r = lastr + rand.nextDouble() * MULTIPLIER;
if (r > MULTIPLIER) r -= MULTIPLIER;
lastr = r;
return r/(double)MULTIPLIER;
}

If I've done this right, it is just like the previous version, except
the addition, comparison, and subtraction are done in long, which can
deal exactly with 54 bit numbers.

Patricia

 

[Patricia Shanahan]

Eric Sosman wrote:
...
...

This problem, and the general difficulties of double rounding, could be
avoided by the rewriting of the problem I have previously suggested.
Think of it as the problem of getting a long in the closed range 0
through 2**53, converting it to double, and dividing it by 2**53.

The nextDouble result multiplied by 2**53 is a integer in the range 0
through 2**53 - 1, exactly convertible to long.

private static final long MULTIPLIER = 1L<<53;

private long lastr = (long)(rand.nextDouble()*MULTIPLIER);

public double nextInclusiveDouble() {
long r = lastr + rand.nextDouble() * MULTIPLIER;

That should be:

long r = lastr + (long)(rand.nextDouble() * MULTIPLIER);

 

[Piotr Kobzda]

Patricia Shanahan wrote:

[...]

That should be:

long r = lastr + (long)(rand.nextDouble() * MULTIPLIER);

And that should be:

if (r > MULTIPLIER) r -= MULTIPLIER + 1;

Without that 1 more, zero is possible only once, just when a
nextDouble() initializing 'lastr' is zero, and is never possible again
later (the same mistake was in my initial double-based approach).


This approach eliminates the Eric's "MAYBE" answer. However, does it
really eliminate his "NO"? I'm not so sure about that...


piotr

 

[Leonard Milcin]

Hi..

I've been searching at the Internet, and in this newsgroup about this
issue, but have not been able to find what I'm searching...

What I need is to generate a lot of random real numbers between 0 and
1, both inclusive (normally it's possible to generate the number
between 0(inclusive) and 1 (exclusive!!!) )

Can anyone tell me how this is done?

/Peter

Are you really sure about the requirements? If we forget limitations of
computational technology the chance to get 1 is infinitesimal. As such
for any practical purposes you should be able to exclude 1 and be
perfectly happy with results.

Even if we take into account the fact that you don't get fine grained
results the chance to get 1 is still very close to 0 for any practical
purposes.

Best Regards,
Leonard Milcin

 

[Eric Sosman]

Piotr Kobzda wrote On 11/06/07 10:45,:

Patricia Shanahan wrote:

[...]

That should be:

long r = lastr + (long)(rand.nextDouble() * MULTIPLIER);

And that should be:

if (r > MULTIPLIER) r -= MULTIPLIER + 1;

Without that 1 more, zero is possible only once, just when a
nextDouble() initializing 'lastr' is zero, and is never possible again
later (the same mistake was in my initial double-based approach).




This approach eliminates the Eric's "MAYBE" answer. However, does it
really eliminate his "NO"? I'm not so sure about that...

Looks like it does. It certainly gets rid of the
rounding problem.

There's still a tiny bias because the amount added each
time is just a bit too small: it's chosen from [0,2**53-1]
instead of from [0,2**53] as it should be. The mean of the
generated distribution will be just a touch less than 0.5,
not 0.5 spot-on. Small values are just a hair more probable
than large values.

(I think. But I'm sometimes rong ...)

The error is too small to pay heed to, but it's really
the same error the O.P. was trying to eliminate, just sort
of "spread out" instead of concentrated on the single value
unity.

 

[Dr J R Stockton]

In comp.lang.java.programmer message <jnasi31porsaqq0nb4u3r6mkdt3otanu9r
@4ax.com>, Sun, 4 Nov 2007 20:28:11, Roedy Green <[email protected]
om.invalid> posted:

What I need is to generate a lot of random real numbers between 0 and
1, both inclusive (normally it's possible to generate the number
between 0(inclusive) and 1 (exclusive!!!) )
Multiply the number by 1+ulp where ulp in the smallest increment in a
double in the vicinity of 1.0

By stretching [0, 1) to cover [0, 1] there is necessarily a hole
created, perhaps at 0.5 .
By creating random R in [0, 1) and doing if (R==0.5) { R = 1 } ;
(apologies; I'm slowly learning to read Java;
writing it properly may follow)
one also gets a hole.
Now replace 0.5 by a random in [0, 1) independently generated and
renewed each time it is hit (so not correlated with the main Random),
and ISTM that all values in [0, 1] should be obtained with virtually
equal probability.

The Java code would be something like:

double nextInclusiveDouble(){
double r = nextDouble();
return r == nextDouble() ? 1 : r;
}

That's not the code for what I suggested. Whether it will do well
enough must depend on the details behind nextDouble. I suggested an
independent generator to maintain the Hole. Pseudo?-JavaSCRIPT :

var Hole = SomeOtherRandomFunction()

function nextInclusiveDouble() { var T
if ((T=Math.random())==Hole) {
T = 1 ; Hole = SomeOtherRandomFunction() }
return T }

Since SomeOtherRandomFunction is called only when Math.random() hits
Hole, it can be slow : something from Knuth implemented in Java[script].

If we were using a true random number generator, this would be close to
perfect. Let N be 2**53, the number of different nextDouble results, and
assume they are uniformly distributed and the results of consecutive
calls are independent.

If the works behind nextDouble were to have 53 bits (I suspect they have
more) consecutive calls would never give the same value.

The two nextDouble calls have N*N possible results.

They have whichever is the less of N**2 and 2**(GeneratorBits), ISTM.

Of those, N give
result 1. N-1 of the pairs result in e.g. 0. N numbers are each chosen
with probability (N-1)/(N*N), and 1.0 is chosen with probability
N/(N*N)=1/N.

This is indeed close enough that I don't think the difference would be
measurable over the probable lifetime of Java.

Does anyone know whether nextDouble produces pairs of consecutive equal
numbers with the correct probability, remembering that it is based on a
pseudo-random number generator?

Is nextDouble based on *A* specific PRNG, or is it based on whichever
PRNG the system author happened to like? In the latter case, your
question may have multiple answers.

 

[Patricia Shanahan]

In comp.lang.java.programmer message <jnasi31porsaqq0nb4u3r6mkdt3otanu9r
@4ax.com>, Sun, 4 Nov 2007 20:28:11, Roedy Green <[email protected]
om.invalid> posted:
What I need is to generate a lot of random real numbers between 0 and
1, both inclusive (normally it's possible to generate the number
between 0(inclusive) and 1 (exclusive!!!) )
Multiply the number by 1+ulp where ulp in the smallest increment in a
double in the vicinity of 1.0
By stretching [0, 1) to cover [0, 1] there is necessarily a hole
created, perhaps at 0.5 .
By creating random R in [0, 1) and doing if (R==0.5) { R = 1 } ;
(apologies; I'm slowly learning to read Java;
writing it properly may follow)
one also gets a hole.
Now replace 0.5 by a random in [0, 1) independently generated and
renewed each time it is hit (so not correlated with the main Random),
and ISTM that all values in [0, 1] should be obtained with virtually
equal probability.

The Java code would be something like:

double nextInclusiveDouble(){
double r = nextDouble();
return r == nextDouble() ? 1 : r;
}

That's not the code for what I suggested. Whether it will do well
enough must depend on the details behind nextDouble. I suggested an
independent generator to maintain the Hole. Pseudo?-JavaSCRIPT :

var Hole = SomeOtherRandomFunction()

function nextInclusiveDouble() { var T
if ((T=Math.random())==Hole) {
T = 1 ; Hole = SomeOtherRandomFunction() }
return T }

Since SomeOtherRandomFunction is called only when Math.random() hits
Hole, it can be slow : something from Knuth implemented in Java[script].

Java has two built-in random number generators, java.util.Random and a
subclass, java.security.SecureRandom. SecureRandom is the class to use
if you want a different, possibly better, but slower generator.

Is nextDouble based on *A* specific PRNG, or is it based on whichever
PRNG the system author happened to like? In the latter case, your
question may have multiple answers.

Random's nextDouble is specified to be based on this method:

synchronized protected int next(int bits) {
seed = (seed * 0x5DEECE66DL + 0xBL) & ((1L << 48) - 1);
return (int)(seed >>> (48 - bits));
}

seed is a long.

Each call to nextDouble uses a next(26) and a next(27).

You can see a lot of this information in the Random API documentation at
http://java.sun.com/j2se/1.5.0/docs/api/java/util/Random.html

Patricia

 

[Piotr Kobzda]

Piotr Kobzda wrote On 11/06/07 10:45,:
[...]

This approach eliminates the Eric's "MAYBE" answer. However, does it
really eliminate his "NO"? I'm not so sure about that...

Looks like it does. It certainly gets rid of the
rounding problem.

I still have some doubts, I don't like "rand.nextDouble()*MILTIPLIER"
part -- I've a feeling it somehow favors small doubles, but I may be
wrong, so never mind...

Anyway, how about reimplementing it as follows:

private static final long maxr = 1L << 53;

private long lastr = nextr();

private long nextr() {
return ((long)next(26) << 27) + next(27) + next(1);
}

public double nextInclusiveDouble() {
long r = lastr + nextr();
if (r > maxr) r -= maxr + 1;
lastr = r;
return r / (double)maxr;
}

?


piotr

 

[Tim Smith]

So, what do you think about the following:

private double lastr = rand.nextDouble();

public double nextInclusiveDouble() {
double r = lastr + rand.nextDouble();
if (r > 1d) r -= 1d;
return lastr = r;
}

The probability of each double in range [0, 1] seams to equal. Is it
correct?

Consider what happens if nextInclusiveDouble() ever returns 0. Then the
next value cannot be 1. Similarly, if it ever returns 1, the next value
cannot be 0. A correct solution would not discriminate against the
sequences 0,1 and 1,0.

 

[Piotr Kobzda]

Anyway, how about reimplementing it as follows:

[...]

Or, like that:

private static final long maxr = 1L << 53;

private long lastr;

public double nextInclusiveDouble() {
final int b1 = next(27), b2 = next(27);
long r = lastr + ((long)b1 << 26) + (b2 >>> 1) + (b2 & 1);
if (r > maxr) r -= maxr + 1;
lastr = r;
return r / (double)maxr;
}

?

Now, because we have a whole "bits domain" possible with each drawing,
initial value of 'lastr' is not important, it now just serves the role
of "error shifter". In addition, that variant is modified to "eat" one
value less from the generator's sequence (pair of 27 bits is enough).


piotr

 

[Piotr Kobzda]

Consider what happens if nextInclusiveDouble() ever returns 0. Then the
next value cannot be 1. Similarly, if it ever returns 1, the next value
cannot be 0. A correct solution would not discriminate against the
sequences 0,1 and 1,0.

Yes. The trivial solution may be calling it twice. However, the recent
versions posted by me are free of that "inconvenience". The problem is,
that they are not proven yet...


piotr

 

[Patricia Shanahan]

Piotr Kobzda wrote On 11/06/07 10:45,:
[...]

This approach eliminates the Eric's "MAYBE" answer. However, does it
really eliminate his "NO"? I'm not so sure about that...

Looks like it does. It certainly gets rid of the
rounding problem.

I still have some doubts, I don't like "rand.nextDouble()*MILTIPLIER"
part -- I've a feeling it somehow favors small doubles, but I may be
wrong, so never mind...

Take a look at the nextDouble code. It works by collecting a total of 53
bits from next(26) and next(27), concatenating them to form a double in
the range 0 through MULTIPLIER-1, and then double precision dividing the
result by MULTIPLIER. The multiplication just recovers the uniformly
distributed [0,MULTIPLIER) long.

Patricia

 

[Piotr Kobzda]

Piotr Kobzda wrote:
[...]

I still have some doubts, I don't like "rand.nextDouble()*MILTIPLIER"
part -- I've a feeling it somehow favors small doubles, but I may be
wrong, so never mind...

Take a look at the nextDouble code. It works by collecting a total of 53
bits from next(26) and next(27), concatenating them to form a double in
the range 0 through MULTIPLIER-1, and then double precision dividing the
result by MULTIPLIER. The multiplication just recovers the uniformly
distributed [0,MULTIPLIER) long.

OK, you right -- my mistake in analyzing... However, there is no need
for that (long) -> (double) -> (long) conversion, which is 1:1 in this
case.

And that's what my recent implementation is avoiding. There is just one
slight difference more (making it almost the same as was your initial
suggestion), additional bit from the generator's sequence is used to
make the result duly ranged for each nextInclusiveDouble() call. I'm
just not sure, if it is really correct now? Simply, not sure if a
"jumping" from last result is enough to prevent us from non-uniform
distribution of the resulting sequence. Do you think it's still correct?


piotr

 

[Eric Sosman]

[...]
Anyway, how about reimplementing it as follows:

private static final long maxr = 1L << 53;

private long lastr = nextr();

private long nextr() {
return ((long)next(26) << 27) + next(27) + next(1);
}

This returns 42 twice as often as it returns 0 or maxr.
The mean of the distribution is correct, but the variance is
a wee bit too small: there is too much "central tendency."

Is there any point in pursuing this further? We've known
since the beginning that there's no practical need for a [0,1]
as opposed to [0,1) distribution, and even if someone wants to
use it for "impractical" reasons a correct solution has already
been posted. The rest is just intellectual niggling about how
many angels can dance on the least-significant bit, and has
ceased to have much to do with Java programming. If the question
still interests you as a theoretical study -- there's no harm
and much good in pursuing such interests -- perhaps it would be
better explored in comp.programming, say, or maybe sci.math.

 

[Patricia Shanahan]

Piotr Kobzda wrote:
[...]

I still have some doubts, I don't like "rand.nextDouble()*MILTIPLIER"
part -- I've a feeling it somehow favors small doubles, but I may be
wrong, so never mind...

Take a look at the nextDouble code. It works by collecting a total of 53
bits from next(26) and next(27), concatenating them to form a double in
the range 0 through MULTIPLIER-1, and then double precision dividing the
result by MULTIPLIER. The multiplication just recovers the uniformly
distributed [0,MULTIPLIER) long.

OK, you right -- my mistake in analyzing... However, there is no need
for that (long) -> (double) -> (long) conversion, which is 1:1 in this
case.

The intent of the conversion was to do the arithmetic in long, rather
than double.

The issue is floating point rounding. Adding two longs in the range
[0,2**53-1] is an exact operation, because the result has no more than
54 significant bits. In double, it can cause rounding, because double
only supports 53 significant bits. After the subtraction, the result is
guaranteed to be exactly representable in double.

I have not yet analyzed your latest solution from the point of view of
floating point rounding. Did you take it into account in designing the
solution?

Patricia

 

[Eric Sosman]

Patricia Shanahan wrote On 11/06/07 20:38,:

Dr J R Stockton wrote:
[...]

Is nextDouble based on *A* specific PRNG, or is it based on whichever
PRNG the system author happened to like? In the latter case, your
question may have multiple answers.

Random's nextDouble is specified to be based on this method:

synchronized protected int next(int bits) {
seed = (seed * 0x5DEECE66DL + 0xBL) & ((1L << 48) - 1);
return (int)(seed >>> (48 - bits));
}

seed is a long.

Each call to nextDouble uses a next(26) and a next(27).

You can see a lot of this information in the Random API documentation at
http://java.sun.com/j2se/1.5.0/docs/api/java/util/Random.html

Note that java.util.Random is not final, and neither
is its next() method. This is (presumably) deliberate:
by extending Random and reimplementing next(), one can
get Random's nextDouble() et al. to use a different PRNG.
(java.util.SecureRandom does exactly this, although it
also overrides and/or adds some other methods, too.)

 

[Piotr Kobzda]

Piotr Kobzda wrote:

[...] However, there is no need

for that (long) -> (double) -> (long) conversion, which is 1:1 in this
case.

The intent of the conversion was to do the arithmetic in long, rather
than double.

Yes, I know that. Note that the first conversion in a chain mentioned
by me, i.e. (long) -> (double), is the one performed by nextDouble(),
second one is a "recovering" it back in implementation.

I have not yet analyzed your latest solution from the point of view of
floating point rounding. Did you take it into account in designing the
solution?

Yes, I did. For that aspect it do not differs from yours
implementation. It's just more efficient approach (because of lack of
unnecessary conversion).

The only significant difference is in generating a longs within a whole
range, i.e. [0, 2^53], not [0, 2^53 - 1].


piotr