Re: [PATCH, fortran/52879] RANDOM_SEED revisited

[Janne Blomqvist <blomqvist dot janne at gmail dot com>]

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..).

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..

> It is possible to extend the patch to allow a skip count such that
> for example seed(2)=10 would use every tenth value in the LCG
> sequence. I haven't pursued this as it appears to be an unneeded
> complication.

Yes, this seems unnecessary.


-- 
Janne Blomqvist