Olaf Seibert <rhialto@polderland.nl> wrote:
> random.S contains the following comment:
> 
>  * Here is a very good random number generator.  This implementation is
>  * based on ``Two Fast Implementations of the "Minimal Standard" Random
>  * Number Generator'', David G. Carta, Communications of the ACM, Jan 1990,
>  * Vol 33 No 1.  Do NOT modify this code unless you have a very thorough
>  * understanding of the algorithm.  It's trickier than you think.  If
>  * you do change it, make sure that its 10,000'th invocation returns
>  * 1043618065.
>  *
>  * Here is easier-to-decipher pseudocode:
>  *
>  * p = (16807*seed)<30:0>       # e.g., the low 31 bits of the product
>  * q = (16807*seed)<62:31>      # e.g., the high 31 bits starting at bit 32
>  * if (p + q < 2^31)
>  *      seed = p + q
>  * else
>  *      seed = ((p + q) & (2^31 - 1)) + 1
>  * return (seed);
>  *
>  * The result is in (0,2^31), e.g., it's always positive.
>
> Given the dire warning above, including "It's trickier than you
> think." I think this is pretty serious. How is one to implement this, or
> verify an implementation is correct?

I can't vouch for the pseudo code, but the algorithm described should
match src/sys/lib/libkern/random.c.

To verify an implementation, start with a seed of 1, call the routine
10,000 times and the last returned value should be 1043618065.

> I wanted to verify this with the machine-independent C version, which I
> expected to be in src/sys/lib/libkern/random.c, but this seems to be a
> verty different algorithm. I don't know if that is supposed to be a
> problem.

It should be the same.  Although it references Communications of the
ACM, vol. 31, no. 10, the same algorithm was used in Communications of
the ACM, Vol 33 No 1.

Cheers,
Dean.


PS.  Here is a sample test program:

int seed=1;

/* acmrand() - returns a positive random integer
 *
 * This implementation is based on ``Two Fast Implementations of
 * the "Minimal Standard" Random Number Generator", David G. Carta,
 * Communications of the ACM, Jan 1990, Vol 33 No 1.
 *
 * 10,000 invocation should return 1043618065. (with seed of 1)
 */

int
acmrand()
{

        seed = 16807 * (seed % 127773) - 2836 * (seed / 127773);
        if (seed <= 0)
                seed += 0x7fffffff;
        return seed;
}

main()
{
        int i, r;
        for(i=1; i<=10000; i++) {
                r = acmrand();
        }
        printf("%d\n", r);
}