Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » sci.math.* » sci.math

Topic: Ensuring the long period of the KISS4691 RNG.
Replies: 14   Last Post: Aug 30, 2010 6:16 AM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
geo

Posts: 38
Registered: 6/19/08
Ensuring the long period of the KISS4691 RNG.
Posted: Aug 22, 2010 2:07 PM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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);
}
---------------------------------------------------------



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.