help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: bad random numbers


From: John Smith
Subject: Re: bad random numbers
Date: Wed, 14 Jul 1999 21:43:13 +0100 (BST)

Much has been said on these random numbers:

Extracting the various bits of code, I think the random number
initialisation and first generation of a random number looks 
roughly like:

( libcruft/ranlib/inrgcm.f)
      m1 = 2147483563  ( = 2e9 ish )
      m2 = 2147483399

( src/rand.cc )
      time (&now);
      tm = localtime (&now);
 
      int hour = tm->tm_hour + 1;
      int minute = tm->tm_min + 1;
      int second = tm->tm_sec + 1;

      int s1 = tm->tm_mday * hour * minute * second;
      int s2 = hour * minute * second;

( ignlgi.f )
      k1 = s1/53668
      s1 = 40014* ( MOD(s1,53668) ) - k1*12211
      IF (s1.LT.0) s1 = s1 + m1

      k2 = s2/52774
      s2 = 40692* ( MOD(s2,52774) ) - k2*3791
      IF (s2.LT.0) s2 = s2 + m2

      return ( (s1-s2) modulo m1 ) / m1

I think the comparatively surprising initial random numbers are caused
largely by the fact that k1 and k2 change rairly; 
k1 changes every few seconds at Jul 14 21:27:50, k2 is always zero.
Anyone working in the early hours after midnight, 
particularly on 1st of the month will probably see no change! 
mod(s1,53668) spends several seconds stepping up,
before wrapping back to zero. s2 is often stuck on on a ramp 
for a whole minute,

We need a seeding system that sets the initial s1 and s2 values so 
that k1, k2, mod(s1,53668) and mod(s2,52774) all seem to be unrelated
from second to second. Every seed should create initial values of
s1 and s2 which have changes by several hundred thousand from the 
previous seed.  Perhaps:

  s1 = (11+tm->tm_mday) * (101+hour) * (103+minute) * (107+second)
  s2 = (41) * (71+hour) * (103+minute) * (157+second)

or as mentioned previously, warming up the generator by doing a few     

    F77_XFCN (dgenunf, DGENUNF, (0.0, 1.0, val));

after the initial

    F77_XFCN (setall, SETALL, (s0, s1));

Just my thoughts.

John Smith






---------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.  To ensure
that development continues, see www.che.wisc.edu/octave/giftform.html
Instructions for unsubscribing: www.che.wisc.edu/octave/archive.html
---------------------------------------------------------------------



reply via email to

[Prev in Thread] Current Thread [Next in Thread]