[prev in list] [next in list] [prev in thread] [next in thread] 

List:       gcc-fortran
Subject:    Re: [PATCH, fortran/52879] RANDOM_SEED revisited
From:       Steve Kargl <sgk () troutmask ! apl ! washington ! edu>
Date:       2014-02-10 16:26:57
Message-ID: 20140210162657.GA97502 () troutmask ! apl ! washington ! edu
[Download RAW message or body]

On Mon, Feb 10, 2014 at 05:23:39PM +0200, Janne Blomqvist wrote:
> On Sun, Feb 9, 2014 at 2:40 AM, Steve Kargl
> <sgk@troutmask.apl.washington.edu> wrote:
> > All,
> >
> > Here is a potentially contentious patch for PR fortran/52879.
> >
> > http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52879
> >
> > As demonstrated in that PR, it is possible to provide RANDOM_SEED
> > with sets of seeds that result in a very poor sequences of PRN.
> > Gfortran's RANDOM_NUMBER uses 3 independent KISS PRNG to generate
> > the bits of the significands of the real PRN.  Each KISS PRNG
> > requires 4 seeds.  This patch removes the need for a user to
> > specify 12 seeds.  It uses a Lehmer linear congruent generator
> > to generate the 12 KISS seeds.  This LCG requires a single seed.
> > I have selected to have seed(1)=0 correspond to the current
> > default set of KISS seeds.  Internally, the LCG declares the
> > seed as a GFC_UINTEGER_8 (ie., uint64_t) type, so for gfortran's
> > default integer kind one has 2**31-1 possible seeds and if one
> > use -fdefault-integer-8 then one has 2**32 possible seeds.
> 
> I like this approach, as it would make the seeding a bit more
> fool-proof and easier to use. Per se, loss of the ability to retrieve
> or update parts of the internal state of the RNG seems not terribly
> relevant, after all users who might be interested in that are most
> likely not going to be satisfied with the black box RNG provided by
> the language anyway.
> 
> That being said, I wonder whether this conforms to the standard. From N1830:
> 
> "The pseudorandom number generator used by RANDOM NUMBER maintains a
> seed that is updated during the execution of RANDOM NUMBER and that
> may be specified or returned by RANDOM SEED."
> 
> and
> 
> "
> For example, after execution of the statements
> 
>      CALL RANDOM_SEED (PUT=SEED1)
>      CALL RANDOM_SEED (GET=SEED2)
>      CALL RANDOM_NUMBER (X1)
>      CALL RANDOM_SEED (PUT=SEED2)
>      CALL RANDOM_NUMBER (X2)
> 
> X2 equals X1.
> "
> 
> That would imply that GET= and PUT= save/restore the complete internal
> state of the PRNG, which seems incompatible with using the simple LCG
> to generate the seed (if I understood this correctly, I agree with
> Nick that the interface is badly designed, but what can you do..).

I suppose at one time I knew about the above requirement.  The
simple approach I proposed won't work as is.  A possible work-around
(and it is ugly) would be to keep an internal count of the number
of times that a random number is drawn.  Call this count CNT.  Then,
on reseeding, random_seed would burn CNT draws and reset CNT to zero.
In thinking about it, we could make seed(1)=<LCG seed> and seed(2)=CNT.
This would restrict CNT to 2**31-1 before rolling over.  If CNT is
internal to random_seed(), it can be a uint64_t giving 2**64 before
rolling over.  This certainly has some performance impact on any code
that calls random_seed numerous times.

> Or maybe it would be possible to get around this requirement by
> requiring the 12 seeds, but then calculate some kind of entropy value
> of them, and if the entropy is too low (say, all values are the same,
> or only the first value is non-zero) then use the Lehmer LCG,
> otherwise set the state directly? Sounds brittle, though..

Interesting idea.  I can't remember.  Is it possible to issue
a runtime warning from within libgfortran?  Is so, we could do

   seed2 = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
   call random_seed(put=seed2)

Runtime Warning: poor quality seeds detected in random_seed().
Seeds reset to PUT=[<LCG generated seeds>].

I'll think about a possible algorithm.

-- 
Steve
[prev in list] [next in list] [prev in thread] [next in thread]