random number generator random() questionable!!!?

[Shane Youl]

I came across the same problem some time ago, when I needed a good random
number generator for some Monte Carlo simulations. The following is the
response to my plea for assistance then. I subsequently received more
information from a guy at CERN and can dig this up for you if you need it.
This latter stuff is in Fortran.
-------

From:   IN%"IVANOVIC%VAXR at circus.llnl.gov"  "Vladimir Ivanovic 415.423.7786"  9-
MAR-1989 10:03
To:     info-vax at kl.sri.COM
Subj:   Random Number Generators: the minimal standard

    Below are the 4 implementations of the minimal standard of Park
    and Miller (see "Random Number Generators: Good Ones are Hard to
    Find", Communications of the ACM, October 1988 v31 #10.)  I quote
    from this highly recommended article:

        There is really no argument in support of ANY of the example
        generators cited in this section [The section's title: A Sampling
        of Inadequate Generators].  Most of them are naive implementations
        a deceptively simple idea.  Even the best of them, those which
        have a full period generator and are correctly implemented,
        are inferior to the minimal standard.
        ...
        If you have need for a random number generator, particularly
        one that will port to a variety of systems, and if you are not
        a specialist in random number generation and do not want to
        become one, use the minimal standard.  It should not be presumed
        that it is easier to write a better one.
        ...
        We also recommend the articles by L'Ecuyer and by Wichmann and
        Hill.  The COMBINED generators proposed in these articles represent
        a logical extension of the minimal standard.  These generators
        appear to yield better statistical properties in those (rare)
        situations when the minimal standard alone may be inadequate.

    So there you have it.  A previous posting giving 16-bit and 32-bit
    versions of the combined generators and here several portable versions
    of the minimal standard.

    Since the code is SO simple, there NO excuse for not using at least
    the minimal standard.  Happy hacking.

................................................................................

program MinimalStandard;

    var
        ISeed : integer;
        RSeed : real;

    function Random1 : real;
    { Integer Version 1 :                                                   }
    { NOT a correct version unless maxint is 2**46 - 1 or larger.           }
        const
            a = 16807;
            m = 2147483647;
    begin
        ISeed := (a * ISeed) mod m;
        Random1 := ISeed / m;
    end;

    function Random2 : real;
    { Integer Version 2:                                                    }
    { Correct on any system for which maxint is 2**31 - 1 or larger.        }
        const
            a = 16807;
            m = 2147483647;
            q = 127773;     { m div a }
            r = 2836;       { m mod a }
        var
            lo, hi, test : integer;
    begin
        hi := ISeed div q;
        lo := ISeed mod q;
        test := a * lo - r * hi;
        if test > 0
            then ISeed := test
            else ISeed := test + m;
        Random2 := ISeed / m;
    end;

    function Random3 : real;
    { Real Version 1:                                                       }
    { Correct if reals are represented with a 46-bit or larger mantissa     }
    { (excluding the sign bit).                                             }
        const
            a = 16807.0;
            m = 2147483647.0;
        var
            temp : real;
    begin
        temp := a * RSeed;
        RSeed := temp - m * Trunc(temp / m);
        Random3 := RSeed / m;
    end;

    function Random4 : real;
    { Real Version 2:                                                       }
    { Correct on any system for which reals are represented with a 32-bit   }
    { or larger mantissa (including the sign bit).                          }
        const
            a = 16807.0;
            m = 2147483647.0;
            q = 127773.0;   { m div a }
            r = 2836.0;     { m mod a }
        var
            lo, hi, test : real;
    begin
        hi := Trunc(RSeed / q);
        lo := RSeed - q * hi;
        test := a * lo - r * hi;
        if test > 0.0
            then RSeed := test
            else RSeed := test + m;
        Random4 := RSeed / m;
    end;
begin
    { Your code here... }
end.


-- 
                                         ____   _____     ____    ____
  Shane Youl                            /    \ /       / /    \  /    \   
  CSIRO Division of Mineral Products   /      /_____  / /_____/ /     /
  PO Box 124   Port Melbourne  3207   /            / / /   \   /     /
  AUSTRALIA                           \____/ _____/ / /     \  \____/
  Internet : sfy at dmp.CSIRO.AU     
  Phone    : +61-3-647-0211            SCIENCE  ADVANCING  AUSTRALIA