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, a1=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 viewersin 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 period1 seed sets, c=0 and each Q[i]=0 c=8192 and each Q[i]=2^321 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); } 

