KISS4691, a potentially top-ranked RNG.
in [Fortran]
|
Prev: solutions book
Next: real kind declaration
From: glen herrmannsfeldt on 26 Jul 2010 16:58 In comp.lang.fortran nmm1(a)cam.ac.uk wrote: > In article <fa9dd141-823a-4179-a80b-18d214a3e846(a)t2g2000yqe.googlegroups.com>, > Harald Anlauf <anlauf.2008(a)arcor.de> wrote: (snip) >>Does this mean that using different seeds will lead to >>streams that are always statistically independent >>(as long as one does not exhaust the RNG's period)? >>Or are there restrictions on the possible combinations >>of seeds? > No. Without checking it more carefully, I can't say definitely, > but it looks as if you would be very unlikely to notice a PRACTICAL > association with two different seeds, provided that you throw > away the first 10,000 or so numbers - and even that qualification > may be unnecessary. But, unless I have some missed some subtlety, > the sequences cannot be guaranteed to be pseudo-independent. My biggest complaint about the current standard RANDOM_SEED is that it doens't provide a way to get a reliably good seed from a default (likely 32 bit) integer. There are many generators with extrememly long periods, and correspondingly long state. As the designers of the RNG are the ones likely to know how to choose a good seed, it would seem they would be the best ones to supply a good seed generator. > The only two methods I know of of guaranteeing pseudo-independence > are using coprime sequences and by choosing them using the spectral > test or equivalent. Even then, there are some qualifications on > what is meant by pseudo-independence. However, in practice, it's > rare to have trouble with high-quality generators. Now, one can supply an array of the appropriate length to RANDOM_SEED(PUT=...), but how to generate such an array from a smaller seed? There is no way to know. (snip) > You need to be very careful to distinguish separate (i.e. disjoint) > sequences from pseudo-independent ones, and FAR too many papers > written by people who ought to know better confuse the two. Doing > that is a common cause of seriously wrong answers in some types of > calculation. -- glen
From: nmm1 on 26 Jul 2010 17:39 In article <de9ee692-e7b7-4b2b-acff-6f8aaefaeec2(a)i28g2000yqa.googlegroups.com>, Harald Anlauf <anlauf.2008(a)arcor.de> wrote: >> >> The only two methods I know of of guaranteeing pseudo-independence >> are using coprime sequences and by choosing them using the spectral >> test or equivalent. =A0Even then, there are some qualifications on >> what is meant by pseudo-independence. =A0However, in practice, it's >> rare to have trouble with high-quality generators. > >I was shown funny experiences with simple generators... =3D:-o Yup. Very common. >> >I am currently using D.E.Knuth's generator from TAOCP, >> >which IIRC allows for 2^30-2 independent streams, and >> >which are asserted to be independent, but being able to >> >extend the "limit" would be nice. >> >> Eh? =A0Which version? =A0There was assuredly nothing with those >> properties in edition 2. =A0The first papers on the topic were later. > >http://www-cs-faculty.stanford.edu/~uno/news02.html#rng Ah. Thanks. It's new, then. I haven't been active in this area ina couple of decades, so haven't been keeping track. >> You need to be very careful to distinguish separate (i.e. disjoint) >> sequences from pseudo-independent ones, and FAR too many papers >> written by people who ought to know better confuse the two. =A0Doing >> that is a common cause of seriously wrong answers in some types of >> calculation. > >For an application which needs to distribute a large problem over many >processors it is very useful to have many (pseudo-)independent >streams, >at least they should be practically indistinguishable from independent >ones. An accidental correlation is likely worse than a short period. Right. That's precisely why I got involved. Regards, Nick Maclaren.
From: robin on 26 Jul 2010 20:21 "Geoff" <geoff(a)invalid.invalid> wrote in message news:qedr46dbfhrlisiig9pk3osbsmiuna7u4s(a)4ax.com... | On Mon, 26 Jul 2010 23:32:27 +1000, "robin" <robin51(a)dodo.com.au> | wrote: | | >"Gib Bogle" <g.bogle(a)auckland.no.spam.ac.nz> wrote in message news:i2ij74$kd6$1(a)speranza.aioe.org... | >| geo wrote: | >| > Thanks very much for the Fortran version. | >| > I made a mistake in the comment on versions | >| > for signed integers. This: | >| > | >| > Languages requiring signed integers should use equivalents | >| > of the same operations, except that the C statement: | >| > c=(t<x)+(x>>19); | >| > can be replaced by that language's version of | >| > if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19) | >| > else c=1-sign(x<<13+c)+(x>>19) | >| > | >| > Sorry for that error. | >| | >| That produces different c values from those generated by the method based on the | >| value of (t<x), therefore it deviates from the C code. This is what I used: | >| | >| m = shiftl(x,13) + c | >| if (sign(1,m) == sign(1,x)) then | >| c = sign(1,x) + shiftr(x,19) | >| else | >| c = 1 - sign(1,m) + shiftr(x,19) | >| endif | > | >Maybe I missed something, | >but isn't this exactly equivalent to what George wrote? | >Just substitute x<<13+c for m in your last two assignments ... | | I am no FORTRAN hacker but I think there's a difference between | sign(x) and sign(1,x). George gave general advice on how to do it. That advice wasn't specific to Fortran. It's necessary to use the appropriate Fortran construct -- and sign(1,x) is the only way to do that.
From: Geoff on 27 Jul 2010 00:40 On Tue, 27 Jul 2010 10:21:34 +1000, "robin" <robin51(a)dodo.com.au> wrote: >"Geoff" <geoff(a)invalid.invalid> wrote in message news:qedr46dbfhrlisiig9pk3osbsmiuna7u4s(a)4ax.com... >| On Mon, 26 Jul 2010 23:32:27 +1000, "robin" <robin51(a)dodo.com.au> >| wrote: >| >| >"Gib Bogle" <g.bogle(a)auckland.no.spam.ac.nz> wrote in message news:i2ij74$kd6$1(a)speranza.aioe.org... >| >| geo wrote: >| >| > Thanks very much for the Fortran version. >| >| > I made a mistake in the comment on versions >| >| > for signed integers. This: >| >| > >| >| > Languages requiring signed integers should use equivalents >| >| > of the same operations, except that the C statement: >| >| > c=(t<x)+(x>>19); >| >| > can be replaced by that language's version of >| >| > if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19) >| >| > else c=1-sign(x<<13+c)+(x>>19) >| >| > >| >| > Sorry for that error. >| >| >| >| That produces different c values from those generated by the method based on the >| >| value of (t<x), therefore it deviates from the C code. This is what I used: >| >| >| >| m = shiftl(x,13) + c >| >| if (sign(1,m) == sign(1,x)) then >| >| c = sign(1,x) + shiftr(x,19) >| >| else >| >| c = 1 - sign(1,m) + shiftr(x,19) >| >| endif >| > >| >Maybe I missed something, >| >but isn't this exactly equivalent to what George wrote? >| >Just substitute x<<13+c for m in your last two assignments ... >| >| I am no FORTRAN hacker but I think there's a difference between >| sign(x) and sign(1,x). > >George gave general advice on how to do it. >That advice wasn't specific to Fortran. >It's necessary to use the appropriate Fortran construct -- >and sign(1,x) is the only way to do that. > Thanks for that clarification. I have not done any Fortran code since 1975 and it looked a whole lot different than it does today.
From: nmm1 on 27 Jul 2010 03:22
In article <knos461nt5r9srat9tnjcvubc27lmpd55o(a)4ax.com>, Geoff <geoff(a)invalid.invalid> wrote: >>| >>| I am no FORTRAN hacker but I think there's a difference between >>| sign(x) and sign(1,x). >> >>George gave general advice on how to do it. >>That advice wasn't specific to Fortran. >>It's necessary to use the appropriate Fortran construct -- >>and sign(1,x) is the only way to do that. >> >Thanks for that clarification. I have not done any Fortran code since >1975 and it looked a whole lot different than it does today. Actually, that aspect hasn't changed since :-) SIGN(X) always was an implementation-dependent extension (now called processor-dependent). Regards, Nick Maclaren. |