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.
The following KISS version, perhaps call it KISS4691, seems to rank at the top in all of those categories. It is my latest, and perhaps my last, as, at age 86, I am slowing down.
Compiling and running the following commented C listing should produce, within about four seconds, 10^9 calls to the principal component MWC(), then 10^9 calls to the KISS combination in another ~7 seconds.
Try it; you may like it.
George Marsaglia
/* The KISS4691 RNG, a Keep-It-Simple-Stupid combination of MWC (Multiply-With-Carry), Congruential and Xorshift sequences. Expressed as a simple MWC() function and C macros #define CNG ( xcng=69069*xcng+123 ) #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) ) #define KISS ( MWC()+CNG+XS ) but easily expressed in other languages, with a slight complication for signed integers.
With its immense period, >10^45000, and speed: MWC()s at around 238 million/sec or full KISSes at around 138 million, there are few RNGs that do as well as this one on tests of randomness and are comparable in even one of the categories: period, speed, simplicity---not to mention comparable in all of them.
The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32 is based on p=8193*b^4691-1, period ~ 2^150124 or 10^45192 For the MWC (Multiply-With-Carry) process, given a current x and c, the new x is (8193*x+c) mod b, the new c is the integer part of (8193*x+c)/b.
The routine MWC() produces, in reverse order, the base-b=2^32 expansion of some j/p with 0<j<p=8193*2^150112-1 with j determined by the 64 bits of seeds xcng and xs, or more generally, by 150112 random bits in the Q[] array. */
static unsigned long xs=521288629,xcng=362436069,Q[4691];
unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/ second*/ {unsigned long t,x,i; static c=0,j=4691; j=(j<4690)? j+1:0; x=Q[j]; t=(x<<13)+c+x; c=(t<x)+(x>>19); return (Q[j]=t);
#include <stdio.h> int main() {unsigned long i,x; for(i=0;i<4691;i++) Q[i]=CNG+XS; for(i=0;i<1000000000;i++) x=MWC(); printf(" MWC result=3740121002 ?\n%22u\n",x); for(i=0;i<1000000000;i++) x=KISS; printf("KISS result=2224631993 ?\n%22u\n",x);
}
/* This RNG uses two seeds to fill the Q[4691] array by means of CNG+XS mod 2^32. Users requiring more than two seeds will need to randomly seed Q[] in main(). By itself, the MWC() routine passes all tests in the Diehard Battery of Tests, but it is probably a good idea to include it in the KISS combination.
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) */
>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.
>The following KISS version, perhaps call it KISS4691, >seems to rank at the top in all of those categories. >It is my latest, and perhaps my last, as, at age 86, >I am slowing down.
>Compiling and running the following commented >C listing should produce, within about four seconds, >10^9 calls to the principal component MWC(), then >10^9 calls to the KISS combination in another ~7 seconds.
>Try it; you may like it.
>George Marsaglia
>/* >The KISS4691 RNG, a Keep-It-Simple-Stupid combination of >MWC (Multiply-With-Carry), Congruential and Xorshift sequences. >Expressed as a simple MWC() function and C macros > #define CNG ( xcng=69069*xcng+123 ) > #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) ) > #define KISS ( MWC()+CNG+XS ) >but easily expressed in other languages, with a slight >complication for signed integers.
>With its immense period, >10^45000, and speed: MWC()s at >around 238 million/sec or full KISSes at around 138 million, >there are few RNGs that do as well as this one >on tests of randomness and are comparable in even one >of the categories: period, speed, simplicity---not to >mention comparable in all of them.
>The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32 >is based on p=8193*b^4691-1, period ~ 2^150124 or 10^45192 >For the MWC (Multiply-With-Carry) process, given a current >x and c, the new x is (8193*x+c) mod b, > the new c is the integer part of (8193*x+c)/b.
>The routine MWC() produces, in reverse order, the base-b=2^32 >expansion of some j/p with 0<j<p=8193*2^150112-1 with j >determined by the 64 bits of seeds xcng and xs, or more >generally, by 150112 random bits in the Q[] array. >*/
>static unsigned long xs=521288629,xcng=362436069,Q[4691];
>unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/ >second*/ >{unsigned long t,x,i; static c=0,j=4691; > j=(j<4690)? j+1:0; > x=Q[j]; > t=(x<<13)+c+x; c=(t<x)+(x>>19); > return (Q[j]=t); >}
>/* >This RNG uses two seeds to fill the Q[4691] array by >means of CNG+XS mod 2^32. Users requiring more than two >seeds will need to randomly seed Q[] in main(). >By itself, the MWC() routine passes all tests in the >Diehard Battery of Tests, but it is probably a good >idea to include it in the KISS combination.
>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) >*/
The compiler also notices that i is defined but unused in MWC().
Purely cosmetic changes for readability and style:
/* * The KISS4691 RNG, a Keep-It-Simple-Stupid combination of * MWC (Multiply-With-Carry), Congruential and Xorshift sequences. * Expressed as a simple MWC() function and C macros * #define CNG ( xcng=69069*xcng+123 ) * #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) ) * #define KISS ( MWC()+CNG+XS ) * but easily expressed in other languages, with a slight * complication for signed integers. * * With its immense period, >10^45000, and speed: MWC()s at * around 238 million/sec or full KISSes at around 138 million, * there are few RNGs that do as well as this one * on tests of randomness and are comparable in even one * of the categories: period, speed, simplicity---not to * mention comparable in all of them. * * The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32 * is based on p=8193*b^4691-1, period ~ 2^150124 or 10^45192 * For the MWC (Multiply-With-Carry) process, given a current * x and c, the new x is (8193*x+c) mod b, * the new c is the integer part of (8193*x+c)/b. * * The routine MWC() produces, in reverse order, the base-b=2^32 * expansion of some j/p with 0<j<p=8193*2^150112-1 with j * determined by the 64 bits of seeds xcng and xs, or more * generally, by 150112 random bits in the Q[] array. * George Marsaglia */
static unsigned long xs = 521288629; static unsigned long xcng = 362436069; static unsigned long Q[4691];
/* takes about 4.2 nanosecs or 238 million/second */ unsigned long MWC(void) { unsigned long t, x, i; static c = 0, j = 4691;
j = (j < 4690) ? j + 1 : 0; x = Q[j]; t = (x<<13) + c + x; c = (t < x) + (x>>19); return (Q[j] = t);
int main() { unsigned long i, x; for(i = 0; i < 4691; i++) Q[i] = CNG + XS; for(i = 0; i < 1000000000; i++) x = MWC(); printf(" MWC result=3740121002 ?\n%22u\n",x); for(i = 0; i < 1000000000; i++) x = KISS; printf("KISS result=2224631993 ?\n%22u\n",x); return 0;
}
/* * This RNG uses two seeds to fill the Q[4691] array by * means of CNG+XS mod 2^32. Users requiring more than two * seeds will need to randomly seed Q[] in main(). * By itself, the MWC() routine passes all tests in the * Diehard Battery of Tests, but it is probably a good * idea to include it in the KISS combination. * * 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) */
> 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.
> The following KISS version, perhaps call it KISS4691, > seems to rank at the top in all of those categories. > It is my latest, and perhaps my last, as, at age 86, > I am slowing down.
> Compiling and running the following commented > C listing should produce, within about four seconds, > 10^9 calls to the principal component MWC(), then > 10^9 calls to the KISS combination in another ~7 seconds.
> Try it; you may like it.
> George Marsaglia
> /* > The KISS4691 RNG, a Keep-It-Simple-Stupid combination of > MWC (Multiply-With-Carry), Congruential and Xorshift sequences. > Expressed as a simple MWC() function and C macros > #define CNG ( xcng=69069*xcng+123 ) > #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) ) > #define KISS ( MWC()+CNG+XS ) > but easily expressed in other languages, with a slight > complication for signed integers.
> With its immense period, >10^45000, and speed: MWC()s at > around 238 million/sec or full KISSes at around 138 million, > there are few RNGs that do as well as this one > on tests of randomness and are comparable in even one > of the categories: period, speed, simplicity---not to > mention comparable in all of them.
> The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32 > is based on p=8193*b^4691-1, period ~ 2^150124 or 10^45192 > For the MWC (Multiply-With-Carry) process, given a current > x and c, the new x is (8193*x+c) mod b, > the new c is the integer part of (8193*x+c)/b.
> The routine MWC() produces, in reverse order, the base-b=2^32 > expansion of some j/p with 0<j<p=8193*2^150112-1 with j > determined by the 64 bits of seeds xcng and xs, or more > generally, by 150112 random bits in the Q[] array. > */
> static unsigned long xs=521288629,xcng=362436069,Q[4691];
> unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/ > second*/ > {unsigned long t,x,i; static c=0,j=4691; > j=(j<4690)? j+1:0; > x=Q[j]; > t=(x<<13)+c+x; c=(t<x)+(x>>19); > return (Q[j]=t); > }
> /* > This RNG uses two seeds to fill the Q[4691] array by > means of CNG+XS mod 2^32. Users requiring more than two > seeds will need to randomly seed Q[] in main(). > By itself, the MWC() routine passes all tests in the > Diehard Battery of Tests, but it is probably a good > idea to include it in the KISS combination.
> 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) > */
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
Compiling with 32 bit machine yields: MWC result=3740121002 ? 3740121002 KISS result=2224631993 ? 2224631993
> 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) > */
Hi George,
I translated this into Fortran, and found that I get different results than with C. I've tracked the difference into MWC(). The following Fortran code, with my laborious comparison of two signed integers treating them as unsigned, gives correct results. If I comment out the line c = tLTx + SHIFTR(x,19) and uncomment the following lines that implement your suggestion above to compute c, I get different results.
integer function MWC() integer :: t, x, i integer, save :: c = 0, j = 4691 integer :: tLTx
if (j < 4690) then j = j + 1 else j = 0 endif x = Q(j) t = SHIFTL(x,13) + c + x if ((t >= 0 .and. x >= 0) .or. (t < 0 .and. x < 0)) then if (t < x) then tLTx = 1 else tLTx = 0 endif elseif (x < 0) then tLTx = 1 elseif (t < 0) then tLTx = 0 endif
c = tLTx + SHIFTR(x,19)
!if (sign(1,SHIFTL(x,13)+c) == sign(1,x)) then ! c = sign(1,x) + SHIFTR(x,19) !else ! c = 1 - sign(1,t) + SHIFTR(x,19) !endif Q(j) = t MWC = t end function
Is it the case that although your suggested workaround gives different results from the C code in this case, it is still equivalent as a RNG?
> > 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) > > */
> Hi George,
> I translated this into Fortran, and found that I get different results than with > C. I've tracked the difference into MWC(). The following Fortran code, with my > laborious comparison of two signed integers treating them as unsigned, gives > correct results. If I comment out the line > c = tLTx + SHIFTR(x,19) > and uncomment the following lines that implement your suggestion above to > compute c, I get different results.
> integer function MWC() > integer :: t, x, i > integer, save :: c = 0, j = 4691 > integer :: tLTx
> if (j < 4690) then > j = j + 1 > else > j = 0 > endif > x = Q(j) > t = SHIFTL(x,13) + c + x > if ((t >= 0 .and. x >= 0) .or. (t < 0 .and. x < 0)) then > if (t < x) then > tLTx = 1 > else > tLTx = 0 > endif > elseif (x < 0) then > tLTx = 1 > elseif (t < 0) then > tLTx = 0 > endif
> c = tLTx + SHIFTR(x,19)
> !if (sign(1,SHIFTL(x,13)+c) == sign(1,x)) then > ! c = sign(1,x) + SHIFTR(x,19) > !else > ! c = 1 - sign(1,t) + SHIFTR(x,19) > !endif > Q(j) = t > MWC = t > end function
> Is it the case that although your suggested workaround gives different results > from the C code in this case, it is still equivalent as a RNG?
> Cheers > Gib
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.
I still like inline functions in Fortan, so would tend to define isign(x)=ishft(x,-31) and m=ishft(x,13)+c if(isign(m).eq.isign(x)) then c=isign(x)+ishft(x,-19) else c=1-isign(m)+ishft(x,-19) and Q[j]=m+x
If calculating the carry c of the MWC operation fails to fix that extra increment properly, then rather than a systematic expansion, in reverse order, 32 bits at a time, of the binary representation of j/(1893*2^150112-1) for some j determined by the random seeds, we will be jumping around in that expansion, and we can't be sure that the period will still be the order of b=2^32 for the prime p=1893*b^4196-1.
> I still like inline functions in Fortan, > so would tend to define > isign(x)=ishft(x,-31)
I don't know anything about random numbers, but "sign" is an intrinsic function in Fortran and compilers should generate near perfect code inline for things like sign(1,x). The generic sign function was added to Fortran 90. In FORTRAN 77, you'd need to use the specific ISIGN(1,X) function.
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
|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. | | The following KISS version, perhaps call it KISS4691, | seems to rank at the top in all of those categories. | It is my latest, and perhaps my last, as, at age 86, | I am slowing down. | | Compiling and running the following commented | C listing should produce, within about four seconds, | 10^9 calls to the principal component MWC(), then | 10^9 calls to the KISS combination in another ~7 seconds. | | Try it; you may like it. | | George Marsaglia
declare (i,x) fixed binary (32) unsigned; /* Initialize: */ do i = 0 to 4691-1; Q(i) = CNG+XXS; end; put skip list (q(0), q(4690)); put skip list ('initialized'); put skip; do i = 0 to 1000000000-1; x=MWC(); end; put skip edit (" MWC result=3740121002 ",x) (a, f(23)); do i = 0 to 1000000000-1; x=KISS; end; put skip edit ("KISS result=2224631993 ",x) (a, f(23));
| > 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 ...
>"Gib Bogle" <g.bo...@auckland.no.spam.ac.nz> wrote in message news:i2ij74$kd6$1@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).
On Jul 24, 3:02 pm, geo <gmarsag...@gmail.com> wrote:
> This RNG uses two seeds to fill the Q[4691] array by > means of CNG+XS mod 2^32. Users requiring more than two > seeds will need to randomly seed Q[] in main(). > By itself, the MWC() routine passes all tests in the > Diehard Battery of Tests, but it is probably a good > idea to include it in the KISS combination.
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?
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.
robin wrote: > "Gib Bogle" <g.bo...@auckland.no.spam.ac.nz> wrote in message news:i2ij74$kd6$1@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 hope so. Maybe I didn't express myself clearly enough. I'll try again. Using my implementation of George's corrected code, I get results from the Fortran code that differ from those generated by his C code.
In article <fa9dd141-823a-4179-a80b-18d214a3e...@t2g2000yqe.googlegroups.com>, Harald Anlauf <anlauf.2...@arcor.de> wrote:
>On Jul 24, 3:02=A0pm, geo <gmarsag...@gmail.com> wrote: >> This RNG uses two seeds to fill the Q[4691] array by >> means of CNG+XS mod 2^32. Users requiring more than two >> seeds will need to randomly seed Q[] in main(). >> By itself, the MWC() routine passes all tests in the >> Diehard Battery of Tests, but it is probably a good >> idea to include it in the KISS combination.
>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.
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.
>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? Which version? There was assuredly nothing with those properties in edition 2. The first papers on the topic were later.
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.
> In article <fa9dd141-823a-4179-a80b-18d214a3e...@t2g2000yqe.googlegroups.com>, > Harald Anlauf <anlauf.2...@arcor.de> wrote:
> >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.
> 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.
I was shown funny experiences with simple generators... =:-o
> >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? Which version? There was assuredly nothing with those > properties in edition 2. The first papers on the topic were later.
> 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.
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.
In comp.lang.fortran n...@cam.ac.uk wrote: > In article <fa9dd141-823a-4179-a80b-18d214a3e...@t2g2000yqe.googlegroups.com>, > Harald Anlauf <anlauf.2...@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.
In article <de9ee692-e7b7-4b2b-acff-6f8aaefae...@i28g2000yqa.googlegroups.com>, Harald Anlauf <anlauf.2...@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.
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.
| >| > 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.
>"Geoff" <ge...@invalid.invalid> wrote in message news:qedr46dbfhrlisiig9pk3osbsmiuna7u4s@4ax.com... >| On Mon, 26 Jul 2010 23:32:27 +1000, "robin" <robi...@dodo.com.au> >| wrote: >| >| >"Gib Bogle" <g.bo...@auckland.no.spam.ac.nz> wrote in message news:i2ij74$kd6$1@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.
In article <knos461nt5r9srat9tnjcvubc27lmpd...@4ax.com>,
Geoff <ge...@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).
|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 :--
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;
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));
| 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
glen herrmannsfeldt wrote: > In comp.lang.fortran n...@cam.ac.uk wrote: >> In article <fa9dd141-823a-4179-a80b-18d214a3e...@t2g2000yqe.googlegroups.com>, >> Harald Anlauf <anlauf.2...@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
In comp.lang.fortran Uno <merrilljen...@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.
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
|
