Portable random numbers

Tue, 28 Oct 1997 13:48:01 +0100 (MET)

>IMHO, this code is still broken.  On a non-DOS platform it tries to
>pass an integer to a function expecting a cardinal and always returns a
>type check error and doesn't compile.  Is there some advanced switch
>that should be used in situations like this to turn off type-checking
>(I've never really looked at all of the directives and didn't see where
>they were in the docs after a quick look) and the type-casting section
>is empty.
>
>The code in question:
>FUNCTION Rand:Integer; C;
>PROCEDURE SRand(Seed:Cardinal); C;
>FUNCTION GetPID:Integer; C;
>

Regarding portability: in my DJGPP C library, there is no `srand' function,
but a `srandom' function which expects an int (=Integer) as parameter.  I am
not sure whether these C functions belong to the ANSI standard for the C
library (I have left my copy of `The C programing language' by KR at home).

>PROCEDURE Randomize;
>BEGIN
>  SRand(GetPID);
>  RandSeed:=Rand
>END;
>

I don't see how this should emulate the Randomize procedure in Borland
Pascal.  According the info docs for my C library, getpid returns the unique
process number for the program running on the system.  If you call getpid
twice in the same program, you get the same integer.  Why not compute a seed
from the system time (somehow), as Borland Pascal does?

That is, if we need a Randomize procedure at all?  My C library (DJGPP)
doesn't have such a beast.  Ok, there is Borland compatibility ....

I use random numbers a great deal in my work (statistical simulations), and
in most cases (if not all) I specify the seed myself, since this is the only
way of gettting reproducible results.  This may not sound like "random
numbers", but there is really no such thing as random numbers on a digital
computer!  What we generate is a completely deterministic sequence (called
`pseudo random') which just looks random;  Whatever constitutes `random' in
this context depends on the particular application!

When I was programming in Turbo Pascal (the pre-GPC era), I *never* used the
Random function in the system unit because I did not have source code and
thus couldn't see how the pseudo-random sequence was generated.  There are a
lot of bad (= non-random) random number generators out there!  Instead, I
wrote my own random number generator, or, rather found one in books/papers on
the topic [there is a huge literature].

The following piece of code represents a widely used procedure for generating
uniform random numbers (reals between 0.0 and 1.0), but it could easily be
changed to generating integers between 0 and 2^31-1 (just modify the last
line in RanUni).

The routine is portable to all systems with a 32-bit integer type.  It does
not rely on controlled overflow which may not be portable.  The technically
inclined can consult the reference given in the comments.  BTW, GPC users
familiar with the book(s) `Numerical Recipes' may be interested in knowing
that the `ran0´ routine in the latest (2nd) edition (C only; not Pascal)
corresponds to the algorithm below.

Anyways; here's the code:

var
  ru_s1 : Integer;   { Global variable }


{ Initialize the seed with value s1 }

procedure Init_RanUni (s1 : Integer);
begin
  ru_s1 := s1;
end;

{ Basic uniform random generator: Minimal Standard in Park and    }
{ Miller (1988): "Random Number Generators: Good Ones Are Hard to }
{ Find", Comm. of the ACM, 31, 1192-1201.                         }
{ Parameters: m = 2^31-1, a=48271.                                }

function RanUni : double;
const
  im = 2147483647;
  ia = 48271;
  iq = 44488;  { m div a }
  ir = 3399;   { m mod a }

var
  k : Integer;

begin
  k := ru_s1 div iq;

  ru_s1 := ia * (ru_s1 - k*iq) - ir * k;
  if ru_s1 < 1 then ru_s1 := ru_s1 + im;

  RanUni := ru_s1 / im;
end;



Regards,

  Jesper Lund

Jesper Lund (jel@hha.dk)