Random numbers between a and b: what's the ultimate answer?
Andrei Alexandrescu
wanna know the answer.
So there's this rand() function that returns me numbers between 0 and
RAND_MAX inclusive. How do I use that function to generate random numbers
between two arbitrary numbers a and b with a < b?
The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My intention,
however, is to not use floating points at all AND to generate a nice random
sequence (avoid nonuniformities caused by overflow, clustering, etc.) AND
support up to 32-bit random numbers with 32-bit math (if possible).
So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX, then
likely the result ain't gonna be very nicely distributed due to math
overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the math
is in 32 bits as well, we ain't gonna see anything but a.
Thanks!
Andrei
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
[ about comp.lang.c++.moderated. First time posters: do this! ]
Francis Glassborow
Alexandrescu <andre...@hotmail.com> writes
>Ok, here's a silly one. I know I'll likely make a fool of myself, but I
>wanna know the answer.
>
>So there's this rand() function that returns me numbers between 0 and
>RAND_MAX inclusive. How do I use that function to generate random
numbers
>between two arbitrary numbers a and b with a < b?
>
>The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
>however, is to not use floating points at all AND to generate a nice
random
>sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
>support up to 32-bit random numbers with 32-bit math (if possible).
>
>So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
>likely the result ain't gonna be very nicely distributed due to math
>overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
math
>is in 32 bits as well, we ain't gonna see anything but a.
Try a simpler problem. How do you do the above when a is zero? Iff you
can solve the problem of generating pseudo-random numbers in the range 0
-> a can you solve the problem for a -> b.
Now there is a good deal of existing practice for this simpler version.
Do you want me to dig it out for you? :-)
--
ACCU Spring Conference 2003 April 2-5
The Conference you cannot afford to miss
Check the details: http://www.accuconference.co.uk/
Francis Glassborow ACCU
Daniel Frey
>
> Ok, here's a silly one. I know I'll likely make a fool of myself, but
I
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
"*b" should be "*(b-a)", which leads to a reduction of the problem: We
only need to find a function to map the random number to 0-n for some n
(here: b-a). I also assume that RAND_MAX > n. The next "obvious"
solution is modulo, which leads to:
template< typename T > T map_random( T rnd, T max ) // map to 0-max
{
return rnd % ( max + 1 ); // remember I assumed: max < RAND_MAX
}
template< typename T > T map_random( T rnd, T min, T max )
{
return min + map_random( rnd, max - min );
}
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
The only problem is probably "clustering" IIUC. You want to make sure
that the lower numbers are not more likely than the higher numbers,
right? (Example: map 0-9 to 0-3 leads to this distribution:
0: 0 4 8
1: 1 5 9
2: 2 6
3: 3 7
The only thing left to do is to "spread" the inbalance, but conceptually
it cannot be eliminated, so there will always be a set of numbers that
is more likely than the other numbers (expect for RAND_MAX % ( max+1 )
== 0 ). At least not with a stateless/pure function, but I am not sure
that this is what you want...
Regards, Daniel
--
Daniel Frey
aixigo AG - financial training, research and technology
Schloß-Rahe-Straße 15, 52072 Aachen, Germany
fon: +49 (0)241 936737-42, fax: +49 (0)241 936737-99
eMail: danie...@aixigo.de, web: http://www.aixigo.de
Carlos Moreno
> Ok, here's a silly one. I know I'll likely make a fool of myself, but
I
> wanna know the answer.
Hey, don't worry! This way, you give us (common mortals) the
opportunity to feel good about ourselves (hey, we know the
answer to one of Alexandrescu's questions!! :-))
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
>
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
> likely the result ain't gonna be very nicely distributed due to math
> overflow and truncation.
The perfect one includes a loop (well, including a loop would
probably make it non-perfect; but the probability that you
loop more than once is extremely small, depending on the
value of b-a).
The idea is that you want to draw random numbers until you
get one number that is less than the highest multiple of
(b-a). That way, you avoid the non-uniformity that would
happen if you use %, or the rounding error if you divide.
Notice that if RAND_MAX is 2^32-1 and (b-a) is, say, a
few hundreds, then the probability that you loop more
than once is in the 10^-7 range. The probability that
you loop N times is in the 10^-(7*N) range.
The code I have (based on this idea) is:
const int divisor = RAND_MAX / (b - a);
const int highest = divisor * (b - a);
int rnd = 0;
do
{
rnd = rand(); // or lrand48(), or read from /dev/urandom
}
while (rnd >= highest);
return a + rnd / divisor;
The above was slightly modified (the code, as I have it, gives
a random number between 0 and N-1, and not between a and b), and
it assumes both b and a positive, and b > a. (hmmm, now that
I think about it, this one gives values between a and b-1,
but anyway, the idea is the same).
You could combine it with a "max_retries" idea, if you are
really paranoid. But again, this should loop only once in
most cases.
HTH,
Carlos
--
John Lloyd
> Ok, here's a silly one. I know I'll likely make a fool of myself, but
I
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
> however, is to not use floating points at all AND to generate a nice
> random sequence (avoid nonuniformities caused by overflow, clustering,
> etc.) AND support up to 32-bit random numbers with 32-bit math (if
> possible).
>
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
Don't you mean a + (rand() * (b - a)) / RAND_MAX ?
> then likely the result ain't gonna be very nicely distributed due to
math
> overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
> math is in 32 bits as well, we ain't gonna see anything but a.
If we assume a and b are integers with a <= b, maybe something like
a + ( ( rand() % RAND_MAX ) / ( RAND_MAX / (b - a + 1) ) )
Try it, or something similar, and have a look at the empirical
distribution.
For better random number generators than rand() you could try the
following:
Random Number Generation and Monte Carlo Methods, J. E. Gentle,
Springer-Verlag.
http://www.amazon.com/exec/obidos/tg/detail/-/0387985220/qid=1045054058/
sr=1-1/ref=sr_1_1/104-9791777-8159107?v=glance&s=books
Thomas Mang
Andrei Alexandrescu schrieb:
> Ok, here's a silly one. I know I'll likely make a fool of myself, but
I
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
>
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
> likely the result ain't gonna be very nicely distributed due to math
> overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
math
> is in 32 bits as well, we ain't gonna see anything but a.
My solution is that simple it may well miss the point, which I apologize
if it
doesn't suit your needs:
2 steps:
a) request a random number in the range [0, b - a]
b) add this number to a
In order to achieve point a), write your own max_rand functions which
looks
something like this (algorithm taken from "Acceleratied C++" by
Koenig/Moo) and
provides a random distribution between 0 and maxValue:
int max_rand(int maxValue)
{
const int buckets = RAND_MAX / maxValue;
int result = 0;
do
result = rand() / buckets;
while(result > maxValue);
return result;
}
hope this helps,
Thomas
Elite Software
I would have written:
range = b - a; // b > a
result = rand() % range;
result += a;
or more succinctly:
result = a + rand() % (b - a);
Does this do what you want or does it fail in some way?
Regards,
Michael Hardwick
michael@hawthylands.+.com replace plus with plus
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
Matt Sable
generators,
I highly recommend reading the chapter in Numerical Recipes on the
topic -
they cover many details comprehensively, including tradeoffs between
efficiency
and "randomness" of the generated numbers, various algorithms, etc...
PDFs
of the chapters are available for free online :
http://lib-www.lanl.gov/numerical/bookcpdf.html See Chapter 7.
Regards,
Matthias
> Ok, here's a silly one. I know I'll likely make a fool of myself, but
I
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
>
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
> likely the result ain't gonna be very nicely distributed due to math
> overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
math
> is in 32 bits as well, we ain't gonna see anything but a.
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
Randy Maddox
news:<b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de>...
> Ok, here's a silly one. I know I'll likely make a fool of myself, but
I
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
>
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
> likely the result ain't gonna be very nicely distributed due to math
> overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
math
> is in 32 bits as well, we ain't gonna see anything but a.
>
> Thanks!
>
> Andrei
>
>
What about:
x = rand();
if((x < a) || (x > b))
{
x %= (b - a); // x now <= (b - a)
x += a; // now a <= x <= b
}
return x;
Not completely sure about the overflow, truncation or clustering
issues, but at least this does generate numbers in the correct range,
and I believe it might actually be correct. I'm definitely interested
to see what others have to say since I ran into a similar need and
solved it using doubles for the obvious calculation you showed above
and converting back to int for the return value. If the calculation
above is correct, then it should certainly be faster.
Any numerical wizards out there with some insight into why this will
or will not work correctly?
Randy.
Andrei Alexandrescu
>
Try a simpler problem. How do you do the above when a is zero? Iff you
> can solve the problem of generating pseudo-random numbers in the range
0
> -> a can you solve the problem for a -> b.
>
> Now there is a good deal of existing practice for this simpler
version.
> Do you want me to dig it out for you? :-)
Yes please. Because I did search before posting here. I found a number
of
functions that I could prove are incorrect (generate nonuniform
sequences)
or slow (use floating point).
Andrei
Thant Tessman
[...]
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
First off, assuming you want a random number in a range that includes
'a' but not 'b', you might mean:
r = a + (double(rand())/RAND_MAX) * (b-a)
But even this is wrong. There is a 1 in RAND_MAX chance that you'll get
a 1 for (rand/rand_max), and what you want is a number between zero and
one that doesn't include one. So I guess you have to do something like
this:
int raw_rand = rand();
r = a + (double(raw_rand?raw_rand-1:0)/RAND_MAX)) * (b-a);
But even this is wrong because it ever so slightly increases the chance
of getting a value of 'a'.
Another approach:
Assuming you have a high-quality rand function (a big assumption), and
speed isn't an issue, the answer might be to call 'rand' multiple times
until you get an integer value <= (b-a). You can narrow down the misses
by masking off the bits higher than the power of two greater than your
target.
Or better yet, consume all the bits from your generator one-by-one as
you need them. (I vaguely remember reading somewhere that often the bits
in random number generators aren't equally random. Applications will use
xor to randomize bits with other bits from different places.)
-thant
Radoslav Getov
news:<b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de>...
I
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
>
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
> likely the result ain't gonna be very nicely distributed due to math
> overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
math
> is in 32 bits as well, we ain't gonna see anything but a.
>
> Thanks!
>
> Andrei
>
It wasn't quite clear, so I assume that the result shold be in range
[a..b).
The base idea is to achieve a flat distribution by throwing away some
rand() values.
assert (b > a);
int range = b - a,
max_rand = RAND_MAX / range * range;
while (true)
{
int r = rand();
if (r >= max_rand)
continue;
return a + r % (b - a);
}
Radoslav Getov
www.ultraheap.com
John Lloyd
>
> If we assume a and b are integers with a <= b, maybe something like
>
> a + ( ( rand() % RAND_MAX ) / ( RAND_MAX / (b - a + 1) ) )
>
Actually, I realised about 5 secs after posting this that it won't give a
uniform 'distribution' - I was trying to get rid of the potential overflow
problems by dividing rather than multiplying and didn't take into account
the truncation that would occur during the integer division.
Francis Glassborow
Alexandrescu <SeeWebsit...@moderncppdesign.com> writes
>Yes please. Because I did search before posting here. I found a number
>of
>functions that I could prove are incorrect (generate nonuniform
>sequences)
>or slow (use floating point).
There was a long debate on this in C Vu several years ago. However if
you want a uniform distribution you just about have to be willing to
discard some high values. The idea is to find (RAND_MAX-1) / n * n where
n is the number of elements in the required range. That give a ceiling,
all values above that must be ignored (note that templates allow you to
determine and apply that statically. Now the guts of the function looks
something like:
do {
int const value(rand());
if (value < ceiling) return value % n;
} while( true);
I am sure the details are wrong, but I am also sure the principle is
right. I am sure you can add your elegance and correct it.
--
ACCU Spring Conference 2003 April 2-5
The Conference you cannot afford to miss
Check the details: http://www.accuconference.co.uk/
Francis Glassborow ACCU
Matthew Collett
"Andrei Alexandrescu" <SeeWebsit...@moderncppdesign.com> wrote:
> I found a number
> of
> functions that I could prove are incorrect (generate nonuniform
> sequences)
> or slow (use floating point).
>
> Andrei
You are undoubtedly right about the incorrectness, but are you sure
about the sloth? It is my understanding that on most modern processors
floating point operations are about the same speed as integer ones; the
exception is division, which is much faster for floating point. (One
set of sample numbers can be found in Section 7.6 of Kernighan and Pike
"The Practice of Programming".) Converting to and fro might be
expensive though.
Best wishes,
Matthew Collett
--
Those who assert that the mathematical sciences have nothing to say
about the good or the beautiful are mistaken. -- Aristotle
Andrew Koenig
Andrei> a number of functions that I could prove are incorrect
Andrei> (generate nonuniform sequences) or slow (use floating point).
I can't resist. From ``Accelerated C++'', page 135:
// return a random integer in the range [0, n)
int nrand(int n)
{
if (n <= 0 || n > RAND_MAX)
throw domain_error("Argument to nrand is out of range");
const int bucket_size = RAND_MAX / n;
int r;
do r = rand() / bucket_size;
while (r >= n);
return r;
}
There is no fixed, finite bound on the worst-case execution time of
this program, but it terminates with probability 1 and its average
performance is O(1). I strongly doubt it is possible to do any better.
--
Andrew Koenig, a...@research.att.com, http://www.research.att.com/info/ark
Andrew Koenig
Elite> result = a + rand() % (b - a);
Elite> Does this do what you want or does it fail in some way?
It fails. To see why, imagine a machine on which RAND_MAX is 32767
and (b - a) is 30000. Then results between 0 and 2767 are twice as
likely as others.
--
Andrew Koenig, a...@research.att.com,
http://www.research.att.com/info/ark
Chris Saam
in
[a, b) where b - a < RAND_MAX.
Let n = b - a;
Devide [0, RAND_MAX] into n subintevals of equal length x plus one
smaller
inteval
[i *x, (i + 1) *x) where i = 0 (1) n and [n * x, RAND_MAX]
Pick a random number, if that number is greater then or equal to n*x
throw
it away and pick another.
Otherwise use the mapping
f(z) = i for all z in [i *x, (i + 1) *x)
The code looks like this
int RandRange(int a, int b)
{
const int nIntSize = RAND_MAX / (b - a);
const int nMax = nIntSize * (b - a);
int nRes = rand();
for (nRes = rand(); nRes >= nMax; nRes = rand());
return a + nRes / nIntSize;
}
"Andrei Alexandrescu" <andre...@hotmail.com> wrote in message
news:b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de...
> Ok, here's a silly one. I know I'll likely make a fool of myself, but
I
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
>
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
> likely the result ain't gonna be very nicely distributed due to math
> overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
math
> is in 32 bits as well, we ain't gonna see anything but a.
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
Carlos Moreno
>
> result = a + rand() % (b - a);
>
> Does this do what you want or does it fail in some way?
This is a bad idea, since it doesn't give a perfectly
flat distribution (as Andrei explicitly said was one of
the requirements).
An example to easily visualize this: imagine that RAND_MAX
is 256 (you get one random byte), and that you want a random
number between 0 and 199.
You would do rand() % 200, right? Well, your random
output will indeed be between 0 and 199 -- but you will
get a 2/256 probability for all the numbers between 0
and 55, and a 1/256 probability for the rest -- not
good, huh? :-(
Of course, when RAND_MAX is 2.1 billion, and you want
a random number between 0 and, say, a few hundreds or
a few thousands, the non-uniformity is not that obvious
(my example was taken to the extreme, so that the
problem becomes evident), but it is still non-uniform.
The other problem is that in general, the lower bits
carry less "randomness" than the higher bits -- when
you take %, you are discarding the higher bits; so,
it is better to always reduce to the required range
by dividing instead of taking modulo. (see other
posts -- including mine -- for a reasonable solution
where you loop discarding values outside the range
that gives an exact multiple of the range you want).
Cheers,
Carlos
--
Jim Melton
news:b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de...
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random numbers
> between two arbitrary numbers a and b with a < b?
Not a problem. I'll make a fool of myself by giving a simple answer. I know
some people take this random number thing *very* seriously. My system
(Solaris 8) has a drand48() function that returns a uniformly distributed
random number between 0 and 1. From there, your formula should be simple:
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My intention,
> however, is to not use floating points at all AND to generate a nice
random
Why the aversion to floating point? Modern hardware can actually be faster
at floating point than integer.
a + (drand48() * b)
--
<disclaimer>
Opinions posted are those of the author.
My company doesn't pay me enough to speak for them.
</disclaimer>
--
Jim Melton
Software Architect, Fusion Programs
Lockheed Martin Astronautics
(303) 971-3846
galathaea
news:<b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de>...
I
> wanna know the answer.
>
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random
numbers
> between two arbitrary numbers a and b with a < b?
>
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b.
I would say the obvious formula is better:
a + (b - a + 1) * int(double(rand()) / (RAND_MAX + 1.0))
for uniformity and correct range, but your intent is clear.
> My intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
>
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
> likely the result ain't gonna be very nicely distributed due to math
> overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
math
> is in 32 bits as well, we ain't gonna see anything but a.
It seems a tricky problem to me trying to provide the correct
distribution and not get stuck with RAND_MAX + 1 problems, but the
ugliest solution I can think of would be something like:
int rangeSize = b - a + 1;
int randExcess = (RAND_MAX % rangeSize + 1) % rangeSize;
int myRandomNumber;
if (randExcess)
{
while ((myRandomNumber = rand()) > (RAND_MAX - randExcess));
}
else
{
myRandomNumber = rand();
}
myRandomNumber = a + myRandomNumber % rangeSize;
This is obviously ugly and slower (I particularly hate the double
modulo for randExcess, but did not come up with anything better -- and
the while is pitiful) than what one would desire (and it looks like
there is a problem with a = 0, b = largest integer -- things get
uglier still if that is the case), but it is the only thing I can work
out to avoid the distribution problems of modulo clipping and RAND_MAX
+ 1 overflow by adapting the rand() using only integer arithmetic...
I sure hope someone else posts a better response.
Shimshon Duvdevan
"Andrei Alexandrescu" <andre...@hotmail.com> wrote in message
news:<b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de>...
> So what is the ultimate formula? If I use a + (rand() * b) / RAND_MAX,
then
> likely the result ain't gonna be very nicely distributed due to math
> overflow and truncation. For example, if RAND_MAX is 2**32 - 1 and the
math
> is in 32 bits as well, we ain't gonna see anything but a.
When writing a Randomizer class for my project, I thought about the
problem of getting uniformly distributed numbers in [0, max), and,
after consulting with several other approaches (see java.util.Random
documentation), came up with this function (the line "*this >> rnd"
puts random bits in rnd, in my implementation I use /dev/urandom).
Note that all enthropy is used, even if the random value is in
"non-uniform" range (all other uniform distribution approaches I saw
just try again). In short, read the comments :).
Compare this to
http://java.sun.com/j2se/1.4/docs/api/java/util/Random.html#nextInt(int)
I remember that when analyzing this function, I came up with expected
number of loops in the worst case being the same as in Sun's approach,
but I think the distribution in my case is better (I actually use the
entropy provided by "bad" random value), and I think that my approach
is cleaner (less tricky, and I think their's wouldn't work on
non-2-complement representations), but it does assume perfect random
source.
Of course, the transition from [0 .. max) to [a .. b] is not hard, and
I think that finding a source of random bits is separate issue.
Shuki.
PS. I wonder whether using the entropy provided by "bad" random values
is genuine, or I was just reinventing the wheel. All comments are
welcome. :)
/**
* Returns random value in [0 .. @a max).
*
* This function is not as effective as operator>>()
* overloaders, since we can't just take an unsigned
* value modulo @a max.
*
* @param max exclusive maximum (must be positive)
*
* @return random value in [0 .. @a max - 1].
* The function may also never return, although
* the probability of this happening is 0.
*/
unsigned
Randomizer::
random(unsigned max)
{
const unsigned TYPE_MAX = std::numeric_limits<unsigned>::max();
// random value accumulator
unsigned acc = 0;
while (true) {
// rnd <- [0 .. UINT_MAX]
unsigned rnd; *this >> rnd;
// [UINT_MAX - @a edge .. UINT_MAX] gets special treatment
const unsigned edge = TYPE_MAX % max;
const unsigned edgeStart = TYPE_MAX - edge;
/**
* If @a rnd is in the @em safe range of
* [0 .. UINT_MAX - @a edge - 1], we just
* return @a rnd modulo @a max
*
* There is also the case when @a edge equals
* @a max - 1, so there is no edge. This
* case can also be taken care of by the
* following @a else clause.
*/
if (rnd < edgeStart || edge == max - 1) {
return acc + rnd % max;
}
/**
* Otherwise, we have the partial random value
* in [UINT_MAX - @a edge .. UINT_MAX] range,
* and we add it with random(@a max - @a edge)
*
* If this range is a singleton, we have the special
* case of @a edge = 0 and @a rnd = UINT_MAX, which
* requires rerunning the algorithm with the same input.
* This case is too rare to devote an @c if clause
* to it, though, and current treatment is faster anyway.
*/
else {
max -= edge;
acc += rnd - edgeStart;
Erik Max Francis
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b.
I think you meant a + ((double(rand())/RAND_MAX)*(b - a) :-).
--
Erik Max Francis / m...@alcyone.com / http://www.alcyone.com/max/
__ San Jose, CA, USA / 37 20 N 121 53 W / &tSftDotIotE
/ \ You cannot step into the same river once.
\__/ Cratylus
Python chess module / http://www.alcyone.com/pyos/chess/
A chess game adjudicator in Python.
Pierre Baillargeon
news:<b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de>...
>
intention,
> however, is to not use floating points at all AND to generate a nice
random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
AND
> support up to 32-bit random numbers with 32-bit math (if possible).
Note that even when using floating points you do not avoid clustering
and non-uniformities: rand() returns a discrete integer value anyway,
so the FP multiplication is just a compact mapping table, but it does
not magically avoids clustering.
The only way to have uniformity is to divide the RAND_MAX range in
equal-sized chunks, discarding unused values at the end of the
[0,RAND_MAX] range. For example, for [a,b] = [1,3] and [0,RAND_MAX] =
[0,15], you would discards the value 15 and divide the rest in thirds.
Note that if the size of [a,b] is prime and larger than [0,RAND_MAX]
(or non-prime but all factors are larger than RAND_MAX) then you will
have to use two or more rand() calls as seeds into a larger
randomizer.
If you require to be sure that your function will run in constant
time, by avoiding discarding values, then it is impossible unless the
size of [a,b] is an exact divisor of the size [0,MAX_RAND].
Andrew Simmons
Andrew Koenig wrote:
> I strongly doubt it is possible to do any better.
>
>
I don't know much about this, but have for some time been relying on Knuth
(TAOCP volume 2, 3rd edition, p185) - "To compute a random integer between 0
and k-1, one should multiply by k and truncate the result." Is he wrong?
Marcelo E. Magallon
> It is my understanding that on most modern processors floating point
> operations are about the same speed as integer ones
It depends on how flexible you want to be with that "about". Integer
operations are certainly not 1000% faster than floating point ones (which
used to be the case -- back then when FP coprocessors weren't commodity),
but the speed difference is large enough to be significant (> 100%).
> the exception is division, which is much faster for floating point.
Uhm, I'm curious what do you mean with division here... I've got code
that says this isn't true. If you mean cast to float, divide and cast
back to int... yeah, perhaps.
--
Email sent to this address will be considered spam.
Google my name if you wish to contact me by email.
John Potter
<moreno_at_mo...@xx.xxx> wrote:
> The other problem is that in general, the lower bits
> carry less "randomness" than the higher bits
Do you have measurements to support this? When I considered
individual bits and measured 0/1 and for pairs 00/01/10/11,
it made no difference which bit was selected.
> -- when
> you take %, you are discarding the higher bits;
That is only true if the divisor is a power of two. In other
cases, remainder is the result of all bits.
> so,
> it is better to always reduce to the required range
> by dividing instead of taking modulo.
True, but nothing in the above supports that.
John
Risto Lankinen
"Matthew Collett" <m.co...@auckland.ac.nz> wrote in message
news:m.collett-53252...@lust.ihug.co.nz...
> You are undoubtedly right about the incorrectness, but are you sure
processors
> floating point operations are about the same speed as integer ones; the
> exception is division, which is much faster for floating point.
Speed is not the concern, but the floating point solution,
too, suffers from a non-uniform distribution. In fact, any
algorithm that tries to map all RAND_MAX outputs to a
range whose size does not evenly divide RAND_MAX
will suffer the same (by pigeonhole principle).
Example: suppose RAND_MAX = 32767 ; generate numbers
between 0..32766 :
0: 0.0/32767.0 * 32766.0 -> to int = 0
1: 1.0/32767.0 * 32766.0 -> to int = 0
2: 2.0/32767.0 * 32766.0 -> to int = 1
3: 3.0/32767.0 * 32766.0 -> to int = 2
4: ...
Already it's obvious that '0' is twice more likely than '1'.
For a truly uniform distribution, some output values of
rand() need to be discarded.
Cheers!
- Risto -
Beta Jin
news:<yu9965rp...@europa.research.att.com>...
> Elite> or more succinctly:
>
> Elite> result = a + rand() % (b - a);
>
> Elite> Does this do what you want or does it fail in some way?
>
> It fails. To see why, imagine a machine on which RAND_MAX is 32767
> and (b - a) is 30000. Then results between 0 and 2767 are twice as
> likely as others.
Then, how about
result = a + (rand()+rand()) % (b - a)
for the case? Anyway, a good rand() implemenation is required.
In the real world, I think we might consider a and b are const values.
So it is okay to specialize every range of [a, b] to meet a good
RAND_MAX_X, where RAND_MAX_X = RAND_MAX * X.
Gerhard Wesp
> The obvious formula is: a + (double(rand()) / RAND_MAX) * b. My
> intention,
> however, is to not use floating points at all AND to generate a nice
> random
> sequence (avoid nonuniformities caused by overflow, clustering, etc.)
> AND
> support up to 32-bit random numbers with 32-bit math (if possible).
Just use the obvious formula. The general advice to avoid premature
optimization applies here as well---Unless a profiler tells you that the
floating point operations really cost you much performance (which I
personally think would be highly unlikely), go ahead and use them.
If you need high quality (or perhaps even cryptographic) random numbers,
don't use rand(). Instead, I recommend some system's /dev/random which
typically yields high-quality non-deterministic random numbers, or some
encryption algorithm for deterministic ones. E.g.
http://www.cosy.sbg.ac.at/~gwesp/sw/rijndael-1.0.tar.gz
-Gerhard
--
| voice: +43 (0)676 6253725 *** web: http://www.cosy.sbg.ac.at/~gwesp/
|
| Passts auf, seid's vuasichdig, und lossds eich nix gfoin!
| -- Dr. Kurt Ostbahn
Matthew Collett
"Marcelo E. Magallon" <hm087...@sneakemail.com> wrote:
>> It is my understanding that on most modern processors floating point
>> operations are about the same speed as integer ones
>
> It depends on how flexible you want to be with that "about". Integer
> operations are certainly not 1000% faster than floating point ones
> (which
> used to be the case -- back then when FP coprocessors weren't
> commodity),
> but the speed difference is large enough to be significant (> 100%).
Obviously this depends on the hardware. I did some quick tests on the
machine I use for reading & posting news, which happens to be a 450MHz
G4 PowerPC. Approximate times per operation:
int+int, int-int 1.3ns
int*int 4.3ns
int/int 48ns
int%int 53ns
float+float, float-float, float*float 2.9ns
float/float 43ns
double+double, double-double, double*double 2.9ns
double/double 77ns
So on this particular machine, integer division (either / or %) is
slower than single precision, but faster than double precision; integer
addition is faster than either size of floating point, but
multiplication is slower. All operations on short (omitted from the
table) take 3-4 ns longer than the corresponding ones on int.
Can anyone point to somewhere on the net that has this sort of
'low-level benchmark' data for a range of hardware?
Best wishes,
Matthew Collett
--
Those who assert that the mathematical sciences have nothing to say
about the good or the beautiful are mistaken. -- Aristotle
Stephen Howe
>
> result = a + (rand()+rand()) % (b - a)
>
> for the case? Anyway, a good rand() implemenation is required.
No that is no good because it is non-uniform. The expression (rand()+rand())
is similar to tossing 2 dice. If you toss 2 dice, 12 can come up 1 way but 7
can up 6 ways.
For (rand()+rand()) there will be a bulge at RAND_MAX and the values 0 and 2
* RAND_MAX will be the most unlikely values. The expression
(rand()+rand()) % (b - a)
will skew the distribution in a odd pattern such that some numbers wil lbe
more likely than others.
Stephen Howe
Andrew Koenig
Andrew> relying on Knuth (TAOCP volume 2, 3rd edition, p185) - "To
Andrew> compute a random integer between 0 and k-1, one should
Andrew> multiply by k and truncate the result." Is he wrong?
I don't have my copy handy, but from the description it sounds like
he's starting with a uniformly distributed real number in [0,1).
That's a different problem entirely.
--
Andrew Koenig, a...@research.att.com, http://www.research.att.com/info/ark
Andrew Koenig
Beta> result = a + (rand()+rand()) % (b - a)
Beta> for the case? Anyway, a good rand() implemenation is required.
If you add two uniformly distributed random numbers, the result is no
longer uniformly distributed. So the approach makes the problem
worse, not better.
--
Andrew Koenig, a...@research.att.com, http://www.research.att.com/info/ark
Carlos Moreno
>>Elite> result = a + rand() % (b - a);
>>
>>Elite> Does this do what you want or does it fail in some way?
>>
>>It fails. To see why, imagine a machine on which RAND_MAX is 32767
>>and (b - a) is 30000. Then results between 0 and 2767 are twice as
>>likely as others.
>
> Then, how about
>
> result = a + (rand()+rand()) % (b - a)
>
> for the case?
It suffers from exactly the same problem. Could you maybe explain
why would you think it's any different? (I might be missing something)
Carlos
--
Carlos Moreno
> On 13 Feb 2003 05:12:04 -0500, Carlos Moreno
> <moreno_at_mo...@xx.xxx> wrote:
>
>
>>The other problem is that in general, the lower bits
>>carry less "randomness" than the higher bits
>
>
> Do you have measurements to support this?
No, but I trust that Donald Knuth does :-) (i.e., he
states that in his book, in the section of linear
congruential random number generators)
>>-- when
>>you take %, you are discarding the higher bits;
>
> That is only true if the divisor is a power of two. In other
> cases, remainder is the result of all bits.
Hmmm, I'm not sure this is entirely true. If the divider
is even, then the lowest bit will have a direct effect
in the remainder -- i.e., if your base number has higher
probability of producing an even number, then your output
after taking the modulo would also have a higher probability
of being even.
Actually, I think he [Knuth] also suggests that before taking
the modulo, the base number should be divided by a number
that is relatively prime to the divider. That way you
avoid the problem I just described.
But then again, why all the trouble, when we have a *perfect*
solution? :-) (i.e., draw numbers until you get one that
is less than the highest multiple of the range you want,
and then divide by the proper number to bring it down to
your range)
Carlos
--
Carlos Moreno
>
> I don't know much about this, but have for some time been relying on Knuth
> (TAOCP volume 2, 3rd edition, p185) - "To compute a random integer between 0
> and k-1, one should multiply by k and truncate the result." Is he wrong?
It's a different problem. That one works given a floating-point
random number in [0,1), but it suffers from the inconvenience of
using floating-point arithmetic, which *might* introduce non-
uniformity in the distribution.
If you start off with an integer random number, the technique
posted by Andrew (and several others) is much better, in that
it gives you a *perfect* (not a figure of speech -- *truly
perfect*) uniform distribution (of course, based on the
distribution of the base random number that you use... But
anyway, you know what I'm saying: the manipulations done to
convert the base number to your output introduce *absolutely
zero* non-uniformity using this technique).
Carlos
--
Carlos Moreno
>>The obvious formula is: a + (double(rand()) / RAND_MAX) * b.
>
> I would say the obvious formula is better:
> a + (b - a + 1) * int(double(rand()) / (RAND_MAX + 1.0))
> for uniformity and correct range
Actually, both of these obvious solutions are incorrect.
They *do not* give uniformity in the required range.
I quote the nice example from Risto Lankinen in this thread:
> Example: suppose RAND_MAX = 32767 ; generate numbers
> between 0..32766 :
>
> 0: 0.0/32767.0 * 32766.0 -> to int = 0
> 1: 1.0/32767.0 * 32766.0 -> to int = 0
> 2: 2.0/32767.0 * 32766.0 -> to int = 1
> 3: 3.0/32767.0 * 32766.0 -> to int = 2
> 4: ...
>
> Already it's obvious that '0' is twice more likely than '1'.
>
> For a truly uniform distribution, some output values of
> rand() need to be discarded.
I think he missed the RAND_MAX + 1.0, but it doesn't change
the result of the example.
Carlos
--
Craig Powers
>
> Andrew Koenig wrote:
>
> > I strongly doubt it is possible to do any better.
>
> I don't know much about this, but have for some time been relying on Knuth
> (TAOCP volume 2, 3rd edition, p185) - "To compute a random integer between 0
> and k-1, one should multiply by k and truncate the result." Is he wrong?
With a suggestion like that, I strongly suspect Knuth is starting from
a floating point distribution between 0 and 1, which is not applicable
to Andrei's question.
Craig Powers
>
> The other problem is that in general, the lower bits
> carry less "randomness" than the higher bits -- when
> you take %, you are discarding the higher bits; so,
> it is better to always reduce to the required range
> by dividing instead of taking modulo. (see other
> posts -- including mine -- for a reasonable solution
> where you loop discarding values outside the range
> that gives an exact multiple of the range you want).
That is a property of a relatively poor PRNG (and I believe I've seen
it reported as an issue with naive implementations of rand() ). It's
not a property of a good one. If I care about the quality of my
random numbers, I'll use a test to make sure that the numbers are
sufficiently random in all of the bits, not just in aggregate.
Carlos Moreno
>
> Just use the obvious formula. The general advice to avoid premature
> optimization applies here as well---Unless a profiler tells you that the
> floating point operations really cost you much performance (which I
> personally think would be highly unlikely), go ahead and use them.
As already said in this thread, the problem with the "obvious"
formula that uses floating-point division is that it is as
wrong as it is obvious :-)
It does not give a uniform distribution, unless RAND_MAX is a
multiple of the range that you require.
> If you need high quality (or perhaps even cryptographic) random numbers,
> don't use rand(). Instead, I recommend some system's /dev/random which
> typically yields high-quality non-deterministic random numbers
Agreed. You could also use /dev/urandom, since /dev/random may
block, and urandom wouldn't -- at the cost of producing random
numbers from a source that has no real "true randomness" left
(still, it is said that from any conceivable point of view,
/dev/urandom would give random numbers impossible to predict)
Carlos
--
David B. Held
news:b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de...
> So there's this rand() function that returns me numbers between 0
> and RAND_MAX inclusive. How do I use that function to generate
> random numbers between two arbitrary numbers a and b with a < b?
Technically, the only way to guarantee a uniform distribution when b - a
is relatively prime to RAND_MAX and not pass over values is to take
b - a random values, and combine them linearly: r_0 +
r_1 * RAND_MAX + r_2 * RAND_MAX^2 + ... +
r_n-1 * RAND_MAX^(n-1), then divide the result by RAND_MAX.
However, that's obviously not very economical, and you won't be able
to do that with 32-bit math for most pairs a, b and most sensible values
of RAND_MAX. However, it is a general solution, and works for any
integer value of RAND_MAX > 1, even when RAND_MAX < b - a. You
get deterministic performance, but you're obviously better off just
discarding out-of-range values.
Dave
Carlos Moreno
> Beta> Then, how about
>
> Beta> result = a + (rand()+rand()) % (b - a)
>
> Beta> for the case? Anyway, a good rand() implemenation is required.
>
> If you add two uniformly distributed random numbers, the result is no
> longer uniformly distributed.
Actually, it depends. If you add two numbers that are
uniformly distributed in the full range of an unsigned
data type (i.e., if rand() RAND_MAX were 2^32-1), then
the result is uniformly distributed.
But I guess in the real case (RAND_MAX is 2^31-1), things
are exactly as you point out: it's making it worse; I
think even worse than you pointed out -- the sum of the
numbers can be negative!!
Carlos
--
Beta Jin
is drifting to discussion of algorithms, and away from discussion of the C++
language. Please try to bring any replies back on-topic. --mod]
Andrew Koenig <a...@research.att.com> wrote in message
news:<yu99wuk3...@europa.research.att.com>...
> Beta> Then, how about
>
> Beta> result = a + (rand()+rand()) % (b - a)
>
> Beta> for the case? Anyway, a good rand() implemenation is required.
>
> If you add two uniformly distributed random numbers, the result is no
> longer uniformly distributed. So the approach makes the problem
> worse, not better.
Ah... yes, I found this issue soon after I posted the comment. Here goes
the new design.
A distributed random sequence at the pace of 1/RAND_MAX in [0, RAND_MAX],
rand() + (1/RAND_MAX)*rand()
So, at the pace of 1 in [0, RAND_MAX*RAND_MAX],
RAND_MAX * (rand() + (1/RAND_MAX)*rand())
namely,
rand() * RAND_MAX + rand()
Express the generic concept in C++,
// struct rand_t
// use rand_t<Number, RandMax, Level>::get() to
// generate a distributed random number in [0, pow(RandMax, Level)-1]
// with Level > 1
template <typename Number, unsigned int RandMax, int Level>
struct rand_t
{
static Number get()
{
Number ret = rand_t<Number, RandMax, Level-1>::get() ;
ret *= RandMax ;
unsigned int plus ;
while((plus=rand()) == RandMax) {}
ret += plus ;
return ret ;
}
} ;
template <typename Number, unsigned int RandMax>
struct rand_t<Number, RandMax, 1>
{
static Number get()
{
return rand() ;
}
} ;
The above facility helps us to create a large RAND_MAX. Then, deal with the
original problem, to generate a distributed random number in [a, b]
with (0 <= a < b). The problem is equal to generating a distibuted random
number in [0, MyRandMax] with (0 < MyRandMax).
template <unsigned int MyRandMax, bool Large=(MyRandMax>RAND_MAX)>
struct my_rand_t
{
// pow(256, sizeof(unsigned int))-1 will be the max unsigned int
enum
{
RandMax = UNSIGNED_INT_MAX,
NearestMax = (RandMax - (RandMax % MyRandMax))
} ;
static unsigned int get()
{
unsigned int ret ;
while((ret = rand_t<unsigned int, 256, sizeof(unsigned int)>::get())
>= NearestMax)
{ }
return ret % RandMax ;
}
} ;
template <unsigned int MyRandMax>
struct my_rant_t<MyRandMax, false>
{
enum { NearestMax = (RAND_MAX - (RAND_MAX % MyRandMax)) } ;
static unsigned int get()
{
unsigned int ret ;
while((ret = rand()) >= NearestMax)
{ }
return ret % RandMax ;
}
} ;
template <>
struct my_rant_t<RAND_MAX, false>
{
static unsigned int get()
{
return rand() ;
}
} ;
Andrew Simmons
Andrew Koenig wrote:
> Andrew> I don't know much about this, but have for some time been
> Andrew> relying on Knuth (TAOCP volume 2, 3rd edition, p185) - "To
> Andrew> compute a random integer between 0 and k-1, one should
> Andrew> multiply by k and truncate the result." Is he wrong?
>
> I don't have my copy handy, but from the description it sounds like
> he's starting with a uniformly distributed real number in [0,1).
> That's a different problem entirely.
>
> --
Knuth was actually talking in the context of the generation of random
integers X between 0 and m-1, using the algorithm
X = (aX + c) mod m
with suitable a, c, and m.
A more extensive quote:
"It is generally best to think of X as a random fraction X/m between 0 and
1, that is, to visualize X with a radix point at its left, rather than to
regard X as a random integer between 0 and m-1. To compute a random integer
between 0 and k -1 , one should multiply by k and truncate the result."
Obviously integer overflow is a potential problem, but in my case, where I
need random integers from 0 to 127, and m, aka RAND_MAX+1, is 32768, Knuth's
method seems to work fine.
Bruce G. Stewart
<moreno_at_mo...@xx.xxx> wrote:
Although this suffers the same malady that you cited above: if the
base number has a excessively high probability of being even than odd,
so likely will the result. The distribution of values over a subrange
isn't likely to any better than distribution over the full range.
Andrei Alexandrescu
(private email appreciated, but newsgroup answers *always* considered
better because they can be further discussed publically and don't cause
duplicated communication).
I understand that the "try until you get a random number in a range
multiple of the upper limit" technique as described in Andrew Koenig's book
(I forgot about that) is the most recommended. The "wrong" generators I had
mentioned were also discussed here - those that just use modulus and no
loops.
I plan to use the algo with the loop, with a little twist and I was
wondering if that twist would improve quality.
I'm thinking of using rand() as a back-end for a function that generates
random bits (not random numbers 0 to RAND_MAX). So imagine we have a
function bool random_bit() that internally uses an accumulator word and
invokes rand() whenever the accumulator is empty (like every 32nd call for
a 32-bit rand()). Effectively, we stream each and every bit from the random
numbers generated by rand().
When asked for a random number in the range [a, b], the algorithm would be
like this:
1. Grab with random_bit() bits to fill a number with as many bits as (b -
a).
2. If the number formed with those bits is greater than (b - a), go to step
1.
3. Return a + grabbed_bits
The advantage of this algorithm is that it calls rand() less times and
therefore (1) has a longer period (2) might work better with random number
generators that don't generate good distributions for all their bits (3)
might be a bit faster if rand() uses a complicated algorithm. The
disadvantage is that random_bit() has to keep state.
So what do you all think? This sounds like a pretty nice solution to me.
Andrei
Carlos Moreno
Stefan Ram wrote:
> Carlos Moreno <moreno_at_mo...@xx.xxx> writes:
> >But then again, why all the trouble, when we have a *perfect*
> >solution? :-) (i.e., draw numbers until you get one that
> >is less than the highest multiple of the range you want,
> >and then divide by the proper number to bring it down to
> >your range)
>
> "The rand function computes a sequence of pseudo-random
> integers in the range 0 to RAND_MAX." (ANSI X3.159-1989).
>
> AFAIK ANSI X3.159-1989 does not specify the meaning of
> "pseudo-random". Are real random numbers (derived from
> a quantum device) also "pseudo-random" numbers?
>
> Does the given specification imply uniformity at all?
What does all this have to do with anything?
What I (and several others) have been trying to imply is that
you need a transformation that maps uniformly distributed random
numbers in [0,RAND_MAX] to uniformly distributed random numbers
in [0, range) -- thus, the *transformation* should be made in
such a way that it introduces no disuniformity (is this even
a real word? :-)).
That way, if you are lucky enough to count of uniformly
distributed random numbers in [0,RAND_MAX], then you will get
your uniformly distributed numbers in your range. If your
transformation introduces disuniformities, then you will never
be able to get u.d. random numbers, no matter how perfect your
"base" random generator may be.
Notice that I normally don't use this technique combined with
rand(), or even with lrand48() -- I use it with data read from
/dev/urandom, which is a rather good random number generator.
> Can one prove that the "perfect" solution will terminate?
It has a probability of terminating = 1 (which theoretically
speaking, does not mean that it necessarily terminates; it
means that from an infinite number of tries, only a finite
number of them it will not terminate).
In practical terms, it always terminates, unless the base
random number generator may fall in a loop that repeats
numbers in the "wrong range" (which is extremely unlikely).
Still, in my actual code (I'm the paranoid kind :-)), I put
a counter to try a maximum number of times -- if it exceeds
it, I simply keep whatever random number I got, and use any
other technique. I doubt that that code has ever executed
the portion corresponding to the maximum number of retries
reached.
> A binary computer can give us (pseudo-)random bits.
> This means, we're fine, if we need a random number
> with ::std::pow( 2., n ) possibilities. A number
> from the set {0, 1, 2} has 1,585 bits of information,
> each outcome has the probability 0,333 (assuming
> uniformity). I can not ask my binary pseudo-random number
> generator to give me exactly that amount of information,
> as it would be possible on a ternary computer.
> It can give me only 1 bit or 2 bits. So I must accept
> an information >= 1,585 and discard the information not
> needed. I can not use all information given, because
> the outcome only has 3 possible states, so it can
> not contain an information of 2 bits, which needs
> four states.
>
> This differs from real numbers in mathematics, because
> each real number has an unlimited amount of information.
> Therefore, a mapping exists to map a uniform distribution
> from the range [0,4) to the range [0,3) uniformly, which
> is even injective, i.e., does not erase any information.
All of the above is true, of course. But it doesn't apply
directly to the discussion, because:
1) Andrei [the OP] explicitly asked for integer random
numbers.
2) Real numbers can not be represented in actual computers;
floating-point numbers can only represent a given number
of different values (2^64 possible values, for an 8-byte
representation -- well, or whatever it is; in any case,
no more than 2^64 different values)
(you probably were not impliying that it applies directly
to the discussion... but just in case)
Daniel Frey
>
> When asked for a random number in the range [a, b], the algorithm
would be
> like this:
>
> 1. Grab with random_bit() bits to fill a number with as many bits as
(b -
> a).
>
> 2. If the number formed with those bits is greater than (b - a), go to
step
> 1.
>
> 3. Return a + grabbed_bits
>
> The advantage of this algorithm is that it calls rand() less times and
> therefore (1) has a longer period (2) might work better with random
number
> generators that don't generate good distributions for all their bits
(3)
> might be a bit faster if rand() uses a complicated algorithm. The
> disadvantage is that random_bit() has to keep state.
>
> So what do you all think? This sounds like a pretty nice solution to
me.
I'm not sure whether it's a good solution for everybody. What about
complexity? O(1) seems desirable to me for such a function, thus a loop
is IMHO inappropriate. I think the key is - as you also noticed - that
you have to keep a state somewhere. That given, what about something
like this:
template< typename T > T map_random( T rnd, T max ) // map to 0-max
{
static T state = 0;
state += rnd;
return state % ( max + 1 );
}
This is of course not complete as you have to add checks for overflows,
etc. but I think you get the idea.
Regards, Daniel
--
Daniel Frey
aixigo AG - financial training, research and technology
Schloß-Rahe-Straße 15, 52072 Aachen, Germany
fon: +49 (0)241 936737-42, fax: +49 (0)241 936737-99
eMail: danie...@aixigo.de, web: http://www.aixigo.de
Louis Lavery
> Beta> Then, how about
>
> Beta> result = a + (rand()+rand()) % (b - a)
>
> Beta> for the case? Anyway, a good rand() implemenation is required.
>
> If you add two uniformly distributed random numbers, the result is no
> longer uniformly distributed. So the approach makes the problem
> worse, not better.
Surely, if r and s are uniformly distributed in [0,n) and t is
their sum modulo n then t is uniformly distributed in [0,n).
Or am I missing something?
Louis.
P.S. Not that this helps Beta, but it doesn't make it worse.
Louis Lavery
David B. Held <dh...@codelogicconsulting.com> wrote in message
news:v4qsvst...@corp.supernews.com...
> "Andrei Alexandrescu" <andre...@hotmail.com> wrote in message
> news:b2d3oj$1b58b2$1...@ID-14036.news.dfncis.de...
> > [...]
> > So there's this rand() function that returns me numbers between 0
> > and RAND_MAX inclusive. How do I use that function to generate
> > random numbers between two arbitrary numbers a and b with a < b?
> > [...]
>
> Technically, the only way to guarantee a uniform distribution when b -
a
> is relatively prime to RAND_MAX and not pass over values is to take
> b - a random values, and combine them linearly: r_0 +
> r_1 * RAND_MAX + r_2 * RAND_MAX^2 + ... +
> r_n-1 * RAND_MAX^(n-1), then divide the result by RAND_MAX.
Could you expand a little on that?
TAI,
Louis.
Andrew Koenig
Carlos> floating-point random number in [0,1), but it suffers from the
Carlos> inconvenience of using floating-point arithmetic, which
Carlos> *might* introduce non- uniformity in the distribution.
It suffers from the additional problem of requiring a random-number
generator that returns a useful floating-point value. If all you have
is a generator that returns a 16-bit integer, you're still stuck.
--
Andrew Koenig, a...@research.att.com,
http://www.research.att.com/info/ark
John Lloyd
> Andrew Koenig wrote:
> > Beta> Then, how about
> >
> > Beta> result = a + (rand()+rand()) % (b - a)
> >
> > Beta> for the case? Anyway, a good rand() implemenation is
required.
> >
> > If you add two uniformly distributed random numbers, the result is
no
> > longer uniformly distributed.
>
> Actually, it depends. If you add two numbers that are
> uniformly distributed in the full range of an unsigned
> data type (i.e., if rand() RAND_MAX were 2^32-1), then
> the result is uniformly distributed.
>
Are you sure? Isn't the distribution of the sum of a pair of discrete
random variables given by the convolution of the two distributions
(resulting in a triangular distribution in this case)?
Carlos Moreno
Andrei Alexandrescu wrote:
If the implementation of rand() is good, then this may be ok.
Basically, what you're proposing is similar to taking rand()
modulo the highest power of 2 in which your range fits (but
then repeating the same for the next bits from your 31 bits,
and thus, not discarding those bits).
The only potential red flag I see is that not necessarily
each and every bit of the result of rand carry as much
"randomness" (or entropy, if you will) as the others, or as
the whole word.
Notice that if you are concerned with efficiency (reducing
the number of calls to rand() because they might be expensive),
this would be the wrong way to go (your method is more
expensive). I know, I know I may be insulting you with
such an obvious comment :-) But I thought I'd mention
it just in case.
If you simply want to increase the period, then just use a
linear congruential generator with higher number of bits,
and thus higher period (e.g., the rand48 family uses 48
bits, and achieves a period of 2^47 -- or was it 2^48??)
That's pretty long, I think. Boost has several random
number generators, and if I remember correctly, they even
have one that uses 64 bits!!!! I don't think you can
ever exhaust such period!
Another clever solution is a trick proposed by Knuth:
use a buffer of N numbers; fill that buffer with the
first N numbers generated by your base RNG (e.g., rand());
then, to draw a random number, use the base generator
to obtain the offset in the buffer. Return the number
contained in that buffer, and replace the contents of
that value of the buffer with another number generated
from the base generator.
I'm not 100% sure if he [Knuth] says that this only
impreoves the "randomness" (i.e., improves the statistical
characteristics of the generated sequence), or if it
increases the period of the generator (I *think* it
does both... You may want to check Knuth's to double
check).
HTH,
Carlos
--
PS: BTW, why not simply using an external source of
random data? (/dev/random and /dev/urandom use
external random events and give you *very* high
quality random numbers).
Carlos Moreno
John Lloyd wrote:
>> > If you add two uniformly distributed random numbers, the result is
>>
> no
>
>> > longer uniformly distributed.
>>
>>Actually, it depends. If you add two numbers that are
>>uniformly distributed in the full range of an unsigned
>>data type (i.e., if rand() RAND_MAX were 2^32-1), then
>>the result is uniformly distributed.
>
> Are you sure?
Yes, I am sure :-)
Notice that I'm not talking about adding as in adding two
integer or real numbers the usual way (in which case, as
you point out, you get a triangular distribution).
I'm talking about taking two numbers in the range [0,n)
and adding them *modulo n*. In that case, the distribution
remains flat: the triangle that would result if you don't
take modulo n, has a peak at n. So, cut that interval
(as the result of the overflow at n) and move the second
half (interval [n,2n) ) and add it to the first one: the
result is a constant value (i.e., uniform distribution)
But anyway, my comment was not necessarily relevant,
since rand() returns a number in an interval different
from the integers range. So, adding them together doesn't
give a nice overflow that brings the distribution back to
flat: it may result in negative values!!! :-(
Carlos
--
Craig Powers
>
> On 15 Feb 2003 04:44:45 -0500, Carlos Moreno
> <moreno_at_mo...@xx.xxx> wrote:
>
> >John Potter wrote:
> >> On 13 Feb 2003 05:12:04 -0500, Carlos Moreno
> >> <moreno_at_mo...@xx.xxx> wrote:
> >>
> >>
> >>>The other problem is that in general, the lower bits
> >>>carry less "randomness" than the higher bits
[...]
> >But then again, why all the trouble, when we have a *perfect*
> >solution? :-) (i.e., draw numbers until you get one that
> >is less than the highest multiple of the range you want,
> >and then divide by the proper number to bring it down to
> >your range)
>
> Although this suffers the same malady that you cited above: if the
> base number has a excessively high probability of being even than odd,
> so likely will the result. The distribution of values over a subrange
> isn't likely to any better than distribution over the full range.
That's not the problem he cited above. The problem is insufficient
randomness in the low-order bits, which isn't the same thing as
non-uniform distribution. Rather, the problem is that the distribution
is a little *too* uniform.
e.g.
1, 0, 1, 0, 1, 0, etc. is a uniform, but non-random, distribution.
The probability of being even or odd is 0.5.
Dividing down from a multiple of the range of interest would tend to
reduce the contribution of the "bad" bit, but there might still be
issues in the affected bucket.
e.g. if it's the lowest-order bit, the first bucket might still see
some non-randomness.
Nom De Plume
> So there's this rand() function that returns me numbers between 0 and
> RAND_MAX inclusive. How do I use that function to generate random numbers
> between two arbitrary numbers a and b with a < b?
Ok, here is my thought off the back of a napkin, just so
I can look foolish too :-)
// Warning, untested code follows, if you use this without
// testing, you deserve what you get!!!!
class RandomNum
{
public:
RandomNum(unsigned a, unsigned b);
~RandomNum(); // Dtor omitted below for brevity
unsigned rand();
private:
unsigned a, * data, *data_end, *current;
}
RandomNum::RandomNum(unsigned a, unsigned b): a(a)
{
if(b<=a) throw range_error(); // Or whatever
data = new unsigned[b-a];
data_end = &data[b-a];
unsigned *datap = data;
unsigned i = 0;
while(datap!=data_end) *datap = i++;
current = data_end;
}
unsigned RandomNum::rand()
{
if(currentp == data_end)
{
random_shuffle(data, data_end);
curentp = data;
}
return (*currentp++) + a;
}
Looking at this code, since random_shuffle is linear, then
this calls std::rand the same number of times as a simplistic
"return rand()%(b-a)". (It does do it in a peaky fashion
rather than smoothly though.) The additional overhead above
that is (b-a) unsigned swaps, one test per call, and one
deref and increment per call. This doesn't seem to be
much of an increase above the simplistic approach.
Of couse, it also uses some memory, and it doesn't meet
Andrei's requirement of up to MAX_INT for b. It is probably
limited to fairly small ranges.
Gerhard Wesp
> It does not give a uniform distribution, unless RAND_MAX is a
> multiple of the range that you require.
OK, but that's a problem with _any_ method which maps rand() to the
desired range. Let A be the range of rand() and B be the desired range.
Ideally, we'd have a uniform probability density function (PDF) on A (in
practice, this is more often than not not the case!). Then, if we have
_any_ mapping A -> B, be it floating point or modulo based or whatever,
we can achieve a uniform PDF on B only if |A| is a multiple of |B|. The
proof is left as an exercise for the reader. :-)
That said, if |A| >> |B| (_much_ larger), then in practice the deviation
from uniformity on B might be tolerable even if |A| is not a multiple of
|B|. In particular, if a linear congruential generator for A is
tolerable, then I venture the guess that the deviation will also be
tolerable.
-Gerhard
--
| voice: +43 (0)676 6253725 *** web: http://www.cosy.sbg.ac.at/~gwesp/
|
| Passts auf, seid's vuasichdig, und lossds eich nix gfoin!
| -- Dr. Kurt Ostbahn
Andrei Alexandrescu
> I'm not sure whether it's a good solution for everybody. What about
> complexity? O(1) seems desirable to me for such a function, thus a
loop
> is IMHO inappropriate. I think the key is - as you also noticed - that
> you have to keep a state somewhere. That given, what about something
> like this:
>
> template< typename T > T map_random( T rnd, T max ) // map to 0-max
> {
> static T state = 0;
> state += rnd;
> return state % ( max + 1 );
> }
The numbers generated with this method are very non-random. Think of
using
this device to generate 32-bit random numbers wit a 16-bit random
generator.
Andrei
Andrei Alexandrescu
news:3E4FFEC0...@xx.xxx...
> If the implementation of rand() is good, then this may be ok.
> Basically, what you're proposing is similar to taking rand()
> modulo the highest power of 2 in which your range fits (but
> then repeating the same for the next bits from your 31 bits,
> and thus, not discarding those bits).
Si.
> The only potential red flag I see is that not necessarily
> each and every bit of the result of rand carry as much
> "randomness" (or entropy, if you will) as the others, or as
> the whole word.
True; however, the method of discarding the upper bits suffers of this
same
problem, and in a worse way because it solely relies on the lower bits.
For
example, if you generate 1-bit numbers with the loop, you simply take
the
LSB. How random is that, I dunno.
> Notice that if you are concerned with efficiency (reducing
> the number of calls to rand() because they might be expensive),
> this would be the wrong way to go (your method is more
> expensive). I know, I know I may be insulting you with
> such an obvious comment :-) But I thought I'd mention
> it just in case.
Insult, insult, because it's a good point. But without the loop or
floating-point math I dunno how to generate quality random numbers.
By the way, now that the loop is the way of doing it in integer math,
doesn't the floating point version start to look more attractive?
Somehow I
feel it's best to stay away from floating point when it comes about
generating random integers.
> If you simply want to increase the period, then just use a
> linear congruential generator with higher number of bits,
> and thus higher period (e.g., the rand48 family uses 48
> bits, and achieves a period of 2^47 -- or was it 2^48??)
>
> That's pretty long, I think. Boost has several random
> number generators, and if I remember correctly, they even
> have one that uses 64 bits!!!! I don't think you can
> ever exhaust such period!
You're so naive :o). Check out the Mersenne twister, with a period of...
2**19937-1. Boggles my mind how you could even compute that bound, I've
never been too good at math.
> Another clever solution is a trick proposed by Knuth:
> use a buffer of N numbers; fill that buffer with the
> first N numbers generated by your base RNG (e.g., rand());
> then, to draw a random number, use the base generator
> to obtain the offset in the buffer. Return the number
> contained in that buffer, and replace the contents of
> that value of the buffer with another number generated
> from the base generator.
But this is orthogonal on scaling the numbers in the range a - b.
> PS: BTW, why not simply using an external source of
> random data? (/dev/random and /dev/urandom use
> external random events and give you *very* high
> quality random numbers).
Bicoz I wanna stay portable and only Unix has'em?
Andrei
James Kanze
news:<3E4EB87B...@monitorbm.co.nz>...
> Andrew Koenig wrote:
> > Andrew> I don't know much about this, but have for some time been
> > Andrew> relying on Knuth (TAOCP volume 2, 3rd edition, p185) - "To
> > Andrew> compute a random integer between 0 and k-1, one should
> > Andrew> multiply by k and truncate the result." Is he wrong?
> > I don't have my copy handy, but from the description it sounds like
> > he's starting with a uniformly distributed real number in [0,1).
> > That's a different problem entirely.
> Knuth was actually talking in the context of the generation of random
> integers X between 0 and m-1, using the algorithm
> X = (aX + c) mod m
> with suitable a, c, and m.
> A more extensive quote:
> "It is generally best to think of X as a random fraction X/m between 0
> and 1, that is, to visualize X with a radix point at its left, rather
> than to regard X as a random integer between 0 and m-1. To compute a
> random integer between 0 and k -1 , one should multiply by k and
> truncate the result."
> Obviously integer overflow is a potential problem, but in my case,
> where I need random integers from 0 to 127, and m, aka RAND_MAX+1, is
> 32768, Knuth's method seems to work fine.
Reread the passage you just quoted. His method is basically:
trunc( (double)rand() / (double)RAND_MAX * (double)k )
I'm not sure, but I believe that recent studies have revealed flaws in
this method. For that matter, even the simplest logic says that there
must be a flaw somewhere -- rand() returns (hopefully) RAND_MAX
different values. Dividing by a constant changes the range, but it
doesn't change the number of distinct values (unless floating point
precision is very low). Ditto multiplying by k. The trunc function
basically distributes these values into k buckets, one for each integral
value. Unless RAND_MAX is perfectly divisible by k, a uniform
distribution is impossible.
--
James Kanze mailto:jka...@caicheuvreux.com
Conseils en informatique orientée objet/
Beratung in objektorientierter Datenverarbeitung
James Kanze
[...]
> The other problem is that in general, the lower bits carry less
> "randomness" than the higher bits
I think you have it backwards. With the usual linear congruent
generators (see "Random Number Generators: Good Ones Are Hard to Find",
Park and Miller, CACM, Oct.@1988, Vol.@31, No.@10, pp.@1192-1201), it is
the high order bits which are the least random; a "small" value will be
inevitably followed by a large one -- if the multiplier is e.g. 48271
and the mod 2147483647 (one of the recommended sets of values), any time
the generator returns a value less than 44488, the next value is
garanteed to be greater than 44488.
If you are seeing little randomness in the lower bits, it is probable
that your generator is using a power of 2 as a multiplier. Throw it
out.
James Kanze
> John Potter wrote:
> > On 13 Feb 2003 05:12:04 -0500, Carlos Moreno
> > <moreno_at_mo...@xx.xxx> wrote:
> >>The other problem is that in general, the lower bits carry less
> >>"randomness" than the higher bits
> > Do you have measurements to support this?
> No, but I trust that Donald Knuth does :-) (i.e., he states that in
> his book, in the section of linear congruential random number
> generators)
Does he state it for ALL linear congruent generators, or just a certain
subset of them. It is pretty well known that using a power of two as
the modulo in the generator doesn't give very good results in general,
and in particular, results in a period of 2 for the low order bit, 4 for
the next lowest, etc.
Thomas Richter
> > It does not give a uniform distribution, unless RAND_MAX is a
> > multiple of the range that you require.
> OK, but that's a problem with _any_ method which maps rand() to the
> desired range. Let A be the range of rand() and B be the desired range.
> Ideally, we'd have a uniform probability density function (PDF) on A (in
> practice, this is more often than not not the case!). Then, if we have
> _any_ mapping A -> B, be it floating point or modulo based or whatever,
> we can achieve a uniform PDF on B only if |A| is a multiple of |B|. The
> proof is left as an exercise for the reader. :-)
And here's the proof that the above statement is incorrect. (-; The simplest
algorithm that does perform such a mapping and does, indeed, return a
uniform probability density is to run rand() as long as it returns a number
within the desired range. It's highly inefficient, though:
do {
x = rand();
} while(x >= B);
Optimizations of this technique are rather obvious. (I.e. discard
some LSBs of x, for example).
> That said, if |A| >> |B| (_much_ larger), then in practice the deviation
> from uniformity on B might be tolerable even if |A| is not a multiple of
> |B|. In particular, if a linear congruential generator for A is
> tolerable, then I venture the guess that the deviation will also be
> tolerable.
That, too.
So long,
Thomas
James Kanze
message news:<b2nc30$1eojjk$1...@ID-14036.news.dfncis.de>...
I can see several possible weaknesses, but the only way to say for sure
is to implement and test:
1) If a long period is the goal, using two linear congruent generators
with different periods is doubtlessly a better (albeit probably
slower) solution. Depending on the number of bits used each time,
your period will be multiplied by some small integer. The period of
a compound generator can easily be n*m, where n and m are the
periods of the individual generators. (For typical good generators,
n and m will be just a little under 2^31.)
2) The solution for generators which don't generate good distributions
for all of their bits is simple. Change the generator. There is
really no excuse for using rand() in C++ anyway, given its potential
(and real) inconvenients -- IMHO, the standard random number
generators for C++ are those in boost/random.
3) I'm not sure about the question of speed. Bit manipulations aren't
that cheap either. And on my machine, the standard linear congruent
generator uses one multiply and one divide per operation -- that's
not much. The only way to really tell is to measure. (On a 32 bit
machine, it's quite likely that you do speed things up some, at
least when small random values are being used.)
On the other hand, it is an excellent solution for generating really
large random numbers.
> So what do you all think? This sounds like a pretty nice solution to
> me.
I think that it's a good solution for very large random numbers. For
the rest, I doubt that it will bring any real improvement over what's
already available from boost.
(Don't feel bad about it. Lately, it seems that every time I've an idea
for some clever new class as well, it's already been done by Boost.)
--
James Kanze mailto:jka...@caicheuvreux.com
Conseils en informatique orientée objet/
Beratung in objektorientierter Datenverarbeitung
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
James Kanze
> If you simply want to increase the period, then just use a linear
> congruential generator with higher number of bits, and thus higher
> period (e.g., the rand48 family uses 48 bits, and achieves a period of
> 2^47 -- or was it 2^48??)
That's an expensive solution in CPU time. It's probably cheaper to use
a compound generator. Note that a good 32 bit generator will require
either something like a 48 bit CPU (or more) or extra operations.
> That's pretty long, I think. Boost has several random number
> generators, and if I remember correctly, they even have one that uses
> 64 bits!!!! I don't think you can ever exhaust such period!
Sure you can. It's finite. You've just got to have patience:-).
> Another clever solution is a trick proposed by Knuth: use a buffer of
> N numbers; fill that buffer with the first N numbers generated by your
> base RNG (e.g., rand()); then, to draw a random number, use the base
> generator to obtain the offset in the buffer. Return the number
> contained in that buffer, and replace the contents of that value of
> the buffer with another number generated from the base generator.
That's basically the way I implement my combined generator. I've two
linear congruent generators with slightly different periods (both prime
numbers). I use one to fill the table, and replace the entries used,
and the other to choose which entry to use.
--
James Kanze mailto:jka...@caicheuvreux.com
Conseils en informatique orientée objet/
Beratung in objektorientierter Datenverarbeitung
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
KIM Seungbeom
Unless we divide by (RAND_MAX+1), we don't get a uniform distribution.
Suppose RAND_MAX is 99 and we want random numbers from 0 to 9(incl.):
the expression inside trunc( ) will be in the range of [0,9], and if
it is truncated, there will be very low possibility that it will give 9.
--
KIM Seungbeom <musi...@bawi.org>
Carlos Moreno
Andrei Alexandrescu wrote:
> For
> example, if you generate 1-bit numbers with the loop, you simply take
> the
> LSB. How random is that, I dunno.
Actually, wouldn't it be better in that case to use the fantastic
compression algorithm (published in the Journal of Useless Algorithms,
in 1996, I believe) that compresses any block of data to 1 bit?
The compression algorithm goes like this:
1) Take the data, and count how many 1's it has;
2) Assign data = binary representation of the count, and
goto step #1 until the count is 1-bit long
Hmmm, now that you ask the question "how random is that", I
think it's trivial to prove that this results in 1 for all
values of the data that have at least one non-zero bit...
After all, it *was* a truly useless algorithm!! (as opposed
to useless as a compression algorithm :-))
> By the way, now that the loop is the way of doing it in integer math,
> doesn't the floating point version start to look more attractive?
> Somehow I
> feel it's best to stay away from floating point when it comes about
> generating random integers.
I'm inclined to believe so. Rounding errors alone may introduce
deviations from the uniform distribution. Unless absolutely
necessary, why bother with the folating-point solution? (and
given that it is not absolutely necessary, gievn that we can
do the same with just integer arithmetic), then why bother?
>>Another clever solution is a trick proposed by Knuth:
>>use a buffer of N numbers; fill that buffer with the
>>first N numbers generated by your base RNG (e.g., rand());
>>then, to draw a random number, use the base generator
>>to obtain the offset in the buffer. Return the number
>>contained in that buffer, and replace the contents of
>>that value of the buffer with another number generated
>>from the base generator.
>>
>
> But this is orthogonal on scaling the numbers in the range a - b.
Hmmm, I'm not sure I follow. I was talking about increasing
the period of the RNG -- I wasn't sure if that is what you
wanted to do with the trick of making use of every single
bit that you get from rand(); that's why this trick proposed
by Knuth came to mind. After all, even if it just increases
the quality of the random numbers, it is worth using it as
your "base" RNG -- after having applied that trick, then
you use the loop to discard the values outside the highest
range divisible by (b-a).
>>PS: BTW, why not simply using an external source of
>> random data? (/dev/random and /dev/urandom use
>> external random events and give you *very* high
>> quality random numbers).
>
> Bicoz I wanna stay portable and only Unix has'em?
I meant using external random data in general (I cited
/dev/random as an example of using external random data).
But I guess you answered my question: any attempt to
use external random data is not only platform-specific;
I guess it is *application* specific, so it may not be
a viable solution in your case.
Carlos
--
Daniel Frey
>
> "Daniel Frey" <danie...@aixigo.de> wrote in message
> > I'm not sure whether it's a good solution for everybody. What about
> > complexity? O(1) seems desirable to me for such a function, thus a
> loop
> > is IMHO inappropriate. I think the key is - as you also noticed - that
> > you have to keep a state somewhere. That given, what about something
> > like this:
> >
> > template< typename T > T map_random( T rnd, T max ) // map to 0-max
> > {
> > static T state = 0;
> > state += rnd;
> > return state % ( max + 1 );
> > }
>
> The numbers generated with this method are very non-random. Think of
> using
> this device to generate 32-bit random numbers wit a 16-bit random
> generator.
Yes, I forgot to mention that the assumption I made in my first answer
to your question (max < RAND_MAX) should still be valid for this one. I
am sure you know how you can handle this limitation in an appropriate
manner :)
Regards, Daniel
--
Daniel Frey
aixigo AG - financial training, research and technology
Schlo-Rahe-Strae 15, 52072 Aachen, Germany
fon: +49 (0)241 936737-42, fax: +49 (0)241 936737-99
eMail: danie...@aixigo.de, web: http://www.aixigo.de
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
galathaea
> wondering if that twist would improve quality.
>
> I'm thinking of using rand() as a back-end for a function that generates
> random bits (not random numbers 0 to RAND_MAX). So imagine we have a
> function bool random_bit() that internally uses an accumulator word and
> invokes rand() whenever the accumulator is empty (like every 32nd call for
> a 32-bit rand()). Effectively, we stream each and every bit from the random
> numbers generated by rand().
>
> When asked for a random number in the range [a, b], the algorithm would be
> like this:
>
> 1. Grab with random_bit() bits to fill a number with as many bits as (b -
> a).
>
> 2. If the number formed with those bits is greater than (b - a), go to step
> 1.
>
> 3. Return a + grabbed_bits
>
> The advantage of this algorithm is that it calls rand() less times and
> therefore (1) has a longer period (2) might work better with random number
> generators that don't generate good distributions for all their bits (3)
> might be a bit faster if rand() uses a complicated algorithm. The
> disadvantage is that random_bit() has to keep state.
>
> So what do you all think? This sounds like a pretty nice solution to me.
I would just mention that point 3) should really be checked out in
comparison to the previous loop approaches. This is because you may
decrease calls to rand(), but you may not as well with this method.
The problem is that you are going to be discarding _more_ values (here
only bits, though, and not the whole word) with this method
_on_average_ because a general routine that generates just enough bits
to fit the range and discards those values outside the range will fail
on average 1 out of 4 times (the probability that the second highest
bit is 1 -- looking at the range max as randomly distributed input --
and the second highest generated bit is 0... assuming you force the
high bit to 1 by the generator to match the always 1 high bit of the
range maximum) when the range has 0 as its minimum. This can be quite
a bit worse than discarding only the small range near RAND_MAX that
makes the other methods useful.
Of course, if you know the range maximum is not going to be used
equally distributed or (as you rightfully mention) you discover that
your particular use actually does call rand() less, then you may find
the technique you mention useful. And I doubt I have to mention to
you to always profile code whose purpose is to optimize, but it makes
me feel better. :)
Also, if your use is involved with cryptography then it probably is
not as effective as using a proven secure stream generator (though
discarding bits _does_ make cryptanalysis more difficult). And if you
are using this for simulations, you will want to check your generator
with DIEHARD to ensure that the stream looks random across various
vector decompositions (though I think here you may find a fairly good
distribution when you check ranges on byte boundaries because of your
points 1 and 2). In other words: build it and see (= verify)!
Nom De Plume
> Carlos Moreno <moreno_at_mo...@xx.xxx> wrote in message
> news:<3E4FFEC0...@xx.xxx>...
> > 64 bits!!!! I don't think you can ever exhaust such period!
>
> Sure you can. It's finite. You've just got to have patience:-).
Patience and life extension therapy. By my calculations, if
you call it one million times a second you'd have had to start
a few ice ages back before you could exhaust this one.
There is a important difference between infinity and huge,
but, in many practical respects, huge can become effectively
infinity. (For example, when the number of pseudo random numbers
you want is massively smaller than the period of your generator.)
I remember the realization I had when learning probability and
the factorial function, that, since there are 52! ways to shuffle
a pack of cards, and that that number is huge (8e67) that in the
whole history of the world only a tiny fraction of possible
shuffles have ever been seen. Which is amazing when you think
about how many games of poker and blackjack have been played over
history.
Louis Lavery
>
> >>Another clever solution is a trick proposed by Knuth:
> >>use a buffer of N numbers; fill that buffer with the
> >>first N numbers generated by your base RNG (e.g., rand());
> >>then, to draw a random number, use the base generator
> >>to obtain the offset in the buffer. Return the number
> >>contained in that buffer, and replace the contents of
> >>that value of the buffer with another number generated
> >>from the base generator.
> >>
> >
> > But this is orthogonal on scaling the numbers in the range a - b.
>
>
> Hmmm, I'm not sure I follow. I was talking about increasing
> the period of the RNG -- I wasn't sure if that is what you
> wanted to do with the trick of making use of every single
> bit that you get from rand(); that's why this trick proposed
> by Knuth came to mind. After all, even if it just increases
> the quality of the random numbers, it is worth using it as
> your "base" RNG -- after having applied that trick, then
> you use the loop to discard the values outside the highest
> range divisible by (b-a).
Yes, I'm trying to grasp this thread.
Is it that you're trying to improve on Andrew Koenig's nrand()?
If so, in what way?
Or have I lost it altogether?
Thanks,
Louis.
James Kanze
news:<3E52DC25...@bawi.org>...
> James Kanze wrote:
> > > A more extensive quote:
Reread what I wrote. Whether you divide by RAND_MAX or RAND_MAX+1 is
irrelevant. You don't get a random distribution unless RAND_MAX + 1 is
a multiple of the number of values in our range (and, of course, all
values in the range [0..RAND_MAX] are equally likely).
To make the point clearer, consider the simple case where RAND_MAX is
9, and we want integral random values in the range [0..6) (to simulate
the throw of a die, for example). The 10 values that rand() can return,
divided by 10 (RAND_MAX+1) give the sequence (0.0, 0.6, 1.2, 1.8, 2.4,
3.0, 3.6, 4.2, 4.8, 5.4) -- if you truncate, it is evident that 3 and 5
will only occur with half the frequency of the others.
The conclusion is simple: no matter what algorithm you use, you need to
throw away some of the values, unless RAND_MAX+1 is a multiple of the
number of values you need. There are three possible solutions: using
division, using modulo, or using floating point, basically as above (but
as you say, dividing by RAND_MAX+1). But all three, if done correctly,
basically start out by throwing out all values greater than the maximum
multiple of the number of values needed. Beyond this, I'd be worried
that rounding in the floating point operations could possibily skew the
results. And if the number of values requires approaches the modulo
divided by the multiplier in a linear congruent generator, I'd be
worried that using division would result in 0 never (or too rarely)
being followed by a second 0. That leaves modulo.
--
James Kanze mailto:jka...@caicheuvreux.com
Conseils en informatique orientée objet/
Beratung in objektorientierter Datenverarbeitung
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
James Kanze
news:<b2t7lh$mpt$1...@mamenchi.zrz.TU-Berlin.DE>...
> > Carlos Moreno <moreno_at_mo...@xx.xxx> wrote:
> > > It does not give a uniform distribution, unless RAND_MAX is a
> > > multiple of the range that you require.
> > OK, but that's a problem with _any_ method which maps rand() to the
> > desired range. Let A be the range of rand() and B be the desired
> > range. Ideally, we'd have a uniform probability density function
> > (PDF) on A (in practice, this is more often than not not the case!).
> > Then, if we have _any_ mapping A -> B, be it floating point or
> > modulo based or whatever, we can achieve a uniform PDF on B only if
> > |A| is a multiple of |B|. The proof is left as an exercise for the
> > reader. :-)
> And here's the proof that the above statement is incorrect. (-; The
> simplest algorithm that does perform such a mapping and does, indeed,
> return a uniform probability density is to run rand() as long as it
> returns a number within the desired range.
The only mapping that that performs is the identity mapping. You're
mapping the values in the range 0-n to 0-n using the identity function,
and throwing out all other values.
What has been proved, and what was claimed to be proved, is that any
method to convert the output of rand() to a given range must throw out
values, unless RAND_MAX+1 is an exact multiple of the target range.
> It's highly inefficient, though:
> do {
> x = rand();
> } while(x >= B);
> Optimizations of this technique are rather obvious. (I.e. discard some
> LSBs of x, for example).
Discarding bits, regardless of at what end, will make the results of a
good random generator non uniform. Discarding bits basically comes down
to dividing by a power of 2, or taking the value modulo a power of 2 --
to reducing the range. As proven above, this reduction can only
preserve uniform distribution if the old range (RAND_MAX + 1) is an
exact multiple of the new range. Since RAND_MAX which is a power of 2
is almost a proof of a poor rand() implementation (at least if a linear
congruent generator is used), so discarding bits will distroy uniformity
of distribution in a good implementation.
The only valid solution is to introduce an artificial cap on the output
of rand(), which is the largest multiple of the number of values we are
interested in that is still less or equal to RAND_MAX. This is what my
preferred implementation does explicitly, and Andy's does implicitly.
--
James Kanze mailto:jka...@caicheuvreux.com
Conseils en informatique orientée objet/
Beratung in objektorientierter Datenverarbeitung
[ Send an empty e-mail to c++-...@netlab.cs.rpi.edu for info ]
Nom De Plume
> like this:
>
> 1. Grab with random_bit() bits to fill a number with as many bits as (b -
> a).
>
> 2. If the number formed with those bits is greater than (b - a), go to step
> 1.
>
> 3. Return a + grabbed_bits
>
> The advantage of this algorithm is that it calls rand() less times and
It is not necessarily true that it calls rand less times.
In fact there are circumstances when it calls rand more
times.
Lets look at the alternative approach:
Andrew Koenig wrote:
>I can't resist. From ``Accelerated C++'', page 135:
>
> // return a random integer in the range [0, n)
> int nrand(int n)
> {
> if (n <= 0 || n > RAND_MAX)
> throw domain_error("Argument to nrand is out of range");
>
> const int bucket_size = RAND_MAX / n;
> int r;
>
> do r = rand() / bucket_size;
> while (r >= n);
>
> return r;
> }
Basically what Andrew has done is split the results of
rand up into a+1 groups [0..n-1], [n, 2n-1], [2n, 3n-1] ...
where a = RAND_MAX/n. However the last bucket is short,
that is to say not n in length, because RAND_MAX does
not divide evenly by n. (Of course if RAND_MAX%n == 0,
then the whole problem is trivial.) If rand chooses a
number in the last bucket, Andrew's code loops to try
again.
This leads to a simple observation, namely that if you
call this function RAND_MAX times, and assume that rand,
during those calls, generates every number from 0 to RAND_MAX
(as would be the case with a good linear congruent
generator), then you have to go round the loop an extra
k times cumulatively, where k is the size of the last
bucket (that is to say, RAND_MAX%n.) So calling Andrew's
function RAND_MAX times causes rand to be called
RAND_MAX + RAND_MAX%n times.
However, with your code, Andrei, you could potentially
have to do worse than that. In the case of your code,
you have to regenerate the number if you go over
your maximum. Lets use a specific example. (For convenience
I'll convert your code to generate a number 0 to n, rather
than a..b, since the transformation from one to the other
is trivial.)
If n is 20000, and RAND_MAX is 65535, then you will shift
in 15 bits from your random_bit function. However, if
the generated number is 20000 to 32767 then you will have
to call rand again. Similarly, if this second call is between
20000 to 32767 you'll have to call it a third time, and
so forth. This means that for RAND_MAX calls, you must call
rand 15/16*RAND_MAX times plus there is a 12767/32767 chance
of a second call, and a (12767/32767)^2 change of a third
call and so forth. That is to say, the total number of calls
is:
15/16 * RAND_MAX + 15/16*RAND_MAX * (sum n->infinity of (12767/32767)^n)
= 15/16 * RAND_MAX * (1/(1-12767/32767))
= 100659
[Remember your high school math? Sum to infinity of a geometric
series is 1/(1-r) -- OK, I admit it, I had to look it up on
the web :-)]
By my calculation, the "extra bucket" at the end using
Andrew's algorithm is 5535 in size, (65535 = 3*20000+5535)
so his algorithm will take RAND_MAX + 5535, which in this
case is 71070, or 29589 less than Andrei's (i.e., Andrew
calls rand on average 1.08 times per call, whereas, Andrei
calls it 1.54 times per call.)
So this is a specific case where your algorithm needs more
calls to rand than Andrew's. You gain because you don't
always use all the bits from every call to rand, but you loose
because you hit your excess more often.
I'm shooting from the hip here, math was never my strong
suit, but I believe this seems a reasonable reflection of
reality.
However, all this leads me to a different approach to this
function. In the case of Andrew's code, what he does is
throw away the extra, short, bucket at the end. However,
a different approach would be to fill in the rest of the
short bucket to a full n. How would that work?
If the short bucket is k in length (k>=0, k<n), then we know
that a full cycle of the random number generator generates
the numbers 0..k-1 one more time than it generates the numbers
k..n-1. (because these numbers don't exist in the last
bucket.) Again I am assuming a linear congruent generator
here with a cycle of RAND_MAX, and that in that cycle it
generates every number between 0 and RAND_MAX.
In that case a trivial solution would be this (using unsigned
for convenience.)
class Random
{
public:
Random(unsigned n)
: n(n), k(RAND_MAX%n), rand_calls(0), overflow(false), counter(0)
{;};
unsigned n, k; // as defined in the text
unsigned rand_calls; // Number of times rand has been
// called this cycle
bool overflow; // have we called it RAND_MAX time?
unsigned counter; // Yes! So emit this extra number to
// balance the distribution
unsigned nrand()
{
unsigned r;
if(!overflow && rand_calls < RAND_MAX)
{
rand_calls++;
return rand();
}
if(overflow) // In overflow, so emit a number to
// fill up the short bucket
{
r = counter++;
if(counter == n)
overflow = false;
}
else // about to go into overflow because
// rand_calls == RAND_MAX
{
r = counter = k;
overflow = true;
rand_calls = 0;
}
return r;
}
};
In essence what this does is for the first RAND_MAX calls
it simply calls and returns rand(). Then for the next
k..n-1 calls it simply returns the numbers k..n-1 in
sequence. (That is to say, filling in the last bucket.)
What this means is that over RAND_MAX + (n-k) calls to
this function, the distribution is even. Of course it
also means that at one point in the random number sequence
you get a little block of sequential numbers. However,
that problem is fairly trivial to fix. You can simply
spread the k..n-1 sequence throughout the rest of the
calls evenly so that they don't all come in one block.
(I'll leave it up to the reader to write the code, my
head hurts after all the math and stuff already.)
Interestingly, the performance of this is actually
*better* than the dumb "return rand()%n" approach,
since rand is called on average less times. It is also
better (in terms of calls to rand, but giving an even
distribution) than both Koenig and Alexandrescu. By
my quick calculations, using the n=20000 and RAND_MAX
of 65565, then my code calls rand an average of
RAND_MAX/(RAND_MAX+RAND_MAX%n), or 0.92 times per
call to my function, which is certainly a weird and
unexpected result.
Andrew Koenig tells us that he strongly doubts that
it is possible to do better than his code. Did I
change your mind Andrew?
Carlos Moreno
Louis Lavery wrote:
> [snip]
>
>> >>Another clever solution is a trick proposed by Knuth:
>> >>use a buffer of N numbers; fill that buffer with the
>> >>first N numbers generated by your base RNG (e.g., rand());
>> >>then, to draw a random number, use the base generator
>> >>to obtain the offset in the buffer. Return the number
>> >>contained in that buffer, and replace the contents of
>> >>that value of the buffer with another number generated
>> >>from the base generator.
>> >>
>> >
>> > But this is orthogonal on scaling the numbers in the range a - b.
>>
>>
>>Hmmm, I'm not sure I follow. I was talking about increasing
>>the period of the RNG -- I wasn't sure if that is what you
>>wanted to do with the trick of making use of every single
>>bit that you get from rand(); that's why this trick proposed
>>by Knuth came to mind. After all, even if it just increases
>>the quality of the random numbers, it is worth using it as
>>your "base" RNG -- after having applied that trick, then
>>you use the loop to discard the values outside the highest
>>range divisible by (b-a).
>>
>
> Yes, I'm trying to grasp this thread.
>
> Is it that you're trying to improve on Andrew Koenig's nrand()?
Not really. Well, it depends on how you look at it, I guess.
The bottom line is that getting a random number between a and b
requires two (independent) steps: obtain random numbers (in the
range determined by the particular random number generator -- e.g.,
0 to RAND_MAX), and step 2 is doing soemthing to obtain a random
number (uniformly distributed, presumably) between a and b given
random numbers in a different range (e.g., 0 to RAND_MAX).
The trick of looping discarding random values outside the range
of interest (the key idea in Andrew's nrand function) addresses
step 2 only -- it is what guarantees that you will have uniform
distribution in the range a-b provided that the base random
numbers are uniformly distributed in their range. But that's
it: it uses a random generator (e.g., rand(), which is the only
one in the standard C++ library) and then does some tricks with
it.
But obtaining the "base" random numbers is independent of the
conversion to a random number between a and b; so, you can
combine the trick with other random number generators (such
as lrand48, or /dev/random or /dev/urandom, on Unix systems),
or some of the boost generators, or with some elaborate trick
that uses rand() yet it improves the quality of the sequence
(which is what the example I cited from Knuth [The Art of
Computer Programming - Vol. 2, Seminumerical Algorithms] does.
So, you could summarize the whole thing with the following
fragment of (untested) pseudo-code:
int better_rand() // according to Knuth
{
const int size = 256;
const int ratio = (RAND_MAX + 1) / size;
static int buf[size];
if (first_time)
fill buf with 256 values from std::rand()
int offset = rand() / ratio;
int rnd = buf[offset];
buf[offset] = rand();
return rnd;
}
int nrand (int n)
{
// code the exact Andrew Koenig's function, but use
// the function better_rand() instead of rand()
}
HTH,
Carlos
--
Andrew Simmons
> Unless we divide by (RAND_MAX+1), we don't get a uniform distribution.
Indeed, which is why I equated Knuth's m to RAND_MAX+1, and why
James's code is not exactly an implementation of Knuth's method. Nor,
in any case, do you necessarily need floating point arithmetic. On my
system, where RAND_MAX=32767, I can use Knuth's method for values of k
up to 32768 using unsigned 32 bit integer arithmetic. In this case I
am effectively regarding the least significant 15 bits as the
fractional part and the rest as the integer part.
And in a desperate attempt to bring this thread back on topic, James,
isn't it politically incorrect to write (double) rather than static
cast<double> these days?
Carlos Moreno
James Kanze wrote:
>
>>It's highly inefficient, though:
>>
>
>>do {
>> x = rand();
>>} while(x >= B);
>>
>
>>Optimizations of this technique are rather obvious. (I.e. discard some
>>LSBs of x, for example).
>
> Discarding bits, regardless of at what end, will make the results of a
> good random generator non uniform. Discarding bits basically comes down
> to dividing by a power of 2, or taking the value modulo a power of 2 --
> to reducing the range. As proven above, this reduction can only
> preserve uniform distribution if the old range (RAND_MAX + 1) is an
> exact multiple of the new range. Since RAND_MAX which is a power of 2
No it's not. RAND_MAX+1 is a power of 2. But now that I think about
it, I don't think the Standard specifies anything about this, which
makes your statement potentially correct (i.e., dividing a number
generated by rand() by a power of 2 *is not guaranteed* to preserve
uniform distribution).
Assuming that RAND_MAX+1 is a power of two, then there is absolutely
nothing wrong with discarding bits -- other than the fact that it will
be much slower than the trick of discarding values outside the range
given by the highest multiple of the required range (e.g., if you
want a random number between 0 and 129, then you would have to discard
all but 8 of the bits, and then discard the resulting numbers that
fall in the interval 130 to 255 -- that means discard half of the
values: the probability of having to loop N times is now (126/256)^N,
instead of (126/(RAND_MAX+1))^N
I guess this goes also as additional comments/feedback for Andrei's
solution of drawing random bits to avoid discarding bits from the
base random generator.
Carlos
--
Andrei Alexandrescu
news:b59f4084.03022...@posting.google.com...
> However, with your code, Andrei, you could potentially
> have to do worse than that.
Aha, now I understand (... why math has never been my stronghold :o)).
Cool. So the idea is that under certain conditions my method would need to
"waste" more values while searching for an eligible random number than
Andy's.
For example, if I wanna generate numbers in the range [0..8] (to use a
Haskell notation :o)), I would need to draw four bits at a time, which would
force me to waste almost half of the numbers I encounter. With Andy's
method, however, using a 15-bit generator, any number between 0 and 32669
inclusive is fair game (given that 32670 is the closest number to 32767 that
divides 9). So only 8 numbers out of 32768 will be wasted, which is so much
better.
However, if the range is a multiple of 2, I feel uneasy about Andy's method
practically considering only the lower bits of the random number. At the
extreme, again, if all I want is a [0..1] range, only the LSB ever matters.
I'm not sure whether that is a valid concern at all, and I'm not qualified
to assess Nom de Plume's solution.
Andrei
Louis Lavery
in
message news:<b2nc30$1eojjk$1...@ID-14036.news.dfncis.de>...
> > When asked for a random number in the range [a, b], the algorithm
would
be
> > like this:
[snip - details of revolutionary new RNG]
> Interestingly, the performance of this is actually
> *better* than the dumb "return rand()%n" approach,
> since rand is called on average less times. It is also
> better (in terms of calls to rand, but giving an even
> distribution) than both Koenig and Alexandrescu. By
> my quick calculations, using the n=20000 and RAND_MAX
> of 65565, then my code calls rand an average of
> RAND_MAX/(RAND_MAX+RAND_MAX%n), or 0.92 times per
> call to my function, which is certainly a weird and
> unexpected result.
Although my knowledge of maths is not up to following all
your logic, and far be it from me to advise one of your
obvious abilities, but you seem to be onto a winner here.
No doubt you've already twigged that the performance of
your nrand() can be improved further by having it throw
away the first number it gets and calling it self!
This second call reduces the average number of calls
from 0.92 to 0.92*0.92 = 0.85. Of course, this recursion
must be stopped at some level as, eventually, the overhead
of the recusive calls will outweigh the savings in the
reduced number of calls to rand().
I'm pretty sure, by posting it here, that you've not
compromised your copyright to the code. So if I were
you I'd knock out the code, profile it to get the
optimum recusion depth, make it available via the Web
(you don't need to charge much as demand will be so
great) and wait for the dosh to roll in.
Good luck,
Louis.
Risto Lankinen
"Andrew Koenig" <a...@research.att.com> wrote in message
news:yu99adh1...@europa.research.att.com...
>
> // return a random integer in the range [0, n)
> int nrand(int n)
> {
> if (n <= 0 || n > RAND_MAX)
> throw domain_error("Argument to nrand is out of range");
>
> const int bucket_size = RAND_MAX / n;
> int r;
>
> do r = rand() / bucket_size;
> while (r >= n);
>
> return r;
> }
>
> this program, but it terminates with probability 1 and its average
> performance is O(1). I strongly doubt it is possible to do any better.
I think it is possible to get a flat distribution even with
a deterministic rejection algorithm:
int nrand( int n )
{
static unsigned long rejector;
++rejector;
int bias = 1 + (RAND_MAX+1) / n;
int pivot = (RAND_MAX+1) % n;
int result = rand() % n;
if( result < pivot && rejector%bias == 0 )
{
result = rand() % n;
}
return result;
}
The 'bias' is the size of biased bucket. The 'pivot' is the
bucket below which the outcomes are biased upward
by one. The algorithm rejects biased outcomes with
probability '1/bias', hence returning 'bias-1' values for
every bucket whether below or above 'pivot'.
One remaining source of uneven distribution is when
'rejector' wraps around ULONG_MAX. Duh, I can
live with that... :-)
Cheers!
- Risto -
Andrea Griffini
<andre...@hotmail.com> wrote:
>However, if the range is a multiple of 2, I feel uneasy about Andy's
method
>practically considering only the lower bits of the random number. At
the
>extreme, again, if all I want is a [0..1] range, only the LSB ever
matters.
>I'm not sure whether that is a valid concern at all, and I'm not
qualified
>to assess Nom de Plume's solution.
After reading about your concerns I investigated a bit on
this problem considering really "precious" the random bits
used to generate the desired output. Note that I personally
don't think random-like bits are really so interesting,
but probably its me... I don't think diamonds or gold are
really that precious either :-).
Anyway... the idea I came to was use Andrew's approach just
using less bits than a full rand() extraction...
template<int N, // Output will be 0..N-1
int M, // Number of random bits to use
int S, // Number of random bits provided by the source
class R> // Source of randomness
class RndGen
{
private:
R r; // Random source
unsigned acc; // Currently available random bits
int acc_sz; // Number of random bits currently in acc
public:
RndGen(const R& r) : r(r),acc(0),acc_sz(0)
{ }
int operator()()
{
const int bucket_size = (1<<M)/N;
int bucket;
do
{
// Make sure we've enough random bits available
while (acc_sz<M)
{
unsigned int s = r() & ((1<<S)-1);
acc |= ( s << acc_sz );
acc_sz += S;
}
// Draw the random number 0..2^M-1
unsigned int draw = acc & ((1<<M)-1);
acc >>= M; acc_sz -= M;
// Compute the bucket number and
// stop if we got a valid result
bucket = draw/bucket_size;
} while (bucket>=N);
return bucket;
}
};
For example
template<unsigned N> struct LOG2
{ enum { value = 1 + LOG2< N>>1 >::value }; };
template<> struct LOG2<1>
{ enum { value = 0 }; };
...
RndGen<10,4,LOG2<1+RAND_MAX>::value,int(*)()> rg(rand);
will create a random generator that uses 4 bits at a time
to generate a number between 0 and 9, using as source the
standard rand() function.
What I think an interesting problem is however considering
what is the optimal number of bits to use (4 in the previous
example) so that the number of bits used per extraction will
be, on average, the minimum.
The program just consider that the toal number of bits
used will be on average
k = b / ( 1 - s )
where b is the number of used bits and s is the
probability of the extraction needing to be repeated
(i.e. ((2^b)%n)/(2^b)). Here are some excerpt from
the resulting table:
n=2, b=1, i=0, k=1
n=3, b=2, i=0, k=2.66666666666667
n=4, b=2, i=0, k=2
n=5, b=4, i=1, k=4.26666666666667
n=6, b=3, i=0, k=4
n=7, b=3, i=0, k=3.42857142857143
n=8, b=3, i=0, k=3
n=9, b=5, i=1, k=5.92592592592593
n=10, b=5, i=1, k=5.33333333333333
n=11, b=4, i=0, k=5.81818181818182
n=12, b=4, i=0, k=5.33333333333333
...
n=337, b=10, i=1, k=10.1285855588526
n=338, b=10, i=1, k=10.0986193293886
n=339, b=10, i=1, k=10.0688298918387
n=340, b=10, i=1, k=10.0392156862745
n=341, b=10, i=1, k=10.0097751710655
n=342, b=12, i=3, k=13.0653907496013
n=343, b=12, i=3, k=13.0272992313809
n=344, b=12, i=3, k=12.9894291754757
n=345, b=12, i=3, k=12.9517786561265
n=346, b=12, i=3, k=12.9143457698371
...
n=370, b=12, i=3, k=12.0766584766585
n=371, b=12, i=3, k=12.0441068365597
n=372, b=12, i=3, k=12.0117302052786
n=373, b=11, i=2, k=12.0793565683646
n=374, b=11, i=2, k=12.0470588235294
n=375, b=11, i=2, k=12.0149333333333
n=376, b=11, i=2, k=11.9829787234043
n=377, b=11, i=2, k=11.9511936339523
where n is the desired output range width, b is
the best number of input bits to use, i is the
number of bits "in excess" used above the minimum
usable number of bits for a single extraction and
k is the number of bits that will be used on average.
For the example it shows that for extracting a number
between 0 and 9 the optimal first two template
parameters are 10 and 5 (not 4).
In a little experiment I found that using 4 to extract
1000 random numbers I used 6465 random bits from the
source, and using 5 bits the total was 5370, an average
of 5.37 bits per generated number that looks to me
close enough to the expected value.
After a quick watch at the table I've the impression
that the number of bits in excess (i) that is worth
using is generally quite low, and shows up often when
the desired range has the form n=2^t+j for small j.
Just as a sensation I don't expect an exact direct
formula for i(n) however being simple, because looking
for when this number of "excess" bits reaches 3 for
the first time I found n=342 and probably chinese
mathematicians reading this message will share with me
the fear for such a coincidence (p=341 is the first
composite number for which p exactly divides 2^p-2:
those numbers are extremely rare and for this property
they gained the name of "pseudo-primes").
For n in the range (2..1000000) these are the first
times i(n) reachs an higher value
n=2, b=1, i=0, k=1
n=5, b=4, i=1, k=4.26666666666667
n=33, b=8, i=2, k=8.86580086580087
n=342, b=12, i=3, k=13.0653907496013
n=699051, b=24, i=4, k=25.0434663192143
The wasted number of bits when using as b the smallest
allowable value (that is the smalles b so that 2^b>=n)
increases as n increase and here is a table in which "e"
shows the difference between the best k and the k for
the lowest allowable b value. I computed the table only
for n=2^m+1, but for all b<27 these are the values in
which e(n) gets to a new peak anyway.
m= 1, b= 2, i=0, e=0, k=2.66666666
m= 2, b= 4, i=1, e=0.53333333, k=4.26666666
m= 3, b= 5, i=1, e=1.18518518, k=5.92592592
m= 4, b= 6, i=1, e=1.88235294, k=7.52941176
m= 5, b= 8, i=2, e=2.77056277, k=8.86580086
m= 6, b= 9, i=2, e=3.65714285, k=10.1274725
m= 7, b=10, i=2, e=4.53599114, k=11.3399778
m= 8, b=11, i=2, e=5.40744858, k=12.5225125
m= 9, b=12, i=2, e=6.27346143, k=13.6875522
m=10, b=13, i=2, e=7.13588850, k=14.8426480
m=11, b=14, i=2, e=7.99609565, k=15.9921913
m=12, b=16, i=3, e=8.93115287, k=17.0625010
m=13, b=17, i=3, e=9.86546238, k=18.1311200
m=14, b=18, i=3, e=10.7993408, k=19.1988281
m=15, b=19, i=3, e=11.7329752, k=20.2660481
m=16, b=20, i=3, e=12.6664733, k=21.3330078
m=17, b=21, i=3, e=13.5998962, k=22.3998291
m=18, b=22, i=3, e=14.5332778, k=23.4665771
m=19, b=23, i=3, e=15.4666371, k=24.5332865
m=20, b=24, i=3, e=16.3999843, k=25.5999755
m=21, b=25, i=3, e=17.3333250, k=26.6666539
m=22, b=26, i=3, e=18.2666623, k=27.7333267
m=23, b=27, i=3, e=19.1999977, k=28.7999965
m=24, b=28, i=3, e=20.1333321, k=29.8666648
m=25, b=29, i=3, e=21.0666660, k=30.9333324
m=26, b=30, i=3, e=21.9999996, k=31.9999995
m=27, b=32, i=4, e=22.9677417, k=33.0322578
m=28, b=33, i=4, e=23.9354837, k=34.0645160
m=29, b=34, i=4, e=24.9032257, k=35.0967741
m=30, b=35, i=4, e=25.8709677, k=36.1290322
The line m=27 for example tells that using 32 random bits
instead of 28 brings on average a saving of almost 23 bits
per computation. Note that you were right in your feeling...
using a full RAND_MAX extraction (supposing that giving 32
random bits) is always a waste of bits if n is below 2^27+1.
Using the smallest number of bits shouldn't be too
complex with a bit of template magic, and that's
not too bad compared to the extraction of always the
maximum number of bits to do the bucket computation
(if we suppose the cost of a random extraction being
high enough to justify the bit manipulation code).
Anyway here is the value for the optimal i(n) (number
of additional bits to read before using the bucket
approach) for all values of n less than 10,000,000.
A row in the table was printed only if the i(n) value
on the row was different from the previous row.
The "e" value is the wasted number of random bits on
average if the minimum number of bits is used.
The table has just 112 entries for representing all
numbers less than 10,000,000 and I wouldn't be surprised
in finding that extending the table for higher values
(say up to 2^32) is still reasonable.
I've no idea if it's simple to bury this table in a
template metaprogram to generate optimal random
number range reduction code (optimal in terms of source
random bits used). Also checking the real gain obtained
by saving random bits is worth investigation, for
example I've seen generators like the Mersenne Twister
that on average really do just a handful of operations
to produce a plethora of bits. "precious" true random
bits could be used to just initialize the internal state
of such a generator, and then it should be simple to
quickly get good enough noise from that.
n=1, b=1, i=0, e=0, k=1
n=5, b=4, i=1, e=0.53333333, k=4.26666666
n=6, b=3, i=0, e=0, k=4
n=9, b=5, i=1, e=1.18518518, k=5.92592592
n=11, b=4, i=0, e=0, k=5.81818181
n=17, b=6, i=1, e=1.88235294, k=7.52941176
n=22, b=5, i=0, e=0, k=7.27272727
n=33, b=8, i=2, e=2.77056277, k=8.86580086
n=37, b=7, i=1, e=2.30630630, k=8.07207207
n=43, b=6, i=0, e=0, k=8.93023255
n=65, b=9, i=2, e=3.65714285, k=10.1274725
n=74, b=8, i=1, e=2.88288288, k=9.22522522
n=86, b=7, i=0, e=0, k=10.4186046
n=129, b=10, i=2, e=4.53599114, k=11.3399778
n=147, b=9, i=1, e=3.48299319, k=10.4489795
n=171, b=8, i=0, e=0, k=11.9766081
n=257, b=11, i=2, e=5.40744858, k=12.5225125
n=293, b=10, i=1, e=4.07736063, k=11.6496018
n=342, b=12, i=3, e=0.40829346, k=13.0653907
n=373, b=11, i=2, e=0.27453083, k=12.0793565
n=410, b=9, i=0, e=0, k=11.2390243
n=513, b=12, i=2, e=6.27346143, k=13.6875522
n=586, b=11, i=1, e=4.65984072, k=12.8145620
n=683, b=13, i=3, e=0.81778251, k=14.1748968
n=745, b=12, i=2, e=0.54979865, k=13.1951677
n=820, b=10, i=0, e=0, k=12.4878048
n=1025, b=13, i=2, e=7.13588850, k=14.8426480
n=1171, b=12, i=1, e=5.24679760, k=13.9914602
n=1366, b=14, i=3, e=1.22667376, k=15.2652735
n=1490, b=13, i=2, e=0.82469798, k=14.2947651
n=1639, b=11, i=0, e=0, k=13.7449664
n=2049, b=14, i=2, e=7.99609565, k=15.9921913
n=2341, b=13, i=1, e=5.83226541, k=15.1638900
n=2731, b=15, i=3, e=1.63616390, k=16.3616390
n=2979, b=14, i=2, e=1.09996643, k=15.3995300
n=3277, b=12, i=0, e=0, k=14.9990845
n=4097, b=16, i=3, e=8.93115287, k=17.0625010
n=4370, b=15, i=2, e=8.30179797, k=16.0679960
n=4682, b=14, i=1, e=6.41549195, k=16.3303431
n=5462, b=16, i=3, e=2.04520488, k=17.4524150
n=5958, b=15, i=2, e=1.37495803, k=16.4994964
n=6554, b=13, i=0, e=0, k=16.2490082
n=8193, b=17, i=3, e=9.86546238, k=18.1311200
n=8739, b=16, i=2, e=9.10623968, k=17.1411570
n=9363, b=15, i=1, e=6.99946598, k=17.4986649
n=10923, b=17, i=3, e=2.45447055, k=18.5448886
n=11916, b=16, i=2, e=1.64994964, k=17.5994629
n=13108, b=14, i=0, e=0, k=17.4989319
n=16385, b=18, i=3, e=10.7993408, k=19.1988281
n=17477, b=17, i=2, e=9.91029843, k=18.2135214
n=18725, b=16, i=1, e=7.58315976, k=18.6662394
n=21846, b=18, i=3, e=2.86354897, k=19.6357644
n=23832, b=17, i=2, e=1.92494125, k=18.6994293
n=26215, b=15, i=0, e=0, k=18.7495708
n=32769, b=19, i=3, e=11.7329752, k=20.2660481
n=34953, b=18, i=2, e=10.7141426, k=19.2854567
n=37450, b=17, i=1, e=8.16647975, k=19.8328793
n=43691, b=19, i=3, e=3.27270230, k=20.7271145
n=47663, b=18, i=2, e=2.19997901, k=19.7998111
n=52429, b=16, i=0, e=0, k=19.9999237
n=65537, b=20, i=3, e=12.6664733, k=21.3330078
n=69906, b=19, i=2, e=11.5177033, k=20.3568710
n=74899, b=18, i=1, e=8.74991655, k=20.9997997
n=87382, b=20, i=3, e=3.68179009, k=21.8180153
n=95326, b=19, i=2, e=2.47497639, k=20.8998006
n=104858, b=17, i=0, e=0, k=21.2499189
n=131073, b=21, i=3, e=13.5998962, k=22.3998291
n=139811, b=20, i=2, e=12.3213521, k=21.4284385
n=149797, b=19, i=1, e=9.33330663, k=22.1666032
n=174763, b=21, i=3, e=4.09090128, k=22.9090472
n=190651, b=20, i=2, e=2.74998819, k=21.9999055
n=209716, b=18, i=0, e=0, k=22.4999141
n=262145, b=22, i=3, e=14.5332778, k=23.4665771
n=279621, b=21, i=2, e=13.1249655, k=22.4999409
n=299594, b=20, i=1, e=9.91663829, k=23.3332665
n=349526, b=22, i=3, e=4.49999141, k=23.9999542
n=381301, b=21, i=2, e=3.02499495, k=23.0999614
n=419431, b=19, i=0, e=0, k=23.7499660
n=524289, b=23, i=3, e=15.4666371, k=24.5332865
n=559241, b=22, i=2, e=13.9285598, k=23.5714089
n=599187, b=21, i=1, e=10.4999874, k=24.4999707
n=699051, b=24, i=4, e=4.95651937, k=25.0434663
n=729445, b=23, i=3, e=4.70454012, k=24.0454273
n=762601, b=22, i=2, e=3.29999881, k=24.1999913
n=838861, b=20, i=0, e=0, k=24.9999940
n=1048577, b=24, i=3, e=16.3999843, k=25.5999755
n=1118482, b=23, i=2, e=14.7321305, k=24.6428365
n=1198373, b=22, i=1, e=11.0833293, k=25.6666574
n=1398102, b=25, i=4, e=5.41304089, k=26.0869440
n=1458889, b=24, i=3, e=5.09658863, k=25.0908978
n=1525202, b=23, i=2, e=3.57499872, k=25.2999909
n=1677722, b=21, i=0, e=0, k=26.2499937
n=2097153, b=25, i=3, e=17.3333250, k=26.6666539
n=2236963, b=24, i=2, e=15.5357082, k=25.7142757
n=2396746, b=23, i=1, e=11.6666624, k=26.8333237
n=2796203, b=26, i=4, e=5.86956451, k=27.1304315
n=2917777, b=25, i=3, e=5.48863579, k=26.1363609
n=3050403, b=24, i=2, e=3.84999988, k=26.3999992
n=3355444, b=26, i=4, e=0.13157891, k=27.3684145
n=3532046, b=22, i=0, e=0, k=26.1249961
n=4194305, b=26, i=3, e=18.2666623, k=27.7333267
n=4473925, b=25, i=2, e=16.3392830, k=26.7857098
n=4793491, b=24, i=1, e=12.2499981, k=27.9999958
n=5592406, b=27, i=4, e=6.32608620, k=28.1739096
n=5835554, b=26, i=3, e=5.88068120, k=27.1818153
n=6100806, b=25, i=2, e=4.12499987, k=27.4999991
n=6710887, b=27, i=4, e=0.32894733, k=28.4210500
n=7064091, b=23, i=0, e=0, k=27.3124997
n=8388609, b=27, i=3, e=19.1999977, k=28.7999965
n=8947849, b=26, i=2, e=17.1428562, k=27.8571414
n=9586981, b=28, i=4, e=12.9629623, k=29.0370357
n=9942054, b=27, i=3, e=12.4615383, k=28.0384613
One last thought is that probably there's a way to
do even better than using b(n) bits... it's just
an impression, but I think some investigation on
the line of arithmetic compressors could save even
more bits. For example a simple quasi-arithmetic
approach that I used for data compression with good
results is just considering the extraction of the
numbers in small groups.
This would mean extracting a random number between 0
and N^R-1 if this fits better, and then use that number
as the source of R random numbers between 0 and N-1
using a repeated modulo/division process.
The code is really straightforward, and much simpler
than a true arithmetic approach, with just a very
tiny extra waste.
HTH
Andrea
PS:
I wrote most of this on a sunday morning, waking
up after just a few sleep hours. If you find this
useful then fine. If it turns out it's just a
pile of rubbish I'll pretend it wasn't me
writing this message :-)
Andrea Griffini
>Anyway... the idea I came to was use Andrew's approach just
>using less bits than a full rand() extraction...
>
>template<int N, // Output will be 0..N-1
> int M, // Number of random bits to use
> int S, // Number of random bits provided by the source
> class R> // Source of randomness
....
While walking to work I realized that my code example
is very bad for at least two reasons. First I think that
no one guarantees that MAX_RAND has the form 2^M-1,
even if to be honest I think this is probably by far
the most common case. The second problem is that I'm using
an unsigned int to hold all the needed bits, and this
easily becomes insufficient when N+S is bigger than the
number of bits in an unsigned int. Even in the what I think
is a common case of having S being the number of bits
in an unsigned int the code code doesn't work correctly.
A possible fix to this big problem is
template<int N, // Output will be 0..N-1
int M, // Number of random bits to use
int S, // Number of random bits provided by the source
class R> // Source of randomness
class RndGen
{
private:
R r; // Random source
unsigned acc; // Currently available random bits
int acc_sz; // Number of random bits currently in acc
public:
RndGen(const R& r) : r(r),acc(0),acc_sz(0),count(0)
{ }
int operator()()
{
const int bucket_size = (1<<M)/N;
int bucket;
do
{
// Get M random bits
unsigned draw=0, draw_sz=0;
do
{
// Refill if there are no more available random bits
if (acc_sz==0)
{
acc = r() & ((1<<S)-1);
acc_sz = S;
}
// Compute the number of needed bits
int bsz = M-draw_sz;
if (bsz>acc_sz) bsz=acc_sz; // We don't have that many
// Move the bits
draw |= (acc & ((1<<bsz)-1)) << draw_sz;
draw_sz += bsz;
acc >>= bsz;
acc_sz -= bsz;
} while (draw_sz<M);
// Compute the bucket number and
// stop if we got a valid result
bucket = draw/bucket_size;
} while (bucket>=N);
return bucket;
}
};
In this version the bits are first extracted and then
inserted in the acc temporary storage, so at least in
that point there is no risk of overflowing.
I'm actually not sure at all that this code works
portably and/or in limit cases (all those shift opertions
look scary)... but the original code was totally broken
and failing almost always (of course except with that
example I used to test the code and BCC55 where
RAND_MAX is just 32767 and ints have 32 bits).
The more I think to the complexity of doing that bit
manipulation, the more I'm wondering if it's a bad
idea trying to save on the consumed noise.
Sorry for the first horrible version ...
Andrea
PS:
Once a time I found a nice quiz on how to extract
a random number between 1 and 10 using a dice.
A proposed solution was to use a dice shaped as
an icosahedron.
Louis Lavery
Risto Lankinen <rlan...@hotmail.com> wrote in message
news:EOq6a.20$7b....@news1.nokia.com...
>
[snip]
> I think it is possible to get a flat distribution even with
> a deterministic rejection algorithm:
>
>
> int nrand( int n )
> {
> static unsigned long rejector;
>
> ++rejector;
>
> int bias = 1 + (RAND_MAX+1) / n;
> int pivot = (RAND_MAX+1) % n;
> int result = rand() % n;
>
> if( result < pivot && rejector%bias == 0 )
> {
> result = rand() % n;
> }
>
> return result;
> }
>
Here's another deterministic distribution flattening algorithm...
int nrand(int n)
{
static int i;
return (rand() + i++)%n;
}
...well, it needs a little work to avoid problems with overflow,
but you get the gist.
The trouble is it suffers from the same malady that any other
deterministic distribution flattening algorithm, that depends on
previous state, does, which is just that - it depends on previous
state.
In the case of your nrand, on the ith call
prob(nrand < pivot) = {if i%bias == 0 prob(rand()%n < pivot)^2
{otherwise prob(rand()%n < pivot)
>
> The 'bias' is the size of biased bucket. The 'pivot' is the
> bucket below which the outcomes are biased upward
> by one. The algorithm rejects biased outcomes with
> probability '1/bias',
No it doesn't.
Louis.
Mary K. Kuhner
Nom De Plume <nde_...@ziplip.com> wrote:
>What this means is that over RAND_MAX + (n-k) calls to
>this function, the distribution is even. Of course it
>also means that at one point in the random number sequence
>you get a little block of sequential numbers. However,
>that problem is fairly trivial to fix. You can simply
>spread the k..n-1 sequence throughout the rest of the
>calls evenly so that they don't all come in one block.
>(I'll leave it up to the reader to write the code, my
>head hurts after all the math and stuff already.)
If your application happens to be periodic with the period
you chose for "spreading the sequence throughout the rest
of the code evenly" then the non-randomness will bite you
bigtime.
For example, in my simulation programs, I may make one major
branch choice for every 9,999 minor branch choices. If
your spreading-out happens to put the high bin numbers every
10,000 numbers, my major branch choices will all be very
high!
This may not be likely, but one wants better out of an
industrial grade RNG.
You can't scatter the high bin numbers aperiodically without needing
a random function to scatter them, which is a chicken and egg
situation that will slow down your code.
Mary Kuhner mkku...@gs.washington.edu