simple PRNG?

[jacob navia]

jacob navia said:


Show us.

This program:
---------------------------------------------------------------
/*
* Pseudo-random number generator. The result is uniform on
* [0, 2^31 - 1].
*/
static unsigned long _randseed = 1;

unsigned long random(void)
{
register long x, hi, lo, t;

/*
* Compute x[n + 1] = (7^5 * x[n]) mod (2^31 - 1).
* From "Random number generators: good ones are hard to find",
* Park and Miller, Communications of the ACM, vol. 31, no. 10,
* October 1988, p. 1195.
*/
x = _randseed;
hi = x / 127773;
lo = x % 127773;
t = 16807 * lo - 2836 * hi;
if (t <= 0)
t += 0x7fffffff;
_randseed = t;
return (t);
}

void srandom(unsigned long seed)
{
_randseed=seed;
}


#include <stdio.h>
int main(void)
{
for (int i=0; i<10;i++) {
printf("%d\n",random()&0xffff);
}
}

-----------------------------------------------------------------------
produces the following output:
16807
15089
44249
3114
46978
56008
36568
2558
12099
1101

(1)
I am of course unable to write a fully compliant C99 compiler as you
say, but I can still read a list of numbers.

Your assertion that the random number generator produces always the
same result is FALSE.


(2) I can't tell you if the properties of those numbers are the same as
the 32 bit numbers. I *know* that this algorithm has been ported to a
16 bit DSP, but I can't tell you its equivalent in standard C.

(3) For a very good discussion about the history of this algorithm,
random numbers, and related topics see

http://www.firstpr.com.au/dsp/rand31/

 

[Richard Tobin]

This needs 32 bit maths as you see. I think it could be ported to 16
bits by making the result mod 2^16 -1

Did you mean mod 2^16, i.e. using the bottom 16 bits.

-- Richard

 

[jacob navia]

Did you mean mod 2^16, i.e. using the bottom 16 bits.

-- Richard

Yes, but I do not know if the results satisfy the required
properties of a good random number generator. The 32 bit
sequence has a long cycle number, but the 16 bit mod of that
sequence could have completely different properties, I do not know.

I know that the algorithm was ported to a 16 bit dsp, and most 16 bit
machines provide a 32 bit result of two 16 bit number multiply...

 

[Walter Roberson]

Walter Roberson wrote:

1) I have no license for this, it is not my algorithm as stated
in the comment. How could I license that?

Your -expression- of the algorithm is copyrighted, and you
routinely license your -copyright- to others, usually requiring
a fee for commercial usage . The original poster requested a PNRG
that was license free. Are you offering your copyrighted code
for this routine as being license free? That is, are you offering
to let the original poster (and everyone else) copy that
routine (which is presently copyrighted by you) without restriction,
including for commercial use?

The original paper by Steve Park and Keith Miller contains code
in Pascal, and code in FORTRAN, but no code in C, so the code
you posted was not -copied- from the paper but rather written
by you (or copied by you from someone else), so the -copyright-
on the code is not resident in the paper you cited.

(Note: in mentioning that you require a fee for commercial usage,
I am in no way criticising you for chosing to charge licensing fees;
I've made my living for many years from organizations that charge
software licensing fees, and I do not subscribe to the notion
that "Software demands to be free!". The point of my mentioning your
commercial fees is solely with respect to the original poster's
requirement about *not* having to worry about licensing (and thus
potentially fees) for use of the desired PNRG.)

3) You can't patent an algorithm, as far as I know.

In the paper you cited, in the About the Authors section,
Steve Park and Keith Miller are both listed as being at College
of William and Mary, in Williamsburg Virginia USA -- and in
the USA, you certainly *can* patent algorithms. See for example,
http://www.bitlaw.com/software-patent/history.html#1990s

But even before that, we can see successful algorithm patents
such as the LZW patent filed in June 1983,
http://patft.uspto.gov/netacgi/nph-Parser?patentnumber=4558302

 

[Walter Roberson]

A short code fragment implementing a published algorithm is not
copyrightable anyway.

