From: "robin" <robin51@dodo.com.au>
Subject: Re: Ensuring the long period of the KISS4691 RNG.
Date: Mon, 23 Aug 2010 17:18:30 +1000
Date: 2010-08-23T17:18:30+10:00 [thread overview]
Message-ID: <4c722971$0$1261$c30e37c6@exi-reader.telstra.net> (raw)
In-Reply-To: 08db9776-d933-44a5-a20f-7a137e81cce5@5g2000yqz.googlegroups.com
"geo" <gmarsaglia@gmail.com> wrote in message news:08db9776-d933-44a5-a20f-7a137e81cce5@5g2000yqz.googlegroups.com...
| 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.
Below is my PL/I version, using unsigned arithmetic,
which I confirm delivers the same test results as before.
(NOSIZE, NOFOFL):
RNG: PROCEDURE OPTIONS (MAIN, REORDER);
declare (xs initial (521288629), xcng initial (362436069),
Q(0:4690) ) static fixed binary (32) unsigned;
MWC: procedure () returns (fixed binary (32) unsigned);
declare (t,x,i) fixed binary (32) unsigned;
declare (c initial (0), j initial (4691) ) fixed binary (32) unsigned static;
if j < hbound(Q,1) then j = j + 1; else j = 0;
x = Q(j); if x = 0 then do; Q(j) = c; c = 0; return(Q(j)); end;
t = isll(x,13)+c+x; c = (t<=x) + isrl(x,19);
Q(j)=t;
return (t);
end MWC;
CNG: procedure returns (fixed binary (32) unsigned);
xcng=bin(69069)*xcng+bin(123);
return (xcng);
end CNG;
XXS: procedure returns (fixed binary (32) unsigned);
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 (32) unsigned);
return ( MWC()+CNG+XXS );
end KISS;
declare (i,x) fixed binary (32) unsigned;
Q = CNG + XXS; /* Initialize: */
do i = 1 to 1000000000; x=MWC(); end;
put skip edit (" MWC result=3740121002",x) (a, skip, f(22));
do i = 1 to 1000000000; x=KISS; end;
put skip edit ("KISS result=2224631993",x) (a, skip, f(22));
end RNG;
next parent reply other threads:[~2010-08-23 7:18 UTC|newest]
Thread overview: 4+ messages / expand[flat|nested] mbox.gz Atom feed top
[not found] <08db9776-d933-44a5-a20f-7a137e81cce5@5g2000yqz.googlegroups.com>
2010-08-23 7:18 ` robin [this message]
2010-08-23 8:01 ` Ensuring the long period of the KISS4691 RNG robin
2010-08-23 8:17 ` robin
2010-08-23 17:31 ` io_x
replies disabled
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox