views:

666

answers:

6

I'd like to get uniform distribution in range [0.0, 1.0)

If possible, please let the implementation make use of random bytes from /dev/urandom.

It would also be nice if your solution was thread-safe. If you're not sure, please indicate that.

See some solution I thought about after reading other answers.

A: 
#include <stdlib.h>
printf("%f\n", drand48());

/dev/random:

double c;
fd = open("/dev/random", O_RDONLY);
unsigned int a, b;
read(fd, &a, sizeof(a));
read(fd, &b, sizeof(b));
if (a > b)
   c = fabs((double)b / (double)a);
else
    c = fabs((double)a / (double)b);

c is your random value

There's a bug in this - if a or b is negative then the comparison doesn't necessarily choose the value with smaller magnitude as the numerator. Also, I'm too dim to easily see that it produces uniformly distributed output.
Steve Jessop
Also, one time in 2^64 you hit double jackpot, and get division by 0 undefined behaviour.
Steve Jessop
Thinking of which, this is never uniform. You get 0 if either value is 0, which is almost 1 in 2^31, which is far more than you should in a double with 52 bits of mantissa precision.
Steve Jessop
+1  A: 

Reading from files is thread-safe AFAIK, so using fopen() to read from /dev/urandom will yield "truly random" bytes.

Although there might be potential gotchas, methinks any set of such bytes accessed as an integer, divided by the maximum integer of that size, will yield a floating point value between 0 and 1 with approximately that distribution.

Eg:

FILE* f = fopen("/dev/urandom", "r");
int32_t int;
fread(&int, sizeof(int32_t), 1, f);
fclose(f);
double theRandomValue = int / (double) (2 ** 32 - 1);
millenomi
A few notes:(1) You don't want the "- 1" after 2 ^ 32 since phjr asked for output 1.0 to be excluded.(2) There is no ** in C++. Only << (use 1 << 32).(3) You probably want an unsigned int here.
Tyler
And you can't use 'int' as a variable because it is a keyword.
Jonathan Leffler
A: 

The trick is you need a 54 bit randomizer that meets your requirements. A few lines of code with a union to stick those 54 bits in the mantissa and you have your number. The trick is not double float the trick is your desired randomizer.

dwelch
Does POSIX mandate IEEE floats?
Steve Jessop
According to java.util.Random you need 53 bits, not 54. Read the comments on the nextDouble() method to see why. Otherwise, you're exactly on the right track.
Chris Jester-Young
Although the Java approach doesn't use a union for that purpose, it just makes a large 53-bit number, then divides it by (1 << 53).
Chris Jester-Young
That should also avoid generating denormalised values 50% of the time, which the union will do. Your floating point unit and libraries may or may not handle denorms gracefully.
Steve Jessop
+2  A: 

Simple: A double has 52 bits of precision assuming IEEE. So generate a 52 bit (or larger) unsigned random integer (for example by reading bytes from dev/urandom), convert it into a double and divide it by 2^(number of bits it was).

This gives a numerically uniform distribution (in that the probability of a value being in a given range is proportional to the range) down to the 52nd binary digit.

Complicated: However, there are a lot of double values in the range [0,1) which the above cannot generate. To be specific, half the values in the range [0,0.5) (the ones that have their least significant bit set) can't occur. Three quarters of the values in the range [0,0.25) (the ones that have either of their least 2 bits set) can't occur, etc, all the way to only one positive value less than 2^-51 being possible, despite a double being capable of representing squillions of such values. So it can't be said to be truly uniform across the specified range to full precision.

Of course we don't want to choose one of those doubles with equal probability, because then the resulting number will on average be too small. We still need the probability of the result being in a given range to be proportional to the range, but with a higher precision on what ranges that works for.

I think the following works. I haven't particularly studied or tested this algorithm (as you can probably tell by the way there's no code), and personally I wouldn't use it without finding proper references indicating it's valid. But here goes:

  • Start the exponent off at 52 and choose a 52-bit random unsigned integer (assuming 52 bits of mantissa).
  • If the most significant bit of the integer is 0, increase the exponent by one, shift the integer left by one, and fill the least significant bit in with a new random bit.
  • Repeat until either you hit a 1 in the most significant place, or else the exponent gets too big for your double (1023. Or possibly 1022).
  • If you found a 1, divide your value by 2^exponent. If you got all zeroes, return 0 (I know, that's not actually a special case, but it bears emphasis how very unlikely a 0 return is [Edit: actually it might be a special case - it depends whether or not you want to generate denorms. If not then once you have enough 0s in a row you discard anything left and return 0. But in practice this is so unlikely as to be negligible, unless the random source isn't random).

I don't know whether there's actually any practical use for such a random double, mind you. Your definition of random should depend to an extent what it's for. But if you can benefit from all 52 of its significant bits being random, this might actually be helpful.

Steve Jessop
I was going to disagree with your complicated approach, but once I read through the whole thing it actually made some sense. Although I'm not an expert at any of this, so my opinion shouldn't count. :-P
Chris Jester-Young
Me neither - I think the best approach is not to design algorithms which require "random floats" unless you're sure you know what that's supposed to mean. I'm not sure I do...
Steve Jessop
A: 

This seems to be pretty good way:

unsigned short int r1, r2, r3;
// let r1, r2 and r3 hold random values
double result = ldexp(r1, -48) + ldexp(r2, -32) + ldexp(r3, -16);

This is based on NetBSD's drand48 implementation.

phjr
A: 

/dev/urandom is not POSIX, and is not generally available.

The standard way of generating a double uniformly in [0,1) is to generate an integer in the range [0,2^N) and divide by 2^N. So pick your favorite random number generator and use it. For simulations, mine is the Mersenne Twister, as it is extremely quick, yet still not well correlated. Actually, it can do this for you, and even has a version that will give more precision for the smaller numbers. Typically you give it a seed to start with, which helps for repeatability for debugging or showing others your results. Of course, you can have your code grab a random number from /dev/urandom as the seed if one isn't specified.

For cryptographic purposes you should use one of the standard cryptographic libraries out there instead, such as openssl) which will indeed use /dev/urandom when available.

As to thread safety, most won't be, at least with the standard interfaces, so you'll need to build a layer on top, or only use them in one thread. The ones that are thread safe have you supply a state that they modify, so that instead you are effectively running multiple non-interacting random number generators, which may not be what you are looking for.

wnoise