That is completely incorrect from a legal standpoint. Copyright
law is about the "original expression of an idea", not about the
idea itself. For example in books and movies, basic plots get
reused over and over again, but each of them is copyrightable
because each of them is a different -expression- of the idea.

Do you have any reason to suppose that the algorithm is patented?

The authors of the publication were in the USA at the time
of publication, so the possibility of a patent on the algorithm
is real. Any offering of specific code that does not *know*
whether the underlying algorithm is patented or not (or at
least specifically consider the point) risks being interpreted
as a definitive statement about the legal status of the algorithm.

Do
you even think it's possible to patent such an algorithm in most of
the world?

Oh yes, definitely; all you have to do is wrap it within the phrasing
of being part of a "computer system" (or equivilent) and that
often gives enough physicallity to qualify for a patent.

In short, your post is mean-minded FUD, a most unworthy response
to someone who is trying to help.

The matter of whether the PNRG was license free or not was
a critical part of the original poster's request. Any referral
to specific code that fails to consider whether the code is
license free or not fails to address this critical part of the
original request.

 

[Walter Roberson]

This program:
---------------------------------------------------------------
/*
* Pseudo-random number generator. The result is uniform on
* [0, 2^31 - 1].
*/

#include <stdio.h>
int main(void)
{
for (int i=0; i<10;i++) {
printf("%d\n",random()&0xffff);
}
}

-----------------------------------------------------------------------
produces the following output:
16807
15089
44249
3114
46978
56008

On a system where int was 16 bit, random()&0xffff
risks exceeding INT_MAX, so printing the result with a %d format
could result in negative numbers, such as would occur for
the results 44249 46978 56008. An unsigned output format
should be used instead.

 

[jacob navia]

Walter Roberson wrote:
[snip]

On a system where int was 16 bit, random()&0xffff
risks exceeding INT_MAX, so printing the result with a %d format
could result in negative numbers, such as would occur for
the results 44249 46978 56008. An unsigned output format
should be used instead.

Granted, but this doesn't mean that the algorithm will
always print the same number as Heathfield supposed.

I should have written "%u" of course.

 

[jacob navia]

Walter Roberson wrote:
[snip legalese]

OK. Mr Roberson you know FAR more legal stuff than I do.

The only thing that I can answer is:

DISCLAIMER:
-----------

The code published by me in this group is released into
the public domain. Anyone can use it for any purpose
whatsoever.

No liabilities can be construed against me. Use at your
own risk.

 

[copx]

Can anyone point me to a simple, fast RRNG function to generate random
ints within a specified range?

To answer my own question:

I cooked this up based on C FAQ entries. I use the portable "minimal
standard" algorithm, and the recommended method to get numbers within a
certain range based (the float free variant):

------- code start ----
#include <stdio.h>
#include <time.h>

static long int seed = 1;


/*
* generate int between lb and ub (inclusive)
*/
int mrand(int lb, int ub)
{
long int hi = seed / 44488;
long int lo = seed % 44488;

seed = 48271 * lo - 3399 * hi;

if (seed <= 0) seed += 2147483647;

return lb + (seed / (2147483646 / (ub + 1 - lb) + 1));
}


/*
* test distribution
*/
void test_it(void)
{
int lim = 11000;
int oc[11] = {0,0,0,0,0,0,0,0,0,0,0};
int i;

while (lim-- > 0) ++oc[mrand(0, 10)];

for (i = 0; i < 11; i++) printf("%d: %d\n", i, oc);
}


int main(void)
{

seed = time(NULL);

test_it();

return 0;
}

--------- code ends here ------

Feel free to point out any problems.

It seems to work, at least on my system. Notice that I have also tested the
algorithm with billions of random numbers, the distribution stays uniform
(or at least close enough to be called "non-biased" - I do not know how to
test for real uniformity - but this is only for a game so it's definitely
good enough).

 

[Walter Roberson]

