geo
Posts:
38
Registered:
6/19/08
|
|
Ensuring the long period of the KISS4691 RNG.
Posted:
Aug 22, 2010 2:07 PM
|
|
My posting for KISS4691 has become so mixed with
diverse comments over various newsgroups that I
thought it better to use a separate posting to
bring suggested changes to interested potential
users via sci.math, comp.lang.c, comp.lang.fortran.
At least two readers have pointed out special cases that
might arise from the code I suggested for the KISS4691.
christain.bau pointed out that (x<<13)+c
can exceed 2^32 when the rightmost 19 bits of x are 1's
and c is its maximum possible, a-1=8192=2^13.
Changing
c=(t<x)+(x>>19); to c=(t<=x)+(x>>19); fixes that.
But io_x points out that this overlooks the case x=c=0,
for which c=(t<=x)+(x>>19) would incorrectly make c=1
rather than 0.
This is perhaps best handled by merely, when x=0,
setting Q[j]=c; c=0; return Q[j];
Below is a revised listing of the C version;
changes should be transparent for versions
kindly provided by interested viewers---in Fortran,
PL/1 and others.
While neither correction is likely to have an
appreciable effect on use of KISS4691 as an RNG,
both corrections are necessary to ensure that
the period of MWC() will be the value (> 10^45191)
called for by the underlying mathematics of MWC,
(assuming the two period-1 seed sets,
c=0 and each Q[i]=0
c=8192 and each Q[i]=2^32-1
are not used).
And the results after 10^9 calls to MWC() then
10^9 call to KISS as MWC()+CNG+XS,
confirmed in various versions, are not changed.
George Marsaglia
----------------------------------------------------------
static unsigned long xs=521288629,xcng=362436069,Q[4691];
unsigned long MWC(void) /*238 million/second*/
{unsigned long t,x; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j]; if(x==0) {Q[j]=c; c=0; return Q[j];}
t=(x<<13)+c+x; c=(t<=x)+(x>>19);
return (Q[j]=t);
}
#define CNG ( xcng=69069*xcng+123 )
#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
#define KISS ( MWC()+CNG+XS ) /*138 million/sec*/
#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);
}
---------------------------------------------------------
|
|
