KISS4691, a potentially top-ranked RNG.
in [Fortran]
|
Prev: solutions book
Next: real kind declaration
From: robin on 27 Jul 2010 01:46 "geo" <gmarsaglia(a)gmail.com> wrote in message news:a82cebe3-cdb9-48af-8080-bca935eeb9b1(a)l14g2000yql.googlegroups.com... |I have been asked to recommend an RNG | (Random Number Generator) that ranks | at or near the top in all of the categories: | performance on tests of randomness, | length of period, simplicity and speed. | The most important measure, of course, is | performance on extensive tests of randomness, and for | those that perform well, selection may well depend | on those other measures. I have already posted a PL/I version using unsigned arithmetic. Here is another version, this time using signed arithmetic :-- (NOSIZE, NOFOFL): RNG: PROCEDURE OPTIONS (MAIN, REORDER); declare (xs initial (521288629), xcng initial (362436069), Q(0:4690) ) static fixed binary (31); MWC: procedure () returns (fixed binary (31)); declare (t,x,i) fixed binary (31); declare (c initial (0), j initial (4691) ) fixed binary (31) static; declare (t1, t2, t3) fixed binary (31); if j < hbound(Q,1) then j = j + 1; else j = 0; x = Q(j); t = isll(x,13)+c+x; t1 = iand(x, 3) - iand(t, 3); t2 = isrl(x, 2) - isrl(t, 2); if t2 = 0 then t2 = t1; if t2 > 0 then t3 = 1; else t3 = 0; c = t3 + isrl(x, 19); Q(j)=t; return (t); end MWC; CNG: procedure returns (fixed binary (31)); xcng=bin(69069)*xcng+bin(123); return (xcng); end CNG; XXS: procedure returns (fixed binary (31)); xs = ieor (xs, isll(xs, 13) ); xs = ieor (xs, isrl(xs, 17) ); xs = ieor (xs, isll(xs, 5) ); return (xs); end XXS; KISS: procedure returns (fixed binary (31)); return ( MWC()+CNG+XXS ); end KISS; declare (i,x) fixed binary (31); declare y fixed decimal (11); Q = CNG+XXS; /* Initialize. */ do i = 1 to 1000000000; x=MWC(); end; put skip edit (" Expected MWC result = 3740121002", 'computed =', x) (a, skip, x(12), a, f(11)); y = iand(x, 2147483647); if x < 0 then y = y + 2147483648; put skip edit (y) (x(11), f(22)); put skip; do i = 1 to 1000000000; x=KISS; end; put skip edit ("Expected KISS result = 2224631993", 'computed =', x) (a, skip, x(12), a, f(11)); y = iand(x, 2147483647); if x < 0 then y = y + 2147483648; put skip edit (y) (x(11), f(22)); end RNG;
From: robin on 27 Jul 2010 06:19 "jacob navia" <jacob(a)spamsink.net> wrote in message news:i2fir2$op4$1(a)speranza.aioe.org... | This doesn't work with systems that have unsigned long as a 64 bit quantity. | | I obtain: | | MWC result=3740121002 ? | 4169348530 | KISS result=2224631993 ? | 1421918629 For a 64-bit machine (using 64-bit integer arithmetic), you'd need to truncate each result to 32 bits. That not only applies to the multiplication, it also applies to addition, etc. On a 32-bit machine, these extra bits are discarded, but in 64-bit arithmetic, they are retained, and unless they are similarly discarded, you won't get the same results. I suggest using IAND(k, 2*2147483647+1) for the truncation. With such modifications in the program, it should then produce the same results on both 32-bit and 64-bit machines. P.S. the product 2*2... is best obtained using ISHFT. | Compiling with 32 bit machine yields: | MWC result=3740121002 ? | 3740121002 | KISS result=2224631993 ? | 2224631993
From: Uno on 29 Jul 2010 15:28 glen herrmannsfeldt wrote: > 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. So for the put= values in fortran, you need a vector of pseudorandom integers, which is as good as it gets without truly random devices, making--one hopes-a period that is large with respect to the interval you're interested in. It doesn't seem like a problem with epistemology as much a mathematical ceiling on how much randomness you can create by a handful of values. -- Uno
From: glen herrmannsfeldt on 29 Jul 2010 17:47 In comp.lang.fortran Uno <merrilljensen(a)q.com> wrote: (snip, I wrote) >> 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. (snip) > So for the put= values in fortran, you need a vector of pseudorandom > integers, which is as good as it gets without truly random devices, > making--one hopes-a period that is large with respect to the interval > you're interested in. In a large fraction of the cases, 2 billion different seeds are enough, but one can still desire the appropriate randomness from those different seeds. > It doesn't seem like a problem with epistemology as much > a mathematical ceiling on how much randomness you can create > by a handful of values. Given a default integer, one might fill an array with that integer and use that as a seed. That might be good for some, not so good for others. Even more, for standard Fortran such should be done without knowing the range of values of an integer variable. R has two seeding functions, one that takes a full length state array, and the other takes a single integer. That makes sense to me. -- glen
From: Uno on 29 Jul 2010 20:38
Gib Bogle wrote: > 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(t)+(x>>19) >> >> should have been: >> >> >> 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 Gib, can you post your entire fortran version? -- Uno |