unsigned long random(void)
{
register long x, hi, lo, t;

/*
* Compute x[n + 1] = (7^5 * x[n]) mod (2^31 - 1).
* From "Random number generators: good ones are hard to find",
* Park and Miller, Communications of the ACM, vol. 31, no. 10,
* October 1988, p. 1195.
*/
x = _randseed;
hi = x / 127773;
lo = x % 127773;
t = 16807 * lo - 2836 * hi;
if (t <= 0)
t += 0x7fffffff;
_randseed = t;
return (t);
}

It would be safer (or at least clearer) if most of your constants
were suffixed with L to avoid their being interpreted (if only
mentally) as overflowing the range of UINT_MAX on 16 bit systems.

 

[Eric Sosman]

Can anyone point me to a simple, fast RRNG function to generate random ints
within a specified range? It is important that each value within the range
has the same probability (uniform distribution).
I do not want to use the unreliable rand() function, but I do not want to
bloat my code with something as complex as MT either. I am just looking for
a short code snippet that I can copy without worrying about licensing.
The function should work on limited platforms (no floating-point math
please, one that works even on platforms where int is only 16 bit would be
perfect).
I did search this group and the web but I could not find anything which
meets the requirements.

MT isn't all that large, but if you seek things even
shorter you might want to search for Park & Miller's "Minimal
Standard Random Number Generator." Also, George Marsaglia
posted a few of his generators to this group a few years ago,
before the ignorant yahoos drove him away. You should be
able to find them in the archives.

 

[Richard Heathfield]

jacob navia said:

This program:

printf("%d\n",random()&0xffff);

That's not mod, which is what you suggested earlier. It's bitwise-AND. mod
2^16-1 gives the result I told you about before. But let's take your
modification and run with it. Doing so produces 32767 on my system, every
time.

produces the following output:
16807
15089
44249
3114
46978
56008
36568
2558
12099
1101

Interesting. I get different results.

(1)
I am of course unable to write a fully compliant C99 compiler as you
say,

I didn't say you can't. I said you haven't. Whether you can or not remains
to be seen.

but I can still read a list of numbers.

So can I. I ran your function 33554432 times, so I'd expect approximately
512 results (give or take) for each possible value in the range 0 to
65535. What I actually get is 33554432 lots of 32767, and *nothing else*.

Now, the possibility remains that my driver has a bug. I've looked at it
very carefully, and I can't find such a bug. (Of course, that might be
because I didn't use a debugger.)

Here, then, is my driver. Feel free to point out the bug.

#include <stdio.h>

#define FREQ_EXPECTED 512

#define ARRAY_SIZE 65536UL

int main(void)
{
unsigned long frequency[ARRAY_SIZE] = {0};
unsigned long counter = 0;

srandom(0);

for(counter = 0; counter < ARRAY_SIZE * FREQ_EXPECTED; counter++)
{
++frequency[random() & 0xFFFF];
}

for(counter = 0; counter < ARRAY_SIZE; counter++)
{
if(frequency[counter] != 0)
{
printf("%lu,%lu\n", counter, frequency[counter]);
}
}
return 0;
}

Output:

65535,33554432

Your assertion that the random number generator produces always the
same result is FALSE.

Not from where I'm sitting - not for your first suggestion (modulo), or
indeed for your second (bitwise-AND).

<snip>

 

[Eric Sosman]

/*
* Pseudo-random number generator. The result is uniform on
* [0, 2^31 - 1].

ITYM (0, 2^31 - 1) or [1, 2^31 - 2].

 

[Walter Roberson]

DISCLAIMER:
-----------

The code published by me in this group is released into
the public domain. Anyone can use it for any purpose
whatsoever.

No liabilities can be construed against me. Use at your
own risk.

Fair enough; I asked you to clarify the license terms of the
code and you did.


For clarity: did you intend your release into the public domain to
apply only to this one routine, or to apply to -all- of your code that
you post in comp.lang.c? If you intend it to apply to -all- of your
code that you post, then it would be better from a legal standpoint to
have the release into the public domain accompany each block of code
that you so release. Under the Berne Convention on Copyright, original
expressions of ideas are automatically copyrighted by their author
(except perhaps "works for hire") unless copyright is specifically
released, so for any of your code that you published, the legal
presumption would be that it is copyrighted unless the copyright
release was "right there".

 

[Walter Roberson]

/*
* generate int between lb and ub (inclusive)
*/
int mrand(int lb, int ub)
{
long int hi = seed / 44488;
long int lo = seed % 44488;

seed = 48271 * lo - 3399 * hi;

if (seed <= 0) seed += 2147483647;

return lb + (seed / (2147483646 / (ub + 1 - lb) + 1));
}

