On Mon, Feb 10, 2014 at 05:23:39PM +0200, Janne Blomqvist wrote: > On Sun, Feb 9, 2014 at 2:40 AM, Steve Kargl > <s...@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