Caution: you are mixing int arithmetic with long arithmetic.
ub + 1 could overflow INT_MAX, and ub + 1 - lb could
certainly overflow INT_MAX.

Also, as I indicated in response to Jacob's code: it would
be safer (or at least clearer) to suffix your constants with
L when you are working with long arithmetic, to avoid
problems (or at least mental confusion) with constants that
exceed INT_MAX.

 

[Ben Bacarisse]

Can anyone point me to a simple, fast RRNG function to generate random ints
within a specified range? It is important that each value within the range
has the same probability (uniform distribution).
I do not want to use the unreliable rand() function, but I do not want to
bloat my code with something as complex as MT either. I am just looking for
a short code snippet that I can copy without worrying about licensing.
The function should work on limited platforms (no floating-point math
please, one that works even on platforms where int is only 16 bit would be
perfect).

Probably the best bit-size independent algorithm is

Robert M. Ziff, "Four-tap shift-register-sequence random-number
generators," Computers in Physics 12(4), Jul/Aug 1998, pp 385-392.

It is very simple if you can spare the space for the state (min is
about 9700 previous numbers). It can be implemented slightly faster
if you can spare the space for 16384 of them (the next power of 2).

The algorithm is just

R[n] = R[n-A] ^ R[n-B] ^ R[n-C] ^ R[n-D]

with "good" choices of A, B, C and D being 471, 1586, 6988 and 9689.
It is fast and has a huge period. The advantage over arithmetic PRNGs
is that it is in essence parallel: each number is result is combining
N bits from N random bit streams (N being the width of each R) so it
works just as well with N=1 (requiring about 9700 bits of state) or
N=64 with each R being a 64-bit int.

Implementation probably no more than 20 lines. GPL version exists in
the GSL library.

 

[jacob navia]

Richard Heathfield wrote:

[snip drivel]

Here, then, is my driver. Feel free to point out the bug.

[snip]


srandom(0);

There

I never set the seed to zero in my code.

Seed must be bigger than zero obviously
And please do not say that you "did not know that"
since I did not use srandom() at all in the example
I sent.

The whole content of Heathfield's messages is this:

Word games, more word games, bad faith, apparently innocent
remarks full of hate... etc.

 

[Eric Sosman]

[...]
unsigned simple_fast_uniform_free_16(void)
{ static unsigned x;
return x = 0x0FFFFu & ++x;
}

is restricted to 16-bit integer values, meets speed requirement
for 8086, and distribution is perfectly uniform. It's free and
there is no licensing to worry about.

Undefined behavior is, I guess, a kind of "randomness."

 

[Morris Dovey]

[...]
unsigned simple_fast_uniform_free_16(void)
{ static unsigned x;
return x = 0x0FFFFu & ++x;
}

is restricted to 16-bit integer values, meets speed requirement
for 8086, and distribution is perfectly uniform. It's free and
there is no licensing to worry about.

Undefined behavior is, I guess, a kind of "randomness."

I decided to give socialized programming a shot. <g>

 

[jacob navia]

[...]
unsigned simple_fast_uniform_free_16(void)
{ static unsigned x;
return x = 0x0FFFFu & ++x;
}

is restricted to 16-bit integer values, meets speed requirement
for 8086, and distribution is perfectly uniform. It's free and
there is no licensing to worry about.

Undefined behavior is, I guess, a kind of "randomness."

(!!!!!!!!!!!!!!!)

Why this nasty habit of laughing at people that
ask questions here?

And then presenting a program that has UB!

I think he was trying to mimic Heathfield with his other
example mocking the OP.